Mercurial > repos > public > sbplib_julia

changeset 824:18b22737fc26 operator_storage_array_of_table

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 21 Dec 2021 15:36:14 +0100
parents 219c9661e700 (diff) 181c3c93463d (current diff)
children 69700b0b1452
files
diffstat 23 files changed, 572 insertions(+), 481 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Tue Dec 21 15:34:51 2021 +0100
+++ b/Notes.md	Tue Dec 21 15:36:14 2021 +0100
@@ -330,3 +330,5 @@
 A different approach would be to include it as a trait for operators so that you can specify what the adjoint for that operator is.
 
 
+## Name of the `VolumeOperator` type for constant stencils
+It seems that the name is too general. The name of the method `volume_operator` makes sense. It should return different types of `TensorMapping` specialized for the grid. A suggetion for a better name is `ConstantStencilVolumeOperator`
--- a/src/SbpOperators/SbpOperators.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -4,11 +4,16 @@
 using Sbplib.LazyTensors
 using Sbplib.Grids
 
+@enum Parity begin
+    odd = -1
+    even = 1
+end
+
 include("stencil.jl")
-include("d2.jl")
 include("readoperator.jl")
 include("volumeops/volume_operator.jl")
-include("volumeops/derivatives/secondderivative.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")
 include("volumeops/inner_products/inverse_inner_product.jl")
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -9,7 +9,7 @@
 On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
 a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
 """
-boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
-boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+boundary_restriction(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
+boundary_restriction(grid::EquidistantGrid{1}, closure_stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
 
 export boundary_restriction
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -9,10 +9,10 @@
 On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
 a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
 """
-function normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+function normal_derivative(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
 end
-normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+normal_derivative(grid::EquidistantGrid{1}, closure_stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
 export normal_derivative
--- a/src/SbpOperators/d2.jl	Tue Dec 21 15:34:51 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-export D2, closuresize
-
-@enum Parity begin
-    odd = -1
-    even = 1
-end
-
-struct D2{T,M}
-    innerStencil::Stencil{T}
-    closureStencils::NTuple{M,Stencil{T}}
-    eClosure::Stencil{T}
-    dClosure::Stencil{T}
-    quadratureClosure::NTuple{M,Stencil{T}}
-    parity::Parity
-end
-
-closuresize(D::D2{T,M}) where {T,M} = M
--- a/src/SbpOperators/operators/standard_diagonal.toml	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Tue Dec 21 15:36:14 2021 +0100
@@ -1,36 +1,42 @@
 [meta]
 authors = "Ken Mattson"
-descripion = "Standard operators for equidistant grids"
+description = "Standard operators for equidistant grids"
 type = "equidistant"
+cite = "A paper a long time ago in a galaxy far far away."
 
-[order2]
-H.inner = ["1"]
+[[stencil_set]]
+
+order = 2
+
+H.inner = "1"
 H.closure = ["1/2"]
 
 D1.inner_stencil = ["-1/2", "0", "1/2"]
 D1.closure_stencils = [
-    ["-1", "1"],
+    {s = ["-1", "1"], c = 1},
 ]
 
 D2.inner_stencil = ["1", "-2", "1"]
 D2.closure_stencils = [
-    ["1", "-2", "1"],
+    {s = ["1", "-2", "1"], c = 1},
 ]
 
 e.closure = ["1"]
-d1.closure = ["-3/2", "2", "-1/2"]
+d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1}
 
-[order4]
-H.inner = ["1"]
+[[stencil_set]]
+
+order = 4
+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"]
 D2.closure_stencils = [
-    [     "2",    "-5",      "4",       "-1",     "0",     "0"],
-    [     "1",    "-2",      "1",        "0",     "0",     "0"],
-    [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"],
-    [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"],
+    {s = [     "2",    "-5",      "4",       "-1",     "0",     "0"], c = 1},
+    {s = [     "1",    "-2",      "1",        "0",     "0",     "0"], c = 2},
+    {s = [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"], c = 3},
+    {s = [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"], c = 4},
 ]
 
 e.closure = ["1"]
-d1.closure = ["-11/6", "3", "-3/2", "1/3"]
+d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1}
--- a/src/SbpOperators/readoperator.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/readoperator.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -1,141 +1,110 @@
 using TOML
 
-export read_D2_operator
-export read_stencil
-export read_stencils
-export read_tuple
-
-export get_stencil
-export get_stencils
-export get_tuple
+export read_stencil_set
+export get_stencil_set
 
-function read_D2_operator(fn; order)
-    operators = TOML.parsefile(fn)["order$order"]
-    D2 = operators["D2"]
-    H = operators["H"]
-    e = operators["e"]
-    d1 = operators["d1"]
+export parse_stencil
+export parse_rational
 
-    # Create inner stencil
-    innerStencil = get_stencil(operators, "D2", "inner_stencil")
+export sbp_operators_path
 
-    # Create boundary stencils
-    boundarySize = length(D2["closure_stencils"])
-    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
-    for i ∈ 1:boundarySize
-        closureStencils = (closureStencils..., get_stencil(operators, "D2", "closure_stencils", i; center=i))
-    end
-    # TODO: Get rid of the padding here. Any padding should be handled by the consturctor accepting the stencils.
-    eClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, e["closure"]), boundarySize)..., center=1)
-    dClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, d1["closure"]), boundarySize)..., center=1)
+# The read_stencil_set and get_stencil_set functions return the freshly parsed
+# toml. The generic code in these functions can't be expected to know anyhting
+# about how to read different stencil sets as they may contain many different
+# kinds of stencils. We should how ever add read_ and get_ functions for all
+# the types of stencils we know about.
+#
+# After getting a stencil set the user can use parse functions to parse what
+# they want from the TOML dict. I.e no more "paths".
+# Functions needed:
+#   * parse stencil
+#   * parse rational
+#
+# maybe there is a better name than parse?
+# Would be nice to be able to control the type in the stencil
 
