Mercurial > repos > public > sbplib_julia

changeset 864:9a2776352c2a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge operator_storage_array_of_table
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 19 Jan 2022 11:08:43 +0100
parents 5088de9b6d65 (current diff) 336e8d1a553c (diff)
children 1784b1c0af3e 313648b01504 dd2ab001a7b6 a378ab959b6f ffde7dad9da5
files
diffstat 22 files changed, 699 insertions(+), 371 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Wed Jan 19 07:24:36 2022 +0100
+++ b/Notes.md	Wed Jan 19 11:08:43 2022 +0100
@@ -53,6 +53,23 @@
 
 * Remove order as a table name and put it as a variable.
 
+### Parsing of stencil sets
+At the moment the only parsing that can be done at the top level is conversion
+from the toml file to a dict of strings. This forces the user to dig through
+the dictionary and apply the correct parsing methods for the different parts,
+e.g. `parse_stencil` or `parse_tuple`. While very flexible there is a tight
+coupling between what is written in the file and what code is run to make data
+in the file usable. While this coupling is hard to avoid it should be made
+explicit. This could be done by putting a reference to a parsing function in
+the operator-storage format or somehow specifying the type of each object.
+This mechanism should be extensible without changing the package. Perhaps
+there could be a way to register parsing functions or object types for the
+toml.
+
+If possible the goal should be for the parsing to get all the way to the
+stencils so that a user calls `read_stencil_set` and gets a
+dictionary-structure containing stencils, tuples, scalars and other types
+ready for input to the methods creating the operators.
 
 ## Variable second derivative
 
