Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 924:12e8e431b43c feature/laplace_opset

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