changeset 867:313648b01504 feature/variable_derivatives

Start implementing nested stencils
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 19 Jan 2022 21:49:34 +0100
parents 9a2776352c2a
children e37ee63bf9ac
files src/SbpOperators/stencil.jl test/SbpOperators/stencil_test.jl
diffstat 2 files changed, 74 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/stencil.jl	Wed Jan 19 11:08:43 2022 +0100
+++ b/src/SbpOperators/stencil.jl	Wed Jan 19 21:49:34 2022 +0100
@@ -1,4 +1,5 @@
 export CenteredStencil
+export CenteredNestedStencil
 
 struct Stencil{T,N}
     range::Tuple{Int,Int}
@@ -78,3 +79,50 @@
     end
     return w
 end
+
+
+struct NestedStencil{T,N}
+    s::Stencil{Stencil{T,N},N}
+end
+
+NestedStencil(s::Vararg{Stencil}; center) = NestedStencil(Stencil(s... ; center))
+
+function NestedStencil(weights::Vararg{Tuple}; center)
+    inner_stencils = map((w,c) -> Stencil(w...,center=c), weights, 1:length(weights))
+    return NestedStencil(Stencil(inner_stencils... ; center))
+end
+
+CenteredNestedStencil(s...) = NestedStencil(CenteredStencil(s...))
+
+Base.eltype(::NestedStencil{T}) where T = T
+
+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_stencil(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_stencil(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_stencil_backwards(ns::NestedStencil, c::AbstractVector, v::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_stencil_backwards(ns,c,i)
+    return apply_stencil_backwards(s, v, i)
+end
--- a/test/SbpOperators/stencil_test.jl	Wed Jan 19 11:08:43 2022 +0100
+++ b/test/SbpOperators/stencil_test.jl	Wed Jan 19 21:49:34 2022 +0100
@@ -1,6 +1,8 @@
 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.))
@@ -29,3 +31,27 @@
         @test convert(Stencil{Rational}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
     end
 end
+
+@testset "NestedStencil" begin
+
+    @testset "Constructors" begin
+        s1 = Stencil(-1, 1, 0; center = 1)
+        s2 = Stencil(-1, 0, 1; center = 2)
+        s3 = Stencil( 0,-1, 1; center = 3)
+
+        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
+
+
+        @testset "Error handling" begin
+
+        end
+    end
+end