Mercurial > repos > public > sbplib_julia

changeset 1221:b3b4d29b46c3 refactor/grids

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 Fri, 10 Feb 2023 08:36:56 +0100
parents 93bba649aea2 (current diff) 7ee258e5289e (diff)
children 5f677cd6f0b6
files
diffstat 25 files changed, 882 insertions(+), 321 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Tue Nov 01 22:44:00 2022 +0100
+++ b/Manifest.toml	Fri Feb 10 08:36:56 2023 +0100
@@ -1,40 +1,40 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.8.2"
+julia_version = "1.8.5"
 manifest_format = "2.0"
 project_hash = "b024d6898b484706c36ee3b2a041918f3a9d2088"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "195c5505521008abea5aee4f96930717958eac6f"
+git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.4.0"
+version = "3.5.0"
 
 [[deps.ArrayInterface]]
-deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "Static"]
-git-tree-sha1 = "d6173480145eb632d6571c148d94b9d3d773820e"
+deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"]
+git-tree-sha1 = "dedc16cbdd1d32bead4617d27572f582216ccf23"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "6.0.23"
+version = "6.0.25"
 
 [[deps.ArrayInterfaceCore]]
-deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "5bb0f8292405a516880a3809954cb832ae7a31c5"
+deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d"
 uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2"
-version = "0.1.20"
+version = "0.1.29"
 
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
 [[deps.Compat]]
 deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "5856d3031cdb1f3b2b6340dfdc66b6d9a149a374"
+git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.2.0"
+version = "4.6.0"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "0.5.2+0"
+version = "1.0.1+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -54,15 +54,21 @@
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "1ea784113a6aa054c5ebd95945fa5e52c2f378e7"
+git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.12.7"
+version = "1.12.9"
 
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 version = "0.3.20+0"
 
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
 [[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -78,15 +84,21 @@
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[deps.SnoopPrecompile]]
+deps = ["Preferences"]
+git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
+uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+version = "1.0.3"
+
 [[deps.SparseArrays]]
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "de4f0a4f049a4c87e4948c04acff37baf1be01a6"
+git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.7.7"
+version = "0.8.3"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
@@ -99,9 +111,9 @@
 
 [[deps.TiledIteration]]
 deps = ["ArrayInterface", "OffsetArrays"]
-git-tree-sha1 = "5e02b75701f1905e55e44fc788bd13caedb5a6e3"
+git-tree-sha1 = "1bf2bb587a7fc99fefac2ff076b18b500128e9c0"
 uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
-version = "0.4.1"
+version = "0.4.2"
 
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/WORKFLOW.md	Fri Feb 10 08:36:56 2023 +0100
@@ -0,0 +1,46 @@
+# Branch model
+
+The `default` branch contains stable code and the goal is that all tests should be passing on this branch. Version of the code are marked with mercurial tags. Changes in the code are developed in named branches which are closed and merged into `default` when the code is ready. During development merging `default` into the development branch is encouraged to avoid complicated conflicts.
+
+Branches are named using slash separated keywords. The first keyword describes the type of change being pursued in the branch. Important type keywords are
+ * feature
+ * bugfix
+ * refactor
+Further keywords may describe where, e.g. what sub package, the change happens. The last keyword should describe the change.
+
+Some examples:
+ * refactor/grids: Branch to refactor the grids module
+ * bugfix/lazy_tensors/lazyfunctionarray: Branch to fix a bug in LazyFunctionArray
+
+## Merging a branch into `default`
+The changes in a branch has been reviewed and deemed ready to merge the branch is closed and then merged.
+
+Before merging a development branch, `default` should be merge into the development branch to make sure the whole state of the code is reviewed and tested before it ends up on `default`.
+
+With the development branch active the following commands can be used to complete the merging of a development branch.
+```shell
+hg merge default
+hg commit --close-branch -m "Close before merge"
+hg update default
+hg merge development/branch
+hg commit -m "Merge development/branch"
+```
+
+# Review
+
+## Checklist for review
+* Push and pull new changes
+* Search and check TODOs
+* Search and check TBDs
+* Search and check REVIEWs
+* Review code
+* Review tests
+* Review docstrings
+* Render Documenter and check docstrings in browser
+* Run full tests
+
+# Special comments
+The following special comments are used:
+* `# TODO: `: Something that should be done at some point.
+* `# TBD: `:  "To be determined", i.e a decision that has to be made.
+* `# REVIEW: `: A review comment. Should only exist on development branches.
--- a/docs/Manifest.toml	Tue Nov 01 22:44:00 2022 +0100
+++ b/docs/Manifest.toml	Fri Feb 10 08:36:56 2023 +0100
@@ -1,7 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.1"
+julia_version = "1.8.5"
 manifest_format = "2.0"
+project_hash = "4f0756199bb5f6739a5f4697152617efc4e0705c"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -10,9 +11,21 @@
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
+git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.3.3"
+version = "3.5.0"
+
+[[deps.ArrayInterface]]
+deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"]
+git-tree-sha1 = "dedc16cbdd1d32bead4617d27572f582216ccf23"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "6.0.25"
+
+[[deps.ArrayInterfaceCore]]
+deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d"
+uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2"
+version = "0.1.29"
 
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -20,9 +33,16 @@
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[deps.Compat]]
+deps = ["Dates", "LinearAlgebra", "UUIDs"]
+git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.6.0"
+
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.0.1+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -30,15 +50,15 @@
 
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
-git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.6"
+version = "0.9.3"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "f425293f7e0acaf9144de6d731772de156676233"
+git-tree-sha1 = "58fea7c536acd71f3eef6be3b21c0df5f3df88fd"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.10"
+version = "0.27.24"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
@@ -46,15 +66,20 @@
 uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
 version = "0.2.2"
 
