Mercurial > repos > public > sbplib_julia

changeset 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 330c39505a94
children 37b05221beda
files TODO.md src/SbpOperators/volumeops/laplace/laplace.jl test/SbpOperators/volumeops/laplace/laplace_test.jl
diffstat 3 files changed, 87 insertions(+), 40 deletions(-) [+]
line wrap: on
line diff
--- a/TODO.md	Sun May 26 18:18:17 2024 -0700
+++ b/TODO.md	Sun May 26 18:19:02 2024 -0700
@@ -21,6 +21,10 @@
        See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing)
  - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe
  - [ ] Move export statements to top of each module
+ - [ ] Implement apply_transpose for
+      - [ ] ElementwiseTensorOperation
+      - [ ] VolumeOperator
+      - [ ] Laplace
 
 
  - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative. Kan vi använda promote istället?
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun May 26 18:18:17 2024 -0700
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun May 26 18:19:02 2024 -0700
@@ -53,12 +53,31 @@
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
-# TODO: Add sat_tensor for Diirichlet condition
+"""
+sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+
+The operators required to construct the SAT for imposing a Dirichlet condition.
+`tuning` specifies the strength of the penalty. See
+
+See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
+"""
+function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+    id = bc.id
+    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)
+    sat_op = H⁻¹∘(d' - B*e')∘Hᵧ
+    return sat_op, e
+end
+BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition) = BoundaryConditions.sat_tensors(Δ, g, bc, (1.,1.))
 
 """
 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
 
-The operators required to construct the SAT for imposing Neumann condition
+The operators required to construct the SAT for imposing a Neumann condition
 
 
 See also: [`sat`,`NeumannCondition`](@ref).
@@ -71,6 +90,19 @@
     e = boundary_restriction(g, set, id)
     d = normal_derivative(g, set, id)
 
-    sat_op = H⁻¹∘e'∘Hᵧ
+    sat_op = -H⁻¹∘e'∘Hᵧ
     return sat_op, d
 end
+
+function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+    pos_prop = positivity_properties(Δ)
+    h = spacing(orthogonal_grid(g, bc.id))
+    θ_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
+end
+
+positivity_properties(Δ::Laplace) = parse_named_tuple(Δ.stencil_set["Positivity"]["D2"])
\ No newline at end of file
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sun May 26 18:18:17 2024 -0700
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sun May 26 18:19:02 2024 -0700
@@ -91,47 +91,58 @@
 
 @testset "sat_tensors" begin
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
-    stencil_set = read_stencil_set(operator_path; order=4)
-    g = equidistant_grid((101,102), (-1.,-1.), (1.,1.))
-    W,E,S,N = boundary_identifiers(g)
+    orders = (2,4)
+    tols = (5e-2,5e-4)
+    sz = (201,401)
+    g = equidistant_grid((0.,0.), (1.,1.), sz...)
     
-    u = eval_on(g, (x,y) -> sin(x+y))
-    uWx = eval_on(boundary_grid(g,W), (x,y) -> -cos(x+y))
-    uEx = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y))
-    uSy = eval_on(boundary_grid(g,S), (x,y) -> -cos(x+y))
-    uNy = eval_on(boundary_grid(g,N), (x,y) -> cos(x+y))
-    
-    
-    v = eval_on(g, (x,y) -> cos(x+y))
-    vW = eval_on(boundary_grid(g,W), (x,y) -> cos(x+y))
-    vE = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y))
-    vS = eval_on(boundary_grid(g,S), (x,y) -> cos(x+y))
-    vN = eval_on(boundary_grid(g,N), (x,y) -> cos(x+y))
+    # Verify implementation of sat_tesnors by testing accuracy and symmetry (TODO) 
+    # of the operator D = Δ + SAT, where SAT is the tensor composition of the 
+    # operators from sat_tensor. Note that SAT*u should approximate 0 for the 
+    # conditions chosen.
 
+    @testset "Dirichlet" begin
+        for (o, tol) ∈ zip(orders,tols)
+            stencil_set = read_stencil_set(operator_path; order=o)
+            Δ = Laplace(g, stencil_set)
+            H = inner_product(g, stencil_set)
+            u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
+            Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
+            D = Δ 
+            for id ∈ boundary_identifiers(g)
+                D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
+            end
+            e = D*u .- Δu
+            # Accuracy
+            @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+            # Symmetry
+            # TODO: # Consider generating the matrices to H and D and test D'H == H'D
+            r = randn(size(u))
+            @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+        end
+    end
 
     @testset "Neumann" begin
-        Δ = Laplace(g, stencil_set)
-        H = inner_product(g, stencil_set)
-        HW = inner_product(boundary_grid(g,W), stencil_set)
-        HE = inner_product(boundary_grid(g,E), stencil_set)
-        HS = inner_product(boundary_grid(g,S), stencil_set)
-        HN = inner_product(boundary_grid(g,N), stencil_set)
-        
-        ncW = NeumannCondition(0., W)
-        ncE = NeumannCondition(0., E)
-        ncS = NeumannCondition(0., S)
-        ncN = NeumannCondition(0., N)
-        
-        SATW = foldl(∘,sat_tensors(Δ, g, ncW))
-        SATE = foldl(∘,sat_tensors(Δ, g, ncE))
-        SATS = foldl(∘,sat_tensors(Δ, g, ncS))
-        SATN = foldl(∘,sat_tensors(Δ, g, ncN))
-        
-
-        @test sum((H*SATW*u).*v) ≈ sum((HW*uWx).*vW) rtol = 1e-6
-        @test sum((H*SATE*u).*v) ≈ sum((HE*uEx).*vE) rtol = 1e-6
-        @test sum((H*SATS*u).*v) ≈ sum((HS*uSy).*vS) rtol = 1e-6
-        @test sum((H*SATN*u).*v) ≈ sum((HN*uNy).*vN) rtol = 1e-6
+        @testset "Dirichlet" begin
+            for (o, tol) ∈ zip(orders,tols)
+                stencil_set = read_stencil_set(operator_path; order=o)
+                Δ = Laplace(g, stencil_set)
+                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 = Δ 
+                for id ∈ boundary_identifiers(g)
+                    op = op + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
+                end
+                e = op*u .- Δu
+                # Accuracy
+                @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+                # Symmetry
+                # TODO: # Consider generating the matrices to H and D and test D'H == H'D
+                r = randn(size(u))
+                @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+            end
+        end
     end
 end