Mercurial > repos > public > sbplib_julia

changeset 980:f885e1de6dc4 feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 15 Mar 2022 21:38:55 +0100
parents bc12be1b1ae5 (current diff) a52bd2f1126c (diff)
children df562695b1b5
files Notes.md TODO.md src/SbpOperators/SbpOperators.jl test/SbpOperators/boundaryops/boundary_operator_test.jl
diffstat 17 files changed, 211 insertions(+), 104 deletions(-) [+]
line wrap: on
line diff
diff -r bc12be1b1ae5 -r f885e1de6dc4 Notes.md
--- a/Notes.md	Mon Mar 14 10:20:39 2022 +0100
+++ b/Notes.md	Tue Mar 15 21:38:55 2022 +0100
@@ -147,6 +147,7 @@
  - [ ] How do we handle mixes of periodic and non-periodic grids? Seems it should be supported on the grid level and on the 1d operator level. Between there it should be transparent.
  - [ ] Can we have a trait to tell if a TensorMapping is transposable?
  - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input?
+ - [ ] Figure out how to treat the borrowing parameters of operators. Include in into the struct? Expose via function dispatched on the operator type and grid?
 
 ## Identifiers for regions
 The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probabily be included in the grid module. This allows new grid types to come with their own regions.
diff -r bc12be1b1ae5 -r f885e1de6dc4 TODO.md
--- a/TODO.md	Mon Mar 14 10:20:39 2022 +0100
+++ b/TODO.md	Tue Mar 15 21:38:55 2022 +0100
@@ -1,6 +1,7 @@
 # TODO
 
 
+## Coding
  - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
  - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?)
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
@@ -21,6 +22,8 @@
  - [ ] Add custom pretty printing to LazyTensors/SbpOperators to enhance readability of e.g error messages.
        See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing)
  - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe
+ - [ ] Move export statements to top of each module
+ - [ ] Add a type StencilSet for easier dispatch
 
 
  - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative. Kan vi använda promote istället?
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/Grids/AbstractGrid.jl
--- a/src/Grids/AbstractGrid.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/Grids/AbstractGrid.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -7,7 +7,7 @@
 
 """
 abstract type AbstractGrid end
-
+export AbstractGrid
 function dimension end
 function points end
 export dimension, points
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/SbpOperators.jl
--- a/src/SbpOperators/SbpOperators.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -1,5 +1,14 @@
 module SbpOperators
 
+export boundary_quadrature
+export boundary_restriction
+export inner_product
+export inverse_inner_product
+export Laplace
+export laplace
+export normal_derivative
+export second_derivative
+
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 using Sbplib.Grids
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/boundaryops/boundary_restriction.jl
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -1,18 +1,27 @@
+# TODO: The type parameter closure_stencil::Stencil is required since there isnt any suitable type
+# for stencil_set. We should consider adding type ::StencilSet and dispatch on that instead.
+# The same goes for other operators
 """
-    boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    boundary_restriction(grid, closure_stencil::Stencil, boundary)
 
-Creates the boundary restriction operator `e` as a `TensorMapping`
+Creates boundary restriction operators `e` as `TensorMapping`s on `boundary`
 
-`e` is the restriction of a grid function to the boundary specified by `boundary` or `region` using some `closure_stencil`.
-`e'` is the prolongation of a grid function on the boundary to the whole grid using the same `closure_stencil`.
+`e` is the restriction of a grid function to `boundary` using a `Stencil` `closure_stencil`.
+`e'` is the prolongation of a grid function on `boundary` to the whole grid using the same `closure_stencil`.
 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.
+a `BoundaryOperator`.
+
+See also: [`boundary_operator`](@ref).
 """
-function boundary_restriction(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
+function boundary_restriction(grid, closure_stencil::Stencil, boundary)
     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)}())
+
+"""
+    boundary_restriction(grid, stencil_set, boundary)
 
-export boundary_restriction
+Creates a `boundary_restriction` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+boundary_restriction(grid, stencil_set, boundary) = boundary_restriction(grid, parse_stencil(stencil_set["e"]["closure"]), boundary)
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/boundaryops/normal_derivative.jl
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -1,18 +1,25 @@
 """
-    normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    normal_derivative(grid, closure_stencil::Stencil, boundary)
 
 Creates the normal derivative boundary operator `d` as a `TensorMapping`
 
-`d` is the normal derivative of a grid function at the boundary specified by `boundary` or `region` using some `closure_stencil`.
+`d` computes the normal derivative of a grid function  on `boundary` a `Stencil` `closure_stencil`.
 `d'` is the prolongation of the normal derivative of a grid function to the whole grid using the same `closure_stencil`.
 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.
+a `BoundaryOperator`.
+
+See also: [`boundary_operator`](@ref).
 """
