Mercurial > repos > public > sbplib_julia
changeset 1091:e3b41d48b5aa feature/variable_derivatives
Merge refactor/grids
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 26 Apr 2022 14:29:08 +0200 |
| parents | de972d825289 (diff) 9b40aeac4269 (current diff) |
| children | 703eaa3e50c4 |
| files | |
| diffstat | 21 files changed, 1149 insertions(+), 26 deletions(-) [+] |
line wrap: on
line diff
--- a/src/LazyTensors/lazy_tensor_operations.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/src/LazyTensors/lazy_tensor_operations.jl Tue Apr 26 14:29:08 2022 +0200 @@ -1,3 +1,5 @@ +using Sbplib.RegionIndices + """ TensorApplication{T,R,D} <: LazyArray{T,R} @@ -268,6 +270,20 @@ LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms) + +""" + inflate(tm, sz, dir) + +Inflate `tm` with identity tensors in all directions `d` for `d != dir`. + +# TODO: Describe when it is useful +""" +function inflate(tm::LazyTensor, sz, dir) + Is = IdentityTensor{eltype(tm)}.(sz) + parts = Base.setindex(Is, tm, dir) + return foldl(⊗, parts) +end + function check_domain_size(tm::LazyTensor, sz) if domain_size(tm) != sz throw(DomainSizeMismatch(tm,sz)) @@ -290,7 +306,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 +315,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/LazyTensors/tuple_manipulation.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/src/LazyTensors/tuple_manipulation.jl Tue Apr 26 14:29:08 2022 +0200 @@ -74,3 +74,23 @@ flatten_tuple(t::NTuple{N, Number} where N) = t flatten_tuple(t::Tuple) = ((flatten_tuple.(t)...)...,) # simplify? flatten_tuple(ts::Vararg) = flatten_tuple(ts) + + +function left_pad_tuple(t, val, N) + if N < length(t) + throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements")) + end + + padding = ntuple(i->val, N-length(t)) + return (padding..., t...) +end + +function right_pad_tuple(t, val, N) + if N < length(t) + throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements")) + end + + padding = ntuple(i->val, N-length(t)) + return (t..., padding...) +end +
--- a/src/RegionIndices/RegionIndices.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/src/RegionIndices/RegionIndices.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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 Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/SbpOperators.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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 @@ -19,6 +20,7 @@ export normal_derivative export first_derivative export second_derivative +export undivided_dissipation using Sbplib.RegionIndices using Sbplib.LazyTensors @@ -29,12 +31,17 @@ even = 1 end +export closure_size + include("stencil.jl") include("stencil_set.jl") include("volumeops/volume_operator.jl") +include("volumeops/stencil_operator_distinct_closures.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/derivatives/dissipation.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 Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Tue Apr 26 14:29:08 2022 +0200 @@ -12,19 +12,14 @@ function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary) #TODO:Check that dim(boundary) <= Dim? - # Create 1D boundary operator - r = region(boundary) d = dim(boundary) - op = BoundaryOperator(restrict(grid, d), closure_stencil, r) + op = BoundaryOperator(restrict(grid, d), closure_stencil, region(boundary)) # Create 1D IdentityTensors for each coordinate direction one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids)) - # Formulate the correct outer product sequence of the identity mappings and - # the boundary operator - parts = Base.setindex(Is, op, d) - return foldl(⊗, parts) + return LazyTensors.inflate(op, size(grid), d) end """ @@ -84,6 +79,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 Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/operators/standard_diagonal.toml Tue Apr 26 14:29:08 2022 +0200 @@ -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 Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/stencil.jl Tue Apr 26 14:29:08 2022 +0200 @@ -85,6 +85,21 @@ return w end +function left_pad(s::Stencil, N) + weights = LazyTensors.left_pad_tuple(s.weights, zero(eltype(s)), N) + range = (first(s.range) - (N - length(s.weights))):last(s.range) + + return Stencil(range, weights) +end + +function right_pad(s::Stencil, N) + weights = LazyTensors.right_pad_tuple(s.weights, zero(eltype(s)), N) + range = first(s.range):(last(s.range) + (N - length(s.weights))) + + return Stencil(range, weights) +end + + struct NestedStencil{T,N,M} s::Stencil{Stencil{T,N},M}
--- a/src/SbpOperators/stencil_set.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/stencil_set.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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 Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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/dissipation.jl Tue Apr 26 14:29:08 2022 +0200 @@ -0,0 +1,87 @@ +function undivided_dissipation(g::EquidistantGrid, p, direction) + T = eltype(g) + interior_weights = T.(dissipation_interior_weights(p)) + + D = stencil_operator_distinct_closures( + g, + dissipation_interior_stencil(interior_weights), + dissipation_lower_closure_stencils(interior_weights), + dissipation_upper_closure_stencils(interior_weights), + direction, + ) + Dᵀ = stencil_operator_distinct_closures( + g, + dissipation_transpose_interior_stencil(interior_weights), + dissipation_transpose_lower_closure_stencils(interior_weights), + dissipation_transpose_upper_closure_stencils(interior_weights), + direction, + ) + + return D, Dᵀ +end + +undivided_dissipation(g::EquidistantGrid{1}, p) = undivided_dissipation(g, p, 1) + +function dissipation_interior_weights(p) + if p == 0 + return (1,) + end + + return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0) +end + +midpoint(weights) = length(weights)÷2 + 1 +midpoint_transpose(weights) = length(weights)+1 - midpoint(weights) + +function dissipation_interior_stencil(weights) + return Stencil(weights..., center=midpoint(weights)) +end +function dissipation_transpose_interior_stencil(weights) + if iseven(length(weights)) + weights = map(-, weights) + end + + return Stencil(weights..., center=midpoint_transpose(weights)) +end + +dissipation_lower_closure_size(weights) = midpoint(weights) - 1 +dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights) + +dissipation_lower_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=i ), dissipation_lower_closure_size(interior_weights)) +dissipation_upper_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=length(interior_weights)-dissipation_upper_closure_size(interior_weights)+i), dissipation_upper_closure_size(interior_weights)) + +function dissipation_transpose_lower_closure_stencils(interior_weights) + closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights)) + + N = maximum(s->length(s.weights), closure) + return right_pad.(closure, N) +end + +function dissipation_transpose_upper_closure_stencils(interior_weights) + closure = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights))) + + N = maximum(s->length(s.weights), closure) + return left_pad.(closure, N) +end + + +function dissipation_transpose_lower_closure_stencil(interior_weights, i) + w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights)) + + for k ∈ i:-1:1 + w = (w..., interior_weights[k]) + end + + return Stencil(w..., center = i) +end + +function dissipation_transpose_upper_closure_stencil(interior_weights, i) + j = length(interior_weights)+1-i + w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights)) + + for k ∈ j:1:length(interior_weights) + w = (interior_weights[k], w...) + end + + return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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 Tue Apr 26 14:29:08 2022 +0200 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl Tue Apr 26 14:29:08 2022 +0200 @@ -0,0 +1,58 @@ +function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction) + op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure) + return LazyTensors.inflate(op, size(grid), direction) +end + +""" + StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1} + +A one dimensional stencil operator with separate closures for the two boundaries. + +`StencilOperatorDistinctClosures` can be contrasted to `VolumeOperator` in +that it has different closure stencils for the upper and lower boundary. +`VolumeOperator` uses the same closure for both boundaries. Having distinct +closures is useful for representing operators with skewed stencils like upwind +operators. + +See also: [`VolumeOperator`](@ref) +""" +struct StencilOperatorDistinctClosures{T,K,N,M,LC<:NTuple{N,Stencil{T,L}} where L, UC<:NTuple{M,Stencil{T,L}} where L} <: LazyTensor{T,1,1} + inner_stencil::Stencil{T,K} + lower_closure::LC + upper_closure::UC + size::Tuple{Int} +end + +function StencilOperatorDistinctClosures(grid::EquidistantGrid{1}, inner_stencil, lower_closure, upper_closure) + return StencilOperatorDistinctClosures(inner_stencil, Tuple(lower_closure), Tuple(upper_closure), size(grid)) +end + +lower_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = N +upper_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = M + +LazyTensors.range_size(op::StencilOperatorDistinctClosures) = op.size +LazyTensors.domain_size(op::StencilOperatorDistinctClosures) = op.size + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Lower}) + return @inbounds apply_stencil(op.lower_closure[Int(i)], v, Int(i)) +end + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Interior}) + return apply_stencil(op.inner_stencil, v, Int(i)) +end + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Upper}) + stencil_index = Int(i) - (op.size[1]-upper_closure_size(op)) + return @inbounds apply_stencil(op.upper_closure[stencil_index], v, Int(i)) +end + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i) + if i <= lower_closure_size(op) + LazyTensors.apply(op, v, Index(i, Lower)) + elseif i > op.size[1]-upper_closure_size(op) + LazyTensors.apply(op, v, Index(i, Upper)) + else + LazyTensors.apply(op, v, Index(i, Interior)) + end +end +# TODO: Move this to LazyTensors when we have the region communication down.
--- a/src/SbpOperators/volumeops/volume_operator.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/src/SbpOperators/volumeops/volume_operator.jl Tue Apr 26 14:29:08 2022 +0200 @@ -12,16 +12,10 @@ function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction) #TODO: Check that direction <= Dim? - # Create 1D volume operator in along coordinate direction op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) - # Create 1D IdentityTensors for each coordinate direction - one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) - Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids)) - # Formulate the correct outer product sequence of the identity mappings and - # the volume operator - parts = Base.setindex(Is, op, direction) - return foldl(⊗, parts) + return LazyTensors.inflate(op, size(grid), direction) end +# TBD: Should the inflation happen here or should we remove this method and do it at the caller instead? """ VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1} @@ -56,6 +50,6 @@ 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 +# TODO: Move this to LazyTensors when we have the region communication down.
--- a/test/LazyTensors/lazy_tensor_operations_test.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -358,3 +358,17 @@ @test I1⊗Ã⊗I2 == InflatedTensor(I1, Ã, I2) end end + +@testset "inflate" begin + I = LazyTensors.inflate(IdentityTensor(),(3,4,5,6), 2) + @test I isa LazyTensor{Float64, 3,3} + @test range_size(I) == (3,5,6) + @test domain_size(I) == (3,5,6) + + @test LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 1) == InflatedTensor(IdentityTensor{Float64}(),ScalingTensor(1., (4,)),IdentityTensor(4,5,6)) + @test LazyTensors.inflate(ScalingTensor(2., (1,)),(3,4,5,6), 2) == InflatedTensor(IdentityTensor(3),ScalingTensor(2., (1,)),IdentityTensor(5,6)) + @test LazyTensors.inflate(ScalingTensor(3., (6,)),(3,4,5,6), 4) == InflatedTensor(IdentityTensor(3,4,5),ScalingTensor(3., (6,)),IdentityTensor{Float64}()) + + @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 0) + @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 5) +end
--- a/test/LazyTensors/tuple_manipulation_test.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/test/LazyTensors/tuple_manipulation_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -60,3 +60,19 @@ @test LazyTensors.flatten_tuple((1,2,(3,(4,5)),6)) == (1,2,3,4,5,6) @test LazyTensors.flatten_tuple(((1,2),(3,4),(5,),6)) == (1,2,3,4,5,6) end + +@testset "left_pad_tuple" begin + @test LazyTensors.left_pad_tuple((1,2), 0, 2) == (1,2) + @test LazyTensors.left_pad_tuple((1,2), 0, 3) == (0,1,2) + @test LazyTensors.left_pad_tuple((3,2), 1, 6) == (1,1,1,1,3,2) + + @test_throws DomainError(0, "Can't pad tuple of length 2 to 0 elements") LazyTensors.left_pad_tuple((1,2), 0, 0) == (1,2) +end + +@testset "right_pad_tuple" begin + @test LazyTensors.right_pad_tuple((1,2), 0, 2) == (1,2) + @test LazyTensors.right_pad_tuple((1,2), 0, 3) == (1,2,0) + @test LazyTensors.right_pad_tuple((3,2), 1, 6) == (3,2,1,1,1,1) + + @test_throws DomainError(0, "Can't pad tuple of length 2 to 0 elements") LazyTensors.right_pad_tuple((1,2), 0, 0) == (1,2) +end
--- a/test/SbpOperators/stencil_set_test.jl Tue Apr 26 14:28:42 2022 +0200 +++ b/test/SbpOperators/stencil_set_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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 Apr 26 14:28:42 2022 +0200 +++ b/test/SbpOperators/stencil_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -62,6 +62,23 @@ end end +@testset "left_pad" begin + @test SbpOperators.left_pad(Stencil(1,1, center = 1), 2) == Stencil(1,1, center=1) + @test SbpOperators.left_pad(Stencil(1,1, center = 1), 3) == Stencil(0,1,1, center=2) + @test SbpOperators.left_pad(Stencil(2,3, center = 2), 4) == Stencil(0,0,2,3, center=4) + + @test SbpOperators.left_pad(Stencil(2.,3., center = 2), 4) == Stencil(0.,0.,2.,3., center=4) +end + +@testset "right_pad" begin + @test SbpOperators.right_pad(Stencil(1,1, center = 1), 2) == Stencil(1,1, center=1) + @test SbpOperators.right_pad(Stencil(1,1, center = 1), 3) == Stencil(1,1,0, center=1) + @test SbpOperators.right_pad(Stencil(2,3, center = 2), 4) == Stencil(2,3,0,0, center=2) + + @test SbpOperators.right_pad(Stencil(2.,3., center = 2), 4) == Stencil(2.,3.,0.,0., center=2) +end + + @testset "NestedStencil" begin @testset "Constructors" begin @@ -170,5 +187,4 @@ @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/dissipation_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -0,0 +1,217 @@ +using Test + +using Sbplib.SbpOperators +using Sbplib.Grids +using Sbplib.LazyTensors + +using Sbplib.SbpOperators: Stencil + +using Sbplib.SbpOperators: dissipation_interior_weights +using Sbplib.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil +using Sbplib.SbpOperators: midpoint, midpoint_transpose +using Sbplib.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size +using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils +using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils + +""" + monomial(x,k) + +Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``. +Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)`` +""" +function monomial(x,k) + if k < 0 + return zero(x) + end + x^k/factorial(k) +end + +@testset "undivided_dissipation" begin + g = EquidistantGrid(20, 0., 11.) + D,Dᵀ = undivided_dissipation(g, 1) + + @test D isa LazyTensor{Float64,1,1} where T + @test Dᵀ isa LazyTensor{Float64,1,1} where T + + @testset "Accuracy conditions" begin + N = 20 + g = EquidistantGrid(N, 0//1,2//1) + h = only(spacing(g)) + @testset "D_$p" for p ∈ [1,2,3,4] + D,Dᵀ = undivided_dissipation(g, p) + + @testset "x^$k" for k ∈ 0:p + v = evalOn(g, x->monomial(x,k)) + vₚₓ = evalOn(g, x->monomial(x,k-p)) + + @test D*v == h^p * vₚₓ + end + end + end + + @testset "tanspose equality" begin + function get_matrix(D) + N = only(range_size(D)) + M = only(domain_size(D)) + + Dmat = zeros(N,M) + e = zeros(M) + for i ∈ 1:M + if i > 1 + e[i-1] = 0. + end + e[i] = 1. + Dmat[:,i] = D*e + end + + return Dmat + end + + g = EquidistantGrid(11, 0., 1.) + @testset "D_$p" for p ∈ [1,2,3,4] + D,Dᵀ = undivided_dissipation(g, p) + + D̄ = get_matrix(D) + D̄ᵀ = get_matrix(Dᵀ) + + @test D̄ == D̄ᵀ' + end + end +end + +@testset "dissipation_interior_weights" begin + @test dissipation_interior_weights(1) == (-1, 1) + @test dissipation_interior_weights(2) == (1,-2, 1) + @test dissipation_interior_weights(3) == (-1, 3,-3, 1) + @test dissipation_interior_weights(4) == (1, -4, 6, -4, 1) +end + +@testset "dissipation_interior_stencil" begin + @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1, 1, center=2) + @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil( 1,-2, 1, center=2) + @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1, 3,-3, 1, center=3) + @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil( 1,-4, 6,-4, 1, center=3) +end + +@testset "dissipation_transpose_interior_stencil" begin + @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(1,-1, center=1) + @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2, 1, center=2) + @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(1,-3, 3,-1, center=2) + @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1,-4, 6,-4, 1, center=3) +end + +@testset "midpoint" begin + @test midpoint((1,1)) == 2 + @test midpoint((1,1,1)) == 2 + @test midpoint((1,1,1,1)) == 3 + @test midpoint((1,1,1,1,1)) == 3 +end + +@testset "midpoint_transpose" begin + @test midpoint_transpose((1,1)) == 1 + @test midpoint_transpose((1,1,1)) == 2 + @test midpoint_transpose((1,1,1,1)) == 2 + @test midpoint_transpose((1,1,1,1,1)) == 3 +end + +@testset "dissipation_lower_closure_size" begin + @test dissipation_lower_closure_size((1,1)) == 1 + @test dissipation_lower_closure_size((1,1,1)) == 1 + @test dissipation_lower_closure_size((1,1,1,1)) == 2 + @test dissipation_lower_closure_size((1,1,1,1,1)) == 2 +end + +@testset "dissipation_upper_closure_size" begin + @test dissipation_upper_closure_size((1,1)) == 0 + @test dissipation_upper_closure_size((1,1,1)) == 1 + @test dissipation_upper_closure_size((1,1,1,1)) == 1 + @test dissipation_upper_closure_size((1,1,1,1,1)) == 2 +end + +@testset "dissipation_lower_closure_stencils" begin + cases = ( + (-1,1) => ( + Stencil(-1, 1, center=1), + ), + (1,-2,1) => ( + Stencil( 1,-2, 1, center=1), + ), + (-1,3,-3,1) => ( + Stencil(-1,3,-3,1, center=1), + Stencil(-1,3,-3,1, center=2), + ), + (1, -4, 6, -4, 1) => ( + Stencil(1, -4, 6, -4, 1, center=1), + Stencil(1, -4, 6, -4, 1, center=2), + ) + ) + @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases + @test dissipation_lower_closure_stencils(w) == closure_stencils + end +end + +@testset "dissipation_upper_closure_stencils" begin + cases = ( + (-1,1) => (), + (1,-2,1) => ( + Stencil( 1,-2, 1, center=3), + ), + (-1,3,-3,1) => ( + Stencil(-1,3,-3,1, center=4), + ), + (1, -4, 6, -4, 1) => ( + Stencil(1, -4, 6, -4, 1, center=4), + Stencil(1, -4, 6, -4, 1, center=5), + ) + ) + @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases + @test dissipation_upper_closure_stencils(w) == closure_stencils + end +end + + +@testset "dissipation_transpose_lower_closure_stencils" begin + cases = ( + (-1,1) => ( + Stencil(-1,-1, 0, center=1), + Stencil( 1, 1,-1, center=2), + ), + (1,-2,1) => ( + Stencil( 1, 1, 0, 0, center=1), + Stencil(-2,-2, 1, 0, center=2), + Stencil( 1, 1,-2, 1, center=3), + ), + (-1,3,-3,1) => ( + Stencil(-1,-1,-1, 0, 0, 0, center=1), + Stencil( 3, 3, 3,-1, 0, 0, center=2), + Stencil(-3,-3,-3, 3,-1, 0, center=3), + Stencil( 1, 1, 1,-3, 3,-1, center=4), + ), + ) + @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases + @test dissipation_transpose_lower_closure_stencils(w) == closure_stencils + end +end + +@testset "dissipation_transpose_upper_closure_stencils" begin + cases = ( + (-1,1) => ( + Stencil( 1,-1, center = 1), + Stencil( 0, 1, center = 2), + ), + (1,-2,1) => ( + Stencil( 1,-2, 1, 1, center=2), + Stencil( 0, 1,-2,-2, center=3), + Stencil( 0, 0, 1, 1, center=4), + ), + (-1,3,-3,1) => ( + Stencil( 1,-3, 3,-1,-1, center=2), + Stencil( 0, 1,-3, 3, 3, center=3), + Stencil( 0, 0, 1,-3,-3, center=4), + Stencil( 0, 0, 0, 1, 1, center=5), + ), + ) + @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases + @test dissipation_transpose_upper_closure_stencils(w) == closure_stencils + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl Tue Apr 26 14:29:08 2022 +0200 @@ -0,0 +1,83 @@ +using Test + +using Sbplib.SbpOperators +using Sbplib.Grids +# using Sbplib.RegionIndices +using Sbplib.LazyTensors + +import Sbplib.SbpOperators.Stencil +import Sbplib.SbpOperators.StencilOperatorDistinctClosures +import Sbplib.SbpOperators.stencil_operator_distinct_closures + +@testset "stencil_operator_distinct_closures" begin + lower_closure = ( + Stencil(-1,1, center=1), + ) + + inner_stencil = Stencil(-2,2, center=1) + + upper_closure = ( + Stencil(-3,3, center=1), + Stencil(-4,4, center=2), + ) + + g₁ = EquidistantGrid(5, 0., 1.) + g₂ = EquidistantGrid((5,5), (0.,0.), (1.,1.)) + h = 1/4 + + A₁ = stencil_operator_distinct_closures(g₁, inner_stencil, lower_closure, upper_closure, 1) + A₂¹ = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 1) + A₂² = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 2) + + v₁ = evalOn(g₁, x->x) + + u = [1., 2., 2., 3., 4.]*h + @test A₁*v₁ == u + + v₂ = evalOn(g₂, (x,y)-> x + 3y) + @test A₂¹*v₂ == repeat(u, 1, 5) + @test A₂²*v₂ == repeat(3u', 5, 1) +end + +@testset "StencilOperatorDistinctClosures" begin + g = EquidistantGrid(11, 0., 1.) + + lower_closure = ( + Stencil(-1,1, center=1), + Stencil(-2,2, center=2), + ) + + inner_stencil = Stencil(-3,3, center=1) + + upper_closure = ( + Stencil(4,-4,4, center=1), + Stencil(5,-5,5, center=2), + Stencil(6,-6,6, center=3), + ) + + A = StencilOperatorDistinctClosures(g, inner_stencil, lower_closure, upper_closure) + @test A isa LazyTensor{T,1,1} where T + + @test SbpOperators.lower_closure_size(A) == 2 + @test SbpOperators.upper_closure_size(A) == 3 + + @test domain_size(A) == (11,) + @test range_size(A) == (11,) + + v = rand(11) + @testset "apply" begin + # Lower closure + @test LazyTensors.apply(A, v, 1) ≈ 1*(-v[1] + v[2]) + @test LazyTensors.apply(A, v, 2) ≈ 2*(-v[1] + v[2]) + + # Interior + @test LazyTensors.apply(A, v, 3) ≈ 3*(-v[3] + v[4]) + @test LazyTensors.apply(A, v, 4) ≈ 3*(-v[4] + v[5]) + @test LazyTensors.apply(A, v, 8) ≈ 3*(-v[8] + v[9]) + + # Upper closure + @test LazyTensors.apply(A, v, 9) ≈ 4*(v[9] - v[10] + v[11]) + @test LazyTensors.apply(A, v, 10) ≈ 5*(v[9] - v[10] + v[11]) + @test LazyTensors.apply(A, v, 11) ≈ 6*(v[9] - v[10] + v[11]) + end +end