+[[deps.IfElse]]
+git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
+uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+version = "0.1.1"
+
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
+git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.2"
+version = "0.21.3"
 
 [[deps.LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
@@ -79,22 +104,30 @@
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
+git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.10.8"
+version = "1.12.9"
 
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.20+0"
 
 [[deps.Parsers]]
-deps = ["Dates"]
-git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918"
+deps = ["Dates", "SnoopPrecompile"]
+git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.1.3"
+version = "2.5.4"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -110,6 +143,7 @@
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
 
 [[deps.Sbplib]]
 deps = ["TOML", "TiledIteration"]
@@ -120,22 +154,47 @@
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[deps.SnoopPrecompile]]
+deps = ["Preferences"]
+git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
+uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+version = "1.0.3"
+
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.Static]]
+deps = ["IfElse"]
+git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c"
+uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+version = "0.8.3"
+
+[[deps.SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
 [[deps.TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
 
 [[deps.Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[deps.TiledIteration]]
-deps = ["OffsetArrays"]
-git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297"
+deps = ["ArrayInterface", "OffsetArrays"]
+git-tree-sha1 = "1bf2bb587a7fc99fefac2ff076b18b500128e9c0"
 uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
-version = "0.3.1"
+version = "0.4.2"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
@@ -143,3 +202,4 @@
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.1.1+0"
--- a/src/LazyTensors/lazy_tensor_operations.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -269,6 +269,29 @@
 LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms)
 
 
+
+"""
+    inflate(tm::LazyTensor, sz, dir)
+
+Inflate `tm` such that it gets the size `sz` in all directions except `dir`.
+Here `sz[dir]` is ignored and replaced with the range and domains size of
+`tm`.
+
+An example of when this operation is useful is when extending a one
+dimensional difference operator `D` to a 2D grid of a ceratin size. In that
+case we could have
+
+```julia
+Dx = inflate(D, (10,10), 1)
+Dy = inflate(D, (10,10), 2)
+```
+"""
+function inflate(tm::LazyTensor, sz, dir)
+    Is = IdentityTensor{eltype(tm)}.(sz)
+    parts = Base.setindex(Is, tm, dir)
+    return foldl(⊗, parts)
+end
+
 function check_domain_size(tm::LazyTensor, sz)
     if domain_size(tm) != sz
         throw(DomainSizeMismatch(tm,sz))
--- a/src/LazyTensors/tuple_manipulation.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/LazyTensors/tuple_manipulation.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -74,3 +74,32 @@
 flatten_tuple(t::NTuple{N, Number} where N) = t
 flatten_tuple(t::Tuple) = ((flatten_tuple.(t)...)...,) # simplify?
 flatten_tuple(ts::Vararg) = flatten_tuple(ts)
+
+"""
+    left_pad_tuple(t, val, N)
+
+Left pad the tuple `t` to length `N` using the value `val`.
+"""
+function left_pad_tuple(t, val, N)
+    if N < length(t)
+        throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements"))
+    end
+
+    padding = ntuple(i->val, N-length(t))
+    return (padding..., t...)
+end
+
+"""
+    right_pad_tuple(t, val, N)
+
+Right pad the tuple `t` to length `N` using the value `val`.
+"""
+function right_pad_tuple(t, val, N)
+    if N < length(t)
+        throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements"))
+    end
+
+    padding = ntuple(i->val, N-length(t))
+    return (t..., padding...)
+end
+
--- a/src/SbpOperators/SbpOperators.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -19,6 +19,7 @@
 export normal_derivative
 export first_derivative
 export second_derivative
+export undivided_skewed04
 
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
@@ -32,9 +33,11 @@
 include("stencil.jl")
 include("stencil_set.jl")
 include("volumeops/volume_operator.jl")
+include("volumeops/stencil_operator_distinct_closures.jl")
 include("volumeops/constant_interior_scaling_operator.jl")
 include("volumeops/derivatives/first_derivative.jl")
 include("volumeops/derivatives/second_derivative.jl")
+include("volumeops/derivatives/dissipation.jl")
 include("volumeops/laplace/laplace.jl")
 include("volumeops/inner_products/inner_product.jl")
 include("volumeops/inner_products/inverse_inner_product.jl")
--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -1,32 +1,3 @@
-"""
-    boundary_operator(grid,closure_stencil,boundary)
-
-Creates a boundary operator on a `Dim`-dimensional grid for the
-specified `boundary`. The action of the operator is determined by `closure_stencil`.
-
-When `Dim=1`, the corresponding `BoundaryOperator` tensor mapping is returned.
-When `Dim>1`, the `BoundaryOperator` `op` is inflated by the outer product
-of `IdentityTensors` in orthogonal coordinate directions, e.g for `Dim=3`,
-the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
-"""
-function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
-    #TODO:Check that dim(boundary) <= Dim?
-
-    # Create 1D boundary operator
-    r = region(boundary)
-    d = dim(boundary)
-    op = BoundaryOperator(restrict(grid, d), closure_stencil, r)
-
-    # Create 1D IdentityTensors for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(dims(grid)))
-    Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids))
-
-    # Formulate the correct outer product sequence of the identity mappings and
-    # the boundary operator
-    parts = Base.setindex(Is, op, d)
-    return foldl(⊗, parts)
-end
-
 """
     BoundaryOperator{T,R,N} <: LazyTensor{T,0,1}
 
@@ -41,8 +12,6 @@
     size::Int
 end
 
-BoundaryOperator{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryOperator{T,R,N}(stencil, size)
-
 """
     BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region)
 
@@ -55,6 +24,7 @@
 
 """
     closure_size(::BoundaryOperator)
+
 The size of the closure stencil.
 """
 closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -8,11 +8,13 @@
 On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
 a `BoundaryOperator`.
 
-See also: [`boundary_operator`](@ref).
+See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function boundary_restriction(grid, closure_stencil, boundary)
     converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil)
-    return SbpOperators.boundary_operator(grid, converted_stencil, boundary)
+
+    op = BoundaryOperator(restrict(grid, dim(boundary)), converted_stencil, region(boundary))
+    return LazyTensors.inflate(op, size(grid), dim(boundary))
 end
 
 """
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -8,12 +8,14 @@
 On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
 a `BoundaryOperator`.
 
-See also: [`boundary_operator`](@ref).
+See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function normal_derivative(grid, closure_stencil, boundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
-    return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
+
+    op = BoundaryOperator(restrict(grid, dim(boundary)), scale(closure_stencil,h_inv), region(boundary))
+    return LazyTensors.inflate(op, size(grid), dim(boundary))
 end
 
 """
--- a/src/SbpOperators/stencil.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/stencil.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -85,6 +85,21 @@
     return w
 end
 