-function normal_derivative(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
+function normal_derivative(grid, closure_stencil::Stencil, boundary)
     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, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
-export normal_derivative
+
+"""
+    normal_derivative(grid, stencil_set, boundary)
+
+Creates a `normal_derivative` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+normal_derivative(grid, stencil_set, boundary) = normal_derivative(grid, parse_stencil(stencil_set["d1"]["closure"]), boundary)
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/volumeops/derivatives/second_derivative.jl
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -9,11 +9,23 @@
 
 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.
+
+See also: [`volume_operator`](@ref).
 """
 function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
     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
+second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
+
+"""
+    second_derivative(grid, stencil_set, direction)
+
+Creates a `second_derivative` operator on `grid` along coordinate dimension `direction` given a parsed TOML
+`stencil_set`.
+"""
+function second_derivative(grid::EquidistantGrid, stencil_set, direction)
+    inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+    second_derivative(grid,inner_stencil,closure_stencils,direction);
+end 
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/volumeops/inner_products/inner_product.jl
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -10,8 +10,10 @@
 
 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`.
+product operators for each coordinate direction. On a 0-dimensional grid,
+`H` is a 0-dimensional `IdentityMapping`.
+
+See also: [`ConstantInteriorScalingOperator`](@ref).
 """
 function inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
     Hs = ()
@@ -22,7 +24,6 @@
 
     return foldl(⊗, Hs)
 end
-export inner_product
 
 function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
     h = spacing(grid)[1]
@@ -32,3 +33,15 @@
 end
 
 inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
+
+"""
+    inner_product(grid, stencil_set)
+
+Creates a `inner_product` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+function inner_product(grid, stencil_set)
+    inner_stencil = parse_scalar(stencil_set["H"]["inner"])
+    closure_stencils = parse_tuple(stencil_set["H"]["closure"])
+    return inner_product(grid, inner_stencil, closure_stencils)
+end
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -8,7 +8,9 @@
 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`.
+`grid`, `H⁻¹` is a 0-dimensional `IdentityMapping`. 
+
+See also: [`ConstantInteriorScalingOperator`](@ref).
 """
 function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
     H⁻¹s = ()
@@ -25,6 +27,17 @@
     H⁻¹ = SbpOperators.ConstantInteriorScalingOperator(grid, h⁻¹*1/interior_weight, h⁻¹./closure_weights)
     return H⁻¹
 end
-export inverse_inner_product
 
 inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
+
+"""
+    inverse_inner_product(grid, stencil_set)
+
+Creates a `inverse_inner_product` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+function inverse_inner_product(grid, stencil_set)
+    inner_stencil = parse_scalar(stencil_set["H"]["inner"])
+    closure_stencils = parse_tuple(stencil_set["H"]["closure"])
+    return inverse_inner_product(grid, inner_stencil, closure_stencils)
+end
diff -r bc12be1b1ae5 -r f885e1de6dc4 src/SbpOperators/volumeops/laplace/laplace.jl
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -1,5 +1,37 @@
 """
-    laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
+    Laplace{T, Dim, TM} <: TensorMapping{T, Dim, Dim}
+
+Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
+`TensorMapping`. Additionally `Laplace` stores the stencil set (parsed from TOML) 
+used to construct the `TensorMapping`.
+"""
+struct Laplace{T, Dim, TM<:TensorMapping{T, Dim, Dim}} <: TensorMapping{T, Dim, Dim}
+    D::TM       # Difference operator
+    stencil_set # Stencil set of the operator
+end
+
+"""
+    Laplace(grid::Equidistant, stencil_set)
+
+Creates the `Laplace` operator `Δ` on `grid` given a parsed TOML
+`stencil_set`. See also [`laplace`](@ref).
+"""
+function Laplace(grid::EquidistantGrid, stencil_set)
+    inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+    Δ = laplace(grid, inner_stencil,closure_stencils)
+    return Laplace(Δ,stencil_set)
+end
+
+LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
+LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
+LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
+
+# TODO: Implement pretty printing of Laplace once pretty printing of TensorMappings is implemented. 
+# Base.show(io::IO, L::Laplace) = ...
+
+"""
+    laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
 
 Creates the Laplace operator operator `Δ` as a `TensorMapping`
 
@@ -10,6 +42,8 @@
 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
 where the sum is carried out lazily.
+
+See also: [`second_derivative`](@ref).
 """
 function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
     Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
@@ -18,4 +52,3 @@
     end
     return Δ
 end
-export laplace
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/boundaryops/boundary_operator_test.jl
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -32,14 +32,8 @@
             @test e_w isa TensorMapping{T,1,2} where T
         end
     end
-
-    op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}())
-    op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}())
-
-    op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}())
-    op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}())
-    op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}())
-    op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}())
+    op_l, op_r = boundary_operator.(Ref(g_1D), Ref(closure_stencil), boundary_identifiers(g_1D))
+    op_w, op_e, op_s, op_n = boundary_operator.(Ref(g_2D), Ref(closure_stencil), boundary_identifiers(g_2D))
 
     @testset "Sizes" begin
         @testset "1D" begin
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/boundaryops/boundary_restriction_test.jl
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -2,9 +2,8 @@
 
 using Sbplib.SbpOperators
 using Sbplib.Grids
