Mercurial > repos > public > sbplib_julia
changeset 912:f800f97b3a45 feature/variable_derivatives
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 15 Feb 2022 15:15:11 +0100 |
| parents | 33c7e266e1a9 (diff) a378ab959b6f (current diff) |
| children | cad3a9f82009 |
| files | |
| diffstat | 11 files changed, 734 insertions(+), 36 deletions(-) [+] |
line wrap: on
line diff
--- a/Notes.md Tue Feb 15 15:14:28 2022 +0100 +++ b/Notes.md Tue Feb 15 15:15:11 2022 +0100 @@ -154,6 +154,43 @@ - [ ] 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 TensorMappings. 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 Feb 15 15:14:28 2022 +0100 +++ b/TODO.md Tue Feb 15 15:15:11 2022 +0100 @@ -23,6 +23,7 @@ - [ ] 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) + - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative ## Repo - [ ] Rename repo to Sbplib.jl
--- a/src/SbpOperators/SbpOperators.jl Tue Feb 15 15:14:28 2022 +0100 +++ b/src/SbpOperators/SbpOperators.jl Tue Feb 15 15:15:11 2022 +0100 @@ -9,11 +9,14 @@ even = 1 end +export closure_size + include("stencil.jl") include("readoperator.jl") include("volumeops/volume_operator.jl") include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/second_derivative.jl") +include("volumeops/derivatives/second_derivative_variable.jl") include("volumeops/laplace/laplace.jl") include("volumeops/inner_products/inner_product.jl") include("volumeops/inner_products/inverse_inner_product.jl")
--- a/src/SbpOperators/operators/standard_diagonal.toml Tue Feb 15 15:14:28 2022 +0100 +++ b/src/SbpOperators/operators/standard_diagonal.toml Tue Feb 15 15:15:11 2022 +0100 @@ -20,6 +20,10 @@ H.inner = "1" H.closure = ["1/2"] +e.closure = ["1"] +d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1} + + D1.inner_stencil = ["-1/2", "0", "1/2"] D1.closure_stencils = [ {s = ["-1", "1"], c = 1}, @@ -30,8 +34,10 @@ {s = ["1", "-2", "1"], c = 1}, ] -e.closure = ["1"] -d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1} +D2variable.inner_stencil = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]] +D2variable.closure_stencils = [ + {s = [["2", "-1", "0"],["-3", "1", "0"],["1","0","0"]], c = 1}, +] [[stencil_set]] @@ -40,6 +46,9 @@ H.inner = "1" H.closure = ["17/48", "59/48", "43/48", "49/48"] +e.closure = ["1"] +d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1} + D1.inner_stencil = ["1/12","-2/3","0","2/3","-1/12"] D1.closure_stencils = [ {s = [ "-24/17", "59/34", "-4/17", "-3/34", "0", "0"], c = 1}, @@ -56,5 +65,79 @@ {s = [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"], c = 4}, ] +D2variable.inner_stencil = [ + ["-1/8", "1/6", "-1/8", "0", "0" ], + [ "1/6", "1/2", "1/2", "1/6", "0" ], + ["-1/24", "-5/6", "-3/4", "-5/6", "-1/24"], + [ "0", "1/6", "1/2", "1/2", "1/6" ], + [ "0", "0", "-1/8", "1/6", "-1/8" ], +] +D2variable.closure_stencils = [ + {c = 1, s = [ + [ "920/289", "-59/68", "-81031200387/366633756146", "-69462376031/733267512292", "0", "0", "0", "0" ], + ["-1740/289", "0", "6025413881/7482321554", "1612249989/7482321554", "0", "0", "0", "0" ], + [ "1128/289", "59/68", "-6251815797/8526366422", "-639954015/17052732844", "0", "0", "0", "0" ], + [ "-308/289", "0", "1244724001/7482321554", "-752806667/7482321554", "0", "0", "0", "0" ], + [ "0", "0", "-148737261/10783345769", "148737261/10783345769", "0", "0", "0", "0" ], + [ "0", "0", "-3/833", "3/833", "0", "0", "0", "0" ], + ]}, + {c = 2, s = [ + [ "12/17", "0", "102125659/440136562", "27326271/440136562", "0", "0", "0", "0" ], + [ "-59/68", "0", "-156920047993625/159775733917868", "-12001237118451/79887866958934", "0", "0", "0", "0" ], + [ "2/17", "0", "1489556735319/1857857371138", "149729180391/1857857371138", "0", "0", "0", "0" ], + [ "3/68", "0", "-13235456910147/159775733917868", "3093263736297/79887866958934", "0", "0", "0", "0" ], + [ "0", "0", "67535018271/2349643145851", "-67535018271/2349643145851", "0", "0", "0", "0" ], + [ "0", "0", "441/181507", "-441/181507", "0", "0", "0", "0" ], + ]}, + {c = 3, s = [ + [ "-96/731", "59/172", "-6251815797/21566691538", "-639954015/43133383076", "0", "0", "0", "0" ], + [ "118/731", "0", "87883847383821/79887866958934", "8834021643069/79887866958934", "0", "0", "0", "0" ], + [ "-16/731", "-59/172", "-1134866646907639536627/727679167377258785038", "-13777050223300597/23487032885926596", "-26254/557679", "0", "0", "0" ], + [ "-6/731", "0", "14509020271326561681/14850595252597118062", "17220493277981/79887866958934", "1500708/7993399", "0", "0", "0" ], + [ "0", "0", "-4841930283098652915/21402328452272317207", "31597236232005/115132514146699", "-26254/185893", "0", "0", "0" ], + [ "0", "0", "-2318724711/1653303156799", "960119/1147305747", "13564/23980197", "0", "0", "0" ], + ]}, + {c = 4, s = [ + [ "-36/833", "0", "1244724001/21566691538", "-752806667/21566691538", "0", "0", "0", "0" ], + [ "177/3332", "0", "-780891957698673/7829010961975532", "3724542049827/79887866958934", "0", "0", "0", "0" ], + [ "-6/833", "0", "14509020271326561681/16922771334354855466", "2460070468283/13005001597966", "1500708/9108757", "0", "0", "0" ], + [ "-9/3332", "0", "-217407431400324796377/207908333536359652868", "-1950062198436997/3914505480987766", "-7476412/9108757", "-2/49", "0", "0" ], + [ "0", "0", "4959271814984644613/21402328452272317207", "47996144728947/115132514146699", "4502124/9108757", "8/49", "0", "0" ], + [ "0", "0", "-2258420001/1653303156799", "-1063649/8893843", "1473580/9108757", "-6/49", "0", "0" ], + ]}, + {c = 5, s = [ + [ "0", "0", "-49579087/10149031312", "49579087/10149031312", "0", "0", "0", "0" ], + [ "0", "0", "1328188692663/37594290333616", "-1328188692663/37594290333616", "0", "0", "0", "0" ], + [ "0", "0", "-1613976761032884305/7963657098519931984", "10532412077335/42840005263888", "-564461/4461432", "0", "0", "0" ], + [ "0", "0", "4959271814984644613/20965546238960637264", "15998714909649/37594290333616", "375177/743572", "1/6", "0", "0" ], + [ "0", "0", "-8386761355510099813/128413970713633903242", "-2224717261773437/2763180339520776", "-280535/371786", "-5/6", "-1/24", "0" ], + [ "0", "0", "13091810925/13226425254392", "35039615/213452232", "1118749/2230716", "1/2", "1/6", "0" ], + [ "0", "0", "0", "0", "-1/8", "1/6", "-1/8", "0" ], + ]}, + {c = 6, s = [ + [ "0", "0", "-1/784", "1/784", "0", "0", "0", "0" ], + [ "0", "0", "8673/2904112", "-8673/2904112", "0", "0", "0", "0" ], + [ "0", "0", "-33235054191/26452850508784", "960119/1280713392", "3391/6692148", "0", "0", "0" ], + [ "0", "0", "-752806667/539854092016", "-1063649/8712336", "368395/2230716", "-1/8", "0", "0" ], + [ "0", "0", "13091810925/13226425254392", "35039615/213452232", "1118749/2230716", "1/2", "1/6", "0" ], + [ "0", "0", "-660204843/13226425254392", "-3290636/80044587", "-5580181/6692148", "-3/4", "-5/6", "-1/24"], + [ "0", "0", "0", "0", "1/6", "1/2", "1/2", "1/6" ], + [ "0", "0", "0", "0", "0", "-1/8", "1/6", "-1/8" ], + ]} +] + + + +[[stencil_set]] + +order = 6 + +H.inner = "1" +H.closure = ["13649/43200", "12013/8640", "2711/4320", "5359/4320", "7877/8640", "43801/43200"] + + + + + e.closure = ["1"] -d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1} +d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"]
--- a/src/SbpOperators/readoperator.jl Tue Feb 15 15:14:28 2022 +0100 +++ b/src/SbpOperators/readoperator.jl Tue Feb 15 15:15:11 2022 +0100 @@ -4,6 +4,7 @@ export get_stencil_set export parse_stencil +export parse_nested_stencil export parse_scalar export parse_tuple @@ -106,6 +107,33 @@ end end + +""" + parse_nested_stencil(parsed_toml) + +Accept parsed TOML and read it as a nested tuple. + +See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref). +""" +function parse_nested_stencil(parsed_toml) + if parsed_toml isa Array + weights = parse_stencil.(parsed_toml) + return CenteredNestedStencil(weights...) + end + + center = parsed_toml["c"] + weights = parse_tuple.(parsed_toml["s"]) + return NestedStencil(weights...; center) +end + +""" + parse_nested_stencil(T, parsed_toml) + +Parse the input as a nested stencil with element type `T`. +""" +parse_nested_stencil(T, parsed_toml) = NestedStencil{T}(parse_nested_stencil(parsed_toml)) + + """ parse_scalar(parsed_toml)
--- a/src/SbpOperators/stencil.jl Tue Feb 15 15:14:28 2022 +0100 +++ b/src/SbpOperators/stencil.jl Tue Feb 15 15:15:11 2022 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Tue Feb 15 15:15:11 2022 +0100 @@ -0,0 +1,130 @@ +export SecondDerivativeVariable + +# """ +# SecondDerivativeVariable(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 `SecondDerivativeVariable` tensor mapping is returned. When `Dim>1`, the +# returned operator is the appropriate outer product of a one-dimensional +# operators and `IdentityMapping`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 = SecondDerivativeVariable(restrict(grid, direction), inner_stencil, closure_stencils, parity) +# # Create 1D IdentityMappings for each coordinate direction +# one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) +# Is = IdentityMapping{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 + +""" + SecondDerivativeVariable{T,N,M,K} <: TensorOperator{T,1} + +Implements the one-dimensional second derivative with variable coefficients. +""" +struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: TensorMapping{T,D,D} + inner_stencil::IStencil + closure_stencils::NTuple{M,CStencil} + size::NTuple{D,Int} + coefficient::TArray + + function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, size::NTuple{D,Int}, coefficient::AbstractArray) where {Dir,T,D,M} + IStencil = typeof(inner_stencil) + CStencil = eltype(closure_stencils) + TArray = typeof(coefficient) + return new{Dir,T,D,M,IStencil,CStencil,TArray}(inner_stencil,closure_stencils,size, coefficient) + end +end + +function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir) + Δxᵢ = spacing(grid)[dir] + scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2) + scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2) + return SecondDerivativeVariable{dir, dimension(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff) +end + +function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils) + return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1) +end + +function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir) + inner_stencil = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"]) + closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"]) + + return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, dir) +end + +derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir + +closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils) + +LazyTensors.range_size(op::SecondDerivativeVariable) = op.size +LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size + + +function derivative_view(op, a, I) + d = derivative_direction(op) + + Iview = Base.setindex(I,:,d) + return @view a[Iview...] + + # D = domain_dim(op) + # Iₗ, _, Iᵣ = split_tuple(I, Val(d-1), Val(1), Val(D-d)) + # return @view a[Iₗ..., :, Iᵣ...] +end + +function apply_lower(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i) +end + +function apply_interior(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + return apply_stencil(op.inner_stencil, c̃, ṽ, i) +end + +function apply_upper(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I::Vararg{Index}) + if I[derivative_direction(op)] isa Index{Lower} + return apply_lower(op, v, Int.(I)...) + elseif I[derivative_direction(op)] isa Index{Upper} + return apply_upper(op, v, Int.(I)...) + else + return apply_interior(op, v, Int.(I)...) + end +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...) + dir = derivative_direction(op) + + i = I[dir] + r = getregion(i, closure_size(op), op.size[dir]) + + I = map(i->Index(i, Interior), I) + I = Base.setindex(I, Index(i, r), dir) + return LazyTensors.apply(op, v, I...) +end + +# TODO: Rename SecondDerivativeVariable -> VariableSecondDerivative
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Tue Feb 15 15:14:28 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Tue Feb 15 15:15:11 2022 +0100 @@ -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/readoperator_test.jl Tue Feb 15 15:14:28 2022 +0100 +++ b/test/SbpOperators/readoperator_test.jl Tue Feb 15 15:15:11 2022 +0100 @@ -4,6 +4,7 @@ using Sbplib.SbpOperators import Sbplib.SbpOperators.Stencil +import Sbplib.SbpOperators.NestedStencil @testset "readoperator" begin toml_str = """ @@ -170,3 +171,18 @@ @test SbpOperators.parse_rational(2) isa Rational @test SbpOperators.parse_rational(2) == 2//1 end + +@testset "parse_nested_stencil" begin + toml = TOML.parse(""" + s1 = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]] + s2 = {s = [[ "2", "-1", "0"],[ "-3", "1", "0"],["1", "0", "0"]], c = 1} + s3 = {s = [[ "2", "-1", "0"],[ "-3", "1", "0"],["1", "0", "0"]], c = 2} + """) + + @test parse_nested_stencil(toml["s1"]) == CenteredNestedStencil((1//2, 1//2, 0//1),( -1//2, -1//1, -1//2),(0//1, 1//2, 1//2)) + @test parse_nested_stencil(toml["s2"]) == NestedStencil((2//1, -1//1, 0//1),( -3//1, 1//1, 0//1),(1//1, 0//1, 0//1), center = 1) + @test parse_nested_stencil(toml["s3"]) == NestedStencil((2//1, -1//1, 0//1),( -3//1, 1//1, 0//1),(1//1, 0//1, 0//1), center = 2) + + @test parse_nested_stencil(Float64, toml["s1"]) == CenteredNestedStencil((1/2, 1/2, 0.),( -1/2, -1., -1/2),(0., 1/2, 1/2)) + @test parse_nested_stencil(Int, toml["s2"]) == NestedStencil((2, -1, 0),( -3, 1, 0),(1, 0, 0), center = 1) +end
--- a/test/SbpOperators/stencil_test.jl Tue Feb 15 15:14:28 2022 +0100 +++ b/test/SbpOperators/stencil_test.jl Tue Feb 15 15:15:11 2022 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl Tue Feb 15 15:15:11 2022 +0100 @@ -0,0 +1,173 @@ +using Test + +using Sbplib.Grids +using Sbplib.LazyTensors +using Sbplib.SbpOperators +using Sbplib.RegionIndices +using Sbplib.SbpOperators: NestedStencil, CenteredNestedStencil + +using LinearAlgebra + +@testset "SecondDerivativeVariable" begin + interior_stencil = CenteredNestedStencil((1/2, 1/2, 0.),(-1/2, -1., -1/2),( 0., 1/2, 1/2)) + closure_stencils = [ + NestedStencil(( 2., -1., 0.),(-3., 1., 0.), (1., 0., 0.), center = 1), + ] + + @testset "1D" begin + g = EquidistantGrid(11, 0., 1.) + c = [ 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66] + @testset "Constructors" begin + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) isa TensorMapping + + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) + @test range_dim(D₂ᶜ) == 1 + @test domain_dim(D₂ᶜ) == 1 + end + + @testset "sizes" begin + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) + @test closure_size(D₂ᶜ) == 1 + @test range_size(D₂ᶜ) == (11,) + @test domain_size(D₂ᶜ) == (11,) + end + + @testset "application" begin + + function apply_to_functions(; v, c) + g = EquidistantGrid(11, 0., 10.) # h = 1 + c̄ = evalOn(g,c) + v̄ = evalOn(g,v) + + D₂ᶜ = SecondDerivativeVariable(g, c̄, interior_stencil, closure_stencils) + return D₂ᶜ*v̄ + end + + @test apply_to_functions(v=x->1., c=x-> -1.) == zeros(11) + @test apply_to_functions(v=x->1., c=x-> -x ) == zeros(11) + @test apply_to_functions(v=x->x, c=x-> 1.) == zeros(11) + @test apply_to_functions(v=x->x, c=x-> -x ) == -ones(11) + @test apply_to_functions(v=x->x^2, c=x-> 1.) == 2ones(11) + end + + @testset "Inferred" begin + g = EquidistantGrid(11, 0., 10.) # h = 1 + c̄ = evalOn(g,x-> -1) + v̄ = evalOn(g,x->1.) + + D₂ᶜ = SecondDerivativeVariable(g, c̄, interior_stencil, closure_stencils) + + @inferred SbpOperators.apply_lower(D₂ᶜ, v̄, 1) + @inferred SbpOperators.apply_interior(D₂ᶜ, v̄, 5) + @inferred SbpOperators.apply_upper(D₂ᶜ, v̄, 11) + @inferred (D₂ᶜ*v̄)[Index(1,Lower)] + end + end + + @testset "2D" begin + g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1 + c = evalOn(g, (x,y)->x+y) + @testset "Constructors" begin + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) isa TensorMapping + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) isa TensorMapping + + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) + @test range_dim(D₂ᶜ) == 2 + @test domain_dim(D₂ᶜ) == 2 + end + + @testset "sizes" begin + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) + @test range_size(D₂ᶜ) == (11,9) + @test domain_size(D₂ᶜ) == (11,9) + @test closure_size(D₂ᶜ) == 1 + + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) + @test range_size(D₂ᶜ) == (11,9) + @test domain_size(D₂ᶜ) == (11,9) + @test closure_size(D₂ᶜ) == 1 + end + + @testset "application" begin + function apply_to_functions(dir; v, c) + g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1 + c̄ = evalOn(g,c) + v̄ = evalOn(g,v) + + D₂ᶜ = SecondDerivativeVariable(g, c̄, interior_stencil, closure_stencils,dir) + return D₂ᶜ*v̄ + end + + # x-direction + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)-> -1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)->- x ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->x, c=(x,y)-> -x ) == -ones(11,9) + @test apply_to_functions(1,v=(x,y)->x^2, c=(x,y)-> 1.) == 2ones(11,9) + + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)->- y ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y, c=(x,y)-> -y ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y^2, c=(x,y)-> 1.) == zeros(11,9) + + # y-direction + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)-> -1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)->- y ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->y, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->y, c=(x,y)-> -y ) == -ones(11,9) + @test apply_to_functions(2,v=(x,y)->y^2, c=(x,y)-> 1.) == 2ones(11,9) + + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)->- x ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> -x ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x^2, c=(x,y)-> 1.) == zeros(11,9) + + + + # TBD: This should be moved somewhere else right? + @testset "real operators" begin + c(x,y) = exp(x) + exp(1.5(1-y)) + v(x,y) = sin(x) + cos(1.5(1-y)) + + Dxv(x,y) = cos(x)*exp(x) - (exp(x) + exp(1.5 - 1.5y))*sin(x) + Dyv(x,y) = -1.5(1.5exp(x) + 1.5exp(1.5 - 1.5y))*cos(1.5 - 1.5y) - 2.25exp(1.5 - 1.5y)*sin(1.5 - 1.5y) + + # g₁ = EquidistantGrid((30,37), (0.,0.), (1.,2.)) + g₁ = EquidistantGrid((60,67), (0.,0.), (1.,2.)) + g₂ = refine(g₁,2) + + c̄₁ = evalOn(g₁, c) + c̄₂ = evalOn(g₂, c) + + v̄₁ = evalOn(g₁, v) + v̄₂ = evalOn(g₂, v) + + + function convergence_rate_estimate(stencil_set, dir, Dv_true) + D₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, dir) + D₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, dir) + + Dv̄₁ = D₁*v̄₁ + Dv̄₂ = D₂*v̄₂ + + Dv₁ = evalOn(g₁,Dv_true) + Dv₂ = evalOn(g₂,Dv_true) + + e₁ = norm(Dv̄₁ - Dv₁)/norm(Dv₁) + e₂ = norm(Dv̄₂ - Dv₂)/norm(Dv₂) + + return log2(e₁/e₂) + end + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 1.5 rtol = 1e-1 + @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 1.5 rtol = 1e-1 + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) + @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 2.5 rtol = 1e-1 + @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 2.5 rtol = 2e-1 + end + end + end +end +