--- a/docs/make.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/docs/make.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -25,6 +25,7 @@
 
 pages = [
     "Home" => "index.md",
+    "operator_file_format.md",
     "Submodules" => [
         "submodules/grids.md",
         "submodules/diff_ops.md",
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/src/operator_file_format.md	Wed Jan 19 11:08:43 2022 +0100
@@ -0,0 +1,72 @@
+# Operator file format
+
+The intention is that Sbplib.jl should be a general and extensible framework
+for working with finite difference methods. It therefore includes a set of
+tools for storing and sharing operator definitions as well as a set of widely
+used operators.
+
+## Using the included operators
+
+Most users will likely access the included operators by simply passing the
+filename of the wanted operator set to the appropriate function.  The location
+of the included stencil sets can be computed using
+[`sbp_operators_path`](@ref).
+```@meta
+# TODO: provide examples of functions to pass the files to
+```
+Advanced user might want to get access to the individual objects of an
+operator file. This can be accomplished using functions such as
+* [`read_stencil_set`](@ref)
+* [`parse_scalar`](@ref)
+* [`parse_stencil`](@ref)
+* [`parse_tuple`](@ref)
+
+When parsing operator objects they are interpreted using `Rational`s and
+possibly have to be converted to a desired type before use. This allows
+preserving maximum accuracy when needed.
+```@meta
+# TBD: "possibly have to be converted to a desired type before use" Is this the case? Can it be fixed?
+```
+
+## File format
+The file format is based on TOML and can be parsed using `TOML.parse`. A file
+can optionally start with a `[meta]` section which can specify things like who
+the author was, a description and how to cite the operators.
+
+After the `[meta]` section one or more stencil sets follow, each one beginning
+with `[[stencil_set]]`. Each stencil set should include descriptors like
+`order`, `name` or `number_of_bondary_points` to make them unique within the
+TOML-file. What descriptors to use are up to the author of the file to decide.
+
+Beyond identifying information the stencil set can contain any valid TOML.
+This data is then parsed by the functions creating specific operators like
+``D_1`` or ``D_2``.
+
+### Numbers
+Number can be represented as regular TOML numbers e.g. `1`, `-0.4` or
+`4.32e-3`. Alternatively they can be represented as strings which allows
+specifying fraction e.g. `"1/2"` or `"0"`.
+
+All numbers are accurately converted to `Rational`s when using the
+[`parse_scalar`](@ref) function.
+
+### Stencils
+Stencils are parsed using [`parse_stencil`](@ref). They can be specified
+either as a simple arrays
+```toml
+stencil = ["-1/2","0", "1/2"]
+```
+which assumes a centered stencil. Or as a TOML inline table
+```toml
+stencil =  {s = ["-24/17", "59/34", "-4/17", "-3/34", "0", "0"], c = 1},
+```
+which allows specifying the center of the stencil using the key `c`.
+
+## Creating your own operator files
+Operator files can be created either to add new variants of existing types of
+operators like ``D_1`` or ``D_2`` or to describe completely new types of
+operators like for example a novel kind of interpolation operator. In the
+second case new parsing functions are also necessary.
+
+The files can then be used to easily test or share different variants of
+operators.
--- a/src/SbpOperators/SbpOperators.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -4,10 +4,15 @@
 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/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/boundaryops/boundary_restriction.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -9,7 +9,10 @@
 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, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
+function boundary_restriction(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
+    converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil)
+    return SbpOperators.boundary_operator(grid, converted_stencil, boundary)
+end
 boundary_restriction(grid::EquidistantGrid{1}, closure_stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
 
 export boundary_restriction
--- a/src/SbpOperators/d2.jl	Wed Jan 19 07:24:36 2022 +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	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Wed Jan 19 11:08:43 2022 +0100
@@ -1,36 +1,60 @@
 [meta]
 authors = "Ken Mattson"
-descripion = "Standard operators for equidistant grids"
+description = "Standard operators for equidistant grids"
 type = "equidistant"
+cite = """
+    Ken Mattsson, Jan Nordström,
+    Summation by parts operators for finite difference approximations of second derivatives,
+    Journal of Computational Physics,
+    Volume 199, Issue 2,
+    2004,
+    Pages 503-540,
+    ISSN 0021-9991,
+    https://doi.org/10.1016/j.jcp.2004.03.001.
+"""
 
-[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}
+
+[[stencil_set]]
+
+order = 4
 
-[order4]
-H.inner = ["1"]
+H.inner = "1"
 H.closure = ["17/48", "59/48", "43/48", "49/48"]
 
+D1.inner_stencil = ["1/12","-2/3","0","2/3","-1/12"]
+D1.closure_stencils = [
+    {s = [ "-24/17",  "59/34",  "-4/17", "-3/34",     "0",     "0"], c = 1},
+    {s = [   "-1/2",      "0",    "1/2",     "0",     "0",     "0"], c = 2},
+    {s = [   "4/43", "-59/86",      "0", "59/86", "-4/43",     "0"], c = 3},
+    {s = [   "3/98",      "0", "-59/98",     "0", "32/49", "-4/49"], c = 4},
+]
+
 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	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/readoperator.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -1,155 +1,160 @@
 using TOML
 
-export read_D2_operator
-export read_stencil
-export read_stencils
-export read_tuple
+export read_stencil_set
+export get_stencil_set
+
+export parse_stencil
+export parse_scalar
+export parse_tuple
+
+export sbp_operators_path
+
 
-export get_stencil
-export get_stencils
-export get_tuple
+"""
+    read_stencil_set(filename; filters)
 
-function read_D2_operator(fn; order)
-    operators = TOML.parsefile(fn)["order$order"]
-    D2 = operators["D2"]
-    H = operators["H"]
-    e = operators["e"]
-    d1 = operators["d1"]
+Picks out a stencil set from a TOML file based on some key-value
+filters. If more than one set matches the filters an error is raised. The
+returned stencil set contains parsed TOML intended for functions like
+`parse_scalar` and `parse_stencil`.
+
+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.
+
+The reason for this is that since stencil sets are intended to be very
+general, and currently do not include any way to specify how to parse a given
+section, the exact parsing is left to the user.
 
-    # Create inner stencil
-    innerStencil = get_stencil(operators, "D2", "inner_stencil")
+For more information see [Operator file format](@ref) in the documentation.
+
+See also [`sbp_operators_path`](@ref), [`get_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref),.
+"""
+read_stencil_set(filename; filters...) = get_stencil_set(TOML.parsefile(filename); filters...)
+
+"""
+    get_stencil_set(parsed_toml; filters...)
+
+Picks out a stencil set from an already parsed TOML based on some key-value
+filters.
 
-    # 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))
+See also [`read_stencil_set`](@ref).
+"""
+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
+
+        return true
     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)
 
-    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))
+    if length(matches) != 1
+        throw(ArgumentError("filters must pick out a single stencil set"))
     end
 
-    d2 = SbpOperators.D2(
-        innerStencil,
-        closureStencils,
-        eClosure,
-        dClosure,
-        quadratureClosure,
-        even
-    )
+    i = matches[1]
+    return parsed_toml["stencil_set"][i]
+end
+
+"""
+    parse_stencil(parsed_toml)
+
+Accepts parsed TOML and reads it as a stencil.
+
+See also [`read_stencil_set`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref).
+"""
+function parse_stencil(parsed_toml)
+    check_stencil_toml(parsed_toml)
+
+    if parsed_toml isa Array
+        weights = parse_rational.(parsed_toml)
+        return CenteredStencil(weights...)
+    end
+
+    weights = parse_rational.(parsed_toml["s"])
+    return Stencil(weights..., center = parsed_toml["c"])
+end
+
+"""
+    parse_stencil(T, parsed_toml)
+
+Parses the input as a stencil with element type `T`.
+"""
+parse_stencil(T, parsed_toml) = Stencil{T}(parse_stencil(parsed_toml))
+
+function check_stencil_toml(parsed_toml)
+    if !(parsed_toml isa Dict || parsed_toml isa Vector{String})
+        throw(ArgumentError("the TOML for a stencil must be a vector of strings or a table."))
+    end
+
+    if parsed_toml isa Vector{String}
+        return
+    end
 
-    return d2
+    if !(haskey(parsed_toml, "s") && haskey(parsed_toml, "c"))
+        throw(ArgumentError("the table form of a stencil must have fields `s` and `c`."))
+    end
+
+    if !(parsed_toml["s"] isa Vector{String})
+        throw(ArgumentError("a stencil must be specified as a vector of strings."))
+    end
+
+    if !(parsed_toml["c"] isa Int)
+        throw(ArgumentError("the center of a stencil must be specified as an integer."))
+    end
+end
+
+"""
+    parse_scalar(parsed_toml)
+
+Parse a scalar, represented as a string or a number in the TOML, and return it as a `Rational`
+
+See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref) [`parse_tuple`](@ref).
+"""
+function parse_scalar(parsed_toml)
+    try
+        return parse_rational(parsed_toml)
+    catch e
+        throw(ArgumentError("must be a number or a string representing a number."))
+    end
+end
+
+"""
+    parse_tuple(parsed_toml)
+
+Parse an array as a tuple of scalars.
+
+See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref).
+"""
+function parse_tuple(parsed_toml)
+    if !(parsed_toml isa Array)
+        throw(ArgumentError("argument must be an array"))
+    end
+    return Tuple(parse_scalar.(parsed_toml))
 end
 
 
 """
-    read_stencil(fn, path...; [center])
-
-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).
+    parse_rational(parsed_toml)
 
-# 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)
-
-"""
-    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)
-
+Parse a string or a number as a rational.
 """
-    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)
-
-Same as [`read_stencil`](@ref)) but takes 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))
-
-    width = length(stencil_weights)
-
-    if isnothing(center)
-        center = div(width,2)+1
+function parse_rational(parsed_toml)
+    if parsed_toml isa String
+        expr = Meta.parse(replace(parsed_toml, "/"=>"//"))
+        return eval(:(Rational($expr)))
+    else
+        return Rational(parsed_toml)
     end
-
-    return Stencil(stencil_weights..., center=center)
 end
 
 """
