Mercurial > repos > public > sbplib_julia

changeset 800:f91495f23604 operator_storage_array_of_table

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 25 Jul 2021 15:39:29 +0200
parents 24df68453890 (diff) 26bf5b2b3e32 (current diff)
children 04e549669b10
files src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl
diffstat 9 files changed, 199 insertions(+), 94 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl	Thu Jul 22 14:26:34 2021 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -8,6 +8,7 @@
 include("d2.jl")
 include("readoperator.jl")
 include("volumeops/volume_operator.jl")
+include("volumeops/constant_interior_scaling_operator.jl")
 include("volumeops/derivatives/second_derivative.jl")
 include("volumeops/laplace/laplace.jl")
 include("volumeops/inner_products/inner_product.jl")
--- a/src/SbpOperators/operators/standard_diagonal.toml	Thu Jul 22 14:26:34 2021 +0200
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Sun Jul 25 15:39:29 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"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -0,0 +1,48 @@
+"""
+    ConstantInteriorScalingOperator{T,N} <: TensorMapping{T,1,1}
+
+A one-dimensional operator scaling a vector. The first and last `N` points are
+scaled with individual weights while all interior points are scaled the same.
+"""
+struct ConstantInteriorScalingOperator{T,N} <: TensorMapping{T,1,1}
+    interior_weight::T
+    closure_weights::NTuple{N,T}
+    size::Int
+
+    function ConstantInteriorScalingOperator(interior_weight::T, closure_weights::NTuple{N,T}, size::Int) where {T,N}
+        if size < 2length(closure_weights)
+            throw(DomainError(size, "size must be larger that two times the closure size."))
+        end
+
+        return new{T,N}(interior_weight, closure_weights, size)
+    end
+end
+
+function ConstantInteriorScalingOperator(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+    return ConstantInteriorScalingOperator(interior_weight, Tuple(closure_weights), size(grid)[1])
+end
+
+closure_size(::ConstantInteriorScalingOperator{T,N}) where {T,N} = N
+
+LazyTensors.range_size(op::ConstantInteriorScalingOperator) = (op.size,)
+LazyTensors.domain_size(op::ConstantInteriorScalingOperator) = (op.size,)
+
+# TBD: @inbounds in apply methods?
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+    return op.closure_weights[Int(i)]*v[Int(i)]
+end
+
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+    return op.interior_weight*v[Int(i)]
+end
+
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T
+    return op.closure_weights[op.size[1]-Int(i)+1]*v[Int(i)]
+end
+
+function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T
+    r = getregion(i, closure_size(op), op.size[1])
+    return LazyTensors.apply(op, v, Index(i, r))
+end
+
+LazyTensors.apply_transpose(op::ConstantInteriorScalingOperator, v, i) = apply(op, v, i)
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Thu Jul 22 14:26:34 2021 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -1,37 +1,34 @@
-# TODO:refactor to take a tuple instead. Convert the tuple to stencil for now. Could probably be refactored using a diagonal operator later.
-# How would a block-ip be built? A method on inner_product taking a stencil collection for the closure which then returns a different type of tensormapping
+"""
+    inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
 
-"""
-    inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
-
-Creates the discrete inner product operator `H` as a `TensorMapping` on an equidistant
-grid, defined as `(u,v)  = u'Hv` for grid functions `u,v`.
+Creates the discrete inner product operator `H` as a `TensorMapping` on an
+equidistant grid, defined as `(u,v)  = u'Hv` for grid functions `u,v`.
 
-`inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)` creates
-`H` on `grid` the using a set of stencils `closure_stencils` for the points in
-the closure regions and the stencil and `inner_stencil` in the interior.
+`inner_product` creates `H` on `grid` using the `interior_weight` for the
+interior points and the `closure_weights` for the points close to the
+boundary.
 
-On a 1-dimensional `grid`, `H` is a `VolumeOperator`. On a N-dimensional
-`grid`, `H` is the outer product of the 1-dimensional inner product operators in
-each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
-`H` is a 0-dimensional `IdentityMapping`.
+On a 1-dimensional grid, `H` is a `ConstantInteriorScalingOperator`. On a
+N-dimensional grid, `H` is the outer product of the 1-dimensional inner
+product operators for each coordinate direction. Also see the documentation of
+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/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Thu Jul 22 14:26:34 2021 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -1,37 +1,30 @@
 """
-    inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil)
+    inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
 
-Creates the inverse inner product operator `H⁻¹` as a `TensorMapping` on an
-equidistant grid. `H⁻¹` is defined implicitly by `H⁻¹∘H = I`, where
-`H` is the corresponding inner product operator and `I` is the `IdentityMapping`.
+Constructs the inverse inner product operator `H⁻¹` as a `TensorMapping` using
+the weights of `H`, `interior_weight`, `closure_weights`. `H⁻¹` is inverse of
+the inner product operator `H`. The weights are the
 
-`inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil)`
-constructs `H⁻¹` using a set of stencils `inv_closure_stencils` for the points
-in the closure regions and the stencil `inv_inner_stencil` in the interior.
-
-On a 1-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a N-dimensional
-`grid`, `H⁻¹` is the outer product of the 1-dimensional inverse inner product
-operators in each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
-`H⁻¹` is a 0-dimensional `IdentityMapping`.
+On a 1-dimensional grid, `H⁻¹` is a `ConstantInteriorScalingOperator`. On an
+N-dimensional grid, `H⁻¹` is the outer product of the 1-dimensional inverse
+inner product operators for each coordinate direction. On a 0-dimensional
+`grid`, `H⁻¹` is a 0-dimensional `IdentityMapping`.
 """
-function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil)
+function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
     H⁻¹s = ()
 
     for i ∈ 1:dimension(grid)
