diff src/SbpOperators/volumeops/laplace/laplace.jl @ 1598:19cdec9c21cb feature/boundary_conditions

Implement and test sat_tensors for Dirichlet and Neumann conditions
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sun, 26 May 2024 18:19:02 -0700
parents 8d60d045c2a2
children 37b05221beda
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun May 26 18:18:17 2024 -0700
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun May 26 18:19:02 2024 -0700
@@ -53,12 +53,31 @@
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
-# TODO: Add sat_tensor for Diirichlet condition
+"""
+sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+
+The operators required to construct the SAT for imposing a Dirichlet condition.
+`tuning` specifies the strength of the penalty. See
+
+See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
+"""
+function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+    id = bc.id
+    set  = Δ.stencil_set
+    H⁻¹ = inverse_inner_product(g,set)
+    Hᵧ = inner_product(boundary_grid(g, id), set)
+    e = boundary_restriction(g, set, id)
+    d = normal_derivative(g, set, id)
+    B = positivity_decomposition(Δ, g, bc, tuning)
+    sat_op = H⁻¹∘(d' - B*e')∘Hᵧ
+    return sat_op, e
+end
+BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition) = BoundaryConditions.sat_tensors(Δ, g, bc, (1.,1.))
 
 """
 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
 
-The operators required to construct the SAT for imposing Neumann condition
+The operators required to construct the SAT for imposing a Neumann condition
 
 
 See also: [`sat`,`NeumannCondition`](@ref).
@@ -71,6 +90,19 @@
     e = boundary_restriction(g, set, id)
     d = normal_derivative(g, set, id)
 
-    sat_op = H⁻¹∘e'∘Hᵧ
+    sat_op = -H⁻¹∘e'∘Hᵧ
     return sat_op, d
 end
+
+function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+    pos_prop = positivity_properties(Δ)
+    h = spacing(orthogonal_grid(g, bc.id))
+    θ_H = pos_prop.theta_H
+    τ_H = tuning[1]*ndims(g)/(h*θ_H)
+    θ_R = pos_prop.theta_R
+    τ_R = tuning[2]/(h*θ_R)
+    B = τ_H + τ_R
+    return B
+end
+
+positivity_properties(Δ::Laplace) = parse_named_tuple(Δ.stencil_set["Positivity"]["D2"])
\ No newline at end of file