Mercurial > repos > public > sbplib_julia

changeset 1089:2278730f9cee refactor/sbpoperators/inflation

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, 10 May 2022 20:24:20 +0200
parents 62f321caa964 (current diff) 74c54996de6a (diff)
children 6f51160c7ca7 157a78959e5d
files src/LazyTensors/lazy_tensor_operations.jl src/SbpOperators/volumeops/volume_operator.jl test/LazyTensors/lazy_tensor_operations_test.jl
diffstat 15 files changed, 440 insertions(+), 69 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Tue May 10 20:21:15 2022 +0200
+++ b/Manifest.toml	Tue May 10 20:24:20 2022 +0200
@@ -1,13 +1,13 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.0"
+julia_version = "1.7.1"
 manifest_format = "2.0"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
+git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.3.1"
+version = "3.3.3"
 
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
--- a/Notes.md	Tue May 10 20:21:15 2022 +0200
+++ b/Notes.md	Tue May 10 20:24:20 2022 +0200
@@ -146,12 +146,52 @@
  - [ ] 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.
+
 ## Regions and tensormappings
 - [ ] Use a trait to indicate if a LazyTensor uses indices with regions.
     The default should be that they do NOT.
         - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle?
  - [ ] Figure out repeated application of regioned LazyTensors. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work?
 
+### Ideas for information sharing functions
+```julia
+using StaticArrays
+
+function regions(op::SecondDerivativeVariable)
+    t = ntuple(i->(Interior(),),range_dim(op))
+    return Base.setindex(t, (Lower(), Interior(), Upper()), derivative_direction(op))
+end
+
+function regionsizes(op::SecondDerivativeVariable)
+    sz = tuple.(range_size(op))
+
+    cl = closuresize(op)
+    return Base.setindex(sz, (cl, n-2cl, cl), derivative_direction(op))
+end
+
+
+g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1
+c = evalOn(g, (x,y)->x+y)
+
+D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1)
+@test regions(D₂ᶜ) == (
+    (Lower(), Interior(), Upper()),
+    (Interior(),),
+)
+@test regionsizes(D₂ᶜ) == ((1,9,1),(9,))
+
+
+D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2)
+@test regions(D₂ᶜ) == (
+    (Interior(),),
+    (Lower(), Interior(), Upper()),
+)
+@test regionsizes(D₂ᶜ) == ((11,),(1,7,1))
+```
+
+
 ## Boundschecking and dimension checking
 Does it make sense to have boundschecking only in getindex methods?
 This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway.
--- a/TODO.md	Tue May 10 20:21:15 2022 +0200
+++ b/TODO.md	Tue May 10 20:24:20 2022 +0200
@@ -1,14 +1,13 @@
 # TODO
 
-## Skämskudde
- - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
- - [ ] Skriv tester
 
 ## 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.
  - [ ] Replace getindex hack for flattening tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`)
  - [ ] Use `@inferred` in a lot of tests.
+ - [ ] Replace `@inferred` tests with a benchmark suite that automatically tests for regressions.
  - [ ] Make sure we are setting tolerances in tests in a consistent way
  - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index
  - [ ] Clean up RegionIndices
@@ -19,12 +18,11 @@
  - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`?
  - [ ] 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
 
-## Repo
- - [ ] Rename repo to Sbplib.jl
 
-# Wrap up tasks
+ - [ ] 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?
  - [ ] Kolla att vi har @inbounds och @propagate_inbounds på rätt ställen
  - [ ] Kolla att vi gör boundschecks överallt och att de är markerade med @boundscheck
  - [ ] Kolla att vi har @inline på rätt ställen
--- a/src/LazyTensors/LazyTensors.jl	Tue May 10 20:21:15 2022 +0200
+++ b/src/LazyTensors/LazyTensors.jl	Tue May 10 20:24:20 2022 +0200
@@ -3,9 +3,10 @@
 export TensorApplication
 export TensorTranspose
 export TensorComposition