-        H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), inv_closure_stencils, inv_inner_stencil))
+        H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), interior_weight, closure_weights))
     end
 
     return foldl(⊗, H⁻¹s)
 end
 
-function inverse_inner_product(grid::EquidistantGrid{1}, inv_closure_stencils, inv_inner_stencil)
-    h⁻¹ = inverse_spacing(grid)
-    H⁻¹ = SbpOperators.volume_operator(grid, scale(inv_inner_stencil, h⁻¹[1]), scale.(inv_closure_stencils, h⁻¹[1]),even,1)
+function inverse_inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+    h⁻¹ = inverse_spacing(grid)[1]
+    H⁻¹ = SbpOperators.ConstantInteriorScalingOperator(grid, h⁻¹*1/interior_weight, h⁻¹./closure_weights)
     return H⁻¹
 end
 export inverse_inner_product
 
-inverse_inner_product(grid::EquidistantGrid{0}, inv_closure_stencils, inv_inner_stencil) = IdentityMapping{eltype(grid)}()
-
-reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
+inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -0,0 +1,36 @@
+using Test
+
+using Sbplib.LazyTensors
+using Sbplib.SbpOperators
+import Sbplib.SbpOperators: ConstantInteriorScalingOperator
+using Sbplib.Grids
+
+@testset "ConstantInteriorScalingOperator" begin
+    @test ConstantInteriorScalingOperator(1, (2,3), 10) isa ConstantInteriorScalingOperator{Int,2}
+    @test ConstantInteriorScalingOperator(1., (2.,3.), 10) isa ConstantInteriorScalingOperator{Float64,2}
+
+    a = ConstantInteriorScalingOperator(4, (2,3), 10)
+    v = ones(Int, 10)
+    @test a*v == [2,3,4,4,4,4,4,4,3,2]
+    @test a'*v == [2,3,4,4,4,4,4,4,3,2]
+
+    @test range_size(a) == (10,)
+    @test domain_size(a) == (10,)
+
+
+    a = ConstantInteriorScalingOperator(.5, (.1,.2), 7)
+    v = ones(7)
+
+    @test a*v == [.1,.2,.5,.5,.5,.2,.1]
+    @test a'*v == [.1,.2,.5,.5,.5,.2,.1]
+
+    @test range_size(a) == (7,)
+    @test domain_size(a) == (7,)
+
+    @test_throws DomainError ConstantInteriorScalingOperator(4,(2,3), 3)
+
+    @testset "Grid constructor" begin
+        g = EquidistantGrid(11, 0., 2.)
+        @test ConstantInteriorScalingOperator(g, 3., (.1,.2)) isa ConstantInteriorScalingOperator{Float64}
+    end
+end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Thu Jul 22 14:26:34 2021 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sun Jul 25 15:39:29 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
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Thu Jul 22 14:26:34 2021 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -12,34 +12,38 @@
     g_1D = EquidistantGrid(77, 0.0, Lx)
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
     @testset "inverse_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
