Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1751:f3d7e2d7a43f feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 11 Sep 2024 16:26:19 +0200
parents 29b96fc75bee b5690ab5f0b8
children 1f42944d4a72
comparison
equal deleted inserted replaced
1731:3684db043add 1751:f3d7e2d7a43f
81 81
82 The operators required to construct the SAT for imposing a Dirichlet 82 The operators required to construct the SAT for imposing a Dirichlet
83 condition. `H_tuning` and `R_tuning` are used to specify the strength of the 83 condition. `H_tuning` and `R_tuning` are used to specify the strength of the
84 penalty. 84 penalty.
85 85
86 See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref). 86 See also: [`sat`](@ref), [`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
87 """ 87 """
88 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.) 88 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
89 id = boundary(bc) 89 id = boundary(bc)
90 set = Δ.stencil_set 90 set = Δ.stencil_set
91 H⁻¹ = inverse_inner_product(g,set) 91 H⁻¹ = inverse_inner_product(g,set)
92 Hᵧ = inner_product(boundary_grid(g, id), set) 92 Hᵧ = inner_product(boundary_grid(g, id), set)
93 e = boundary_restriction(g, set, id) 93 e = boundary_restriction(g, set, id)
94 d = normal_derivative(g, set, id) 94 d = normal_derivative(g, set, id)
95 B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning) 95 B = positivity_decomposition(Δ, g, boundary(bc); H_tuning, R_tuning)
96 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ 96 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
97 return penalty_tensor, e 97 return penalty_tensor, e
98 end 98 end
99 99
100 """ 100 """
115 penalty_tensor = -H⁻¹∘e'∘Hᵧ 115 penalty_tensor = -H⁻¹∘e'∘Hᵧ
116 return penalty_tensor, d 116 return penalty_tensor, d
117 end 117 end
118 118
119 """ 119 """
120 positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) 120 positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
121 121
122 Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive 122 Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive
123 definite with respect to the boundary quadrature. Here `d` is the normal 123 definite with respect to the boundary quadrature. Here `d` is the normal
124 derivative and `e` is the boundary restriction operator. `B` can then be used 124 derivative and `e` is the boundary restriction operator. `B` can then be used
125 to form a symmetric and energy stable penalty for a Dirichlet condition. The 125 to form a symmetric and energy stable penalty for a Dirichlet condition. The
126 parameters `H_tuning` and `R_tuning` are used to specify the strength of the 126 parameters `H_tuning` and `R_tuning` are used to specify the strength of the
127 penalty and must be greater than 1. For details we refer to 127 penalty and must be greater than 1. For details we refer to
128 https://doi.org/10.1016/j.jcp.2020.109294 128 <https://doi.org/10.1016/j.jcp.2020.109294>
129 """ 129 """
130 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) 130 function positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
131 @assert(H_tuning ≥ 1.) 131 @assert(H_tuning ≥ 1.)
132 @assert(R_tuning ≥ 1.) 132 @assert(R_tuning ≥ 1.)
133 Nτ_H, τ_R = positivity_limits(Δ,g,bc) 133 Nτ_H, τ_R = positivity_limits(Δ,g,b)
134 return H_tuning*Nτ_H + R_tuning*τ_R 134 return H_tuning*Nτ_H + R_tuning*τ_R
135 end 135 end
136 136
137 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then 137 function positivity_limits(Δ::Laplace, g::EquidistantGrid, b::BoundaryIdentifier)
138 # change bc::BoundaryCondition to id::BoundaryIdentifier
139 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
140 h = spacing(g) 138 h = spacing(g)
141 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1]) 139 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
142 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"]) 140 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
143 141
144 τ_H = 1/(h*θ_H) 142 τ_H = one(eltype(Δ))/(h*θ_H)
145 τ_R = 1/(h*θ_R) 143 τ_R = one(eltype(Δ))/(h*θ_R)
146 return τ_H, τ_R 144 return τ_H, τ_R
147 end 145 end
148 146
149 function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition) 147 function positivity_limits(Δ::Laplace, g::TensorGrid, b::BoundaryIdentifier)
150 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc) 148 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(b)], b)
151 return τ_H*ndims(g), τ_R 149 return τ_H*ndims(g), τ_R
152 end 150 end