+using Sbplib.LazyTensors
 using Sbplib.RegionIndices
-using Sbplib.LazyTensors
-
 import Sbplib.SbpOperators.BoundaryOperator
 
 @testset "boundary_restriction" begin
@@ -15,14 +14,14 @@
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D,e_closure,Lower())
-            @test e_l == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
+            e_l = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
+            @test e_l == boundary_restriction(g_1D,stencil_set,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,e_closure,Upper())
-            @test e_r == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
+            e_r = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
+            @test e_r == boundary_restriction(g_1D,stencil_set,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
@@ -30,6 +29,7 @@
 
         @testset "2D" begin
             e_w = boundary_restriction(g_2D,e_closure,CartesianBoundary{1,Upper}())
+            @test e_w == boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}())
             @test e_w isa InflatedTensorMapping
             @test e_w isa TensorMapping{T,1,2} where T
         end
@@ -37,9 +37,7 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Lower}())
-            e_r = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Upper}())
-
+            e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(e_closure), boundary_identifiers(g_1D))
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
 
@@ -49,11 +47,7 @@
         end
 
         @testset "2D" begin
-            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}())
-
+            e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(e_closure), boundary_identifiers(g_2D))
             v = rand(11, 15)
             u = fill(3.124)
 
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/boundaryops/normal_derivative_test.jl
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -2,9 +2,8 @@
 
 using Sbplib.SbpOperators
 using Sbplib.Grids
+using Sbplib.LazyTensors
 using Sbplib.RegionIndices
-using Sbplib.LazyTensors
-
 import Sbplib.SbpOperators.BoundaryOperator
 
 @testset "normal_derivative" begin
@@ -14,8 +13,8 @@
     	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, d_closure, Lower())
-            @test d_l == normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}())
+            @test d_l == normal_derivative(g_1D, stencil_set, CartesianBoundary{1,Lower}())
             @test d_l isa BoundaryOperator{T,Lower} where T
             @test d_l isa TensorMapping{T,0,1} where T
         end
@@ -24,8 +23,9 @@
             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),d_closure,Lower())
-            d_r = normal_derivative(restrict(g_2D,2),d_closure,Upper())
+            d_l = normal_derivative(restrict(g_2D,1),d_closure,CartesianBoundary{1,Lower}())
+            d_r = normal_derivative(restrict(g_2D,2),d_closure,CartesianBoundary{1,Upper}())
+            @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa TensorMapping{T,1,2} where T
@@ -40,10 +40,7 @@
         @testset "2nd order" begin
         	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}())
+            d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(d_closure), boundary_identifiers(g_2D))
 
             @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ v∂x[end,:] atol = 1e-13
@@ -54,11 +51,8 @@
         @testset "4th order" begin
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         	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}())
-
+            d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(d_closure), boundary_identifiers(g_2D))
+            
             @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ v∂x[end,:] atol = 1e-13
             @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/volumeops/derivatives/second_derivative_test.jl
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -7,7 +7,8 @@
 import Sbplib.SbpOperators.VolumeOperator
 
 @testset "SecondDerivative" begin
-    stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    stencil_set = read_stencil_set(operator_path; order=4)
     inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
     closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
     Lx = 3.5
@@ -17,8 +18,9 @@
 
     @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)
+            Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils,1)
+            @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils)
+            @test Dₓₓ == second_derivative(g_1D,stencil_set,1)
             @test Dₓₓ isa VolumeOperator
         end
         @testset "2D" begin
@@ -26,6 +28,7 @@
             D2 = second_derivative(g_1D,inner_stencil,closure_stencils)
             I = IdentityMapping{Float64}(size(g_2D)[2])
             @test Dₓₓ == D2⊗I
+            @test Dₓₓ == second_derivative(g_2D,stencil_set,1)
             @test Dₓₓ isa TensorMapping{T,2,2} where T
         end
     end
@@ -47,7 +50,7 @@
             # 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)
+                stencil_set = read_stencil_set(operator_path; 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)
@@ -60,7 +63,7 @@
             # 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)
+                stencil_set = read_stencil_set(operator_path; 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)
@@ -88,7 +91,7 @@
             # 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)