-            Hi = inverse_inner_product(EquidistantGrid{Float64}(),SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
+            Hi = inverse_inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
             @test Hi == IdentityMapping{Float64}()
             @test Hi isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
-            Hi = inverse_inner_product(g_1D, SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.));
+            Hi = inverse_inner_product(g_1D,  quadrature_interior, quadrature_closure)
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D,op.quadratureClosure, CenteredStencil(1.))
-            Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure, CenteredStencil(1.))
-            Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure, CenteredStencil(1.))
+            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
+            Hi_x = inverse_inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure)
+            Hi_y = inverse_inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure)
             @test Hi == Hi_x⊗Hi_y
             @test Hi 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
-            Hi = inverse_inner_product(g_1D,op.quadratureClosure, CenteredStencil(1.))
+            Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
             @test domain_size(Hi) == size(g_1D)
             @test range_size(Hi) == size(g_1D)
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D,op.quadratureClosure, CenteredStencil(1.))
+            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -50,16 +54,20 @@
             v = evalOn(g_1D,x->sin(x))
             u = evalOn(g_1D,x->x^3-x^2+1)
             @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.))
-                Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(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)
+                Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             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.))
-                Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(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)
+                Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -68,16 +76,20 @@
             v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(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.))
-                Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(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)
+                Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             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.))
-                Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(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)
+                Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Thu Jul 22 14:26:34 2021 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sun Jul 25 15:39:29 2021 +0200
@@ -8,18 +8,20 @@
     g_1D = EquidistantGrid(101, 0.0, 1.)
     g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
     @testset "Constructors" 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)
+        inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+        closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
         @testset "1D" begin
-            L = laplace(g_1D, op.innerStencil, op.closureStencils)
-            @test L == second_derivative(g_1D, op.innerStencil, op.closureStencils)
+            L = laplace(g_1D, inner_stencil, closure_stencils)
+            @test L == second_derivative(g_1D, inner_stencil, closure_stencils)
             @test L isa TensorMapping{T,1,1}  where T
         end
         @testset "3D" begin
-            L = laplace(g_3D, op.innerStencil, op.closureStencils)
+            L = laplace(g_3D, inner_stencil, closure_stencils)
             @test L isa TensorMapping{T,3,3} where T
-            Dxx = second_derivative(g_3D, op.innerStencil, op.closureStencils,1)
-            Dyy = second_derivative(g_3D, op.innerStencil, op.closureStencils,2)
-            Dzz = second_derivative(g_3D, op.innerStencil, op.closureStencils,3)
+            Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
+            Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
+            Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
             @test L == Dxx + Dyy + Dzz
         end
     end
@@ -40,8 +42,10 @@
         # 2nd order interior stencil, 1st order boundary stencil,
         # implies that L*v should be exact for binomials up to order 2.
         @testset "2nd order" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-            L = laplace(g_3D,op.innerStencil,op.closureStencils)
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+            closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+            L = laplace(g_3D, inner_stencil, closure_stencils)
             @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
             @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
             @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
@@ -51,8 +55,10 @@
         # 4th order interior stencil, 2nd order boundary stencil,
         # implies that L*v should be exact for binomials up to order 3.
         @testset "4th order" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            L = laplace(g_3D,op.innerStencil,op.closureStencils)
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+            closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+            L = laplace(g_3D, inner_stencil, closure_stencils)
             # NOTE: high tolerances for checking the "exact" differentiation
             # due to accumulation of round-off errors/cancellation errors?
             @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9