Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/laplace/laplace.jl @ 312:c1fcc35e19cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Laplace consist of tuples of SecondDerivatives. Begin cleanup of functionallity which should be moved to separate files. WIP
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 09 Sep 2020 21:18:28 +0200
parents 6fa2ba769ae3
children 9cc5d1498b2d
comparison
equal deleted inserted replaced
311:f2d6ec89dfc5 312:c1fcc35e19cb
1 export Laplace
1 """ 2 """
2 Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim} 3 Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
3 4
4 Implements the Laplace operator `L` in Dim dimensions as a tensor operator 5 Implements the Laplace operator `L` in Dim dimensions as a tensor operator
5 The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative 6 The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
6 tensor operators. 7 tensor operators.
7 """ 8 """
8 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim} 9 #export quadrature, inverse_quadrature, boundary_quadrature, boundary_value, normal_derivative
9 D2::NTuple(Dim,SecondDerivative{T,N,M,K}) 10 struct Laplace{Dim,T,N,M,K} <: TensorOperator{T,Dim}
11 D2::NTuple{Dim,SecondDerivative{T,N,M,K}}
10 #TODO: Write a good constructor 12 #TODO: Write a good constructor
11 end 13 end
12 export Laplace
13 14
14 LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) = range_size 15 LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
15 16
16 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} 17 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
17 error("not implemented") 18 error("not implemented")
18 end 19 end
19 20
20 function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} = LazyTensors.apply(L, v, I)
21
22 # u = L*v 21 # u = L*v
23 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T 22 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::Index) where T
24 @inbounds u = apply(L.D2[1],v,I) 23 @inbounds u = LazyTensors.apply(L.D2[1],v,I)
25 return u 24 return u
26 end 25 end
27 26
28 27 # TODO: Fix dispatch on tuples!
29 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T 28 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
30 # 2nd x-derivative 29 # 2nd x-derivative
31 @inbounds vx = view(v, :, Int(I[2])) 30 @inbounds vx = view(v, :, Int(J))
32 @inbounds uᵢ = apply(L.D2[1], vx , I[1]) 31 @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , I)
33 32
34 # 2nd y-derivative 33 # 2nd y-derivative
35 @inbounds vy = view(v, Int(I[1]), :) 34 @inbounds vy = view(v, Int(I), :)
36 @inbounds uᵢ += apply(L.D2[2], vy , I[2]) 35 @inbounds uᵢ += LazyTensors.apply(L.D2[2], vy , J)
37 36
38 return uᵢ 37 return uᵢ
39 end 38 end
40 39
41 quadrature(L::Laplace) = Quadrature(L.op, L.grid) 40 @inline function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
42 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid) 41 return LazyTensors.apply(L, v, I)
43 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
44 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
45 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
46 export quadrature
47
48 """
49 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
50
51 Implements the boundary operator `e` as a TensorMapping
52 """
53 struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
54 op::D2{T,N,M,K}
55 grid::EquidistantGrid{2}
56 bId::CartesianBoundary
57 end
58 export BoundaryValue
59
60 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
61 # Can we give special treatment to TensorMappings that go to a higher dim?
62 function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T
63 if dim(e.bId) == 1
64 return (UnknownDim, domain_size[1])
65 elseif dim(e.bId) == 2
66 return (domain_size[1], UnknownDim)
67 end
68 end
69 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
70 # TODO: Make a nicer solution for 3-dim(e.bId)
71
72 # TODO: Make this independent of dimension
73 function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
74 i = I[dim(e.bId)]
75 j = I[3-dim(e.bId)]
76 N_i = size(e.grid)[dim(e.bId)]
77 return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
78 end 42 end
79 43
80 function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T 44 # quadrature(L::Laplace) = Quadrature(L.op, L.grid)
81 u = selectdim(v,3-dim(e.bId),Int(I[1])) 45 # inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
82 return apply_boundary_value_transpose(e.op, u, region(e.bId)) 46 # boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
83 end 47 # normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
84 48 # boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
85 """ 49 # export NormalDerivative
86 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 50 # """
87 51 # NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
88 Implements the boundary operator `d` as a TensorMapping 52 #
89 """ 53 # Implements the boundary operator `d` as a TensorMapping
90 struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 54 # """
91 op::D2{T,N,M,K} 55 # struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
92 grid::EquidistantGrid{2} 56 # op::D2{T,N,M,K}
93 bId::CartesianBoundary 57 # grid::EquidistantGrid{2}
94 end 58 # bId::CartesianBoundary
95 export NormalDerivative 59 # end
96 60 #
97 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue? 61 # # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
98 # Can we give special treatment to TensorMappings that go to a higher dim? 62 # # Can we give special treatment to TensorMappings that go to a higher dim?
99 function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer}) 63 # function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer})
100 if dim(e.bId) == 1 64 # if dim(e.bId) == 1
101 return (UnknownDim, domain_size[1]) 65 # return (UnknownDim, domain_size[1])
102 elseif dim(e.bId) == 2 66 # elseif dim(e.bId) == 2
103 return (domain_size[1], UnknownDim) 67 # return (domain_size[1], UnknownDim)
104 end 68 # end
105 end 69 # end
106 LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],) 70 # LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],)
107 71 #
108 # TODO: Not type stable D:< 72 # # TODO: Not type stable D:<
109 # TODO: Make this independent of dimension 73 # # TODO: Make this independent of dimension
110 function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T 74 # function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
111 i = I[dim(d.bId)] 75 # i = I[dim(d.bId)]
112 j = I[3-dim(d.bId)] 76 # j = I[3-dim(d.bId)]
113 N_i = size(d.grid)[dim(d.bId)] 77 # N_i = size(d.grid)[dim(d.bId)]
114 h_inv = inverse_spacing(d.grid)[dim(d.bId)] 78 # h_inv = inverse_spacing(d.grid)[dim(d.bId)]
115 return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId)) 79 # return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
116 end 80 # end
117 81 #
118 function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T 82 # function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
119 u = selectdim(v,3-dim(d.bId),Int(I[1])) 83 # u = selectdim(v,3-dim(d.bId),Int(I[1]))
120 return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId)) 84 # return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
121 end 85 # end
122 86 #
123 """ 87 # """
124 BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} 88 # BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
125 89 #
126 Implements the boundary operator `q` as a TensorOperator 90 # Implements the boundary operator `q` as a TensorOperator
127 """ 91 # """
128 struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} 92 # export BoundaryQuadrature
129 op::D2{T,N,M,K} 93 # struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
130 grid::EquidistantGrid{2} 94 # op::D2{T,N,M,K}
131 bId::CartesianBoundary 95 # grid::EquidistantGrid{2}
132 end 96 # bId::CartesianBoundary
133 export BoundaryQuadrature 97 # end
134 98 #
135 # TODO: Make this independent of dimension 99 #
136 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T 100 # # TODO: Make this independent of dimension
137 h = spacing(q.grid)[3-dim(q.bId)] 101 # function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T
138 N = size(v) 102 # h = spacing(q.grid)[3-dim(q.bId)]
139 return apply_quadrature(q.op, h, v[I[1]], I[1], N[1]) 103 # N = size(v)
140 end 104 # return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
141 105 # end
142 LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I) 106 #
143 107 # LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I)
144 108 #
145 109 #
146 110 #
147 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 111 #
148 112 # struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
149 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} 113 #
150 e = boundary_value(L, Bid()) 114 # function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
151 d = normal_derivative(L, Bid()) 115 # e = boundary_value(L, Bid())
152 Hᵧ = boundary_quadrature(L, Bid()) 116 # d = normal_derivative(L, Bid())
153 H⁻¹ = inverse_quadrature(L) 117 # Hᵧ = boundary_quadrature(L, Bid())
154 return (-H⁻¹*e*Hᵧ*(d'*v - g))[I] 118 # H⁻¹ = inverse_quadrature(L)
155 end 119 # return (-H⁻¹*e*Hᵧ*(d'*v - g))[I]
156 120 # end
157 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition 121 #
158 tau::Float64 122 # struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
159 end 123 # tau::Float64
160 124 # end
161 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} 125 #
162 e = boundary_value(L, Bid()) 126 # function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
163 d = normal_derivative(L, Bid()) 127 # e = boundary_value(L, Bid())
164 Hᵧ = boundary_quadrature(L, Bid()) 128 # d = normal_derivative(L, Bid())
165 H⁻¹ = inverse_quadrature(L) 129 # Hᵧ = boundary_quadrature(L, Bid())
166 return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I] 130 # H⁻¹ = inverse_quadrature(L)
167 # Need to handle scalar multiplication and addition of TensorMapping 131 # return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I]
168 end 132 # # Need to handle scalar multiplication and addition of TensorMapping
133 # end
169 134
170 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D 135 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
171 # return apply(s.L, v, i) + 136 # return apply(s.L, v, i) +
172 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + 137 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
173 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + 138 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +