Mercurial > repos > public > sbplib_julia
changeset 1049:3bb94ce74697 feature/variable_derivatives
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 23 Mar 2022 12:54:45 +0100 |
| parents | 86aa69ad3304 (diff) c16116e403e2 (current diff) |
| children | eeecdf135912 |
| files | Notes.md TODO.md src/LazyTensors/lazy_tensor_operations.jl src/SbpOperators/SbpOperators.jl src/SbpOperators/boundaryops/boundary_operator.jl src/SbpOperators/readoperator.jl src/SbpOperators/stencil_set.jl src/SbpOperators/volumeops/constant_interior_scaling_operator.jl src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl src/SbpOperators/volumeops/volume_operator.jl test/SbpOperators/boundaryops/boundary_operator_test.jl test/SbpOperators/readoperator_test.jl test/SbpOperators/stencil_set_test.jl test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl |
| diffstat | 17 files changed, 904 insertions(+), 49 deletions(-) [+] |
line wrap: on
line diff
--- a/Notes.md Tue Mar 22 14:33:13 2022 +0100 +++ b/Notes.md Wed Mar 23 12:54:45 2022 +0100 @@ -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 Mar 22 14:33:13 2022 +0100 +++ b/TODO.md Wed Mar 23 12:54:45 2022 +0100 @@ -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/lazy_tensor_operations.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Wed Mar 23 12:54:45 2022 +0100 @@ -1,3 +1,5 @@ +using Sbplib.RegionIndices + """ TensorApplication{T,R,D} <: LazyArray{T,R} @@ -290,7 +292,6 @@ print(io, "domain size $(domain_size(err.tm)) of LazyTensor not matching size $(err.sz)") end - struct RangeSizeMismatch <: Exception tm::LazyTensor sz @@ -300,3 +301,29 @@ print(io, "RangeSizeMismatch: ") print(io, "range size $(range_size(err.tm)) of LazyTensor 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 \ No newline at end of file
--- a/src/RegionIndices/RegionIndices.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/src/RegionIndices/RegionIndices.jl Wed Mar 23 12:54:45 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 certain cases, for example H*H*v export getregion
--- a/src/SbpOperators/SbpOperators.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/SbpOperators.jl Wed Mar 23 12:54:45 2022 +0100 @@ -5,6 +5,7 @@ export read_stencil_set export get_stencil_set export parse_stencil +export parse_nested_stencil export parse_scalar export parse_tuple export sbp_operators_path @@ -29,12 +30,15 @@ even = 1 end +export closure_size + include("stencil.jl") include("stencil_set.jl") include("volumeops/volume_operator.jl") include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/first_derivative.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 Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Wed Mar 23 12:54:45 2022 +0100 @@ -84,6 +84,5 @@ end function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i) - 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 Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/operators/standard_diagonal.toml Wed Mar 23 12:54:45 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/stencil.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/stencil.jl Wed Mar 23 12:54:45 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/stencil_set.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/stencil_set.jl Wed Mar 23 12:54:45 2022 +0100 @@ -110,6 +110,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/volumeops/constant_interior_scaling_operator.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Wed Mar 23 12:54:45 2022 +0100 @@ -41,8 +41,7 @@ end function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i) - 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 Wed Mar 23 12:54:45 2022 +0100 @@ -0,0 +1,195 @@ +export SecondDerivativeVariable + +# REVIEW: Fixa docs +""" + SecondDerivativeVariable{Dir,T,D,...} <: LazyTensor{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} <: LazyTensor{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 `LazyTensor` 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 Wed Mar 23 12:54:45 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 Tue Mar 22 14:33:13 2022 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Wed Mar 23 12:54:45 2022 +0100 @@ -56,6 +56,5 @@ end function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i) - 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 Tue Mar 22 14:33:13 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Wed Mar 23 12:54:45 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/stencil_set_test.jl Tue Mar 22 14:33:13 2022 +0100 +++ b/test/SbpOperators/stencil_set_test.jl Wed Mar 23 12:54:45 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 Mar 22 14:33:13 2022 +0100 +++ b/test/SbpOperators/stencil_test.jl Wed Mar 23 12:54:45 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 Wed Mar 23 12:54:45 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 LazyTensor + + 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 LazyTensor + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) isa LazyTensor + + 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 +