+function left_pad(s::Stencil, N)
+    weights = LazyTensors.left_pad_tuple(s.weights, zero(eltype(s)), N)
+    range = (first(s.range) - (N - length(s.weights))):last(s.range)
+
+    return Stencil(range, weights)
+end
+
+function right_pad(s::Stencil, N)
+    weights = LazyTensors.right_pad_tuple(s.weights, zero(eltype(s)), N)
+    range = first(s.range):(last(s.range) + (N - length(s.weights)))
+
+    return Stencil(range, weights)
+end
+
+
 
 struct NestedStencil{T,N,M}
     s::Stencil{Stencil{T,N},M}
--- a/src/SbpOperators/stencil_set.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/stencil_set.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -14,7 +14,7 @@
 
 
 """
-read_stencil_set(filename; filters)
+    read_stencil_set(filename; filters)
 
 Creates a `StencilSet` from a TOML file based on some key-value
 filters. If more than one set matches the filters an error is raised. The
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -0,0 +1,107 @@
+"""
+    undivided_skewed04(g::EquidistantGrid, p, direction)
+
+Create undivided difference operators approximating the `p`th derivative. The
+operators do not satisfy any SBP-property and are meant to be used for
+building artificial dissipation terms.
+
+The operators and how they are used to create accurate artifical dissipation
+is described in "K. Mattsson, M. Svärd, and J. Nordström, “Stable and Accurate
+Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp.
+57–79, Aug. 2004"
+"""
+function undivided_skewed04(g::EquidistantGrid, p, direction)
+    T = eltype(g)
+    interior_weights = T.(dissipation_interior_weights(p))
+
+    D  = stencil_operator_distinct_closures(
+        g,
+        dissipation_interior_stencil(interior_weights),
+        dissipation_lower_closure_stencils(interior_weights),
+        dissipation_upper_closure_stencils(interior_weights),
+        direction,
+    )
+    Dᵀ = stencil_operator_distinct_closures(
+        g,
+        dissipation_transpose_interior_stencil(interior_weights),
+        dissipation_transpose_lower_closure_stencils(interior_weights),
+        dissipation_transpose_upper_closure_stencils(interior_weights),
+        direction,
+    )
+
+    return D, Dᵀ
+end
+
+undivided_skewed04(g::EquidistantGrid{1}, p) = undivided_skewed04(g, p, 1)
+
+function dissipation_interior_weights(p)
+   if p == 0
+       return (1,)
+   end
+
+   return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0)
+end
+
+midpoint(weights) = length(weights)÷2 + 1
+midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
+
+function dissipation_interior_stencil(weights)
+    return Stencil(weights..., center=midpoint(weights))
+end
+function dissipation_transpose_interior_stencil(weights)
+    if iseven(length(weights))
+        weights = map(-, weights)
+    end
+
+    return Stencil(weights..., center=midpoint_transpose(weights))
+end
+
+dissipation_lower_closure_size(weights) = midpoint(weights) - 1
+dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
+
+function dissipation_lower_closure_stencils(interior_weights)
+    stencil(i) = Stencil(interior_weights..., center=i)
+    return ntuple(i->stencil(i), dissipation_lower_closure_size(interior_weights))
+end
+
+function dissipation_upper_closure_stencils(interior_weights)
+    center(i) = length(interior_weights) - dissipation_upper_closure_size(interior_weights) + i
+    stencil(i) = Stencil(interior_weights..., center=center(i))
+    return ntuple(i->stencil(i), dissipation_upper_closure_size(interior_weights))
+end
+
+function dissipation_transpose_lower_closure_stencils(interior_weights)
+    closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
+
+    N = maximum(s->length(s.weights), closure)
+    return right_pad.(closure, N)
+end
+
+function dissipation_transpose_upper_closure_stencils(interior_weights)
+    closure = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights)))
+
+    N = maximum(s->length(s.weights), closure)
+    return left_pad.(closure, N)
+end
+
+
+function dissipation_transpose_lower_closure_stencil(interior_weights, i)
+    w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights))
+
+    for k ∈ i:-1:1
+        w = (w..., interior_weights[k])
+    end
+
+    return Stencil(w..., center = i)
+end
+
+function dissipation_transpose_upper_closure_stencil(interior_weights, i)
+    j = length(interior_weights)+1-i
+    w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights))
+
+    for k ∈ j:1:length(interior_weights)
+        w = (interior_weights[k], w...)
+    end
+
+    return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1)
+end
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -7,14 +7,16 @@
 `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
 for the points in the closure regions.
 
-On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the outer product of the
-one-dimensional operator with the `IdentityTensor`s in orthogonal coordinate dirrections.
+On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the inflation of
+a `VolumeOperator`.
 
-See also: [`volume_operator`](@ref).
+See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
     h_inv = inverse_spacing(grid)[direction]
-    return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction)
+
+    D₁ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd)
+    return LazyTensors.inflate(D₁, size(grid), direction)
 end
 
 
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -7,14 +7,16 @@
 `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
 for the points in the closure regions.
 
-On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the outer product of the
-one-dimensional operator with the `IdentityTensor`s in orthogonal coordinate dirrections.
+On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the inflation of
+a `VolumeOperator`.
 
-See also: [`volume_operator`](@ref).
+See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@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)
+
+    D₂ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even)
+    return LazyTensors.inflate(D₂, size(grid), direction)
 end
 
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -0,0 +1,72 @@
+"""
+    stencil_operator_distinct_closures(
+        grid::EquidistantGrid,
+        inner_stencil,
+        lower_closure,
+        upper_closure,
+        direction
+    )
+
+Creates a multi-dimensional `StencilOperatorDistinctClosures` acting on grid
+functions of `grid`.
+
+See also: [`StencilOperatorDistinctClosures`](@ref)
+"""
+function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction)
+    op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure)
+    return LazyTensors.inflate(op, size(grid), direction)
+end
+
+"""
+    StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1}
+
+A one dimensional stencil operator with separate closures for the two boundaries.
+
+`StencilOperatorDistinctClosures` can be contrasted to `VolumeOperator` in
+that it has different closure stencils for the upper and lower boundary.
+`VolumeOperator` uses the same closure for both boundaries. Having distinct
+closures is useful for representing operators with skewed stencils like upwind
+operators.
+
+See also: [`VolumeOperator`](@ref), [`stencil_operator_distinct_closures`](@ref)
+"""
+struct StencilOperatorDistinctClosures{T,K,N,M,LC<:NTuple{N,Stencil{T,L}} where L, UC<:NTuple{M,Stencil{T,L}} where L} <: LazyTensor{T,1,1}
+    inner_stencil::Stencil{T,K}
+    lower_closure::LC
+    upper_closure::UC
+    size::Tuple{Int}
+end
+
+function StencilOperatorDistinctClosures(grid::EquidistantGrid{1}, inner_stencil, lower_closure, upper_closure)
+    return StencilOperatorDistinctClosures(inner_stencil, Tuple(lower_closure), Tuple(upper_closure), size(grid))
+end
+
+lower_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = N
+upper_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = M
+
+LazyTensors.range_size(op::StencilOperatorDistinctClosures) = op.size
+LazyTensors.domain_size(op::StencilOperatorDistinctClosures) = op.size
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Lower})
+    return @inbounds apply_stencil(op.lower_closure[Int(i)], v, Int(i))
+end
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Interior})
+    return apply_stencil(op.inner_stencil, v, Int(i))
+end
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Upper})
+    stencil_index = Int(i) - (op.size[1]-upper_closure_size(op))
+    return @inbounds apply_stencil(op.upper_closure[stencil_index], v, Int(i))
+end
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i)
+    if i <= lower_closure_size(op)
+        LazyTensors.apply(op, v, Index(i, Lower))
+    elseif i > op.size[1]-upper_closure_size(op)
+        LazyTensors.apply(op, v, Index(i, Upper))
+    else
+        LazyTensors.apply(op, v, Index(i, Interior))
+    end
+end
+# TODO: Move this to LazyTensors when we have the region communication down.
--- a/src/SbpOperators/volumeops/volume_operator.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -1,30 +1,6 @@
-"""
-    volume_operator(grid, inner_stencil, closure_stencils, parity, direction)
-
-Creates a volume operator on a `Dim`-dimensional grid acting along the
-specified coordinate `direction`. The action of the operator is determined by
-the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the
-corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the
-returned operator is the appropriate outer product of a one-dimensional
-operators and `IdentityTensor`s, e.g for `Dim=3` the volume operator in the
-y-direction is `I⊗op⊗I`.
-"""
-function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction)
-    #TODO: Check that direction <= Dim?
-
-    # Create 1D volume operator in along coordinate direction
-    op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity)
-    # Create 1D IdentityTensors for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(dims(grid)))
-    Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids))
-    # Formulate the correct outer product sequence of the identity mappings and
-    # the volume operator
-    parts = Base.setindex(Is, op, direction)
-    return foldl(⊗, parts)
-end
-
 """
     VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
