Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1673:f28c92ec843c refactor/grids/boundary_identifiers_1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Address TODO in laplace.jl + minor fixes
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sun, 07 Jul 2024 15:19:33 -0700
parents 1937be9502a7
children b5690ab5f0b8
comparison
equal deleted inserted replaced
1672:3714a391545a 1673:f28c92ec843c
58 58
59 The operators required to construct the SAT for imposing a Dirichlet 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 60 condition. `H_tuning` and `R_tuning` are used to specify the strength of the
61 penalty. 61 penalty.
62 62
63 See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref). 63 See also: [`sat`](@ref), [`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
64 """ 64 """
65 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.) 65 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
66 id = boundary(bc) 66 id = boundary(bc)
67 set = Δ.stencil_set 67 set = Δ.stencil_set
68 H⁻¹ = inverse_inner_product(g,set) 68 H⁻¹ = inverse_inner_product(g,set)
69 Hᵧ = inner_product(boundary_grid(g, id), set) 69 Hᵧ = inner_product(boundary_grid(g, id), set)
70 e = boundary_restriction(g, set, id) 70 e = boundary_restriction(g, set, id)
71 d = normal_derivative(g, set, id) 71 d = normal_derivative(g, set, id)
72 B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning) 72 B = positivity_decomposition(Δ, g, boundary(bc); H_tuning, R_tuning)
73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ 73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
74 return penalty_tensor, e 74 return penalty_tensor, e
75 end 75 end
76 76
77 """ 77 """
92 penalty_tensor = -H⁻¹∘e'∘Hᵧ 92 penalty_tensor = -H⁻¹∘e'∘Hᵧ
93 return penalty_tensor, d 93 return penalty_tensor, d
94 end 94 end
95 95
96 """ 96 """
97 positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) 97 positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
98 98
99 Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive 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 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 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 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 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 104 penalty and must be greater than 1. For details we refer to
105 https://doi.org/10.1016/j.jcp.2020.109294 105 <https://doi.org/10.1016/j.jcp.2020.109294>
106 """ 106 """
107 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) 107 function positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning)
108 @assert(H_tuning ≥ 1.) 108 @assert(H_tuning ≥ 1.)
109 @assert(R_tuning ≥ 1.) 109 @assert(R_tuning ≥ 1.)
110 Nτ_H, τ_R = positivity_limits(Δ,g,bc) 110 Nτ_H, τ_R = positivity_limits(Δ,g,b)
111 return H_tuning*Nτ_H + R_tuning*τ_R 111 return H_tuning*Nτ_H + R_tuning*τ_R
112 end 112 end
113 113
114 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then 114 function positivity_limits(Δ::Laplace, g::EquidistantGrid, b::BoundaryIdentifier)
115 # change bc::BoundaryCondition to id::BoundaryIdentifier
116 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
117 h = spacing(g) 115 h = spacing(g)
118 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1]) 116 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
119 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"]) 117 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
120 118
121 τ_H = 1/(h*θ_H) 119 τ_H = one(eltype(Δ))/(h*θ_H)
122 τ_R = 1/(h*θ_R) 120 τ_R = one(eltype(Δ))/(h*θ_R)
123 return τ_H, τ_R 121 return τ_H, τ_R
124 end 122 end
125 123
126 function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition) 124 function positivity_limits(Δ::Laplace, g::TensorGrid, b::BoundaryIdentifier)
127 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc) 125 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(b)], b)
128 return τ_H*ndims(g), τ_R 126 return τ_H*ndims(g), τ_R
129 end 127 end