-    get_stencils(parsed_toml, path...; centers)
-
-Same as [`read_stencils`](@ref)) but takes already parsed toml.
-"""
-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)
+    sbp_operators_path()
 
-    stencils = ()
-    for i ∈ 1:length(parsed_toml)
-        stencil = get_stencil(parsed_toml[i], center = centers[i])
-        stencils = (stencils..., stencil)
-    end
+Calculate the path for the operators folder with included stencil sets.
 
-    return stencils
-end
-
-"""
-    get_tuple(parsed_toml, path...)
-
-Same as [`read_tuple`](@ref)) but takes already parsed toml.
+See also [`read_stencil_set`](@ref)
 """
-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
-
-# 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)))
-end
-
-function parse_rational(str)
-    expr = Meta.parse(replace(str, "/"=>"//"))
-    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	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/stencil.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -22,6 +22,12 @@
     return Stencil(range, weights)
 end
 
+function Stencil{T}(s::Stencil) where T
+    return Stencil(s.range, T.(s.weights))
+end
+
+Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil)
+
 function CenteredStencil(weights::Vararg)
     if iseven(length(weights))
         throw(ArgumentError("a centered stencil must have an odd number of weights."))
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Wed Jan 19 11:08:43 2022 +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 < 2*length(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	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -1,34 +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.
+`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	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -1,37 +1,30 @@
 """
