Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 875:067a322e4f73 laplace_benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/laplace_opset
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jan 2022 10:55:08 +0100
parents 9929c99754fb 86776d06b883
children
comparison
equal deleted inserted replaced
874:7e9ebd572deb 875:067a322e4f73
1 export Laplace
2 export laplace
3 # REVIEW: Makes more sense to me to have the exports at the top of the file.
4 # Might as well start fixing that.
5
6 # REVIEW: The style of name `Laplace` might clash with other concepts. When
7 # thinking about implementing the variable second derivative I think I will
8 # have to create it as a full TM for the full dimensional problem instead of
9 # building it as a 1D operator and then use that with outer products. The
10 # natural name there would be `VariableSecondDerivative` (or something
11 # similar). But the similarity of the two names would suggest that `Laplace`
12 # and `VariableSecondDerivative` are the same kind of thing, which they
13 # wouldn't be.
14 #
15 # How do we distinguish the kind of type we are implementing here and what we
16 # could potentially do for the variable second derivative?
17 #
18 # I see two ways out:
19 # * Come up with a name for these sets of operators and change `Laplace` accordingly.
20 # * Come up with a name for the bare operators and change `VariableSecondDerivative` accordingly.
21
1 """ 22 """
2 Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim} 23 Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim}
3 Laplace(grid::AbstractGrid, fn; order) 24 Laplace(grid, filename; order)
4 25
5 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a 26 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
6 `TensorMapping`. Additionally, `Laplace` stores the inner product and boundary 27 `TensorMapping`. Additionally, `Laplace` stores the inner product and boundary
7 operators relevant for constructing a SBP finite difference scheme as `TensorMapping`s. 28 operators relevant for constructing a SBP finite difference scheme as a `TensorMapping`.
8 29
9 Laplace(grid, fn; order) creates the Laplace operator defined on `grid`, 30 `Laplace(grid, filename; order)` creates the Laplace operator defined on `grid`,
10 where the operators are read from TOML. The differential operator is created 31 where the operators are read from TOML. The differential operator is created
11 using `laplace(grid::AbstractGrid,...)`. 32 using `laplace(grid,...)`.
12 33
13 Note that all properties of Laplace, excluding the Differential operator `D`, are 34 Note that all properties of Laplace, excluding the differential operator `Laplace.D`, are
14 abstract types. For performance reasons, they should therefore be 35 abstract types. For performance reasons, they should therefore be
15 accessed via the provided getter functions (e.g `inner_product(::Laplace)`). 36 accessed via the provided getter functions (e.g `inner_product(::Laplace)`).
16
17 """ 37 """
18 struct Laplace{T, Dim, TMdiffop<:TensorMapping{T,Dim,Dim}} <: TensorMapping{T,Dim,Dim} 38 struct Laplace{T, Dim, TMdiffop<:TensorMapping{T,Dim,Dim}} <: TensorMapping{T,Dim,Dim}
19 D::TMdiffop # Differential operator 39 D::TMdiffop # Differential operator
20 H::TensorMapping # Inner product operator 40 H::TensorMapping # Inner product operator
21 H_inv::TensorMapping # Inverse inner product operator 41 H_inv::TensorMapping # Inverse inner product operator
22 e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators. 42 e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators.
23 d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators 43 d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators
24 H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product? 44 H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators
25 order::Int
26 end 45 end
27 export Laplace
28 46
29 function Laplace(grid::AbstractGrid, fn; order) 47 function Laplace(grid, filename; order)
48
49 # Read stencils
50 stencil_set = read_stencil_set(filename; order)
30 # TODO: Removed once we can construct the volume and 51 # TODO: Removed once we can construct the volume and
31 # boundary operators by op(grid, fn; order,...). 52 # boundary operators by op(grid, read_stencil_set(fn; order,...)).
32 # Read stencils 53 D_inner_stecil = parse_stencil(stencil_set["D2"]["inner_stencil"])
33 op = read_D2_operator(fn; order) 54 D_closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
34 D_inner_stecil = op.innerStencil 55 H_inner_stencils = parse_scalar(stencil_set["H"]["inner"])
35 D_closure_stencils = op.closureStencils 56 H_closure_stencils = parse_tuple(stencil_set["H"]["closure"])
36 H_closure_stencils = op.quadratureClosure 57 e_closure_stencil = parse_stencil(stencil_set["e"]["closure"])
37 e_closure_stencil = op.eClosure 58 d_closure_stencil = parse_stencil(stencil_set["d1"]["closure"])
38 d_closure_stencil = op.dClosure 59 # REVIEW: Do we add the methods to get rid of this in this branch or a new one?
39 60
40 # Volume operators 61 # Volume operators
41 Δ = laplace(grid, D_inner_stecil, D_closure_stencils) 62 Δ = laplace(grid, D_inner_stecil, D_closure_stencils)
42 H = inner_product(grid, H_closure_stencils) 63 H = inner_product(grid, H_inner_stencils, H_closure_stencils)
43 H⁻¹ = inverse_inner_product(grid, H_closure_stencils) 64 H⁻¹ = inverse_inner_product(grid, H_inner_stencils, H_closure_stencils)
44 65
45 # Boundary operator - id pairs 66 # Boundary operator - id pairs
46 ids = boundary_identifiers(grid) 67 ids = boundary_identifiers(grid)
47 n_ids = length(ids) 68 # REVIEW: Change suggestion: Seems more readable to me but pretty subjective so feel free to ignore
48 e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid,e_closure_stencil,ids[i]),n_ids) 69 e_pairs = map(id -> Pair(id, boundary_restriction(grid, e_closure_stencil, id)), ids)
49 d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids) 70 d_pairs = map(id -> Pair(id, normal_derivative(grid, d_closure_stencil, id)), ids)
50 Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids) 71 Hᵧ_pairs = map(id -> Pair(id, inner_product(boundary_grid(grid, id), H_inner_stencils, H_closure_stencils)), ids)
51 72
52 return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs), order) 73 return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs), order)
53 end 74 end
54 75
55 # TODO: Consider pretty printing of the following form 76 # TODO: Consider pretty printing of the following form
56 # Base.show(io::IO, L::Laplace{T, Dim}) where {T,Dim,TM} = print(io, "Laplace{$T, $Dim, $TM}(", L.D, L.H, L.H_inv, L.e, L.d, L.H_boundary, ")") 77 # Base.show(io::IO, L::Laplace{T, Dim}) where {T,Dim,TM} = print(io, "Laplace{$T, $Dim, $TM}(", L.D, L.H, L.H_inv, L.e, L.d, L.H_boundary, ")")
78 # REVIEW: Should leave a todo here to update this once we have some pretty printing for tensor mappings in general.
57 79
58 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D) 80 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
59 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D) 81 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
60 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...) 82 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
61 83
67 """ 89 """
68 closure_size(L::Laplace) = closure_size(L.D) 90 closure_size(L::Laplace) = closure_size(L.D)
69 export closure_size 91 export closure_size
70 92
71 """ 93 """
72 inner_product(L::Lapalace) 94 inner_product(L::Laplace)
73 95
74 Returns the inner product operator associated with `L` 96 Returns the inner product operator associated with `L`
75
76 """ 97 """
77 inner_product(L::Laplace) = L.H 98 inner_product(L::Laplace) = L.H
78 export inner_product
79 99
80 100
81 """ 101 """
82 inverse_inner_product(L::Lapalace) 102 inverse_inner_product(L::Laplace)
83 103
84 Returns the inverse of the inner product operator associated with `L` 104 Returns the inverse of the inner product operator associated with `L`
85
86 """ 105 """
87 inverse_inner_product(L::Laplace) = L.H_inv 106 inverse_inner_product(L::Laplace) = L.H_inv
88 export inverse_inner_product
89 107
90 108
91 """ 109 """
92 boundary_restriction(L::Lapalace,id::BoundaryIdentifier) 110 boundary_restriction(L::Laplace, id::BoundaryIdentifier)
93 boundary_restriction(L::Lapalace,ids::NTuple{N,BoundaryIdentifier}) 111 boundary_restriction(L::Laplace, ids::Tuple)
94 boundary_restriction(L::Lapalace,ids...) 112 boundary_restriction(L::Laplace, ids...)
95 113
96 Returns boundary restriction operator(s) associated with `L` for the boundary(s) 114 Returns boundary restriction operator(s) associated with `L` for the boundary(s)
97 identified by id(s). 115 identified by id(s).
116 """
117 boundary_restriction(L::Laplace, id::BoundaryIdentifier) = L.e[id]
118 boundary_restriction(L::Laplace, ids::Tuple) = map(id-> L.e[id], ids)
119 boundary_restriction(L::Laplace, ids...) = boundary_restriction(L, ids)
120 # REVIEW: I propose changing the following implementations according to the
121 # above. There are some things we're missing with regards to
122 # `BoundaryIdentifier`, for example we should be able to handle groups of
123 # boundaries as a single `BoundaryIdentifier`. I don't know if we can figure
124 # out the interface for that now or if we save it for the future but either
125 # way these methods will be affected.
98 126
99 """
100 boundary_restriction(L::Laplace,id::BoundaryIdentifier) = L.e[id]
101 boundary_restriction(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N)
102 boundary_restriction(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N)
103 export boundary_restriction
104 127
105 128
106 """ 129 """
107 normal_derivative(L::Lapalace,id::BoundaryIdentifier) 130 normal_derivative(L::Laplace, id::BoundaryIdentifier)
108 normal_derivative(L::Lapalace,ids::NTuple{N,BoundaryIdentifier}) 131 normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
109 normal_derivative(L::Lapalace,ids...) 132 normal_derivative(L::Laplace, ids...)
110 133
111 Returns normal derivative operator(s) associated with `L` for the boundary(s) 134 Returns normal derivative operator(s) associated with `L` for the boundary(s)
112 identified by id(s). 135 identified by id(s).
136 """
137 normal_derivative(L::Laplace, id::BoundaryIdentifier) = L.d[id]
138 normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
139 normal_derivative(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
140
113 141
114 """ 142 """
115 normal_derivative(L::Laplace,id::BoundaryIdentifier) = L.d[id] 143 boundary_quadrature(L::Laplace, id::BoundaryIdentifier)
116 normal_derivative(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N) 144 boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
117 normal_derivative(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N) 145 boundary_quadrature(L::Laplace, ids...)
118 export normal_derivative
119
120
121 # TODO: boundary_inner_product?
122 """
123 boundary_quadrature(L::Lapalace,id::BoundaryIdentifier)
124 boundary_quadrature(L::Lapalace,ids::NTuple{N,BoundaryIdentifier})
125 boundary_quadrature(L::Lapalace,ids...)
126 146
127 Returns boundary quadrature operator(s) associated with `L` for the boundary(s) 147 Returns boundary quadrature operator(s) associated with `L` for the boundary(s)
128 identified by id(s). 148 identified by id(s).
129
130 """ 149 """
131 boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id] 150 boundary_quadrature(L::Laplace, id::BoundaryIdentifier) = L.H_boundary[id]
132 boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N) 151 boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
133 boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N) 152 boundary_quadrature(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
134 export boundary_quadrature
135 153
136 abstract type BoundaryConditionType end 154 abstract type BoundaryConditionType end
137 struct NeumannBC <: BoundaryConditionType end 155 struct NeumannBC <: BoundaryConditionType end
138 struct DirichletBC <: BoundaryConditionType end 156 struct DirichletBC <: BoundaryConditionType end
139 export NeumannBC 157 export NeumannBC
151 closure = tau∘d 169 closure = tau∘d
152 # TODO: Return penalty once we have implemented scalar scaling of the operators. 170 # TODO: Return penalty once we have implemented scalar scaling of the operators.
153 return closure 171 return closure
154 end 172 end
155 173
156 function dirichlet_bc(L::Laplace, id::BoundaryIdentifier)
157 error("Not implemented")
158 # H_inv = inverse_inner_product(L)
159 # e = boundary_restriction(L,id)
160 # d = normal_derivative(L,id)
161 # H_b = boundary_quadrature(L,id)
162 # gamma = borrowing_parameter(L)
163 # tuning = 1.2
164 # S = ScalingOperator(tuning * -1.0/gamma)
165 # tau = H_inv∘(S∘e' + d')∘H_b
166 # closure = tau∘e
167 # penalty = ScalingOperator(-1)∘tau
168 # return (closure, penalty)
169 end
170
171 # function borrowing_parameter(L)
172 # if L.order == 2
173 # return 0.4
174 # elseif L.order == 4
175 # return 0.2508
176 # elseif L.order == 6
177 # return 0.1878
178 # elseif L.order == 8
179 # return 0.0015
180 # elseif L.order == 10
181 # return 0.0351
182 # end
183 # end
184
185 """ 174 """
186 laplace(grid, inner_stencil, closure_stencils) 175 laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
187 176
188 Creates the Laplace operator operator `Δ` as a `TensorMapping` 177 Creates the Laplace operator operator `Δ` as a `TensorMapping`
189 178
190 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,N on the N-dimensional 179 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using
191 `grid`, using the stencil `inner_stencil` in the interior and a set of stencils 180 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
192 `closure_stencils` for the points in the closure regions. 181 for the points in the closure regions.
193 182
194 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a 183 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
195 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s 184 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
196 where the sum is carried out lazily. 185 where the sum is carried out lazily.
197 """ 186 """
198 function laplace(grid::AbstractGrid, inner_stencil, closure_stencils) 187 function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
199 Δ = second_derivative(grid, inner_stencil, closure_stencils, 1) 188 Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
200 for d = 2:dimension(grid) 189 for d = 2:dimension(grid)
201 Δ += second_derivative(grid, inner_stencil, closure_stencils, d) 190 Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
202 end 191 end
203 return Δ 192 return Δ
204 end 193 end
205 export laplace