-    q_tuple = pad_tuple(toml_string_array_to_tuple(Float64, H["closure"]), boundarySize)
-    quadratureClosure = Vector{typeof(innerStencil)}()
-    for i ∈ 1:boundarySize
-        quadratureClosure = (quadratureClosure..., Stencil(q_tuple[i], center=1))
-    end
+# TODO: Control type for the stencil
+# TODO: Think about naming and terminology around freshly parsed TOML.
+# Vidar: What about get_stencil instead of parse_stencil for an already parsed
+# toml. It matches get_stencil_set.
 
-    d2 = SbpOperators.D2(
-        innerStencil,
-        closureStencils,
-        eClosure,
-        dClosure,
-        quadratureClosure,
-        even
-    )
-
-    return d2
-end
-
+# TODO: Docs for readoperator.jl
+    # Parsing as rationals is intentional, allows preserving exactness, which can be lowered using converts or promotions later.
+# TODO: readoperator.jl file name?
+# TODO: Remove references to toml for dict-input arguments
 
 """
-    read_stencil(fn, path...; [center])
+    read_stencil_set(fn; filters)
 
-Read a stencil at `path` from the file with name `fn`.
-If a center is specified the given element of the stecil is set as the center.
-
-See also: [`read_stencils`](@ref), [`read_tuple`](@ref), [`get_stencil`](@ref).
+Picks out a stencil set from the given toml file based on some filters.
+If more than one set matches the filters an error is raised.
 
-# Examples
-```
-read_stencil(sbp_operators_path()*"standard_diagonal.toml", "order2", "D2", "inner_stencil")
-read_stencil(sbp_operators_path()*"standard_diagonal.toml", "order2", "d1", "closure"; center=1)
-```
-"""
-read_stencil(fn, path...; center=nothing) = get_stencil(TOML.parsefile(fn), path...; center=center)
-
+The stencil set is not parsed beyond the inital toml parse. To get usable
+stencils use the `parse_stencil` functions on the fields of the stencil set.
 """
-    read_stencils(fn, path...; centers)
-
-Read stencils at `path` from the file `fn`.
-Centers of the stencils are specified as a tuple or array in `centers`.
-
-See also: [`read_stencil`](@ref), [`read_tuple`](@ref), [`get_stencils`](@ref).
-"""
-read_stencils(fn, path...; centers) = get_stencils(TOML.parsefile(fn), path...; centers=centers)
+read_stencil_set(fn; filters...) = get_stencil_set(TOML.parsefile(fn); filters...)
 
 """
-    read_tuple(fn, path...)
-
-Read tuple at `path` from the file `fn`.
-
-See also: [`read_stencil`](@ref), [`read_stencils`](@ref), [`get_tuple`](@ref).
-"""
-read_tuple(fn, path...) = get_tuple(TOML.parsefile(fn), path...)
-
-"""
-    get_stencil(parsed_toml, path...; center=nothing)
+    get_stencil_set(parsed_toml; filters...)
 
-Same as [`read_stencil`](@ref)) but takes already parsed toml.
+Same as `read_stencil_set` but works on already parsed TOML.
 """
