Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1606:93b86625fcfd feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
REVIEW: Suggest split of tuning tuple. Please help with names!
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 08 Jun 2024 23:47:23 +0200
parents fca4a01d60c9
children 7216448d0c5a
comparison
equal deleted inserted replaced
1605:1388149b54ad 1606:93b86625fcfd
60 The operators required to construct the SAT for imposing a Dirichlet condition. 60 The operators required to construct the SAT for imposing a Dirichlet condition.
61 `tuning` specifies the strength of the penalty. See 61 `tuning` specifies the strength of the penalty. See
62 62
63 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref). 63 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
64 """ 64 """
65 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning = (1., 1.)) 65 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
66 id = boundary(bc) 66 id = boundary(bc)
67 set = Δ.stencil_set 67 set = Δ.stencil_set
68 H⁻¹ = inverse_inner_product(g,set) 68 H⁻¹ = inverse_inner_product(g,set)
69 Hᵧ = inner_product(boundary_grid(g, id), set) 69 Hᵧ = inner_product(boundary_grid(g, id), set)
70 e = boundary_restriction(g, set, id) 70 e = boundary_restriction(g, set, id)
71 d = normal_derivative(g, set, id) 71 d = normal_derivative(g, set, id)
72 B = positivity_decomposition(Δ, g, bc, tuning) 72 B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning)
73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ 73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
74 return penalty_tensor, e 74 return penalty_tensor, e
75 end 75 end
76 76
77 """ 77 """
95 end 95 end
96 96
97 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then 97 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
98 # change bc::BoundaryCondition to id::BoundaryIdentifier 98 # change bc::BoundaryCondition to id::BoundaryIdentifier
99 99
100 function positivity_decomposition(Δ::Laplace, g::EquidistantGrid, bc::BoundaryCondition, tuning) 100 function positivity_decomposition(Δ::Laplace, g::EquidistantGrid, bc::BoundaryCondition; H_tuning, R_tuning)
101 pos_prop = positivity_properties(Δ) 101 pos_prop = positivity_properties(Δ)
102 h = spacing(g) 102 h = spacing(g)
103 θ_H = pos_prop.theta_H 103 θ_H = pos_prop.theta_H
104 τ_H = tuning[1]*ndims(g)/(h*θ_H) 104 τ_H = H_tuning*ndims(g)/(h*θ_H)
105 θ_R = pos_prop.theta_R 105 θ_R = pos_prop.theta_R
106 τ_R = tuning[2]/(h*θ_R) 106 τ_R = R_tuning/(h*θ_R)
107 B = τ_H + τ_R 107 B = τ_H + τ_R
108 return B 108 return B
109 end 109 end
110 110
111 function positivity_decomposition(Δ::Laplace, g::TensorGrid, bc::BoundaryCondition, tuning) 111 function positivity_decomposition(Δ::Laplace, g::TensorGrid, bc::BoundaryCondition; H_tuning, R_tuning)
112 pos_prop = positivity_properties(Δ) 112 pos_prop = positivity_properties(Δ)
113 h = spacing(g.grids[grid_id(boundary(bc))]) # grid spacing of the 1D grid normal to the boundary 113 h = spacing(g.grids[grid_id(boundary(bc))]) # grid spacing of the 1D grid normal to the boundary
114 θ_H = pos_prop.theta_H 114 θ_H = pos_prop.theta_H
115 τ_H = tuning[1]*ndims(g)/(h*θ_H) 115 τ_H = H_tuning*ndims(g)/(h*θ_H)
116 θ_R = pos_prop.theta_R 116 θ_R = pos_prop.theta_R
117 τ_R = tuning[2]/(h*θ_R) 117 τ_R = R_tuning/(h*θ_R)
118 B = τ_H + τ_R 118 B = τ_H + τ_R
119 return B 119 return B
120 end 120 end
121 121
122 function positivity_properties(Δ::Laplace) 122 function positivity_properties(Δ::Laplace)