-export DenseTensor
 export IdentityTensor
 export ScalingTensor
+export DiagonalTensor
+export DenseTensor
 export InflatedTensor
 export LazyOuterProduct
 export ⊗
--- a/src/LazyTensors/lazy_tensor_operations.jl	Tue May 10 20:21:15 2022 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Tue May 10 20:24:20 2022 +0200
@@ -98,7 +98,6 @@
     apply_transpose(c.t2, c.t1'*v, I...)
 end
 
-
 """
     TensorComposition(tm, tmi::IdentityTensor)
     TensorComposition(tmi::IdentityTensor, tm)
@@ -120,6 +119,8 @@
     return tmi
 end
 
+Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm)
+Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm
 
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}
--- a/src/LazyTensors/tensor_types.jl	Tue May 10 20:21:15 2022 +0200
+++ b/src/LazyTensors/tensor_types.jl	Tue May 10 20:24:20 2022 +0200
@@ -37,6 +37,24 @@
 
 
 """
+    DiagonalTensor{T,D,...} <: LazyTensor{T,D,D}
+    DiagonalTensor(a::AbstractArray)
+
+A lazy tensor with diagonal `a`.
+"""
+struct DiagonalTensor{T,D,AT<:AbstractArray{T,D}} <: LazyTensor{T,D,D}
+    diagonal::AT
+end
+
+range_size(tm::DiagonalTensor) = size(tm.diagonal)
+domain_size(tm::DiagonalTensor) = size(tm.diagonal)
+
+
+LazyTensors.apply(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...]
+LazyTensors.apply_transpose(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...]
+
+
+"""
     DenseTensor{T,R,D,...}(A, range_indicies, domain_indicies)
 
 LazyTensor defined by the AbstractArray A. `range_indicies` and `domain_indicies` define which indicies of A should
--- a/src/SbpOperators/stencil.jl	Tue May 10 20:21:15 2022 +0200
+++ b/src/SbpOperators/stencil.jl	Tue May 10 20:24:20 2022 +0200
@@ -1,11 +1,12 @@
 export CenteredStencil
+export CenteredNestedStencil
 
 struct Stencil{T,N}
-    range::Tuple{Int,Int}
+    range::UnitRange{Int64}
     weights::NTuple{N,T}
 
-    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N}
-        @assert range[2]-range[1]+1 == N
+    function Stencil(range::UnitRange,weights::NTuple{N,T}) where {T, N}
+        @assert length(range) == N
         new{T,N}(range,weights)
     end
 end
@@ -15,27 +16,30 @@
 
 Create a stencil with the given weights with element `center` as the center of the stencil.
 """
-function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error
+function Stencil(weights...; center::Int)
+    weights = promote(weights...)
     N = length(weights)
-    range = (1, N) .- center
+    range = (1:N) .- center
 
     return Stencil(range, weights)
 end
 
-function Stencil{T}(s::Stencil) where T
-    return Stencil(s.range, T.(s.weights))
-end
+Stencil{T,N}(s::Stencil{S,N}) where {T,S,N} = Stencil(s.range, T.(s.weights))
+Stencil{T}(s::Stencil) where T = Stencil{T,length(s)}(s)
 
-Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil)
+Base.convert(::Type{Stencil{T1,N}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s)
+Base.convert(::Type{Stencil{T1}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s)
 
-function CenteredStencil(weights::Vararg)
+Base.promote_rule(::Type{Stencil{T1,N}}, ::Type{Stencil{T2,N}}) where {T1,T2,N} = Stencil{promote_type(T1,T2),N}
+
+function CenteredStencil(weights...)
     if iseven(length(weights))
         throw(ArgumentError("a centered stencil must have an odd number of weights."))
     end
 
     r = length(weights) ÷ 2
 
-    return Stencil((-r, r), weights)
+    return Stencil(-r:r, weights)
 end
 
 
@@ -48,7 +52,8 @@
     return Stencil(s.range, a.*s.weights)
 end
 
-Base.eltype(::Stencil{T}) where T = T
+Base.eltype(::Stencil{T,N}) where {T,N} = T
+Base.length(::Stencil{T,N}) where {T,N} = N
 
 function flip(s::Stencil)
     range = (-s.range[2], -s.range[1])
@@ -57,24 +62,103 @@
 
 # Provides index into the Stencil based on offset for the root element
 @inline function Base.getindex(s::Stencil, i::Int)
-    @boundscheck if i < s.range[1] || s.range[2] < i
+    @boundscheck if i ∉ s.range
         return zero(eltype(s))
     end
     return s.weights[1 + i - s.range[1]]
 end
 
-Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = s.weights[1]*v[i + s.range[1]]
-    @simd for k ∈ 2:N
-        w += s.weights[k]*v[i + s.range[1] + k-1]
+Base.@propagate_inbounds @inline function apply_stencil(s::Stencil, v::AbstractVector, i::Int)
+    w = zero(promote_type(eltype(s),eltype(v)))
+    @simd for k ∈ 1:length(s)
+        w += s.weights[k]*v[i + s.range[k]]
+    end
+
+    return w
+end
+
+Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil, v::AbstractVector, i::Int)
+    w = zero(promote_type(eltype(s),eltype(v)))
+    @simd for k ∈ length(s):-1:1
+        w += s.weights[k]*v[i - s.range[k]]
     end
     return w
 end
 
-Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = s.weights[N]*v[i - s.range[2]]
-    @simd for k ∈ N-1:-1:1
-        w += s.weights[k]*v[i - s.range[1] - k + 1]
-    end
-    return w
+
+struct NestedStencil{T,N,M}
+    s::Stencil{Stencil{T,N},M}
+end
+
+# Stencil input
+NestedStencil(s::Vararg{Stencil}; center) = NestedStencil(Stencil(s... ; center))
+CenteredNestedStencil(s::Vararg{Stencil}) = NestedStencil(CenteredStencil(s...))
+
+# Tuple input
+function NestedStencil(weights::Vararg{NTuple{N,Any}}; center) where N
+    inner_stencils = map(w -> Stencil(w...; center), weights)
+    return NestedStencil(Stencil(inner_stencils... ; center))
+end
+function CenteredNestedStencil(weights::Vararg{NTuple{N,Any}}) where N
+    inner_stencils = map(w->CenteredStencil(w...), weights)
+    return CenteredNestedStencil(inner_stencils...)
+end
+
+
+# Conversion
+function NestedStencil{T,N,M}(ns::NestedStencil{S,N,M}) where {T,S,N,M}
+    return NestedStencil(Stencil{Stencil{T}}(ns.s))
+end
+
+function NestedStencil{T}(ns::NestedStencil{S,N,M}) where {T,S,N,M}
+    NestedStencil{T,N,M}(ns)
+end
+
+function Base.convert(::Type{NestedStencil{T,N,M}}, s::NestedStencil{S,N,M}) where {T,S,N,M}
+    return NestedStencil{T,N,M}(s)
+end
+Base.convert(::Type{NestedStencil{T}}, stencil) where T = NestedStencil{T}(stencil)
+
+function Base.promote_rule(::Type{NestedStencil{T,N,M}}, ::Type{NestedStencil{S,N,M}}) where {T,S,N,M}
+    return NestedStencil{promote_type(T,S),N,M}
 end
+
+Base.eltype(::NestedStencil{T}) where T = T
+
+function scale(ns::NestedStencil, a)
+    range = ns.s.range
+    weights = ns.s.weights
+
+    return NestedStencil(Stencil(range, scale.(weights,a)))
+end
+
+function flip(ns::NestedStencil)
+    s_flip = flip(ns.s)
+    return NestedStencil(Stencil(s_flip.range, flip.(s_flip.weights)))
+end
+
+Base.getindex(ns::NestedStencil, i::Int) = ns.s[i]
+
+"Apply inner stencils to `c` and get a concrete stencil"
+Base.@propagate_inbounds function apply_inner_stencils(ns::NestedStencil, c::AbstractVector, i::Int)
+    weights = apply_stencil.(ns.s.weights, Ref(c), i)
+    return Stencil(ns.s.range, weights)
+end
+
+"Apply the whole nested stencil"
+Base.@propagate_inbounds function apply_stencil(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int)
+    s = apply_inner_stencils(ns,c,i)
+    return apply_stencil(s, v, i)
+end
+
+"Apply inner stencils backwards to `c` and get a concrete stencil"
+Base.@propagate_inbounds @inline function apply_inner_stencils_backwards(ns::NestedStencil, c::AbstractVector, i::Int)
+    weights = apply_stencil_backwards.(ns.s.weights, Ref(c), i)
+    return Stencil(ns.s.range, weights)
+end
+
+"Apply the whole nested stencil backwards"
+Base.@propagate_inbounds @inline function apply_stencil_backwards(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int)
+    s = apply_inner_stencils_backwards(ns,c,i)
+    return apply_stencil_backwards(s, v, i)
+end
--- a/src/SbpOperators/volumeops/volume_operator.jl	Tue May 10 20:21:15 2022 +0200
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Tue May 10 20:24:20 2022 +0200
@@ -18,7 +18,7 @@
 # TBD: Should the inflation happen here or should we remove this method and do it at the caller instead?
 
 """
-    VolumeOperator{T,N,M,K} <: TensorOperator{T,1}
+    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}
--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Tue May 10 20:21:15 2022 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Tue May 10 20:24:20 2022 +0200
@@ -181,6 +181,14 @@
 
     @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64
     @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64
+
+    a = 2.
+    v = rand(3)
+    @test a*Ã isa TensorComposition
+    @test a*Ã == Ã*a
+    @test range_size(a*Ã) == range_size(Ã)
+    @test domain_size(a*Ã) == domain_size(Ã)
+    @test a*Ã*v == a.*A*v
 end
 
 
--- a/test/LazyTensors/tensor_types_test.jl	Tue May 10 20:21:15 2022 +0200
+++ b/test/LazyTensors/tensor_types_test.jl	Tue May 10 20:24:20 2022 +0200
@@ -1,5 +1,6 @@
 using Test
 using Sbplib.LazyTensors
+using BenchmarkTools
 
 @testset "IdentityTensor" begin
     @test IdentityTensor{Float64}((4,5)) isa IdentityTensor{T,2} where T
@@ -58,6 +59,47 @@
     @inferred (st'*v)[2,2]
 end
 
+@testset "DiagonalTensor" begin
+    @test DiagonalTensor([1,2,3,4]) isa LazyTensor{Int,1,1}
+    @test DiagonalTensor([1 2 3; 4 5 6]) isa LazyTensor{Int,2,2}
+    @test DiagonalTensor([1. 2. 3.; 4. 5. 6.]) isa LazyTensor{Float64,2,2}
+
+    @test range_size(DiagonalTensor([1,2,3,4])) == (4,)
+    @test domain_size(DiagonalTensor([1,2,3,4])) == (4,)
+
+    @test range_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3)
+    @test domain_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3)
+
+    @testset "apply size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)]
+        diag = rand(sz...)
+        tm = DiagonalTensor(diag)
+
+        v = rand(sz...)
+
+        @test tm*v == diag.*v
+        @test tm'*v == diag.*v
+    end
+
+
+    @testset "allocations size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)]
+        diag = rand(sz...)
+        tm = DiagonalTensor(diag)
+
+        v = rand(sz...)
+
+        @test tm*v == diag.*v
+        @test tm'*v == diag.*v
+    end
+
+    sz = (3,2)
+    diag = rand(sz...)
+    tm = DiagonalTensor(diag)
+
+    v = rand(sz...)
+    LazyTensors.apply(tm,v, 2,1)
+    @test (@ballocated LazyTensors.apply($tm,$v, 2,1)) == 0
+end
+
 
 @testset "DenseTensor" begin
     # Test a standard matrix-vector product