-get_stencil(parsed_toml, path...; center=nothing) = get_stencil(parsed_toml[path[1]], path[2:end]...; center=center)
-function get_stencil(parsed_toml; center=nothing)
-    @assert parsed_toml isa Vector{String}
-    stencil_weights = Float64.(parse_rational.(parsed_toml))
+function get_stencil_set(parsed_toml; filters...)
+    matches = findall(parsed_toml["stencil_set"]) do set
+        for (key, val) ∈ filters
+            if set[string(key)] != val
+                return false
+            end
+        end
 
-    width = length(stencil_weights)
-
-    if isnothing(center)
-        center = div(width,2)+1
+        return true
     end
 
-    return Stencil(stencil_weights..., center=center)
+    if length(matches) != 1
+        throw(ArgumentError("filters must pick out a single stencil set"))
+    end
+
+    i = matches[1]
+    return parsed_toml["stencil_set"][i]
 end
 
 """
-    get_stencils(parsed_toml, path...; centers)
+    parse_stencil(toml)
 
-Same as [`read_stencils`](@ref)) but takes already parsed toml.
+Accepts parsed toml and reads it as a stencil
 """
-get_stencils(parsed_toml, path...; centers) = get_stencils(parsed_toml[path[1]], path[2:end]...; centers=centers)
-function get_stencils(parsed_toml; centers)
-    @assert parsed_toml isa Vector{Vector{String}}
-    @assert length(centers) == length(parsed_toml)
+function parse_stencil(toml)
+    check_stencil_toml(toml)
 
-    stencils = ()
-    for i ∈ 1:length(parsed_toml)
-        stencil = get_stencil(parsed_toml[i], center = centers[i])
-        stencils = (stencils..., stencil)
+    if toml isa Array
+        weights = parse_rational.(toml)
+        return CenteredStencil(weights...)
     end
 
-    return stencils
+    weights = parse_rational.(toml["s"])
+    return Stencil(weights..., center = toml["c"])
 end
 
-"""
-    get_tuple(parsed_toml, path...)
+function check_stencil_toml(toml)
+    if !(toml isa Dict || toml isa Vector{String})
+        throw(ArgumentError("the TOML for a stencil must be a vector of strings or a table."))
+    end
+
+    if toml isa Vector{String}
+        return
+    end
 
-Same as [`read_tuple`](@ref)) but takes already parsed toml.
-"""
-get_tuple(parsed_toml, path...) = get_tuple(parsed_toml[path[1]], path[2:end]...)
-function get_tuple(parsed_toml)
-    @assert parsed_toml isa Vector{String}
-    t = Tuple(Float64.(parse_rational.(parsed_toml)))
-    return t
-end
+    if !(haskey(toml, "s") && haskey(toml, "c"))
+        throw(ArgumentError("the table form of a stencil must have fields `s` and `c`."))
+    end
 
-# TODO: Probably should be deleted once we have gotten rid of read_D2_operator()
-function toml_string_array_to_tuple(::Type{T}, arr::AbstractVector{String}) where T
-    return Tuple(T.(parse_rational.(arr)))
+    if !(toml["s"] isa Vector{String})
+        throw(ArgumentError("a stencil must be specified as a vector of strings."))
+    end
+
+    if !(toml["c"] isa Int)
+        throw(ArgumentError("the center of a stencil must be specified as an integer."))
+    end
 end
 
 function parse_rational(str)
@@ -143,13 +112,4 @@
     return eval(:(Rational($expr)))
 end
 
-function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
-    if N >= n
-        return t
-    else
-        return pad_tuple((t..., zero(T)), n)
-    end
-end
-
 sbp_operators_path() = (@__DIR__) * "/operators/"
-export sbp_operators_path
--- a/src/SbpOperators/stencil.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/stencil.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -1,10 +1,10 @@
 export CenteredStencil
 
-struct Stencil{T<:Real,N}
+struct Stencil{T,N}
     range::Tuple{Int,Int}
     weights::NTuple{N,T}
 
-    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
+    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N}
         @assert range[2]-range[1]+1 == N
         new{T,N}(range,weights)
     end
@@ -15,7 +15,7 @@
 
 Create a stencil with the given weights with element `center` as the center of the stencil.
 """
-function Stencil(weights::Vararg{Number}; center::Int)
+function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error
     N = length(weights)
     range = (1, N) .- center
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -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)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -0,0 +1,20 @@
+"""
+    second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
+    second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
+
+Creates the second-derivative operator `D2` as a `TensorMapping`
+
+`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by
+`direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
+for the points in the closure regions.
+
+On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the outer product of the
+one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
+Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
+"""
+function second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
+    h_inv = inverse_spacing(grid)[direction]
+    return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
+end
+second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
+export second_derivative
--- a/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Tue Dec 21 15:34:51 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,20 +0,0 @@
-"""
-    second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
-    second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
-
-Creates the second-derivative operator `D2` as a `TensorMapping`
-
-`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by
-`direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions.
-
-On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the outer product of the
-one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
-Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
-"""
-function second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
-    h_inv = inverse_spacing(grid)[direction]
-    return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
-end
-second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
-export second_derivative
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -1,29 +1,34 @@
 """
-    inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
+    inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
 
-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. If
-`inner_stencil` is omitted a central interior stencil with weight 1 is used.
+`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 = CenteredStencil(one(eltype(grid))))
-    h = spacing(grid)
-    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
-    for i ∈ 2:dimension(grid)
-        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
-        H = H∘Hᵢ
+function inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
+    Hs = ()
+
+    for i ∈ 1:dimension(grid)
+        Hs = (Hs..., inner_product(restrict(grid, i), interior_weight, closure_weights))
     end
-    return H
+
+    return foldl(⊗, Hs)
 end
 export inner_product
 
-inner_product(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()
+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}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -1,43 +1,30 @@
 """
-    inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
-    inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})
+    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`.
-
-`inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)`
-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. If
-`inv_closure_stencils` is omitted, a central interior stencil with weight 1 is used.
+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, closure_stencils::NTuple{M,Stencil{T,1}})`
-constructs a diagonal inverse inner product operator where `closure_stencils` are the
-closure stencils of `H` (not `H⁻¹`!).
+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, interior_weight, closure_weights)
+    H⁻¹s = ()
 
-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`.
-"""
-function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid))))
-    h⁻¹ = inverse_spacing(grid)
-    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
-    for i ∈ 2:dimension(grid)
-        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
-        H⁻¹ = H⁻¹∘Hᵢ⁻¹
+    for i ∈ 1:dimension(grid)
+        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}, 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)}()
-
-function inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T}
-     inv_closure_stencils = reciprocal_stencil.(closure_stencils)
-     inv_inner_stencil = CenteredStencil(one(T))
-     return inverse_inner_product(grid, inv_closure_stencils, inv_inner_stencil)
-end
-
-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)}()
--- a/src/SbpOperators/volumeops/volume_operator.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -1,14 +1,15 @@
 """
-    volume_operator(grid,inner_stencil,closure_stencils,parity,direction)
+    volume_operator(grid, inner_stencil, closure_stencils, parity, direction)
+
 Creates a volume operator on a `Dim`-dimensional grid acting along the
-specified coordinate `direction`. The action of the operator is determined by the
-stencils `inner_stencil` and `closure_stencils`.
-When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned.
-When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product
-of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
-the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
+specified coordinate `direction`. The action of the operator is determined by
+the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the
+corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the
+returned operator is the appropriate outer product of a one-dimensional
+operators and `IdentityMapping`s, e.g for `Dim=3` the volume operator in the
+y-direction is `I⊗op⊗I`.
 """
-function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M}
+function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil, closure_stencils, parity, direction) where {Dim,T}
     #TODO: Check that direction <= Dim?
 
     # Create 1D volume operator in along coordinate direction
@@ -34,7 +35,7 @@
 end
 
 function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity)
-    return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity)
+    return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity)
 end
 
 closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -8,27 +8,28 @@
 import Sbplib.SbpOperators.BoundaryOperator
 
 @testset "boundary_restriction" 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)
+	e_closure = parse_stencil(stencil_set["e"]["closure"])
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D,op.eClosure,Lower())
-            @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
-            @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
+            e_l = boundary_restriction(g_1D,e_closure,Lower())
+            @test e_l == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
+            @test e_l == BoundaryOperator(g_1D,e_closure,Lower())
             @test e_l isa BoundaryOperator{T,Lower} where T
             @test e_l isa TensorMapping{T,0,1} where T
 
-            e_r = boundary_restriction(g_1D,op.eClosure,Upper())
-            @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
-            @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper())
+            e_r = boundary_restriction(g_1D,e_closure,Upper())
+            @test e_r == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
+            @test e_r == BoundaryOperator(g_1D,e_closure,Upper())
             @test e_r isa BoundaryOperator{T,Upper} where T
             @test e_r isa TensorMapping{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            e_w = boundary_restriction(g_2D,e_closure,CartesianBoundary{1,Upper}())
             @test e_w isa InflatedTensorMapping
             @test e_w isa TensorMapping{T,1,2} where T
         end
@@ -36,8 +37,8 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+            e_l = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Lower}())
+            e_r = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Upper}())
 
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
@@ -48,10 +49,10 @@
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
-            e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
-            e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
+            e_w = boundary_restriction(g_2D, e_closure, CartesianBoundary{1,Lower}())
+            e_e = boundary_restriction(g_2D, e_closure, CartesianBoundary{1,Upper}())
+            e_s = boundary_restriction(g_2D, e_closure, CartesianBoundary{2,Lower}())
+            e_n = boundary_restriction(g_2D, e_closure, CartesianBoundary{2,Upper}())
 
             v = rand(11, 15)
             u = fill(3.124)
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -11,21 +11,21 @@
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
     @testset "normal_derivative" 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)
+    	d_closure = parse_stencil(stencil_set["d1"]["closure"])
         @testset "1D" begin
-            d_l = normal_derivative(g_1D, op.dClosure, Lower())
-            @test d_l == normal_derivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_1D, d_closure, Lower())
+            @test d_l == normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}())
             @test d_l isa BoundaryOperator{T,Lower} where T
             @test d_l isa TensorMapping{T,0,1} where T
         end
         @testset "2D" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
+            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
             Ix = IdentityMapping{Float64}((size(g_2D)[1],))
             Iy = IdentityMapping{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(restrict(g_2D,1),op.dClosure,Lower())
-            d_r = normal_derivative(restrict(g_2D,2),op.dClosure,Upper())
+            d_l = normal_derivative(restrict(g_2D,1),d_closure,Lower())
+            d_r = normal_derivative(restrict(g_2D,2),d_closure,Upper())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa TensorMapping{T,1,2} where T
@@ -38,11 +38,12 @@
         v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x)
         # TODO: Test for higher order polynomials?
         @testset "2nd order" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
-            d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
-            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+        	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+        	d_closure = parse_stencil(stencil_set["d1"]["closure"])
+            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
@@ -51,11 +52,12 @@
         end
 
         @testset "4th order" begin
-            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
-            d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
-            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+        	d_closure = parse_stencil(stencil_set["d1"]["closure"])
+            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
--- a/test/SbpOperators/readoperator_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/readoperator_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -18,76 +18,98 @@
 @testset "readoperator" begin
     toml_str = """
         [meta]