-    inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
+    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`.
 
-`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.
-
-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)}()
--- a/src/StaticDicts/StaticDicts.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/src/StaticDicts/StaticDicts.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -10,9 +10,9 @@
 
 The immutable nature means that `StaticDict` can be compared with `===`, in
 constrast to regular `Dict` or `ImmutableDict` which can not. (See
-<https://github.com/JuliaLang/julia/issues/4648> for details) One important
+<https://github.com/JuliaLang/julia/issues/4648> for details.) One important
 aspect of this is that `StaticDict` can be used in a struct while still
-allowing the struct to be comared using the default implementation of `==` for
+allowing the struct to be compared using the default implementation of `==` for
 structs.
 
 Lookups are done by linear search.
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Wed Jan 19 11:08:43 2022 +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,Stencil{Float64}(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,Stencil{Float64}(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	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Jan 19 11:08:43 2022 +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	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/readoperator_test.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -5,6 +5,153 @@
 
 import Sbplib.SbpOperators.Stencil
 
+@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."
+
+        [[stencil_set]]
+
+        order = 2
+        test = 2
+
+        H.inner = ["1"]
+        H.closure = ["1/2"]
+
+        D1.inner_stencil = ["-1/2", "0", "1/2"]
+        D1.closure_stencils = [
+            {s = ["-1", "1"], c = 1},
+        ]
+
+        D2.inner_stencil = ["1", "-2", "1"]
+        D2.closure_stencils = [
+            {s = ["1", "-2", "1"], c = 1},
+        ]
+
+        e.closure = ["1"]
+        d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1}
+
+        [[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 = [
+            {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_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 ArgumentError get_stencil_set(parsed_toml; test = 2)
+        @test_throws ArgumentError get_stencil_set(parsed_toml; order = 4)
+    end
+
+    @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 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_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"])
+
+        stencil_set = get_stencil_set(parsed_toml; order = 4, test = 1)
+
+        @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),
+        ]
+
+
+        @test parse_stencil(Float64, TOML.parse(toml)["s1"]) == CenteredStencil(-1/12, 4/3, -5/2, 4/3, -1/12)
+        @test parse_stencil(Float64, TOML.parse(toml)["s2"]) == Stencil(2/1, -5/1, 4/1, -1/1, 0/1, 0/1; center=1)
+        @test parse_stencil(Float64, TOML.parse(toml)["s3"]) == Stencil(1/1, -2/1, 1/1, 0/1, 0/1, 0/1; center=2)
+    end
+
+    @testset "parse_scalar" begin
+        toml = TOML.parse("""
+            a1 = 1
+            a2 = 1.5
+            a3 = 1.0
+            a4 = 10
+            a5 = "1/2"
+            a6 = "1.5"
+
+            e1 = [1,2,3]
+            e2 = "a string value"
+        """)
+
+        @test parse_scalar(toml["a1"]) == 1//1
+        @test parse_scalar(toml["a2"]) == 3//2
+        @test parse_scalar(toml["a3"]) == 1//1
+        @test parse_scalar(toml["a4"]) == 10//1
+        @test parse_scalar(toml["a5"]) == 1//2
+        @test parse_scalar(toml["a6"]) == 3//2
+
+        @test_throws ArgumentError parse_scalar(toml["e1"])
+        @test_throws ArgumentError parse_scalar(toml["e2"])
+    end
+
+    @testset "parse_tuple" begin
+        toml = TOML.parse("""
+            t1 = [1,3,4]
+            t2 = ["1/2","3/4","2/1"]
+
+            e1 = "not a tuple"
+            e2.a="1"
+            e3 = 1
+            e4 = ["1/2","3/4","not a number"]
+        """)
+
+        @test parse_tuple(toml["t1"]) == (1//1,3//1,4//1)
+        @test parse_tuple(toml["t2"]) == (1//2,3//4,2//1)
+
+        @test_throws ArgumentError parse_tuple(toml["e1"])
+        @test_throws ArgumentError parse_tuple(toml["e2"])
+        @test_throws ArgumentError parse_tuple(toml["e3"])
+        @test_throws ArgumentError parse_tuple(toml["e4"])
+    end
+end
 
 @testset "parse_rational" begin
     @test SbpOperators.parse_rational("1") isa Rational
@@ -13,81 +160,13 @@
     @test SbpOperators.parse_rational("1/2") == 1//2
     @test SbpOperators.parse_rational("37/13") isa Rational
     @test SbpOperators.parse_rational("37/13") == 37//13
-end
 
-@testset "readoperator" begin
-    toml_str = """
-        [meta]
-        type = "equidistant"
-
-        [order2]
-        H.inner = ["1"]
-
-        D1.inner_stencil = ["-1/2", "0", "1/2"]
-        D1.closure_stencils = [
-            ["-1", "1"],
-        ]
-
-        d1.closure = ["-3/2", "2", "-1/2"]
-
-        [order4]
-        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"],
-        ]
-    """
-
-    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)
+    @test SbpOperators.parse_rational(0.5) isa Rational
+    @test SbpOperators.parse_rational(0.5) == 1//2
 
-        @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
-        @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
-    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),
-        )
+    @test SbpOperators.parse_rational("0.5") isa Rational
+    @test SbpOperators.parse_rational("0.5") == 1//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 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 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
-
-    @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")
-    end
+    @test SbpOperators.parse_rational(2) isa Rational
+    @test SbpOperators.parse_rational(2) == 2//1
 end
--- a/test/SbpOperators/stencil_test.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/stencil_test.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -15,4 +15,17 @@
 
     @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
     @test_throws ArgumentError CenteredStencil(1,2,3,4)
+
+    # Changing the type of the weights
+    @test Stencil{Float64}(Stencil(1,2,3,4,5; center=2)) == Stencil(1.,2.,3.,4.,5.; center=2)
+    @test Stencil{Float64}(CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
+    @test Stencil{Int}(Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
+    @test Stencil{Rational}(Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+
+    @testset "convert" begin
+        @test convert(Stencil{Float64}, Stencil(1,2,3,4,5; center=2)) == Stencil(1.,2.,3.,4.,5.; center=2)
+        @test convert(Stencil{Float64}, CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
+        @test convert(Stencil{Int}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
+        @test convert(Stencil{Rational}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+    end
 end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Wed Jan 19 11:08:43 2022 +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
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -7,7 +7,9 @@
 import Sbplib.SbpOperators.VolumeOperator
 
 @testset "SecondDerivative" 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"])
     Lx = 3.5
     Ly = 3.
     g_1D = EquidistantGrid(121, 0.0, Lx)
@@ -15,13 +17,13 @@
 
     @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)
+            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,op.innerStencil,op.closureStencils,1)
-            D2 = second_derivative(g_1D,op.innerStencil,op.closureStencils)
+            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
@@ -45,8 +47,10 @@
             # 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)
+                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
@@ -56,8 +60,10 @@
             # 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)
+                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
@@ -82,8 +88,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)
-                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
+                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
@@ -93,8 +101,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)
-                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
+                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
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -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_scalar(stencil_set["H"]["inner"])
+        quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+        quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Wed Jan 19 11:08:43 2022 +0100
@@ -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_scalar(stencil_set["H"]["inner"])
+        quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+        quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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_scalar(stencil_set["H"]["inner"])
+                quadrature_closure = parse_tuple(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	Wed Jan 19 07:24:36 2022 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Jan 19 11:08:43 2022 +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