+                stencil_set = read_stencil_set(operator_path; 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)
@@ -101,7 +104,7 @@
             # 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)
+                stencil_set = read_stencil_set(operator_path; 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)
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/volumeops/inner_products/inner_product_test.jl
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -4,6 +4,7 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
+import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
 @testset "Diagonal-stencil inner_product" begin
     Lx = π/2.
@@ -19,18 +20,21 @@
         quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
         @testset "0D" begin
             H = inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
+            @test H == inner_product(EquidistantGrid{Float64}(), stencil_set)
             @test H == IdentityMapping{Float64}()
             @test H isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
             H = inner_product(g_1D, quadrature_interior, quadrature_closure)
-            @test H == inner_product(g_1D, quadrature_interior, quadrature_closure)
+            @test H == inner_product(g_1D, stencil_set)
+            @test H isa ConstantInteriorScalingOperator
             @test H isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
             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 == inner_product(g_2D, stencil_set)
             @test H == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -4,7 +4,7 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
-import Sbplib.SbpOperators.Stencil
+import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
 @testset "Diagonal-stencil inverse_inner_product" begin
     Lx = π/2.
@@ -17,17 +17,21 @@
         quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
         @testset "0D" begin
             Hi = inverse_inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
+            @test Hi == inverse_inner_product(EquidistantGrid{Float64}(), stencil_set)
             @test Hi == IdentityMapping{Float64}()
             @test Hi isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
             Hi = inverse_inner_product(g_1D,  quadrature_interior, quadrature_closure)
+            @test Hi == inverse_inner_product(g_1D, stencil_set)
+            @test Hi isa ConstantInteriorScalingOperator
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
             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 == inverse_inner_product(g_2D, stencil_set)
             @test Hi == Hi_x⊗Hi_y
             @test Hi isa TensorMapping{T,2,2} where T
         end
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/SbpOperators/volumeops/laplace/laplace_test.jl
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Mon Mar 14 10:20:39 2022 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Mar 15 21:38:55 2022 +0100
@@ -4,25 +4,25 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
+# Default stencils (4th order)
+operator_path = sbp_operators_path()*"standard_diagonal.toml"
+stencil_set = read_stencil_set(operator_path; order=4)
+inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+g_1D = EquidistantGrid(101, 0.0, 1.)
+g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+
 @testset "Laplace" begin
-    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
-        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, inner_stencil, closure_stencils)
-            @test L == second_derivative(g_1D, inner_stencil, closure_stencils)
-            @test L isa TensorMapping{T,1,1}  where T
+            Δ = laplace(g_1D, inner_stencil, closure_stencils)
+            @test Laplace(g_1D, stencil_set) == Laplace(Δ, stencil_set)
+            @test Laplace(g_1D, stencil_set) isa TensorMapping{T,1,1}  where T
         end
         @testset "3D" begin
-            L = laplace(g_3D, inner_stencil, closure_stencils)
-            @test L isa TensorMapping{T,3,3} where T
-            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
+            Δ = laplace(g_3D, inner_stencil, closure_stencils)
+            @test Laplace(g_3D, stencil_set) == Laplace(Δ,stencil_set)
+            @test Laplace(g_3D, stencil_set) isa TensorMapping{T,3,3} where T
         end
     end
 
@@ -42,30 +42,44 @@
         # 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"])
-            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
-            @test L*v ≈ Δv rtol = 5e-2 norm = l2
+            stencil_set = read_stencil_set(operator_path; order=2)
+            Δ = Laplace(g_3D, stencil_set)
+            @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
+            @test Δ*v ≈ Δv 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"])
-            L = laplace(g_3D, inner_stencil, closure_stencils)
+            stencil_set = read_stencil_set(operator_path; order=4)
+            Δ = Laplace(g_3D, stencil_set)
             # 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
-            @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
-            @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
-            @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9
-            @test L*v ≈ Δv rtol = 5e-4 norm = l2
+            @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
+            @test Δ*polynomials[4] ≈ polynomials[2] atol = 5e-9
+            @test Δ*v ≈ Δv rtol = 5e-4 norm = l2
         end
     end
 end
+
+@testset "laplace" begin
+    @testset "1D" begin
+        Δ = laplace(g_1D, inner_stencil, closure_stencils)
+        @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils)
+        @test Δ isa TensorMapping{T,1,1}  where T
+    end
+    @testset "3D" begin
+        Δ = laplace(g_3D, inner_stencil, closure_stencils)
+        @test Δ isa TensorMapping{T,3,3} where T
+        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 Δ == Dxx + Dyy + Dzz
+        @test Δ isa TensorMapping{T,3,3} where T
+    end
+end
+
diff -r bc12be1b1ae5 -r f885e1de6dc4 test/runtests.jl