Mercurial > repos > public > sbplib_julia

changeset 796:f682e4fe3cef operator_storage_array_of_table

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix tests for inner products
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 25 Jul 2021 14:39:41 +0200
parents 4ff2dad2dd7f
children 92fafe5980dd
files src/SbpOperators/operators/standard_diagonal.toml src/SbpOperators/volumeops/inner_products/inner_product.jl test/SbpOperators/volumeops/inner_products/inner_product_test.jl
diffstat 3 files changed, 39 insertions(+), 26 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/operators/standard_diagonal.toml	Sun Jul 25 14:38:04 2021 +0200
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Sun Jul 25 14:39:41 2021 +0200
@@ -8,7 +8,7 @@
 
 order = 2
 
-H.inner = ["1"]
+H.inner = "1"
 H.closure = ["1/2"]
 
 D1.inner_stencil = ["-1/2", "0", "1/2"]
@@ -27,7 +27,7 @@
 [[stencil_set]]
 
 order = 4
-H.inner = ["1"]
+H.inner = "1"
 H.closure = ["17/48", "59/48", "43/48", "49/48"]
 
 D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Sun Jul 25 14:38:04 2021 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Sun Jul 25 14:39:41 2021 +0200
@@ -17,21 +17,22 @@
 `SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
 `H` is a 0-dimensional `IdentityMapping`.
 """
-function inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
+function inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
     Hs = ()
 
     for i ∈ 1:dimension(grid)
-        Hs = (Hs..., inner_product(restrict(grid, i), closure_stencils, inner_stencil))
+        Hs = (Hs..., inner_product(restrict(grid, i), interior_weight, closure_weights))
     end
 
     return foldl(⊗, Hs)
 end
 export inner_product
 
-function inner_product(grid::EquidistantGrid{1}, closure_stencils, inner_stencil)
-    h = spacing(grid)
-    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
+function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+    h = spacing(grid)[1]
+
+    H = SbpOperators.ConstantInteriorScalingOperator(grid, h*interior_weight, h.*closure_weights)
     return H
 end
 
-inner_product(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()
+inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sun Jul 25 14:38:04 2021 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sun Jul 25 14:39:41 2021 +0200
@@ -14,35 +14,39 @@
     g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
     integral(H,v) = sum(H*v)
     @testset "inner_product" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+        quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
         @testset "0D" begin
-            H = inner_product(EquidistantGrid{Float64}(), op.quadratureClosure, CenteredStencil(1.))
+            H = inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
             @test H == IdentityMapping{Float64}()
             @test H isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
-            H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
-            @test H == inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+            @test H == inner_product(g_1D, quadrature_interior, quadrature_closure)
             @test H isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
-            H_x = inner_product(restrict( g_2D,1),op.quadratureClosure, CenteredStencil(1.))
-            H_y = inner_product(restrict( g_2D,2),op.quadratureClosure, CenteredStencil(1.))
+            H = inner_product(g_2D, quadrature_interior, quadrature_closure)
+            H_x = inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure)
+            H_y = inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure)
             @test H == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
     end
 
     @testset "Sizes" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+        quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
         @testset "1D" begin
-            H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+            H = inner_product(g_1D, quadrature_interior, quadrature_closure)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
+            H = inner_product(g_2D, quadrature_interior, quadrature_closure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -58,8 +62,10 @@
             u = evalOn(g_1D,x->sin(x))
 
             @testset "2nd order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -67,8 +73,10 @@
             end
 
             @testset "4th order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -81,14 +89,18 @@
             v = b*ones(Float64, size(g_2D))
             u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             @testset "2nd order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-4
             end
             @testset "4th order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-8
             end