Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 291:0f94dc29c4bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in branch boundary_conditions
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 22 Jun 2020 21:43:05 +0200
parents 7247e85dc1e8
children
comparison
equal deleted inserted replaced
231:fbabfd4e8f20 291:0f94dc29c4bf
1 struct NormalDerivative{N,M,K} 1 #TODO: move to sbpoperators.jl
2 op::D2{Float64,N,M,K} 2 """
3 grid::EquidistantGrid 3 Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
4 bId::CartesianBoundary 4
5 end 5 Implements the Laplace operator `L` in Dim dimensions as a tensor operator
6 6 The multi-dimensional tensor operator simply consists of a tuple of the 1D
7 function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer) 7 Laplace tensor operator as defined by ConstantLaplaceOperator.
8 u = selectdim(v,3-dim(d.bId),I) 8 """
9 return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId)) 9 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
10 end 10 tensorOps::NTuple(Dim,ConstantLaplaceOperator{T,N,M,K})
11 11 #TODO: Write a good constructor
12 # Not correct abstraction level 12 end
13 # TODO: Not type stable D:< 13 export Laplace
14 function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer}) 14
15 i = I[dim(d.bId)] 15 LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) = range_size
16 j = I[3-dim(d.bId)] 16
17 N_i = d.grid.size[dim(d.bId)] 17 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
18
19 r = getregion(i, closureSize(d.op), N_i)
20
21 if r != region(d.bId)
22 return 0
23 end
24
25 if r == Lower
26 # Note, closures are indexed by offset. Fix this D:<
27 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
28 elseif r == Upper
29 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-j]*v[j]
30 end
31 end
32
33 struct BoundaryValue{N,M,K}
34 op::D2{Float64,N,M,K}
35 grid::EquidistantGrid
36 bId::CartesianBoundary
37 end
38
39 function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer})
40 i = I[dim(e.bId)]
41 j = I[3-dim(e.bId)]
42 N_i = e.grid.size[dim(e.bId)]
43
44 r = getregion(i, closureSize(e.op), N_i)
45
46 if r != region(e.bId)
47 return 0
48 end
49
50 if r == Lower
51 # Note, closures are indexed by offset. Fix this D:<
52 return e.op.eClosure[i-1]*v[j]
53 elseif r == Upper
54 return e.op.eClosure[N_i-j]*v[j]
55 end
56 end
57
58 function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer)
59 u = selectdim(v,3-dim(e.bId),I)
60 return apply_e(e.op, u, region(e.bId))
61 end
62
63 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
64 grid::EquidistantGrid{Dim,T}
65 a::T
66 op::D2{Float64,N,M,K}
67 # e::BoundaryValue
68 # d::NormalDerivative
69 end
70
71 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
72 error("not implemented") 18 error("not implemented")
73 end 19 end
74 20
75 # u = L*v 21 # u = L*v
76 function apply(L::Laplace{1}, v::AbstractVector, i::Int) 22 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
77 uᵢ = L.a * SbpOperators.apply(L.op, L.grid.spacing[1], v, i) 23 return apply(L.tensorOps[1],v,I)
78 return uᵢ 24 end
79 end 25
80 26
81 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2} 27 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
82 # 2nd x-derivative 28 # 2nd x-derivative
83 @inbounds vx = view(v, :, Int(I[2])) 29 @inbounds vx = view(v, :, Int(I[2]))
84 @inbounds uᵢ = L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[1], vx , I[1]) 30 @inbounds uᵢ = apply(L.tensorOps[1], vx , (I[1],)) #Tuple conversion here is ugly. How to do it? Should we use indexing here?
31
85 # 2nd y-derivative 32 # 2nd y-derivative
86 @inbounds vy = view(v, Int(I[1]), :) 33 @inbounds vy = view(v, Int(I[1]), :)
87 @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2]) 34 @inbounds uᵢ += apply(L.tensorOps[2], vy , (I[2],)) #Tuple conversion here is ugly. How to do it?
88 # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors 35
89 return uᵢ 36 return uᵢ
90 end 37 end
91 38
92 # Slow but maybe convenient? 39 quadrature(L::Laplace) = Quadrature(L.op, L.grid)
93 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2}) 40 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
94 I = Index{Unknown}.(Tuple(i)) 41 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
95 apply(L, v, I) 42 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
96 end 43 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
44 export quadrature
45
46 # At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
47 """
48 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
49
50 Implements the quadrature operator `H` of Dim dimension as a TensorMapping
51 """
52 struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
53 op::D2{T,N,M,K}
54 grid::EquidistantGrid{Dim,T}
55 end
56 export Quadrature
57
58 LazyTensors.domain_size(H::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
59
60 @inline function LazyTensors.apply(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
61 N = size(H.grid)
62 # Quadrature in x direction
63 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1])
64 # Quadrature in y-direction
65 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2])
66 return q
67 end
68
69 LazyTensors.apply_transpose(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H,v,I)
70
71
72 """
73 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
74
75 Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
76 """
77 struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
78 op::D2{T,N,M,K}
79 grid::EquidistantGrid{Dim,T}
80 end
81 export InverseQuadrature
82
83 LazyTensors.domain_size(H_inv::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
84
85 @inline function LazyTensors.apply(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
86 N = size(H_inv.grid)
87 # Inverse quadrature in x direction
88 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1])
89 # Inverse quadrature in y-direction
90 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[2], q_inv, I[2], N[2])
91 return q_inv
92 end
93
94 LazyTensors.apply_transpose(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H_inv,v,I)
95
96 """
97 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
98
99 Implements the boundary operator `e` as a TensorMapping
100 """
101 struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
102 op::D2{T,N,M,K}
103 grid::EquidistantGrid{2}
104 bId::CartesianBoundary
105 end
106 export BoundaryValue
107
108 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
109 # Can we give special treatment to TensorMappings that go to a higher dim?
110 function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T
111 if dim(e.bId) == 1
112 return (UnknownDim, domain_size[1])
113 elseif dim(e.bId) == 2
114 return (domain_size[1], UnknownDim)
115 end
116 end
117 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
118 # TODO: Make a nicer solution for 3-dim(e.bId)
119
120 # TODO: Make this independent of dimension
121 function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
122 i = I[dim(e.bId)]
123 j = I[3-dim(e.bId)]
124 N_i = size(e.grid)[dim(e.bId)]
125 return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
126 end
127
128 function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
129 u = selectdim(v,3-dim(e.bId),Int(I[1]))
130 return apply_boundary_value_transpose(e.op, u, region(e.bId))
131 end
132
133 """
134 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
135
136 Implements the boundary operator `d` as a TensorMapping
137 """
138 struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
139 op::D2{T,N,M,K}
140 grid::EquidistantGrid{2}
141 bId::CartesianBoundary
142 end
143 export NormalDerivative
144
145 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
146 # Can we give special treatment to TensorMappings that go to a higher dim?
147 function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer})
148 if dim(e.bId) == 1
149 return (UnknownDim, domain_size[1])
150 elseif dim(e.bId) == 2
151 return (domain_size[1], UnknownDim)
152 end
153 end
154 LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],)
155
156 # TODO: Not type stable D:<
157 # TODO: Make this independent of dimension
158 function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
159 i = I[dim(d.bId)]
160 j = I[3-dim(d.bId)]
161 N_i = size(d.grid)[dim(d.bId)]
162 h_inv = inverse_spacing(d.grid)[dim(d.bId)]
163 return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
164 end
165
166 function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
167 u = selectdim(v,3-dim(d.bId),Int(I[1]))
168 return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
169 end
170
171 """
172 BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
173
174 Implements the boundary operator `q` as a TensorOperator
175 """
176 struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
177 op::D2{T,N,M,K}
178 grid::EquidistantGrid{2}
179 bId::CartesianBoundary
180 end
181 export BoundaryQuadrature
182
183 # TODO: Make this independent of dimension
184 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T
185 h = spacing(q.grid)[3-dim(q.bId)]
186 N = size(v)
187 return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
188 end
189
190 LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I)
191
192
97 193
98 194
99 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 195 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
100 196
101 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} 197 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
102 e = BoundaryValue(L.op, L.grid, Bid()) 198 e = boundary_value(L, Bid())
103 d = NormalDerivative(L.op, L.grid, Bid()) 199 d = normal_derivative(L, Bid())
104 Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid()) 200 Hᵧ = boundary_quadrature(L, Bid())
105 # TODO: Implement BoundaryQuadrature method 201 H⁻¹ = inverse_quadrature(L)
106 202 return (-H⁻¹*e*Hᵧ*(d'*v - g))[I]
107 return -L.Hi*e*Hᵧ*(d'*v - g)
108 # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on
109 end 203 end
110 204
111 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition 205 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
112 tau::Float64 206 tau::Float64
113 end 207 end
114 208
115 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} 209 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
116 e = BoundaryValue(L.op, L.grid, Bid()) 210 e = boundary_value(L, Bid())
117 d = NormalDerivative(L.op, L.grid, Bid()) 211 d = normal_derivative(L, Bid())
118 Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid()) 212 Hᵧ = boundary_quadrature(L, Bid())
119 # TODO: Implement BoundaryQuadrature method 213 H⁻¹ = inverse_quadrature(L)
120 214 return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I]
121 return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g)
122 # Need to handle scalar multiplication and addition of TensorMapping 215 # Need to handle scalar multiplication and addition of TensorMapping
123 end 216 end
124 217
125 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D 218 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
126 # return apply(s.L, v, i) + 219 # return apply(s.L, v, i) +
127 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + 220 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
128 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + 221 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
129 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + 222 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
130 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) 223 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i)
131 # end 224 # end