Mercurial > repos > public > sbplib_julia

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
1597:330c39505a94 1598:19cdec9c21cb
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 55
56 # TODO: Add sat_tensor for Diirichlet condition 56 """
57 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
58
59 The operators required to construct the SAT for imposing a Dirichlet condition.
60 `tuning` specifies the strength of the penalty. See
61
62 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
63 """
64 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
65 id = bc.id
66 set = Δ.stencil_set
67 H⁻¹ = inverse_inner_product(g,set)
68 Hᵧ = inner_product(boundary_grid(g, id), set)
69 e = boundary_restriction(g, set, id)
70 d = normal_derivative(g, set, id)
71 B = positivity_decomposition(Δ, g, bc, tuning)
72 sat_op = H⁻¹∘(d' - B*e')∘Hᵧ
73 return sat_op, e
74 end
75 BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition) = BoundaryConditions.sat_tensors(Δ, g, bc, (1.,1.))
57 76
58 """ 77 """
59 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) 78 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
60 79
61 The operators required to construct the SAT for imposing Neumann condition 80 The operators required to construct the SAT for imposing a Neumann condition
62 81
63 82
64 See also: [`sat`,`NeumannCondition`](@ref). 83 See also: [`sat`,`NeumannCondition`](@ref).
65 """ 84 """
66 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) 85 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
69 H⁻¹ = inverse_inner_product(g,set) 88 H⁻¹ = inverse_inner_product(g,set)
70 Hᵧ = inner_product(boundary_grid(g, id), set) 89 Hᵧ = inner_product(boundary_grid(g, id), set)
71 e = boundary_restriction(g, set, id) 90 e = boundary_restriction(g, set, id)
72 d = normal_derivative(g, set, id) 91 d = normal_derivative(g, set, id)
73 92
74 sat_op = H⁻¹∘e'∘Hᵧ 93 sat_op = -H⁻¹∘e'∘Hᵧ
75 return sat_op, d 94 return sat_op, d
76 end 95 end
96
97 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
98 pos_prop = positivity_properties(Δ)
99 h = spacing(orthogonal_grid(g, bc.id))
100 θ_H = pos_prop.theta_H
101 τ_H = tuning[1]*ndims(g)/(h*θ_H)
102 θ_R = pos_prop.theta_R
103 τ_R = tuning[2]/(h*θ_R)
104 B = τ_H + τ_R
105 return B
106 end
107
108 positivity_properties(Δ::Laplace) = parse_named_tuple(Δ.stencil_set["Positivity"]["D2"])