Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1651:707fc9761c2b 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, 26 Jun 2024 12:47:26 +0200
parents 4f6f5e5daa35 1937be9502a7
children 29b96fc75bee
comparison
equal deleted inserted replaced
1648:b7dcd3dd3181 1651:707fc9761c2b
49 for d = 2:ndims(g) 49 for d = 2:ndims(g)
50 Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d) 50 Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d)
51 end 51 end
52 return Δ 52 return Δ
53 end 53 end
54
54 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set) 55 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
55
56 56
57 function laplace(grid::MappedGrid, stencil_set) 57 function laplace(grid::MappedGrid, stencil_set)
58 J = jacobian_determinant(grid) 58 J = jacobian_determinant(grid)
59 J⁻¹ = DiagonalTensor(map(inv, J)) 59 J⁻¹ = DiagonalTensor(map(inv, J))
60 60
72 Dⱼ = first_derivative(lg, stencil_set, j) 72 Dⱼ = first_derivative(lg, stencil_set, j)
73 J⁻¹∘Dᵢ∘DiagonalTensor(Jgⁱʲ)∘Dⱼ 73 J⁻¹∘Dᵢ∘DiagonalTensor(Jgⁱʲ)∘Dⱼ
74 end 74 end
75 end 75 end
76 end 76 end
77
78
79 """
80 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
81
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
84 penalty.
85
86 See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
87 """
88 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
89 id = boundary(bc)
90 set = Δ.stencil_set
91 H⁻¹ = inverse_inner_product(g,set)
92 Hᵧ = inner_product(boundary_grid(g, id), set)
93 e = boundary_restriction(g, set, id)
94 d = normal_derivative(g, set, id)
95 B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning)
96 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
97 return penalty_tensor, e
98 end
99
100 """
101 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
102
103 The operators required to construct the SAT for imposing a Neumann condition.
104
105 See also: [`sat`](@ref), [`NeumannCondition`](@ref).
106 """
107 function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
108 id = boundary(bc)
109 set = Δ.stencil_set
110 H⁻¹ = inverse_inner_product(g,set)
111 Hᵧ = inner_product(boundary_grid(g, id), set)
112 e = boundary_restriction(g, set, id)
113 d = normal_derivative(g, set, id)
114
115 penalty_tensor = -H⁻¹∘e'∘Hᵧ
116 return penalty_tensor, d
117 end
118
119 """
120 positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
121
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
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
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
128 https://doi.org/10.1016/j.jcp.2020.109294
129 """
130 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
131 @assert(H_tuning ≥ 1.)
132 @assert(R_tuning ≥ 1.)
133 Nτ_H, τ_R = positivity_limits(Δ,g,bc)
134 return H_tuning*Nτ_H + R_tuning*τ_R
135 end
136
137 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
138 # change bc::BoundaryCondition to id::BoundaryIdentifier
139 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
140 h = spacing(g)
141 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
142 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
143
144 τ_H = 1/(h*θ_H)
145 τ_R = 1/(h*θ_R)
146 return τ_H, τ_R
147 end
148
149 function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition)
150 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc)
151 return τ_H*ndims(g), τ_R
152 end