Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1858:4a9be96f2569 feature/documenter_logo

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 12 Jan 2025 21:18:44 +0100
parents b5690ab5f0b8
children f3d7e2d7a43f
comparison
equal deleted inserted replaced
1857:ffde7dad9da5 1858:4a9be96f2569
1 """ 1 """
2 laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) 2 Laplace{T, Dim, TM} <: LazyTensor{T, Dim, Dim}
3 3
4 Creates the Laplace operator operator `Δ` as a `TensorMapping` 4 The Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
5 `LazyTensor`.
6 """
7 struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim}
8 D::TM # Difference operator
9 stencil_set::StencilSet # Stencil set of the operator
10 end
5 11
6 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using 12 """
7 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` 13 Laplace(g::Grid, stencil_set::StencilSet)
8 for the points in the closure regions.
9 14
10 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a 15 Creates the `Laplace` operator `Δ` on `g` given `stencil_set`.
11 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s 16
12 where the sum is carried out lazily. 17 See also [`laplace`](@ref).
13 """ 18 """
14 function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils) 19 function Laplace(g::Grid, stencil_set::StencilSet)
15 Δ = second_derivative(grid, inner_stencil, closure_stencils, 1) 20 Δ = laplace(g, stencil_set)
16 for d = 2:dimension(grid) 21 return Laplace(Δ, stencil_set)
17 Δ += second_derivative(grid, inner_stencil, closure_stencils, d) 22 end
23
24 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
25 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
26 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
27
28 # TODO: Implement pretty printing of Laplace once pretty printing of LazyTensors is implemented.
29 # Base.show(io::IO, L::Laplace) = ...
30
31 """
32 laplace(g::Grid, stencil_set)
33
34 Creates the Laplace operator operator `Δ` as a `LazyTensor` on `g`.
35
36 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `g`. The
37 approximation depends on the type of grid and the stencil set.
38
39 See also: [`second_derivative`](@ref).
40 """
41 function laplace end
42 function laplace(g::TensorGrid, stencil_set)
43 # return mapreduce(+, enumerate(g.grids)) do (i, gᵢ)
44 # Δᵢ = laplace(gᵢ, stencil_set)
45 # LazyTensors.inflate(Δᵢ, size(g), i)
46 # end
47
48 Δ = LazyTensors.inflate(laplace(g.grids[1], stencil_set), size(g), 1)
49 for d = 2:ndims(g)
50 Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d)
18 end 51 end
19 return Δ 52 return Δ
20 end 53 end
21 export laplace 54 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
55
56 """
57 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
58
59 The operators required to construct the SAT for imposing a Dirichlet
60 condition. `H_tuning` and `R_tuning` are used to specify the strength of the
61 penalty.
62
63 See also: [`sat`](@ref), [`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
64 """
65 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
66 id = boundary(bc)
67 set = Δ.stencil_set
68 H⁻¹ = inverse_inner_product(g,set)
69 Hᵧ = inner_product(boundary_grid(g, id), set)
70 e = boundary_restriction(g, set, id)
71 d = normal_derivative(g, set, id)
72 B = positivity_decomposition(Δ, g, boundary(bc); H_tuning, R_tuning)
73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
74 return penalty_tensor, e
75 end
76
77 """
78 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
79
80 The operators required to construct the SAT for imposing a Neumann condition.
81
82 See also: [`sat`](@ref), [`NeumannCondition`](@ref).
83 """
84 function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
85 id = boundary(bc)
86 set = Δ.stencil_set
87 H⁻¹ = inverse_inner_product(g,set)
88 Hᵧ = inner_product(boundary_grid(g, id), set)
89 e = boundary_restriction(g, set, id)
90 d = normal_derivative(g, set, id)
91
92 penalty_tensor = -H⁻¹∘e'∘Hᵧ
93 return penalty_tensor, d
94 end
95
96 """
97 positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
98
99 Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive
100 definite with respect to the boundary quadrature. Here `d` is the normal
101 derivative and `e` is the boundary restriction operator. `B` can then be used
102 to form a symmetric and energy stable penalty for a Dirichlet condition. The
103 parameters `H_tuning` and `R_tuning` are used to specify the strength of the
104 penalty and must be greater than 1. For details we refer to
105 <https://doi.org/10.1016/j.jcp.2020.109294>
106 """
107 function positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
108 @assert(H_tuning ≥ 1.)
109 @assert(R_tuning ≥ 1.)
110 Nτ_H, τ_R = positivity_limits(Δ,g,b)
111 return H_tuning*Nτ_H + R_tuning*τ_R
112 end
113
114 function positivity_limits(Δ::Laplace, g::EquidistantGrid, b::BoundaryIdentifier)
115 h = spacing(g)
116 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
117 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
118
119 τ_H = one(eltype(Δ))/(h*θ_H)
120 τ_R = one(eltype(Δ))/(h*θ_R)
121 return τ_H, τ_R
122 end
123
124 function positivity_limits(Δ::Laplace, g::TensorGrid, b::BoundaryIdentifier)
125 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(b)], b)
126 return τ_H*ndims(g), τ_R
127 end