+        authors = "Ken Mattson"
+        description = "Standard operators for equidistant grids"
         type = "equidistant"
+        cite = "A paper a long time ago in a galaxy far far away."
 
-        [order2]
+        [[stencil_set]]
+
+        order = 2
+        test = 2
+
         H.inner = ["1"]
+        H.closure = ["1/2"]
 
         D1.inner_stencil = ["-1/2", "0", "1/2"]
         D1.closure_stencils = [
-            ["-1", "1"],
+            {s = ["-1", "1"], c = 1},
+        ]
+
+        D2.inner_stencil = ["1", "-2", "1"]
+        D2.closure_stencils = [
+            {s = ["1", "-2", "1"], c = 1},
         ]
 
-        d1.closure = ["-3/2", "2", "-1/2"]
+        e.closure = ["1"]
+        d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1}
 
-        [order4]
+        [[stencil_set]]
+
+        order = 4
+        test = 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"]
         D2.closure_stencils = [
-            [     "2",    "-5",      "4",       "-1",     "0",     "0"],
-            [     "1",    "-2",      "1",        "0",     "0",     "0"],
-            [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"],
-            [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"],
+            {s = [     "2",    "-5",      "4",       "-1",     "0",     "0"], c = 1},
+            {s = [     "1",    "-2",      "1",        "0",     "0",     "0"], c = 2},
+            {s = [ "-4/43", "59/43", "-110/43",   "59/43", "-4/43",     "0"], c = 3},
+            {s = [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"], c = 4},
         ]
+
+        e.closure = ["1"]
+        d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1}
+
+        [[stencil_set]]
+        order = 4
+        test = 2
+
+        H.closure = ["-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"]
     """
 
     parsed_toml = TOML.parse(toml_str)
-    @testset "get_stencil" begin
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil(-1/2, 0., 1/2, center=2)
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil(-1/2, 0., 1/2; center=1)
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil(-1/2, 0., 1/2; center=3)
 
-        @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil(1.; center=1)
+    @testset "get_stencil_set" begin
+        @test get_stencil_set(parsed_toml; order = 2) isa Dict
+        @test get_stencil_set(parsed_toml; order = 2) == parsed_toml["stencil_set"][1]
+        @test get_stencil_set(parsed_toml; test = 1) == parsed_toml["stencil_set"][2]
+        @test get_stencil_set(parsed_toml; order = 4, test = 2) == parsed_toml["stencil_set"][3]
 
-        @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
-        @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
+        @test_throws ArgumentError get_stencil_set(parsed_toml; test = 2)
+        @test_throws ArgumentError get_stencil_set(parsed_toml; order = 4)
     end
 
-    @testset "get_stencils" begin
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil(-1., 1., center=1),)
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil(-1., 1., center=2),)
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil(-1., 1., center=2),)
-
-        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
-            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=1),
-            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=1),
-            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=1),
-            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=1),
-        )
+    @testset "parse_stencil" begin
+        toml = """
+            s1 = ["-1/12","4/3","-5/2","4/3","-1/12"]
+            s2 = {s = ["2", "-5", "4", "-1", "0", "0"], c = 1}
+            s3 = {s = ["1", "-2", "1", "0", "0", "0"], c = 2}
+            s4 = "not a stencil"
+            s5 = [-1, 4, 3]
+            s6 = {k = ["1", "-2", "1", "0", "0", "0"], c = 2}
+            s7 = {s = [-1, 4, 3], c = 2}
+            s8 = {s = ["1", "-2", "1", "0", "0", "0"], c = [2,2]}
+        """
 
-        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
-            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=4),
-            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=2),
-            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=3),
-            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=1),
-        )
+        @test parse_stencil(TOML.parse(toml)["s1"]) == CenteredStencil(-1//12, 4//3, -5//2, 4//3, -1//12)
+        @test parse_stencil(TOML.parse(toml)["s2"]) == Stencil(2//1, -5//1, 4//1, -1//1, 0//1, 0//1; center=1)
+        @test parse_stencil(TOML.parse(toml)["s3"]) == Stencil(1//1, -2//1, 1//1, 0//1, 0//1, 0//1; center=2)
 
-        @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
-            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=1),
-            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=2),
-            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=3),
-            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=4),
-        )
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s4"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s5"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s6"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s7"])
+        @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s8"])
 
-        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
-        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4))
-        @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2))
-    end
+        stencil_set = get_stencil_set(parsed_toml; order = 4, test = 1)
 
-    @testset "get_tuple" begin
-        @test get_tuple(parsed_toml, "order2", "d1", "closure") == (-3/2, 2, -1/2)
-
-        @test_throws AssertionError get_tuple(parsed_toml, "meta", "type")
+        @test parse_stencil.(stencil_set["D2"]["closure_stencils"]) == [
+            Stencil(  2//1,  -5//1,     4//1,    -1//1,   0//1,   0//1; center=1),
+            Stencil(  1//1,  -2//1,     1//1,     0//1,   0//1,   0//1; center=2),
+            Stencil(-4//43, 59//43, -110//43,   59//43, -4//43,   0//1; center=3),
+            Stencil(-1//49,   0//1,   59//49, -118//49, 64//49, -4//49; center=4),
+        ]
     end
 end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -0,0 +1,118 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+import Sbplib.SbpOperators.VolumeOperator
+
+@testset "SecondDerivative" begin
+    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"])
+    Lx = 3.5
+    Ly = 3.
+    g_1D = EquidistantGrid(121, 0.0, Lx)
+    g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
+
+    @testset "Constructors" begin
+        @testset "1D" begin
+            Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+            @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils,1)
+            @test Dₓₓ isa VolumeOperator
+        end
+        @testset "2D" begin
+            Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1)
+            D2 = second_derivative(g_1D,inner_stencil,closure_stencils)
+            I = IdentityMapping{Float64}(size(g_2D)[2])
+            @test Dₓₓ == D2⊗I
+            @test Dₓₓ isa TensorMapping{T,2,2} where T
+        end
+    end
+
+    # Exact differentiation is measured point-wise. In other cases
+    # the error is measured in the l2-norm.
+    @testset "Accuracy" begin
+        @testset "1D" begin
+            l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
+            monomials = ()
+            maxOrder = 4;
+            for i = 0:maxOrder-1
+                f_i(x) = 1/factorial(i)*x^i
+                monomials = (monomials...,evalOn(g_1D,f_i))
+            end
+            v = evalOn(g_1D,x -> sin(x))
+            vₓₓ = evalOn(g_1D,x -> -sin(x))
+
+            # 2nd order interior stencil, 1nd order boundary stencil,
+            # implies that L*v should be exact for monomials up to order 2.
+            @testset "2nd order" begin
+                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"])
+                Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+                @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
+                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
+            end
+
+            # 4th order interior stencil, 2nd order boundary stencil,
+            # implies that L*v should be exact for monomials up to order 3.
+            @testset "4th order" begin
+                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"])
+                Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+                # NOTE: high tolerances for checking the "exact" differentiation
+                # due to accumulation of round-off errors/cancellation errors?
+                @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
+                @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
+                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
+            end
+        end
+
+        @testset "2D" begin
+            l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
+            binomials = ()
+            maxOrder = 4;
+            for i = 0:maxOrder-1
+                f_i(x,y) = 1/factorial(i)*y^i + x^i
+                binomials = (binomials...,evalOn(g_2D,f_i))
+            end
+            v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
+            v_yy = evalOn(g_2D,(x,y) -> -cos(y))
+
+            # 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
+                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"])
+                Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2)
+                @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
+                @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
+            end
+
+            # 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
+                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"])
+                Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2)
+                # NOTE: high tolerances for checking the "exact" differentiation
+                # due to accumulation of round-off errors/cancellation errors?
+                @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
+                @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
+                @test Dyy*v ≈ v_yy rtol = 5e-4 norm = l2
+            end
+        end
+    end
+end
--- a/test/SbpOperators/volumeops/derivatives/secondderivative_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,108 +0,0 @@
-using Test
-
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
-
-import Sbplib.SbpOperators.VolumeOperator
-
-@testset "SecondDerivative" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    Lx = 3.5
-    Ly = 3.
-    g_1D = EquidistantGrid(121, 0.0, Lx)
-    g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
-
-    @testset "Constructors" begin
-        @testset "1D" begin
-            Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
-            @test Dₓₓ == second_derivative(g_1D,op.innerStencil,op.closureStencils,1)
-            @test Dₓₓ isa VolumeOperator
-        end
-        @testset "2D" begin
-            Dₓₓ = second_derivative(g_2D,op.innerStencil,op.closureStencils,1)
-            D2 = second_derivative(g_1D,op.innerStencil,op.closureStencils)
-            I = IdentityMapping{Float64}(size(g_2D)[2])
-            @test Dₓₓ == D2⊗I
-            @test Dₓₓ isa TensorMapping{T,2,2} where T
-        end
-    end
-
-    # Exact differentiation is measured point-wise. In other cases
-    # the error is measured in the l2-norm.
-    @testset "Accuracy" begin
-        @testset "1D" begin
-            l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
-            monomials = ()
-            maxOrder = 4;
-            for i = 0:maxOrder-1
-                f_i(x) = 1/factorial(i)*x^i
-                monomials = (monomials...,evalOn(g_1D,f_i))
-            end
-            v = evalOn(g_1D,x -> sin(x))
-            vₓₓ = evalOn(g_1D,x -> -sin(x))
-
-            # 2nd order interior stencil, 1nd order boundary stencil,
-            # implies that L*v should be exact for monomials up to order 2.
-            @testset "2nd order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
-                @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
-                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
-                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
-                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
-            end
-
-            # 4th order interior stencil, 2nd order boundary stencil,
-            # implies that L*v should be exact for monomials up to order 3.
-            @testset "4th order" begin
-                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
-                # NOTE: high tolerances for checking the "exact" differentiation
-                # due to accumulation of round-off errors/cancellation errors?
-                @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
-                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
-                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
-                @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
-                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
-            end
-        end
-
-        @testset "2D" begin
-            l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
-            binomials = ()
-            maxOrder = 4;
-            for i = 0:maxOrder-1
-                f_i(x,y) = 1/factorial(i)*y^i + x^i
-                binomials = (binomials...,evalOn(g_2D,f_i))
-            end
-            v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
-            v_yy = evalOn(g_2D,(x,y) -> -cos(y))
-
-            # 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)
-                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
-                @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
-                @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
-                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
-                @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
-            end
-
-            # 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)
-                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
-                # NOTE: high tolerances for checking the "exact" differentiation
-                # due to accumulation of round-off errors/cancellation errors?
-                @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
-                @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
-                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
-                @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
-                @test Dyy*v ≈ v_yy rtol = 5e-4 norm = l2
-            end
-        end
-    end
-end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -14,36 +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)
+            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)
-            inner_stencil = CenteredStencil(1.)
-            @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
+            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)
-            H_x = inner_product(restrict(g_2D,1),op.quadratureClosure)
-            H_y = inner_product(restrict(g_2D,2),op.quadratureClosure)
+            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)
+            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)
+            H = inner_product(g_2D, quadrature_interior, quadrature_closure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -59,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)
+                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
@@ -68,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)
+                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
@@ -82,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)
+                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)
+                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	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -12,40 +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}(),op.quadratureClosure)
+            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, op.quadratureClosure);
-            inner_stencil = CenteredStencil(1.)
-            closures = ()
-            for i = 1:length(op.quadratureClosure)
-                closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
-            end
-            @test Hi == inverse_inner_product(g_1D,closures,inner_stencil)
+            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)
-            Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure)
-            Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure)
+            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)
+            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)
+            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
@@ -56,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)
-                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+                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)
-                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+                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
@@ -74,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)
-                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+                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)
-                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+                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	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Dec 21 15:36:14 2021 +0100
@@ -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