Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1649:b02917bcd7d5 feature/grids/curvilinear

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 Wed, 26 Jun 2024 12:41:03 +0200
parents 1937be9502a7
children 707fc9761c2b f28c92ec843c 84aed3abab94
comparison
equal deleted inserted replaced
1645:64452a678e7a 1649:b02917bcd7d5
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 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set) 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, 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, bc::DirichletCondition; 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, bc::DirichletCondition; H_tuning, R_tuning)
108 @assert(H_tuning ≥ 1.)
109 @assert(R_tuning ≥ 1.)
110 Nτ_H, τ_R = positivity_limits(Δ,g,bc)
111 return H_tuning*Nτ_H + R_tuning*τ_R
112 end
113
114 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
115 # change bc::BoundaryCondition to id::BoundaryIdentifier
116 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
117 h = spacing(g)
118 θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
119 θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
120
121 τ_H = 1/(h*θ_H)
122 τ_R = 1/(h*θ_R)
123 return τ_H, τ_R
124 end
125
126 function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition)
127 τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc)
128 return τ_H*ndims(g), τ_R
129 end