Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 667:f3a0d1f7d842 feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Laplace a type storing relevant operators used when writing a scheme, e.g. quadratures, normal derivatives etc.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 31 Jan 2021 21:04:02 +0100
parents d6edde60909b
children 2d56a53a1646
comparison
equal deleted inserted replaced
666:94fe0761e6f9 667:f3a0d1f7d842
1 """ 1 """
2 Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) 2 Laplace{T,Dim,...}()
3 Laplace(grid::EquidistantGrid, fn; order)
4
5 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
6 `TensorMapping`. Additionally, `Laplace` stores the quadrature, and boundary
7 operators relevant for constructing a SBP finite difference scheme as `TensorMapping`s.
8 """
9 struct Laplace{T, Dim, Rb, TMdiffop<:TensorMapping{T,Dim,Dim}, # Differential operator tensor mapping
10 TMqop<:TensorMapping{T,Dim,Dim}, # Quadrature operator tensor mapping
11 TMbop<:TensorMapping{T,Rb,Dim}, # Boundary operator tensor mapping
12 TMbqop<:TensorMapping{T,Rb,Rb}, # Boundary quadrature tensor mapping
13 BID<:BoundaryIdentifier} <: TensorMapping{T,Dim,Dim}
14 D::TMdiffop # Difference operator
15 H::TMqop # Quadrature (norm) operator
16 H_inv::TMqop # Inverse quadrature (norm) operator
17 e::Dict{BID,TMbop} # Boundary restriction operators
18 d::Dict{BID,TMbop} # Normal derivative operators
19 H_boundary::Dict{BID,TMbqop} # Boundary quadrature operators
20 end
21 export Laplace
22
23 function Laplace(grid::EquidistantGrid, fn; order)
24 # TODO: Removed once we can construct the volume and
25 # boundary operators by op(grid, fn; order,...).
26 # Read stencils
27 op = read_D2_operator(fn; order)
28 D_inner_stecil = op.innerStencil
29 D_closure_stencils = op.closureStencils
30 H_closure_stencils = op.quadratureClosure
31 e_closure_stencil = op.eClosure
32 d_closure_stencil = op.dClosure
33
34 # Volume operators
35 Δ = laplace(grid, D_inner_stecil, D_closure_stencils)
36 H = DiagonalQuadrature(grid, H_closure_stencils)
37 H⁻¹ = InverseDiagonalQuadrature(grid, H_closure_stencils)
38
39 # Pair boundary operators and boundary quadratures with the boundary ids
40 e_pairs = ()
41 d_pairs = ()
42 Hᵧ_pairs = ()
43 dims = collect(1:dimension(grid))
44 for id ∈ boundary_identifiers(grid)
45 # Boundary operators
46 e_pairs = (e_pairs...,Pair(id,BoundaryRestriction(grid,e_closure_stencil,id)))
47 d_pairs = (d_pairs...,Pair(id,NormalDerivative(grid,d_closure_stencil,id)))
48 # Boundary quadratures
49 # Construct these on the lower-dimensional grid in the
50 # coordinite directions orthogonal to dim(id)
51 orth_dims = dims[dims .!= dim(id)]
52 orth_grid = restrict(grid,orth_dims)
53 Hᵧ_pairs = (Hᵧ_pairs...,Pair(id,DiagonalQuadrature(orth_grid,H_closure_stencils)))
54 end
55 return Laplace(Δ, H, H⁻¹, Dict(e_pairs), Dict(d_pairs), Dict(Hᵧ_pairs))
56 end
57
58 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
59 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
60 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
61
62 quadrature(L::Laplace) = L.H
63 export quadrature
64 inverse_quadrature(L::Laplace) = L.H_inv
65 export inverse_quadrature
66 boundary_restriction(L::Laplace,bid::BoundaryIdentifier) = L.e[bid]
67 export boundary_restriction
68 normal_derivative(L::Laplace,bid::BoundaryIdentifier) = L.d[bid]
69 export normal_derivative
70 boundary_quadrature(L::Laplace,bid::BoundaryIdentifier) = L.H_boundary[bid]
71 export boundary_quadrature
72
73 """
74 laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
3 75
4 Creates the Laplace operator operator `Δ` as a `TensorMapping` 76 Creates the Laplace operator operator `Δ` as a `TensorMapping`
5 77
6 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using 78 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using
7 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` 79 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
8 for the points in the closure regions. 80 for the points in the closure regions.
9 81
10 On a one-dimensional `grid`, `Δ` is a `SecondDerivative`. On a multi-dimensional `grid`, `Δ` is the sum of 82 On a one-dimensional `grid`, `Δ` is a `SecondDerivative`. On a multi-dimensional `grid`, `Δ` is the sum of
11 multi-dimensional `SecondDerivative`s where the sum is carried out lazily. 83 multi-dimensional `SecondDerivative`s where the sum is carried out lazily.
12 """ 84 """
13 function Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim 85 function laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
14 Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1) 86 Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1)
15 for d = 2:Dim 87 for d = 2:Dim
16 Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d) 88 Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d)
17 end 89 end
18 return Δ 90 return Δ
19 end 91 end
20 export Laplace 92 export laplace