--- a/test/Manifest.toml	Tue May 10 20:21:15 2022 +0200
+++ b/test/Manifest.toml	Tue May 10 20:24:20 2022 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.0"
+julia_version = "1.7.1"
 manifest_format = "2.0"
 
 [[deps.ArgTools]]
@@ -12,11 +12,17 @@
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[deps.BenchmarkTools]]
+deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "1.3.1"
+
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47"
+git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.11.2"
+version = "1.14.0"
 
 [[deps.ChangesOfVariables]]
 deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
@@ -26,9 +32,9 @@
 
 [[deps.Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
+git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "3.41.0"
+version = "3.42.0"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -48,10 +54,10 @@
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
 [[deps.DiffRules]]
-deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "9bc5dac3c8b6706b58ad5ce24cffd9861f07c94f"
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
@@ -78,9 +84,9 @@
 
 [[deps.InverseFunctions]]
 deps = ["Test"]
-git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
+git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.2"
+version = "0.1.3"
 
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
@@ -89,9 +95,15 @@
 
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
-git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.3.0"
+version = "1.4.1"
+
+[[deps.JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.3"
 
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
@@ -118,9 +130,9 @@
 
 [[deps.LogExpFunctions]]
 deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
+git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.6"
+version = "0.3.10"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -140,9 +152,9 @@
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 
 [[deps.NaNMath]]
-git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283"
+git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "0.3.6"
+version = "1.0.0"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -161,20 +173,30 @@
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
+[[deps.Parsers]]
+deps = ["Dates"]
+git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "2.2.3"
+
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
+git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.2.2"
+version = "1.2.5"
 
 [[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[deps.Profile]]
+deps = ["Printf"]
+uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -185,9 +207,9 @@
 
 [[deps.Requires]]
 deps = ["UUIDs"]
-git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.2.0"
+version = "1.3.0"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -208,9 +230,9 @@
 
 [[deps.SpecialFunctions]]
 deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e"
+git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.0.0"
+version = "2.1.4"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -236,9 +258,9 @@
 
 [[deps.Tullio]]
 deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"]
-git-tree-sha1 = "0288b7a395fc412952baf756fac94e4f28bfec65"
+git-tree-sha1 = "7830c974acc69437a3fee35dd7b510a74cbc862d"
 uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
-version = "0.3.2"
+version = "0.3.3"
 
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
--- a/test/Project.toml	Tue May 10 20:21:15 2022 +0200
+++ b/test/Project.toml	Tue May 10 20:24:20 2022 +0200
@@ -1,4 +1,5 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Tue May 10 20:21:15 2022 +0200
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Tue May 10 20:24:20 2022 +0200
@@ -9,7 +9,7 @@
 import Sbplib.SbpOperators.boundary_operator
 
 @testset "BoundaryOperator" begin
-    closure_stencil = Stencil((0,2), (2.,1.,3.))
+    closure_stencil = Stencil(2.,1.,3.; center = 1)
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
--- a/test/SbpOperators/stencil_test.jl	Tue May 10 20:21:15 2022 +0200
+++ b/test/SbpOperators/stencil_test.jl	Tue May 10 20:24:20 2022 +0200
@@ -1,19 +1,25 @@
 using Test
 using Sbplib.SbpOperators
 import Sbplib.SbpOperators.Stencil
+import Sbplib.SbpOperators.NestedStencil
+import Sbplib.SbpOperators.scale
 
 @testset "Stencil" begin
-    s = Stencil((-2,2), (1.,2.,2.,3.,4.))
+    s = Stencil(-2:2, (1.,2.,2.,3.,4.))
     @test s isa Stencil{Float64, 5}
 
     @test eltype(s) == Float64
-    @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
+
+    @test length(s) == 5
+    @test length(Stencil(-1:2, (1,2,3,4))) == 4
+
+    @test SbpOperators.scale(s, 2) == Stencil(-2:2, (2.,4.,4.,6.,8.))
 
-    @test Stencil(1,2,3,4; center=1) == Stencil((0, 3),(1,2,3,4))
-    @test Stencil(1,2,3,4; center=2) == Stencil((-1, 2),(1,2,3,4))
-    @test Stencil(1,2,3,4; center=4) == Stencil((-3, 0),(1,2,3,4))
+    @test Stencil(1,2,3,4; center=1) == Stencil(0:3,(1,2,3,4))
+    @test Stencil(1,2,3,4; center=2) == Stencil(-1:2,(1,2,3,4))
+    @test Stencil(1,2,3,4; center=4) == Stencil(-3:0,(1,2,3,4))
 
-    @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
+    @test CenteredStencil(1,2,3,4,5) == Stencil(-2:2, (1,2,3,4,5))
     @test_throws ArgumentError CenteredStencil(1,2,3,4)
 
     # Changing the type of the weights
@@ -24,8 +30,145 @@
 
     @testset "convert" begin
         @test convert(Stencil{Float64}, Stencil(1,2,3,4,5; center=2)) == Stencil(1.,2.,3.,4.,5.; center=2)
-        @test convert(Stencil{Float64}, CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
-        @test convert(Stencil{Int}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
-        @test convert(Stencil{Rational}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+        @test convert(Stencil{Float64,5}, CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
+        @test convert(Stencil{Int,5}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
+        @test convert(Stencil{Rational,5}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+    end
+
+    @testset "promotion of weights" begin
+        @test Stencil(1.,2; center = 1) isa Stencil{Float64, 2}
+        @test Stencil(1,2//2; center = 1) isa Stencil{Rational{Int64}, 2}
+    end
+
+    @testset "promotion" begin
+        @test promote(Stencil(1,1;center=1), Stencil(2.,2.;center=2)) == (Stencil(1.,1.;center=1), Stencil(2.,2.;center=2))
+    end
+
+    @testset "type stability" begin
+        s_int = CenteredStencil(1,2,3)
+        s_float = CenteredStencil(1.,2.,3.)
+        v_int = rand(1:10,10);
+        v_float = rand(10);
+
+        @inferred SbpOperators.apply_stencil(s_int, v_int, 2)
+        @inferred SbpOperators.apply_stencil(s_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_int,  v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_float, v_int, 2)
+
+        @inferred SbpOperators.apply_stencil_backwards(s_int, v_int, 5)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, v_float, 5)
+        @inferred SbpOperators.apply_stencil_backwards(s_int,  v_float, 5)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, v_int, 5)
     end
 end
+
+@testset "NestedStencil" begin
+
+    @testset "Constructors" begin
+        s1 = CenteredStencil(-1, 1, 0)
+        s2 = CenteredStencil(-1, 0, 1)
+        s3 = CenteredStencil( 0,-1, 1)
+
+        ns = NestedStencil(CenteredStencil(s1,s2,s3))
+        @test ns isa NestedStencil{Int,3}
+
+        @test CenteredNestedStencil(s1,s2,s3) == ns
+
+        @test NestedStencil(s1,s2,s3, center = 2) == ns
+        @test NestedStencil(s1,s2,s3, center = 1) == NestedStencil(Stencil(s1,s2,s3, center=1))
+
+        @test NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=2) == ns
+        @test CenteredNestedStencil((-1,1,0),(-1,0,1),(0,-1,1)) == ns
+        @test NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=1) == NestedStencil(Stencil(
+            Stencil(-1, 1, 0; center=1),
+            Stencil(-1, 0, 1; center=1),
+            Stencil( 0,-1, 1; center=1);
+            center=1
+        ))
+
+        @testset "Error handling" begin
+        end
+    end
+
+    @testset "scale" begin
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=2)
+        @test SbpOperators.scale(ns, 2) == NestedStencil((-2,2,0),(-2,0,2),(0,-2,2), center=2)
+    end
+
+    @testset "conversion" begin
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=2)
+        @test NestedStencil{Float64}(ns) == NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.), center=2)
+        @test NestedStencil{Rational}(ns) == NestedStencil((-1//1,1//1,0//1),(-1//1,0//1,1//1),(0//1,-1//1,1//1), center=2)
+
+        @test convert(NestedStencil{Float64}, ns) == NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.), center=2)
+        @test convert(NestedStencil{Rational}, ns) == NestedStencil((-1//1,1//1,0//1),(-1//1,0//1,1//1),(0//1,-1//1,1//1), center=2)
+    end
+
+    @testset "promotion of weights" begin
+        @test NestedStencil((-1,1,0),(-1.,0.,1.),(0,-1,1), center=2) isa NestedStencil{Float64,3,3}
+        @test NestedStencil((-1,1,0),(-1,0,1),(0//1,-1,1), center=2) isa NestedStencil{Rational{Int64},3,3}
+    end
+
+    @testset "promotion" begin
+        promote(
+            CenteredNestedStencil((-1,1,0),(-1,0,1),(0,-1,1)),
+            CenteredNestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.))
+        ) == (
+            CenteredNestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.)),
+            CenteredNestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.))
+        )
+    end
+
+    @testset "apply" begin
+        c = [  1,  3,  6, 10, 15, 21, 28, 36, 45, 55]
+        v = [  2,  3,  5,  7, 11, 13, 17, 19, 23, 29]
+
+        # Centered
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-2,2), center=2)
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) == Stencil(4,9,10; center=2)
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) == Stencil(-5,-9,-8; center=2)
+
+        @test SbpOperators.apply_stencil(ns, c, v, 4) == 4*5 + 9*7 + 10*11
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -8*5 - 9*7 - 5*11
+
+        # Non-centered
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=1)
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) == Stencil(5,11,6; center=1)
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) == Stencil(-4,-7,-3; center=1)
+
+        @test SbpOperators.apply_stencil(ns, c, v, 4) == 5*7 + 11*11 + 6*13
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -3*3 - 7*5 - 4*7
+    end
+
+    @testset "type stability" begin
+        s_int = CenteredNestedStencil((1,2,3),(1,2,3),(1,2,3))
+        s_float = CenteredNestedStencil((1.,2.,3.),(1.,2.,3.),(1.,2.,3.))
+
+        v_int = rand(1:10,10);
+        v_float = rand(10);
+
+        c_int = rand(1:10,10);
+        c_float = rand(10);
+
+        @inferred SbpOperators.apply_stencil(s_int,   c_int, v_int,   2)
+        @inferred SbpOperators.apply_stencil(s_float, c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_int,   c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_float, c_int, v_int,   2)
+
+        @inferred SbpOperators.apply_stencil(s_int,   c_float, v_int,   2)
+        @inferred SbpOperators.apply_stencil(s_float, c_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_int,   c_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_float, c_float, v_int,   2)
+
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_int, v_int,   2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_int, v_int,   2)
+
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_float, v_int,   2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_float, v_float, 2)
+        @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
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue May 10 20:21:15 2022 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue May 10 20:24:20 2022 +0200
@@ -71,6 +71,7 @@
     end
 
     @testset "Accuracy on function" begin
+        # 1D
         g = EquidistantGrid(30, 0.,1.)
         v = evalOn(g, x->exp(x))
         @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)]
@@ -79,6 +80,18 @@
 
             @test D₁*v ≈ v rtol=tol
         end
+
+        # 2D
+        g = EquidistantGrid((30,60), (0.,0.),(1.,2.))
+        v = evalOn(g, (x,y)->exp(0.8x+1.2*y))
+        @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)]
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
+            Dx = first_derivative(g, stencil_set, 1)
+            Dy = first_derivative(g, stencil_set, 2)
+
+            @test Dx*v ≈ 0.8v rtol=tol
+            @test Dy*v ≈ 1.2v rtol=tol
+        end
     end
 end