Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 866:1784b1c0af3e feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 19 Jan 2022 14:44:24 +0100
parents 1970ebceabe4 b4acd25943f4
children 4bd35ba8f34a
comparison
equal deleted inserted replaced
865:545a6c1a0a0e 866:1784b1c0af3e
1 """ 1 """
2 Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim} 2 Laplace{T, Dim, TMdiffop} <: TensorMapping{T,Dim,Dim}
3 Laplace(grid::AbstractGrid, fn; order) 3 Laplace(grid, filename; order)
4 4
5 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a 5 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
6 `TensorMapping`. Additionally, `Laplace` stores the inner product and boundary 6 `TensorMapping`. Additionally, `Laplace` stores the inner product and boundary
7 operators relevant for constructing a SBP finite difference scheme as `TensorMapping`s. 7 operators relevant for constructing a SBP finite difference scheme as a `TensorMapping`.
8 8
9 Laplace(grid, fn; order) creates the Laplace operator defined on `grid`, 9 `Laplace(grid, filename; order)` creates the Laplace operator defined on `grid`,
10 where the operators are read from TOML. The differential operator is created 10 where the operators are read from TOML. The differential operator is created
11 using `laplace(grid::AbstractGrid,...)`. 11 using `laplace(grid,...)`.
12 12
13 Note that all properties of Laplace, excluding the Differential operator `D`, are 13 Note that all properties of Laplace, excluding the differential operator `Laplace.D`, are
14 abstract types. For performance reasons, they should therefore be 14 abstract types. For performance reasons, they should therefore be
15 accessed via the provided getter functions (e.g `inner_product(::Laplace)`). 15 accessed via the provided getter functions (e.g `inner_product(::Laplace)`).
16 16
17 """ 17 """
18 struct Laplace{T, Dim, TMdiffop<:TensorMapping{T,Dim,Dim}} <: TensorMapping{T,Dim,Dim} 18 struct Laplace{T, Dim, TMdiffop<:TensorMapping{T,Dim,Dim}} <: TensorMapping{T,Dim,Dim}
19 D::TMdiffop # Differential operator 19 D::TMdiffop # Differential operator
20 H::TensorMapping # Inner product operator 20 H::TensorMapping # Inner product operator
21 H_inv::TensorMapping # Inverse inner product operator 21 H_inv::TensorMapping # Inverse inner product operator
22 e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators. 22 e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators.
23 d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators 23 d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators
24 H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product? 24 H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators
25 end 25 end
26 export Laplace 26 export Laplace
27 27
28 function Laplace(grid::AbstractGrid, fn; order) 28 function Laplace(grid, filename; order)
29
30 # Read stencils
31 stencil_set = read_stencil_set(filename; order)
29 # TODO: Removed once we can construct the volume and 32 # TODO: Removed once we can construct the volume and
30 # boundary operators by op(grid, fn; order,...). 33 # boundary operators by op(grid, read_stencil_set(fn; order,...)).
31 # Read stencils 34 D_inner_stecil = parse_stencil(stencil_set["D2"]["inner_stencil"])
32 op = read_D2_operator(fn; order) 35 D_closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
33 D_inner_stecil = op.innerStencil 36 H_inner_stencils = parse_scalar(stencil_set["H"]["inner"])
34 D_closure_stencils = op.closureStencils 37 H_closure_stencils = parse_tuple(stencil_set["H"]["closure"])
35 H_closure_stencils = op.quadratureClosure 38 e_closure_stencil = parse_stencil(stencil_set["e"]["closure"])
36 e_closure_stencil = op.eClosure 39 d_closure_stencil = parse_stencil(stencil_set["d1"]["closure"])
37 d_closure_stencil = op.dClosure
38 40
39 # Volume operators 41 # Volume operators
40 Δ = laplace(grid, D_inner_stecil, D_closure_stencils) 42 Δ = laplace(grid, D_inner_stecil, D_closure_stencils)
41 H = inner_product(grid, H_closure_stencils) 43 H = inner_product(grid, H_inner_stencils, H_closure_stencils)
42 H⁻¹ = inverse_inner_product(grid, H_closure_stencils) 44 H⁻¹ = inverse_inner_product(grid, H_inner_stencils, H_closure_stencils)
43 45
44 # Boundary operator - id pairs 46 # Boundary operator - id pairs
45 ids = boundary_identifiers(grid) 47 ids = boundary_identifiers(grid)
46 n_ids = length(ids) 48 n_ids = length(ids)
47 e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid,e_closure_stencil,ids[i]),n_ids) 49 e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid, e_closure_stencil, ids[i]), n_ids)
48 d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids) 50 d_pairs = ntuple(i -> ids[i] => normal_derivative(grid, d_closure_stencil, ids[i]), n_ids)
49 Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids) 51 Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid, ids[i]), H_inner_stencils, H_closure_stencils), n_ids)
50 52
51 return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs)) 53 return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs))
52 end 54 end
53 55
54 # TODO: Consider pretty printing of the following form 56 # TODO: Consider pretty printing of the following form
58 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D) 60 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
59 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...) 61 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
60 62
61 63
62 """ 64 """
63 inner_product(L::Lapalace) 65 inner_product(L::Laplace)
64 66
65 Returns the inner product operator associated with `L` 67 Returns the inner product operator associated with `L`
66 68
67 """ 69 """
68 inner_product(L::Laplace) = L.H 70 inner_product(L::Laplace) = L.H
69 export inner_product 71 export inner_product
70 72
71 73
72 """ 74 """
73 inverse_inner_product(L::Lapalace) 75 inverse_inner_product(L::Laplace)
74 76
75 Returns the inverse of the inner product operator associated with `L` 77 Returns the inverse of the inner product operator associated with `L`
76 78
77 """ 79 """
78 inverse_inner_product(L::Laplace) = L.H_inv 80 inverse_inner_product(L::Laplace) = L.H_inv
79 export inverse_inner_product 81 export inverse_inner_product
80 82
81 83
82 """ 84 """
83 boundary_restriction(L::Lapalace,id::BoundaryIdentifier) 85 boundary_restriction(L::Laplace, id::BoundaryIdentifier)
84 boundary_restriction(L::Lapalace,ids::NTuple{N,BoundaryIdentifier}) 86 boundary_restriction(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
85 boundary_restriction(L::Lapalace,ids...) 87 boundary_restriction(L::Laplace, ids...)
86 88
87 Returns boundary restriction operator(s) associated with `L` for the boundary(s) 89 Returns boundary restriction operator(s) associated with `L` for the boundary(s)
88 identified by id(s). 90 identified by id(s).
89 91
90 """ 92 """
91 boundary_restriction(L::Laplace,id::BoundaryIdentifier) = L.e[id] 93 boundary_restriction(L::Laplace, id::BoundaryIdentifier) = L.e[id]
92 boundary_restriction(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N) 94 boundary_restriction(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.e[ids[i]],N)
93 boundary_restriction(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N) 95 boundary_restriction(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.e[ids[i]],N)
94 export boundary_restriction 96 export boundary_restriction
95 97
96 98
97 """ 99 """
98 normal_derivative(L::Lapalace,id::BoundaryIdentifier) 100 normal_derivative(L::Laplace, id::BoundaryIdentifier)
99 normal_derivative(L::Lapalace,ids::NTuple{N,BoundaryIdentifier}) 101 normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
100 normal_derivative(L::Lapalace,ids...) 102 normal_derivative(L::Laplace, ids...)
101 103
102 Returns normal derivative operator(s) associated with `L` for the boundary(s) 104 Returns normal derivative operator(s) associated with `L` for the boundary(s)
103 identified by id(s). 105 identified by id(s).
104 106
105 """ 107 """
106 normal_derivative(L::Laplace,id::BoundaryIdentifier) = L.d[id] 108 normal_derivative(L::Laplace, id::BoundaryIdentifier) = L.d[id]
107 normal_derivative(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N) 109 normal_derivative(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.d[ids[i]],N)
108 normal_derivative(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N) 110 normal_derivative(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.d[ids[i]],N)
109 export normal_derivative 111 export normal_derivative
110 112
111 113
112 # TODO: boundary_inner_product?
113 """ 114 """
114 boundary_quadrature(L::Lapalace,id::BoundaryIdentifier) 115 boundary_quadrature(L::Laplace, id::BoundaryIdentifier)
115 boundary_quadrature(L::Lapalace,ids::NTuple{N,BoundaryIdentifier}) 116 boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier})
116 boundary_quadrature(L::Lapalace,ids...) 117 boundary_quadrature(L::Laplace, ids...)
117 118
118 Returns boundary quadrature operator(s) associated with `L` for the boundary(s) 119 Returns boundary quadrature operator(s) associated with `L` for the boundary(s)
119 identified by id(s). 120 identified by id(s).
120 121
121 """ 122 """
122 boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id] 123 boundary_quadrature(L::Laplace, id::BoundaryIdentifier) = L.H_boundary[id]
123 boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N) 124 boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
124 boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N) 125 boundary_quadrature(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
125 export boundary_quadrature 126 export boundary_quadrature
126 127
127 128
128 """ 129 """
129 laplace(grid, inner_stencil, closure_stencils) 130 laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
130 131
131 Creates the Laplace operator operator `Δ` as a `TensorMapping` 132 Creates the Laplace operator operator `Δ` as a `TensorMapping`
132 133
133 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,N on the N-dimensional 134 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using
134 `grid`, using the stencil `inner_stencil` in the interior and a set of stencils 135 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
135 `closure_stencils` for the points in the closure regions. 136 for the points in the closure regions.
136 137
137 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a 138 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
138 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s 139 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
139 where the sum is carried out lazily. 140 where the sum is carried out lazily.
140 """ 141 """
141 function laplace(grid::AbstractGrid, inner_stencil, closure_stencils) 142 function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
142 Δ = second_derivative(grid, inner_stencil, closure_stencils, 1) 143 Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
143 for d = 2:dimension(grid) 144 for d = 2:dimension(grid)
144 Δ += second_derivative(grid, inner_stencil, closure_stencils, d) 145 Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
145 end 146 end
146 return Δ 147 return Δ