+
 Implements a one-dimensional constant coefficients volume operator
 """
 struct VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
@@ -59,3 +35,4 @@
     r = getregion(i, closure_size(op), op.size[1])
     return LazyTensors.apply(op, v, Index(i, r))
 end
+# TODO: Move this to LazyTensors when we have the region communication down.
--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -366,3 +366,17 @@
         @test I1⊗Ã⊗I2 == InflatedTensor(I1, Ã, I2)
     end
 end
+
+@testset "inflate" begin
+    I = LazyTensors.inflate(IdentityTensor(),(3,4,5,6), 2)
+    @test I isa LazyTensor{Float64, 3,3}
+    @test range_size(I) == (3,5,6)
+    @test domain_size(I) == (3,5,6)
+
+    @test LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 1) == InflatedTensor(IdentityTensor{Float64}(),ScalingTensor(1., (4,)),IdentityTensor(4,5,6))
+    @test LazyTensors.inflate(ScalingTensor(2., (1,)),(3,4,5,6), 2) == InflatedTensor(IdentityTensor(3),ScalingTensor(2., (1,)),IdentityTensor(5,6))
+    @test LazyTensors.inflate(ScalingTensor(3., (6,)),(3,4,5,6), 4) == InflatedTensor(IdentityTensor(3,4,5),ScalingTensor(3., (6,)),IdentityTensor{Float64}())
+
+    @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 0)
+    @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 5)
+end
--- a/test/LazyTensors/tuple_manipulation_test.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/LazyTensors/tuple_manipulation_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -60,3 +60,19 @@
     @test LazyTensors.flatten_tuple((1,2,(3,(4,5)),6)) == (1,2,3,4,5,6)
     @test LazyTensors.flatten_tuple(((1,2),(3,4),(5,),6)) == (1,2,3,4,5,6)
 end
+
+@testset "left_pad_tuple" begin
+    @test LazyTensors.left_pad_tuple((1,2), 0, 2) == (1,2)
+    @test LazyTensors.left_pad_tuple((1,2), 0, 3) == (0,1,2)
+    @test LazyTensors.left_pad_tuple((3,2), 1, 6) == (1,1,1,1,3,2)
+
+    @test_throws DomainError(0, "Can't pad tuple of length 2 to 0 elements") LazyTensors.left_pad_tuple((1,2), 0, 0) == (1,2)
+end
+
+@testset "right_pad_tuple" begin
+    @test LazyTensors.right_pad_tuple((1,2), 0, 2) == (1,2)
+    @test LazyTensors.right_pad_tuple((1,2), 0, 3) == (1,2,0)
+    @test LazyTensors.right_pad_tuple((3,2), 1, 6) == (3,2,1,1,1,1)
+
+    @test_throws DomainError(0, "Can't pad tuple of length 2 to 0 elements") LazyTensors.right_pad_tuple((1,2), 0, 0) == (1,2)
+end
--- a/test/Manifest.toml	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/Manifest.toml	Fri Feb 10 08:36:56 2023 +0100
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.8.2"
+julia_version = "1.8.5"
 manifest_format = "2.0"
 project_hash = "23260eda65ade7d11fffed313a68520d0bc053fc"
 
@@ -16,32 +16,32 @@
 
 [[deps.BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
+git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.3.1"
+version = "1.3.2"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "e7ff6cadf743c098e08fca25c91103ee4303c9bb"
+git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.15.6"
+version = "1.15.7"
 
 [[deps.ChangesOfVariables]]
 deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
-git-tree-sha1 = "38f7a08f19d8810338d4f5085211c7dfa5d5bdd8"
+git-tree-sha1 = "844b061c104c408b24537482469400af6075aae4"
 uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
-version = "0.1.4"
+version = "0.1.5"
 
 [[deps.Compat]]
 deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "5856d3031cdb1f3b2b6340dfdc66b6d9a149a374"
+git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.2.0"
+version = "4.6.0"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "0.5.2+0"
+version = "1.0.1+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -54,9 +54,9 @@
 
 [[deps.DiffRules]]
 deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "992a23afdb109d0d2f8802a30cf5ae4b1fe7ea68"
+git-tree-sha1 = "c5b6685d53f933c11404a3ae9822afe30d522494"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.11.1"
+version = "1.12.2"
 
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
@@ -64,9 +64,9 @@
 
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
-git-tree-sha1 = "5158c2b41018c5f7eb1470d558127ac274eca0c9"
+git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.9.1"
+version = "0.9.3"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -136,9 +136,9 @@
 
 [[deps.LogExpFunctions]]
 deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "94d9c52ca447e23eac0c0f074effbcd38830deb5"
+git-tree-sha1 = "45b288af6956e67e621c5cbb2d75a261ab58300b"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.18"
+version = "0.3.20"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -186,10 +186,10 @@
 version = "0.5.5+0"
 
 [[deps.Parsers]]
-deps = ["Dates"]
-git-tree-sha1 = "3d5bf43e3e8b412656404ed9466f1dcbf7c50269"
+deps = ["Dates", "SnoopPrecompile"]
+git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.4.0"
+version = "2.5.4"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -231,6 +231,12 @@
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[deps.SnoopPrecompile]]
+deps = ["Preferences"]
+git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
+uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+version = "1.0.3"
+
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -6,7 +6,7 @@
 using Sbplib.RegionIndices
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.BoundaryOperator
-import Sbplib.SbpOperators.boundary_operator
+
 
 @testset "BoundaryOperator" begin
     closure_stencil = Stencil(2.,1.,3.; center = 1)
@@ -14,120 +14,50 @@
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "Constructors" begin
-        @testset "1D" begin
-            op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1])
-            @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower())
-            @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}())
-            @test op_l isa LazyTensor{T,0,1} where T
+        @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T
+        @test BoundaryOperator(g_1D, closure_stencil, Upper()) isa LazyTensor{T,0,1} where T
+    end
 
-            op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
-            @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper())
-            @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
-            @test op_r isa LazyTensor{T,0,1} where T
-        end
-
-        @testset "2D" begin
-            e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}())
-            @test e_w isa InflatedTensor
-            @test e_w isa LazyTensor{T,1,2} where T
-        end
-    end
-    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))
+    op_l = BoundaryOperator(g_1D, closure_stencil, Lower())
+    op_r = BoundaryOperator(g_1D, closure_stencil, Upper())
 
     @testset "Sizes" begin
-        @testset "1D" begin
-            @test domain_size(op_l) == (11,)
-            @test domain_size(op_r) == (11,)
-
-            @test range_size(op_l) == ()
-            @test range_size(op_r) == ()
-        end
+        @test domain_size(op_l) == (11,)
+        @test domain_size(op_r) == (11,)
 
-        @testset "2D" begin
-            @test domain_size(op_w) == (11,15)
-            @test domain_size(op_e) == (11,15)
-            @test domain_size(op_s) == (11,15)
-            @test domain_size(op_n) == (11,15)
-
-            @test range_size(op_w) == (15,)
-            @test range_size(op_e) == (15,)
-            @test range_size(op_s) == (11,)
-            @test range_size(op_n) == (11,)
-        end
+        @test range_size(op_l) == ()
+        @test range_size(op_r) == ()
     end
 
     @testset "Application" begin
-        @testset "1D" begin
-            v = evalOn(g_1D,x->1+x^2)
-            u = fill(3.124)
-            @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3]
-            @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2]
-            @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2]
-            @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
-            @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
-
-            v = evalOn(g_1D, x->1. +x*im)
-            @test (op_l*v)[] isa ComplexF64
+        v = evalOn(g_1D,x->1+x^2)
+        u = fill(3.124)
+        @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3]
+        @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2]
+        @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2]
+        @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
+        @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
 
-            u = fill(1. +im)
-            @test (op_l'*u)[1] isa ComplexF64
-            @test (op_l'*u)[5] isa ComplexF64
-            @test (op_l'*u)[11] isa ComplexF64
-        end
-
-        @testset "2D" begin
-            v = rand(size(g_2D)...)
-            u = fill(3.124)
-            @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14
-            @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14
-            @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14
-            @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14
-
-
-            g_x = rand(size(g_2D)[1])
-            g_y = rand(size(g_2D)[2])
-
-            G_w = zeros(Float64, size(g_2D)...)
-            G_w[1,:] = 2*g_y
-            G_w[2,:] = g_y
-            G_w[3,:] = 3*g_y
+        v = evalOn(g_1D, x->1. +x*im)
+        @test (op_l*v)[] isa ComplexF64
 
-            G_e = zeros(Float64, size(g_2D)...)
-            G_e[end,:] = 2*g_y
-            G_e[end-1,:] = g_y
-            G_e[end-2,:] = 3*g_y
-
-            G_s = zeros(Float64, size(g_2D)...)
-            G_s[:,1] = 2*g_x
-            G_s[:,2] = g_x
-            G_s[:,3] = 3*g_x
-
-            G_n = zeros(Float64, size(g_2D)...)
-            G_n[:,end] = 2*g_x
-            G_n[:,end-1] = g_x
-            G_n[:,end-2] = 3*g_x
+        u = fill(1. +im)
+        @test (op_l'*u)[1] isa ComplexF64
+        @test (op_l'*u)[5] isa ComplexF64
+        @test (op_l'*u)[11] isa ComplexF64
 
-            @test op_w'*g_y == G_w
-            @test op_e'*g_y == G_e
-            @test op_s'*g_x == G_s
-            @test op_n'*g_x == G_n
-       end
+        u = fill(3.124)
+        @test (op_l'*u)[Index(1,Lower)] == 2*u[]
+        @test (op_l'*u)[Index(2,Lower)] == u[]
+        @test (op_l'*u)[Index(6,Interior)] == 0
+        @test (op_l'*u)[Index(10,Upper)] == 0
+        @test (op_l'*u)[Index(11,Upper)] == 0
 
-       @testset "Regions" begin
-            u = fill(3.124)
-            @test (op_l'*u)[Index(1,Lower)] == 2*u[]
-            @test (op_l'*u)[Index(2,Lower)] == u[]
-            @test (op_l'*u)[Index(6,Interior)] == 0
-            @test (op_l'*u)[Index(10,Upper)] == 0
-            @test (op_l'*u)[Index(11,Upper)] == 0
-
-            @test (op_r'*u)[Index(1,Lower)] == 0
-            @test (op_r'*u)[Index(2,Lower)] == 0
-            @test (op_r'*u)[Index(6,Interior)] == 0
-            @test (op_r'*u)[Index(10,Upper)] == u[]
-            @test (op_r'*u)[Index(11,Upper)] == 2*u[]
-       end
+        @test (op_r'*u)[Index(1,Lower)] == 0
+        @test (op_r'*u)[Index(2,Lower)] == 0
+        @test (op_r'*u)[Index(6,Interior)] == 0
+        @test (op_r'*u)[Index(10,Upper)] == u[]
+        @test (op_r'*u)[Index(11,Upper)] == 2*u[]
     end
 
     @testset "Inferred" begin
--- a/test/SbpOperators/stencil_test.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/SbpOperators/stencil_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -62,6 +62,23 @@
     end
 end
 
+@testset "left_pad" begin
+    @test SbpOperators.left_pad(Stencil(1,1, center = 1), 2) == Stencil(1,1, center=1)
+    @test SbpOperators.left_pad(Stencil(1,1, center = 1), 3) == Stencil(0,1,1, center=2)
+    @test SbpOperators.left_pad(Stencil(2,3, center = 2), 4) == Stencil(0,0,2,3, center=4)
+
+    @test SbpOperators.left_pad(Stencil(2.,3., center = 2), 4) == Stencil(0.,0.,2.,3., center=4)
+end
+
+@testset "right_pad" begin
+    @test SbpOperators.right_pad(Stencil(1,1, center = 1), 2) == Stencil(1,1, center=1)
+    @test SbpOperators.right_pad(Stencil(1,1, center = 1), 3) == Stencil(1,1,0, center=1)
+    @test SbpOperators.right_pad(Stencil(2,3, center = 2), 4) == Stencil(2,3,0,0, center=2)
+
+    @test SbpOperators.right_pad(Stencil(2.,3., center = 2), 4) == Stencil(2.,3.,0.,0., center=2)
+end
+
+
 @testset "NestedStencil" begin
 
     @testset "Constructors" begin
@@ -170,5 +187,4 @@
         @inferred SbpOperators.apply_stencil_backwards(s_int,   c_float, v_float, 2)
         @inferred SbpOperators.apply_stencil_backwards(s_float, c_float, v_int,   2)
     end
-
 end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -0,0 +1,217 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+using Sbplib.SbpOperators: Stencil
+
+using Sbplib.SbpOperators: dissipation_interior_weights
+using Sbplib.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil
+using Sbplib.SbpOperators: midpoint, midpoint_transpose
+using Sbplib.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size
+using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
+using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
+
+"""
+    monomial(x,k)
+
+Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
+Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
+"""
+function monomial(x,k)
+    if k < 0
+        return zero(x)
+    end
+    x^k/factorial(k)
+end
+
+@testset "undivided_skewed04" begin
+    g = EquidistantGrid(20, 0., 11.)
+    D,Dᵀ = undivided_skewed04(g, 1)
+
+    @test D isa LazyTensor{Float64,1,1}
+    @test Dᵀ isa LazyTensor{Float64,1,1}
+
+     @testset "Accuracy conditions" begin
+        N = 20
+        g = EquidistantGrid(N, 0//1,2//1)
+        h = only(spacing(g))
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = undivided_skewed04(g, p)
+
+            @testset "x^$k" for k ∈ 0:p
+                v  = evalOn(g, x->monomial(x,k))
+                vₚₓ = evalOn(g, x->monomial(x,k-p))
+
+                @test D*v == h^p * vₚₓ
+            end
+        end
+    end
+
+    @testset "transpose equality" begin
+        function get_matrix(D)
+            N = only(range_size(D))
+            M = only(domain_size(D))
+
+            Dmat = zeros(N,M)
+            e = zeros(M)
+            for i ∈ 1:M
+                if i > 1
+                    e[i-1] = 0.
+                end
+                e[i] = 1.
+                Dmat[:,i] = D*e
+            end
+
+            return Dmat
+        end
+
+        g = EquidistantGrid(11, 0., 1.)
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = undivided_skewed04(g, p)
+
+            D̄  = get_matrix(D)
+            D̄ᵀ = get_matrix(Dᵀ)
+
+            @test D̄ == D̄ᵀ'
+        end
+    end
+end
+
+@testset "dissipation_interior_weights" begin
+    @test dissipation_interior_weights(1) == (-1, 1)
+    @test dissipation_interior_weights(2) == (1,-2, 1)
+    @test dissipation_interior_weights(3) == (-1, 3,-3, 1)
+    @test dissipation_interior_weights(4) == (1, -4, 6, -4, 1)
+end
+
+@testset "dissipation_interior_stencil" begin
+    @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1, 1, center=2)
+    @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil( 1,-2, 1, center=2)
+    @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1, 3,-3, 1, center=3)
+    @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil( 1,-4, 6,-4, 1, center=3)
+end
+
+@testset "dissipation_transpose_interior_stencil" begin
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(1,-1, center=1)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2, 1, center=2)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(1,-3, 3,-1, center=2)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1,-4, 6,-4, 1, center=3)
+end
+
+@testset "midpoint" begin
+    @test midpoint((1,1)) == 2
+    @test midpoint((1,1,1)) == 2
+    @test midpoint((1,1,1,1)) == 3
+    @test midpoint((1,1,1,1,1)) == 3
+end
+
+@testset "midpoint_transpose" begin
+    @test midpoint_transpose((1,1)) == 1
+    @test midpoint_transpose((1,1,1)) == 2
+    @test midpoint_transpose((1,1,1,1)) == 2
+    @test midpoint_transpose((1,1,1,1,1)) == 3
+end
+
+@testset "dissipation_lower_closure_size" begin
+    @test dissipation_lower_closure_size((1,1)) == 1
+    @test dissipation_lower_closure_size((1,1,1)) == 1
+    @test dissipation_lower_closure_size((1,1,1,1)) == 2
+    @test dissipation_lower_closure_size((1,1,1,1,1)) == 2
+end
+
+@testset "dissipation_upper_closure_size" begin
+    @test dissipation_upper_closure_size((1,1)) == 0
+    @test dissipation_upper_closure_size((1,1,1)) == 1
+    @test dissipation_upper_closure_size((1,1,1,1)) == 1
+    @test dissipation_upper_closure_size((1,1,1,1,1)) == 2
+end
+
+@testset "dissipation_lower_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil(-1, 1, center=1),
+        ),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, center=1),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,3,-3,1, center=1),
+            Stencil(-1,3,-3,1, center=2),
+        ),
+        (1, -4, 6, -4, 1) => (
+            Stencil(1, -4, 6, -4, 1, center=1),
+            Stencil(1, -4, 6, -4, 1, center=2),
+        )
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_lower_closure_stencils(w) == closure_stencils
+    end
+end
+
+@testset "dissipation_upper_closure_stencils" begin
+    cases = (
+        (-1,1) => (),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, center=3),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,3,-3,1, center=4),
+        ),
+        (1, -4, 6, -4, 1) => (
+            Stencil(1, -4, 6, -4, 1, center=4),
+            Stencil(1, -4, 6, -4, 1, center=5),
+        )
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_upper_closure_stencils(w) == closure_stencils
+    end
+end
+
+
+@testset "dissipation_transpose_lower_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil(-1,-1, 0, center=1),
+            Stencil( 1, 1,-1, center=2),
+        ),
+        (1,-2,1) => (
+            Stencil( 1, 1, 0, 0, center=1),
+            Stencil(-2,-2, 1, 0, center=2),
+            Stencil( 1, 1,-2, 1, center=3),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,-1,-1, 0, 0, 0, center=1),
+            Stencil( 3, 3, 3,-1, 0, 0, center=2),
+            Stencil(-3,-3,-3, 3,-1, 0, center=3),
+            Stencil( 1, 1, 1,-3, 3,-1, center=4),
+        ),
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_transpose_lower_closure_stencils(w) == closure_stencils
+    end
+end
+
+@testset "dissipation_transpose_upper_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil( 1,-1, center = 1),
+            Stencil( 0, 1, center = 2),
+        ),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, 1, center=2),
+            Stencil( 0, 1,-2,-2, center=3),
+            Stencil( 0, 0, 1, 1, center=4),
+        ),
+        (-1,3,-3,1) => (
+            Stencil( 1,-3, 3,-1,-1, center=2),
+            Stencil( 0, 1,-3, 3, 3, center=3),
+            Stencil( 0, 0, 1,-3,-3, center=4),
+            Stencil( 0, 0, 0, 1, 1, center=5),
+        ),
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_transpose_upper_closure_stencils(w) == closure_stencils
+    end
+end
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -6,6 +6,8 @@
 
 import Sbplib.SbpOperators.VolumeOperator
 
+# TODO: Refactor these test to look more like the tests in first_derivative_test.jl.
+
 @testset "SecondDerivative" begin
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
     stencil_set = read_stencil_set(operator_path; order=4)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -0,0 +1,82 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+import Sbplib.SbpOperators.Stencil
+import Sbplib.SbpOperators.StencilOperatorDistinctClosures
+import Sbplib.SbpOperators.stencil_operator_distinct_closures
+
+@testset "stencil_operator_distinct_closures" begin
+    lower_closure = (
+        Stencil(-1,1, center=1),
+    )
+
+    inner_stencil = Stencil(-2,2, center=1)
+
+    upper_closure = (
+        Stencil(-3,3, center=1),
+        Stencil(-4,4, center=2),
+    )
+
+    g₁ = EquidistantGrid(5, 0., 1.)
+    g₂ = EquidistantGrid((5,5), (0.,0.), (1.,1.))
+    h = 1/4
+
+    A₁  = stencil_operator_distinct_closures(g₁, inner_stencil, lower_closure, upper_closure, 1)
+    A₂¹ = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 1)
+    A₂² = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 2)
+
+    v₁ = evalOn(g₁, x->x)
+
+    u = [1., 2., 2., 3., 4.]*h
+    @test A₁*v₁ == u
+
+    v₂ = evalOn(g₂, (x,y)-> x + 3y)
+    @test A₂¹*v₂ == repeat(u, 1, 5)
+    @test A₂²*v₂ == repeat(3u', 5, 1)
+end
+
+@testset "StencilOperatorDistinctClosures" begin
+    g = EquidistantGrid(11, 0., 1.)
+
+    lower_closure = (
+        Stencil(-1,1, center=1),
+        Stencil(-2,2, center=2),
+    )
+
+    inner_stencil = Stencil(-3,3, center=1)
+
+    upper_closure = (
+        Stencil(4,-4,4, center=1),
+        Stencil(5,-5,5, center=2),
+        Stencil(6,-6,6, center=3),
+    )
+
+    A = StencilOperatorDistinctClosures(g, inner_stencil, lower_closure, upper_closure)
+    @test A isa LazyTensor{T,1,1} where T
+
+    @test SbpOperators.lower_closure_size(A) == 2
+    @test SbpOperators.upper_closure_size(A) == 3
+
+    @test domain_size(A) == (11,)
+    @test range_size(A) == (11,)
+
+    v = rand(11)
+    @testset "apply" begin
+        # Lower closure
+        @test LazyTensors.apply(A, v, 1) ≈ 1*(-v[1] + v[2])
+        @test LazyTensors.apply(A, v, 2) ≈ 2*(-v[1] + v[2])
+
+        # Interior
+        @test LazyTensors.apply(A, v, 3) ≈ 3*(-v[3] + v[4])
+        @test LazyTensors.apply(A, v, 4) ≈ 3*(-v[4] + v[5])
+        @test LazyTensors.apply(A, v, 8) ≈ 3*(-v[8] + v[9])
+
+        # Upper closure
+        @test LazyTensors.apply(A, v,  9) ≈ 4*(v[9] - v[10] + v[11])
+        @test LazyTensors.apply(A, v, 10) ≈ 5*(v[9] - v[10] + v[11])
+        @test LazyTensors.apply(A, v, 11) ≈ 6*(v[9] - v[10] + v[11])
+    end
+end
--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Tue Nov 01 22:44:00 2022 +0100
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Fri Feb 10 08:36:56 2023 +0100
@@ -7,120 +7,76 @@
 
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.VolumeOperator
-import Sbplib.SbpOperators.volume_operator
 import Sbplib.SbpOperators.odd
 import Sbplib.SbpOperators.even
 
+
 @testset "VolumeOperator" begin
     inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
-    closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2))
-    g_1D = EquidistantGrid(11,0.,1.)
-    g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
-    g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
+    closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2))
+    g = EquidistantGrid(11,0.,1.)
+
     @testset "Constructors" begin
-        @testset "1D" begin
-            op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)
-            @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even)
-            @test op == volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
-            @test op isa LazyTensor{T,1,1} where T
-        end
-        @testset "2D" begin
-            op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
-            Ix = IdentityTensor{Float64}((11,))
-            Iy = IdentityTensor{Float64}((12,))
-            @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy
-            @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)
-            @test op_x isa LazyTensor{T,2,2} where T
-            @test op_y isa LazyTensor{T,2,2} where T
-        end
-        @testset "3D" begin
-            op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
-            op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
-            Ix = IdentityTensor{Float64}((11,))
-            Iy = IdentityTensor{Float64}((12,))
-            Iz = IdentityTensor{Float64}((10,))
-            @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy⊗Iz
-            @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)⊗Iz
-            @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),even)
-            @test op_x isa LazyTensor{T,3,3} where T
-            @test op_y isa LazyTensor{T,3,3} where T
-            @test op_z isa LazyTensor{T,3,3} where T
-        end
+        op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)
+        @test op == VolumeOperator(g,inner_stencil,closure_stencils,even)
+        @test op isa LazyTensor{T,1,1} where T
     end
 
     @testset "Sizes" begin
-        @testset "1D" begin
-            op = volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
-            @test range_size(op) == domain_size(op) == size(g_1D)
-        end
-
-        @testset "2D" begin
-            op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
-            @test range_size(op_y) == domain_size(op_y) ==
-                  range_size(op_x) == domain_size(op_x) == size(g_2D)
-        end
-        @testset "3D" begin
-            op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
-            op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
-            @test range_size(op_z) == domain_size(op_z) ==
-                  range_size(op_y) == domain_size(op_y) ==
-                  range_size(op_x) == domain_size(op_x) == size(g_3D)
-        end
+        op = VolumeOperator(g,inner_stencil,closure_stencils,even)
+        @test range_size(op) == domain_size(op) == size(g)
     end
 
-    op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-    op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2)
-    v = zeros(size(g_2D))
-    Nx = size(g_2D)[1]
-    Ny = size(g_2D)[2]
-    for i = 1:Nx
-        v[i,:] .= i
+
+    op_even = VolumeOperator(g, inner_stencil, closure_stencils, even)
+    op_odd =  VolumeOperator(g, inner_stencil, closure_stencils, odd)
+
+    N = size(g)[1]
+    v = rand(N)
+
+    r_even = copy(v)
+    r_odd  = copy(v)
+
+    r_even[1] = (v[1] + v[2])/2
+    r_odd[1]  = (v[1] + v[2])/2
+
+    r_even[2] = 2v[1] + v[2]
+    r_odd[2]  = 2v[1] + v[2]
+
+    for i ∈ 3:N-2
+        r_even[i] = (v[i-1] + 2v[i] + v[i+1])/4
+        r_odd[i]  = (v[i-1] + 2v[i] + v[i+1])/4
     end
-    rx = copy(v)
-    rx[1,:] .= 1.5
-    rx[Nx,:] .= (2*Nx-1)/2
-    ry = copy(v)
-    ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny]
+
+    r_even[N-1] =  v[N-1] + 2v[N]
+    r_odd[N-1]  = -v[N-1] - 2v[N]
+
+    r_even[N] =  (v[N-1] + v[N])/2
+    r_odd[N]  = -(v[N-1] + v[N])/2
+
 
     @testset "Application" begin
-        @test op_x*v ≈ rx rtol = 1e-14
-        @test op_y*v ≈ ry rtol = 1e-14
+        @test op_even*v ≈ r_even
+        @test op_odd*v  ≈ r_odd
 
-        @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64
+        @test (op_even*rand(ComplexF64,size(g)))[2] isa ComplexF64
     end
 
     @testset "Regions" begin
-        @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14
-        @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14
-        @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14
-        @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14
-        @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14
-
-        @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)]
-        @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)]
+        @test (op_even*v)[Index(1,Lower)]    ≈ r_even[1]
+        @test (op_even*v)[Index(2,Lower)]    ≈ r_even[2]
+        @test (op_even*v)[Index(6,Interior)] ≈ r_even[6]
+        @test (op_even*v)[Index(10,Upper)]   ≈ r_even[10]
+        @test (op_even*v)[Index(11,Upper)]   ≈ r_even[11]
 
-        @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14
-
-        @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)]
-        @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)]
+        @test_throws BoundsError (op_even*v)[Index(3,Lower)]
+        @test_throws BoundsError (op_even*v)[Index(9,Upper)]
     end
 
     @testset "Inferred" begin
-        @test_skip @inferred apply(op_x, v,1,1)
-        @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower))
-        @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower))
-        @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower))
-        @test_skip @inferred apply(op_y, v,1,1)
-        @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower))
-        @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior))
-        @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper))
+        @inferred apply(op_even, v, 1)
+        @inferred apply(op_even, v, Index(1,Lower))
+        @inferred apply(op_even, v, Index(6,Interior))
+        @inferred apply(op_even, v, Index(11,Upper))
     end
 end