Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1602:3e7438e2a033 feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Address review comments (1 left to be discussed)
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sat, 01 Jun 2024 17:39:54 -0700
parents b2496b001297
children fca4a01d60c9
comparison
equal deleted inserted replaced
1601:fad18896d20a 1602:3e7438e2a033
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 # REVIEW: I think the handling of tuning parameters below should be through kwargs instead.
57 56
58 """ 57 """
59 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning) 58 sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning)
60 59
61 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.
62 `tuning` specifies the strength of the penalty. See 61 `tuning` specifies the strength of the penalty. See
63 62
64 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref). 63 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
65 """ 64 """
66 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning) 65 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning = (1., 1.))
67 id = bc.id 66 id = boundary(bc)
68 set = Δ.stencil_set 67 set = Δ.stencil_set
69 H⁻¹ = inverse_inner_product(g,set) 68 H⁻¹ = inverse_inner_product(g,set)
70 Hᵧ = inner_product(boundary_grid(g, id), set) 69 Hᵧ = inner_product(boundary_grid(g, id), set)
71 e = boundary_restriction(g, set, id) 70 e = boundary_restriction(g, set, id)
72 d = normal_derivative(g, set, id) 71 d = normal_derivative(g, set, id)
73 B = positivity_decomposition(Δ, g, bc, tuning) 72 B = positivity_decomposition(Δ, g, bc, tuning)
74 sat_op = H⁻¹∘(d' - B*e')∘Hᵧ 73 penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
75 return sat_op, e 74 return penalty_tensor, e
76 end 75 end
77 BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition) = BoundaryConditions.sat_tensors(Δ, g, bc, (1.,1.)) # REVIEW: Should be possible to replace this with argument default values.
78 76
79 """ 77 """
80 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) 78 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
81 79
82 The operators required to construct the SAT for imposing a Neumann condition 80 The operators required to construct the SAT for imposing a Neumann condition
83 81
84 82
85 See also: [`sat`,`NeumannCondition`](@ref). 83 See also: [`sat`,`NeumannCondition`](@ref).
86 """ 84 """
87 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) 85 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
88 id = bc.id 86 id = boundary(bc)
89 set = Δ.stencil_set 87 set = Δ.stencil_set
90 H⁻¹ = inverse_inner_product(g,set) 88 H⁻¹ = inverse_inner_product(g,set)
91 Hᵧ = inner_product(boundary_grid(g, id), set) 89 Hᵧ = inner_product(boundary_grid(g, id), set)
92 e = boundary_restriction(g, set, id) 90 e = boundary_restriction(g, set, id)
93 d = normal_derivative(g, set, id) 91 d = normal_derivative(g, set, id)
94 92
95 sat_op = -H⁻¹∘e'∘Hᵧ 93 penalty_tensor = -H⁻¹∘e'∘Hᵧ
96 return sat_op, d 94 return penalty_tensor, d
97 end 95 end
98 96
99 # REVIEW: This function assumes a TensorGrid right? In that case there should probably be a type annotation to get clearer error messages. 97 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
100 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning) 98 # change bc::BoundaryCondition to id::BoundaryIdentifier
99
100 function positivity_decomposition(Δ::Laplace, g::EquidistantGrid, bc::BoundaryCondition, tuning)
101 pos_prop = positivity_properties(Δ) 101 pos_prop = positivity_properties(Δ)
102 h = spacing(orthogonal_grid(g, bc.id)) 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 = tuning[1]*ndims(g)/(h*θ_H)
105 θ_R = pos_prop.theta_R 105 θ_R = pos_prop.theta_R
106 τ_R = tuning[2]/(h*θ_R) 106 τ_R = tuning[2]/(h*θ_R)
107 B = τ_H + τ_R 107 B = τ_H + τ_R
108 return B 108 return B
109 end 109 end
110 110
111 positivity_properties(Δ::Laplace) = parse_named_tuple(Δ.stencil_set["Positivity"]["D2"]) # REVIEW: Can this function extract theta_H from the inner product instead of storing it twice in the TOML? 111 function positivity_decomposition(Δ::Laplace, g::TensorGrid, bc::BoundaryCondition, tuning)
112 pos_prop = positivity_properties(Δ)
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
115 τ_H = tuning[1]*ndims(g)/(h*θ_H)
116 θ_R = pos_prop.theta_R
117 τ_R = tuning[2]/(h*θ_R)
118 B = τ_H + τ_R
119 return B
120 end
121
122 function positivity_properties(Δ::Laplace)
123 D2_pos_prop = parse_named_tuple(Δ.stencil_set["D2"]["positivity"])
124 H_closure = parse_tuple(Δ.stencil_set["H"]["closure"])
125 return merge(D2_pos_prop, (theta_H = H_closure[1],))
126 end