Mercurial > repos > public > sbplib_julia
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