Mercurial > repos > public > sbplib_julia
changeset 958:8b0ff2fddc32 feature/variable_derivatives
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 14 Mar 2022 09:08:23 +0100 |
| parents | b15f39ae1643 (diff) acce83640536 (current diff) |
| children | e9752c1e92f8 |
| files | src/LazyTensors/lazy_tensor_operations.jl |
| diffstat | 17 files changed, 901 insertions(+), 50 deletions(-) [+] |
line wrap: on
line diff
--- a/Notes.md Mon Mar 14 08:37:28 2022 +0100 +++ b/Notes.md Mon Mar 14 09:08:23 2022 +0100 @@ -148,12 +148,52 @@ - [ ] Can we have a trait to tell if a TensorMapping is transposable? - [ ] 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? +## 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 TensorMapping 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 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 Mon Mar 14 08:37:28 2022 +0100 +++ b/TODO.md Mon Mar 14 09:08:23 2022 +0100 @@ -1,15 +1,12 @@ # TODO -## Skämskudde + - [ ] Ändra namn på variabler och funktioner så att det följer style-guide - - [ ] Skriv tester - -## Coding - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?) - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation. - - [ ] Add a quick and simple way of running all tests for all subpackages. - [ ] 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 - [ ] Add check for correct domain sizes to lazy tensor operations using SizeMismatch - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index @@ -23,11 +20,10 @@ - [ ] 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 -## 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/lazy_tensor_operations.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Mon Mar 14 09:08:23 2022 +0100 @@ -431,3 +431,28 @@ print(io, "SizeMismatch: ") print(io, "domain size $(domain_size(err.tm)) of TensorMapping not matching size $(err.sz)") end + +function apply_with_region(op, v, boundary_width::Integer, dim_size::Integer, i) + if 0 < i <= boundary_width + return LazyTensors.apply(op,v,Index(i,Lower)) + elseif boundary_width < i <= dim_size-boundary_width + return LazyTensors.apply(op,v,Index(i,Interior)) + elseif dim_size-boundary_width < i <= dim_size + return LazyTensors.apply(op,v,Index(i,Upper)) + else + error("Bounds error") # TODO: Make this more standard + end +end +# TBD: Can these methods be merge by having a function as an arguement instead? +# TODO: Add inference test that show where things break and how this rewrite fixes it. +function apply_transpose_with_region(op, v, boundary_width::Integer, dim_size::Integer, i) + if 0 < i <= boundary_width + return LazyTensors.apply_transpose(op,v,Index(i,Lower)) + elseif boundary_width < i <= dim_size-boundary_width + return LazyTensors.apply_transpose(op,v,Index(i,Interior)) + elseif dim_size-boundary_width < i <= dim_size + return LazyTensors.apply_transpose(op,v,Index(i,Upper)) + else + error("Bounds error") # TODO: Make this more standard + end +end
--- a/src/RegionIndices/RegionIndices.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/src/RegionIndices/RegionIndices.jl Mon Mar 14 09:08:23 2022 +0100 @@ -65,6 +65,7 @@ error("Bounds error") # TODO: Make this more standard end end +# 2022-02-21: Using the return values of getregion cause type inference to give up in ceratin cases for example H*H*v export getregion
--- a/src/SbpOperators/SbpOperators.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/SbpOperators.jl Mon Mar 14 09:08:23 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/boundaryops/boundary_operator.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Mon Mar 14 09:08:23 2022 +0100 @@ -84,6 +84,5 @@ end function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T - r = getregion(i, closure_size(op), op.size) - apply_transpose(op, v, Index(i,r)) + return LazyTensors.apply_transpose_with_region(op, v, closure_size(op), op.size[1], i) end
--- a/src/SbpOperators/operators/standard_diagonal.toml Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/operators/standard_diagonal.toml Mon Mar 14 09:08:23 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 Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/readoperator.jl Mon Mar 14 09:08:23 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 Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/stencil.jl Mon Mar 14 09:08:23 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
--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Mon Mar 14 09:08:23 2022 +0100 @@ -41,8 +41,7 @@ end function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T - r = getregion(i, closure_size(op), op.size[1]) - return LazyTensors.apply(op, v, Index(i, r)) + return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i) end LazyTensors.apply_transpose(op::ConstantInteriorScalingOperator, v, i) = apply(op, v, i)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Mon Mar 14 09:08:23 2022 +0100 @@ -0,0 +1,195 @@ +export SecondDerivativeVariable + +# REVIEW: Fixa docs +""" + SecondDerivativeVariable{Dir,T,D,...} <: TensorMapping{T,D,D} + +A second derivative operator in direction `Dir` with a variable coefficient. +""" +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) + check_coefficient(grid, coeff) + + Δ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::NestedStencil, closure_stencils) + return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1) +end + +@doc raw""" + SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir) + +Create a `TensorMapping` for the second derivative with a variable coefficient +`coeff` on `grid` from the stencils in `stencil_set`. The direction is +determined by `dir`. + +`coeff` is a grid function on `grid`. + +# Example +With +``` +D = SecondDerivativeVariable(g, c, stencil_set, 2) +``` +then `D*u` approximates +```math +\frac{\partial}{\partial y} c(x,y) \frac{\partial u}{\partial y}, +``` +on ``(0,1)⨯(0,1)`` represented by `g`. +""" +function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir::Int) + 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 + +function check_coefficient(grid, coeff) + if dimension(grid) != ndims(coeff) + throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(dimension(grid))")) + end + + if size(grid) != size(coeff) + throw(DimensionMismatch("the size $(size(coeff)) of the coefficient does not match the size $(size(grid)) of the grid")) + end +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...] +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)...) + elseif I[derivative_direction(op)] isa Index{Interior} + return apply_interior(op, v, Int.(I)...) + else + error("Invalid region") + end +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...) + dir = derivative_direction(op) + + i = I[dir] + + I = map(i->Index(i, Interior), I) + if 0 < i <= closure_size(op) + I = Base.setindex(I, Index(i, Lower), dir) + return LazyTensors.apply(op, v, I...) + elseif closure_size(op) < i <= op.size[dir]-closure_size(op) + I = Base.setindex(I, Index(i, Interior), dir) + return LazyTensors.apply(op, v, I...) + elseif op.size[dir]-closure_size(op) < i <= op.size[dir] + I = Base.setindex(I, Index(i, Upper), dir) + return LazyTensors.apply(op, v, I...) + else + error("Bounds error") # TODO: Make this more standard + end +end + + +## 2D Specific implementations to avoid instability +## TODO: Should really be solved by fixing the general methods instead + + +## x-direction +function apply_lower(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i) +end + +function apply_interior(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, i) +end + +function apply_upper(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) +end + + +## y-direction +function apply_lower(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + return @inbounds apply_stencil(op.closure_stencils[j], c̃, ṽ, j) +end + +function apply_interior(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, j) +end + +function apply_upper(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + stencil = op.closure_stencils[op.size[derivative_direction(op)]-j+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, j) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/inference_trouble.txt Mon Mar 14 09:08:23 2022 +0100 @@ -0,0 +1,53 @@ +Innan ändringarna på den här branchen: + +begin + using Sbplib + using Sbplib.Grids + using Sbplib.SbpOperators + + g = EquidistantGrid((10,10),(0.,0.), (1.,1.)) + v = evalOn(g, (x,y)->x^2+y^2+1) + H = inner_product(g, 1., [1/2]) +end + +# Not type stable +LazyTensors.apply(H.t1, H.t2*v, 1,2) +@code_warntype LazyTensors.apply(H.t1, H.t2*v, 1,2) +@code_warntype LazyTensors.apply(H.t1.tm, view(H.t2*v,:,1), 2) + +# Nedan är halvdåliga +@code_warntype LazyTensors.apply(H.t1.tm, view(H.t2*v,1,:), 2) +@code_warntype LazyTensors.apply(H.t1.tm, view(v,1,:), 2) +@code_warntype LazyTensors.apply(H.t1.tm, view(v,:,1), 2) +@code_warntype LazyTensors.apply(H.t1.tm, v[:,1], 2) + + + + + + + + + +begin + using Sbplib + using Sbplib.Grids + using Sbplib.SbpOperators + import Sbplib.SbpOperators: Stencil + using Sbplib.RegionIndices + + g = EquidistantGrid(10,0., 1.) + v = evalOn(g, (x)->x^2+1) + H = inner_product(g, 1., [1/2]) + V = SbpOperators.volume_operator(g, Stencil(1.,center=1), (Stencil(1/2,center=1),), SbpOperators.even,1) + b = SbpOperators.boundary_operator(g, Stencil(1/2,center=1), CartesianBoundary{1,Lower}()) +end + +@code_warntype LazyTensors.apply(H, H*v, 2) +@code_warntype LazyTensors.apply(V, V*v, 2) +@code_warntype LazyTensors.apply(b, b*v, 2) + + + +begin +end
--- a/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 14 09:08:23 2022 +0100 @@ -56,6 +56,5 @@ end function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T - r = getregion(i, closure_size(op), op.size[1]) - return LazyTensors.apply(op, v, Index(i, r)) + return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i) end
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 14 08:37:28 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 14 09:08:23 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 Mon Mar 14 08:37:28 2022 +0100 +++ b/test/SbpOperators/readoperator_test.jl Mon Mar 14 09:08:23 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 Mon Mar 14 08:37:28 2022 +0100 +++ b/test/SbpOperators/stencil_test.jl Mon Mar 14 09:08:23 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 Mon Mar 14 09:08:23 2022 +0100 @@ -0,0 +1,187 @@ +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 + + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test SecondDerivativeVariable(g, c, stencil_set, 1) isa SecondDerivativeVariable + + @testset "checking c" begin + c_short = rand(5) + c_long = rand(16) + c_higher_dimension = rand(11,11) + + @test_throws DimensionMismatch("the size (5,) of the coefficient does not match the size (11,) of the grid") SecondDerivativeVariable(g, c_short, interior_stencil, closure_stencils) + @test_throws DimensionMismatch("the size (16,) of the coefficient does not match the size (11,) of the grid") SecondDerivativeVariable(g, c_long, interior_stencil, closure_stencils) + @test_throws ArgumentError("The coefficient has dimension 2 while the grid is dimension 1") SecondDerivativeVariable(g, c_higher_dimension, interior_stencil, closure_stencils,1) + end + 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 "type stability" 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 + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test SecondDerivativeVariable(g, c, stencil_set, 1) isa SecondDerivativeVariable + 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) + + + @testset "standard diagonal 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((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 +