Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1616:e41eddc640f3 feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Docs
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sun, 09 Jun 2024 17:51:57 -0700
parents 8315c456e3b4
children 1937be9502a7
comparison
equal deleted inserted replaced
1615:b74e1a21265f 1616:e41eddc640f3
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
57 """ 56 """
58 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning) 57 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
59 58
60 The operators required to construct the SAT for imposing a Dirichlet condition. 59 The operators required to construct the SAT for imposing a Dirichlet condition.
61 `tuning` specifies the strength of the penalty. See 60 `H_tuning` and `R_tuning` are used to specify the strength of the penalty.
62 61
63 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref). 62 See also: [`sat`](@ref),[`DirichletCondition`](@ref),[`positivity_decomposition`](@ref).
64 """ 63 """
65 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.) 64 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
66 id = boundary(bc) 65 id = boundary(bc)
67 set = Δ.stencil_set 66 set = Δ.stencil_set
68 H⁻¹ = inverse_inner_product(g,set) 67 H⁻¹ = inverse_inner_product(g,set)
73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ 72 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
74 return penalty_tensor, e 73 return penalty_tensor, e
75 end 74 end
76 75
77 """ 76 """
78 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) 77 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
79 78
80 The operators required to construct the SAT for imposing a Neumann condition 79 The operators required to construct the SAT for imposing a Neumann condition
81 80
82 81 See also: [`sat`](@ref),[`NeumannCondition`](@ref).
83 See also: [`sat`,`NeumannCondition`](@ref).
84 """ 82 """
85 function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) 83 function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
86 id = boundary(bc) 84 id = boundary(bc)
87 set = Δ.stencil_set 85 set = Δ.stencil_set
88 H⁻¹ = inverse_inner_product(g,set) 86 H⁻¹ = inverse_inner_product(g,set)
92 90
93 penalty_tensor = -H⁻¹∘e'∘Hᵧ 91 penalty_tensor = -H⁻¹∘e'∘Hᵧ
94 return penalty_tensor, d 92 return penalty_tensor, d
95 end 93 end
96 94
95 """
96 positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
97 97
98 Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive definite with respect to
99 the boundary quadrature. Here `d` is the normal derivative and `e` is the boundary restriction operator.
100 `B` can then be used to form a symmetric and energy stable penalty for a Dirichlet condition.
101 The parameters `H_tuning` and `R_tuning` are used to specify the strength of the penalty and
102 must be greater than 1. For details we refer to https://doi.org/10.1016/j.jcp.2020.109294
103 """
98 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) 104 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
105 @assert(H_tuning ≥ 1.)
106 @assert(R_tuning ≥ 1.)
99 Nτ_H, τ_R = positivity_limits(Δ,g,bc) 107 Nτ_H, τ_R = positivity_limits(Δ,g,bc)
100 return H_tuning*Nτ_H + R_tuning*τ_R 108 return H_tuning*Nτ_H + R_tuning*τ_R
101 end 109 end
102
103 110
104 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then 111 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
105 # change bc::BoundaryCondition to id::BoundaryIdentifier 112 # change bc::BoundaryCondition to id::BoundaryIdentifier
106 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition) 113 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
107 h = spacing(g) 114 h = spacing(g)