Mercurial > repos > public > sbplib_julia

changeset 1614:13063028d604 feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge heads
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sun, 09 Jun 2024 16:53:03 -0700
parents 15488c889a50 (diff) df1856b0e2f0 (current diff)
children b74e1a21265f
files
diffstat 6 files changed, 22 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Sun Jun 09 23:03:04 2024 +0200
+++ b/Notes.md	Sun Jun 09 16:53:03 2024 -0700
@@ -1,23 +1,5 @@
 # Notes
 
-## Boundary Conditions and SATs
-
-Types for boundary conditions:
-
- * abstract type `BoundaryData`
- * abstract type `BoundaryCondition{T<:BoundaryData}`
- * concrete types `ConstantBoundaryData <: BoundaryData` and similar
- * concrete types `NeumannCondition{BD<:BoundaryData} <: BoundaryCondition{BD}` and similar
-The concrete `BoundaryData` subtypes are "thin types" wrapping the boundary data, and are used to indicate how the boundary data should be used in e.g. sat routines. The concrete `BoundaryCondition{BD}` subtypes are used for assembling the tensors used to construct e.g. a SAT.
-
-SAT methods:
-There are multiple options for what the SAT methods could return.
-* (Current) a function which returns a `LazyTensorApplication`, e.g. `f = sat(grid,op,bc)`. The the resulting `LazyTensorApplication` can then be added to scheme i.e. `scheme = op*u + f(u)`.  Depdending on the type of data in the BC, e.g. time-depdendent etc one can return f(u,t).
-* `LazyTensor`s `closure, penalty = sat(grid,op,bc)` like in the matlab version. Probably the most general one. Up to the user to make use of the returned `LazyTensor`s. One can for example then easily include the closures to the operator and have eg. `D = (op + closure)*u`.
-* A `LazyTensor` for closure, and a `LazyArray` for `penalty*data`. Mix of the above.
-* Same as first but of  the  form sat = `sat_op*(L*u-g)`. This is how one typically would write the SAT in the litterature. The function `sat_tensors` would return `sat_op` and `L`. Need to get compositions working before we can implement this approach.
-
-
 ## Reading operators
 
 Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and
--- a/src/SbpOperators/SbpOperators.jl	Sun Jun 09 23:03:04 2024 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Sun Jun 09 16:53:03 2024 -0700
@@ -8,11 +8,9 @@
 export parse_nested_stencil
 export parse_scalar
 export parse_tuple
-export parse_named_tuple
 export sbp_operators_path
 
 # Operators
-export boundary_quadrature
 export boundary_restriction
 export inner_product
 export inverse_inner_product
--- a/src/SbpOperators/stencil_set.jl	Sun Jun 09 23:03:04 2024 +0200
+++ b/src/SbpOperators/stencil_set.jl	Sun Jun 09 16:53:03 2024 -0700
@@ -166,17 +166,6 @@
     return Tuple(parse_scalar.(parsed_toml))
 end
 
-"""
-    parse_named_tuple(parsed_toml)
-
-Parse the keys (names) and values (scalars) into a named tuple of rationals.
-
-See also [`parse_scalar`](@ref).
-"""
-function parse_named_tuple(parsed_toml)
-    NamedTuple(Symbol(key) => parse_scalar(val) for (key, val) in parsed_toml)
-end
-
 
 """
     parse_rational(parsed_toml)
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 09 23:03:04 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 09 16:53:03 2024 -0700
@@ -62,14 +62,14 @@
 
 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
 """
-function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning = (1., 1.))
+function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
     id = boundary(bc)
     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)
+    B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning)
     penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
     return penalty_tensor, e
 end
@@ -94,33 +94,26 @@
     return penalty_tensor, d
 end
 
-# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
-# change bc::BoundaryCondition to id::BoundaryIdentifier
 
-function positivity_decomposition(Δ::Laplace, g::EquidistantGrid, bc::BoundaryCondition, tuning)
-    pos_prop = positivity_properties(Δ)
-    h = spacing(g)
-    θ_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
+function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+    Nτ_H, τ_R = positivity_limits(Δ,g,bc)
+    return H_tuning*Nτ_H + R_tuning*τ_R
 end
 
-function positivity_decomposition(Δ::Laplace, g::TensorGrid, bc::BoundaryCondition, tuning)
-    pos_prop = positivity_properties(Δ)
-    h = spacing(g.grids[grid_id(boundary(bc))]) # grid spacing of the 1D grid normal to the boundary
-    θ_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
+
+# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
+# change bc::BoundaryCondition to id::BoundaryIdentifier
+function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
+    h = spacing(g)
+    θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
+    θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
+
+    τ_H = 1/(h*θ_H)
+    τ_R = 1/(h*θ_R)
+    return τ_H, τ_R
 end
 
-function positivity_properties(Δ::Laplace)
-    D2_pos_prop = parse_named_tuple(Δ.stencil_set["D2"]["positivity"])
-    H_closure = parse_tuple(Δ.stencil_set["H"]["closure"])
-    return merge(D2_pos_prop, (theta_H = H_closure[1],))
+function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition)
+    τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc)
+    return τ_H*ndims(g), τ_R
 end
--- a/test/SbpOperators/stencil_set_test.jl	Sun Jun 09 23:03:04 2024 +0200
+++ b/test/SbpOperators/stencil_set_test.jl	Sun Jun 09 16:53:03 2024 -0700
@@ -152,14 +152,6 @@
         @test_throws ArgumentError parse_tuple(toml["e3"])
         @test_throws ArgumentError parse_tuple(toml["e4"])
     end
-
-    @testset "parse_named_tuple" begin
-        toml = TOML.parse("""
-            NamedTuple = {a = "17/48", b = "0.2505765857"}
-        """)
-        @test parse_named_tuple(toml["NamedTuple"]) == (b = Rational(0.2505765857), a = 17//48)
-    end
-
 end
 
 @testset "parse_rational" begin
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sun Jun 09 23:03:04 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sun Jun 09 16:53:03 2024 -0700
@@ -129,11 +129,11 @@
                 H = inner_product(g, stencil_set)
                 u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
                 Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
-                op = Δ 
+                D = Δ 
                 for id ∈ boundary_identifiers(g)
-                    op = op + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
+                    D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
                 end
-                e = op*u .- Δu
+                e = D*u .- Δu
                 # Accuracy
                 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
                 # Symmetry