Mercurial > repos > public > sbplib_julia
changeset 982:2a4f36aca2ea feature/variable_derivatives
Merge feature/variable_derivatives
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 15 Mar 2022 21:42:52 +0100 |
| parents | df562695b1b5 (diff) b90446eb5f27 (current diff) |
| children | 86aa69ad3304 |
| files | src/SbpOperators/SbpOperators.jl |
| diffstat | 22 files changed, 1025 insertions(+), 95 deletions(-) [+] |
line wrap: on
line diff
--- a/Notes.md Tue Mar 15 21:30:56 2022 +0100 +++ b/Notes.md Tue Mar 15 21:42:52 2022 +0100 @@ -149,12 +149,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 TensorMapping uses indices with regions. The default should be that they do NOT. - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle? - [ ] Figure out repeated application of regioned TensorMappings. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work? +### Ideas for information sharing functions +```julia +using StaticArrays + +function regions(op::SecondDerivativeVariable) + t = ntuple(i->(Interior(),),range_dim(op)) + return Base.setindex(t, (Lower(), Interior(), Upper()), derivative_direction(op)) +end + +function regionsizes(op::SecondDerivativeVariable) + sz = tuple.(range_size(op)) + + cl = closuresize(op) + return Base.setindex(sz, (cl, n-2cl, cl), derivative_direction(op)) +end + + +g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1 +c = evalOn(g, (x,y)->x+y) + +D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) +@test regions(D₂ᶜ) == ( + (Lower(), Interior(), Upper()), + (Interior(),), +) +@test regionsizes(D₂ᶜ) == ((1,9,1),(9,)) + + +D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) +@test regions(D₂ᶜ) == ( + (Interior(),), + (Lower(), Interior(), Upper()), +) +@test regionsizes(D₂ᶜ) == ((11,),(1,7,1)) +``` + + ## Boundschecking and dimension checking Does it make sense to have boundschecking only in getindex methods? This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway.
--- a/TODO.md Tue Mar 15 21:30:56 2022 +0100 +++ b/TODO.md Tue Mar 15 21:42:52 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 - [ ] Add check for correct domain sizes to lazy tensor operations using SizeMismatch - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index @@ -22,13 +21,12 @@ - [ ] 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 - [ ] Add a type StencilSet for easier dispatch -## 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/LazyTensors.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/LazyTensors/LazyTensors.jl Tue Mar 15 21:42:52 2022 +0100 @@ -1,5 +1,15 @@ module LazyTensors +export LazyTensorMappingApplication +export LazyTensorMappingTranspose +export TensorMappingComposition +export LazyLinearMap +export IdentityMapping +export InflatedTensorMapping +export LazyOuterProduct +export ⊗ +export SizeMismatch + include("tensor_mapping.jl") include("lazy_array.jl") include("lazy_tensor_operations.jl")
--- a/src/LazyTensors/lazy_tensor_operations.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Tue Mar 15 21:42:52 2022 +0100 @@ -1,3 +1,5 @@ +using Sbplib.RegionIndices + """ LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R} @@ -7,12 +9,17 @@ With a mapping `m` and a vector `v` the LazyTensorMappingApplication object can be created by `m*v`. The actual result will be calcualted when indexing into `m*v`. """ -struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{T,R,D}, AA<:AbstractArray{T,D}} <: LazyArray{T,R} +struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{<:Any,R,D}, AA<:AbstractArray{<:Any,D}} <: LazyArray{T,R} t::TM o::AA + + function LazyTensorMappingApplication(t::TensorMapping{<:Any,R,D}, o::AbstractArray{<:Any,D}) where {R,D} + I = ntuple(i->1, range_dim(t)) + T = typeof(apply(t,o,I...)) + return new{T,R,D,typeof(t), typeof(o)}(t,o) + end end # TODO: Do boundschecking on creation! -export LazyTensorMappingApplication Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Any,R}) where {T,R} = apply(ta.t, ta.o, I...) Base.getindex(ta::LazyTensorMappingApplication{T,1}, I::CartesianIndex{1}) where {T} = apply(ta.t, ta.o, I.I...) # Would otherwise be caught in the previous method. @@ -41,15 +48,14 @@ struct LazyTensorMappingTranspose{T,R,D, TM<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R} tm::TM end -export LazyTensorMappingTranspose # # TBD: Should this be implemented on a type by type basis or through a trait to provide earlier errors? # Jonatan 2020-09-25: Is the problem that you can take the transpose of any TensorMapping even if it doesn't implement `apply_transpose`? Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm) Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm -apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) -apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...) +apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) +apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...) range_size(tmt::LazyTensorMappingTranspose) = domain_size(tmt.tm) domain_size(tmt::LazyTensorMappingTranspose) = range_size(tmt.tm) @@ -65,8 +71,8 @@ end # TODO: Boundschecking in constructor. -apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...) -apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...) +apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...) +apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...) range_size(tmBinOp::LazyTensorMappingBinaryOperation) = range_size(tmBinOp.tm1) domain_size(tmBinOp::LazyTensorMappingBinaryOperation) = domain_size(tmBinOp.tm1) @@ -88,16 +94,15 @@ return new{T,R,K,D, typeof(t1), typeof(t2)}(t1,t2) end end -export TensorMappingComposition range_size(tm::TensorMappingComposition) = range_size(tm.t1) domain_size(tm::TensorMappingComposition) = domain_size(tm.t2) -function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,K,D} +function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,K,D} apply(c.t1, c.t2*v, I...) end -function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,K,D} +function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,K,D} apply_transpose(c.t2, c.t1'*v, I...) end @@ -125,12 +130,11 @@ return new{T,R,D,RD,AA}(A,range_indicies,domain_indicies) end end -export LazyLinearMap range_size(llm::LazyLinearMap) = size(llm.A)[[llm.range_indicies...]] domain_size(llm::LazyLinearMap) = size(llm.A)[[llm.domain_indicies...]] -function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} +function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} view_index = ntuple(i->:,ndims(llm.A)) for i ∈ 1:R view_index = Base.setindex(view_index, Int(I[i]), llm.range_indicies[i]) @@ -139,7 +143,7 @@ return sum(A_view.*v) end -function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} +function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} apply(LazyLinearMap(llm.A, llm.domain_indicies, llm.range_indicies), v, I...) end @@ -153,7 +157,6 @@ struct IdentityMapping{T,D} <: TensorMapping{T,D,D} size::NTuple{D,Int} end -export IdentityMapping IdentityMapping{T}(size::NTuple{D,Int}) where {T,D} = IdentityMapping{T,D}(size) IdentityMapping{T}(size::Vararg{Int,D}) where {T,D} = IdentityMapping{T,D}(size) @@ -162,8 +165,8 @@ range_size(tmi::IdentityMapping) = tmi.size domain_size(tmi::IdentityMapping) = tmi.size -apply(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...] -apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...] +apply(tmi::IdentityMapping{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = v[I...] +apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = v[I...] """ Base.:∘(tm, tmi) @@ -210,7 +213,6 @@ return new{T,R,D,D_before,R_middle,D_middle,D_after, typeof(tm)}(before, tm, after) end end -export InflatedTensorMapping """ InflatedTensorMapping(before, tm, after) InflatedTensorMapping(before,tm) @@ -256,7 +258,7 @@ ) end -function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} +function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} dim_before = range_dim(itm.before) dim_domain = domain_dim(itm.tm) dim_range = range_dim(itm.tm) @@ -268,7 +270,7 @@ return apply(itm.tm, v_inner, inner_index...) end -function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} +function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} dim_before = range_dim(itm.before) dim_domain = domain_dim(itm.tm) dim_range = range_dim(itm.tm) @@ -396,7 +398,6 @@ ``` """ function LazyOuterProduct end -export LazyOuterProduct function LazyOuterProduct(tm1::TensorMapping{T}, tm2::TensorMapping{T}) where T itm1 = InflatedTensorMapping(tm1, IdentityMapping{T}(range_size(tm2))) @@ -412,7 +413,6 @@ LazyOuterProduct(tms::Vararg{TensorMapping}) = foldl(LazyOuterProduct, tms) ⊗(a::TensorMapping, b::TensorMapping) = LazyOuterProduct(a,b) -export ⊗ function check_domain_size(tm::TensorMapping, sz) @@ -425,9 +425,33 @@ tm::TensorMapping sz end -export SizeMismatch function Base.showerror(io::IO, err::SizeMismatch) print(io, "SizeMismatch: ") print(io, "domain size $(domain_size(err.tm)) of TensorMapping not matching size $(err.sz)") end + +function apply_with_region(op, v, boundary_width::Integer, dim_size::Integer, i) + if 0 < i <= boundary_width + return LazyTensors.apply(op,v,Index(i,Lower)) + elseif boundary_width < i <= dim_size-boundary_width + return LazyTensors.apply(op,v,Index(i,Interior)) + elseif dim_size-boundary_width < i <= dim_size + return LazyTensors.apply(op,v,Index(i,Upper)) + else + error("Bounds error") # TODO: Make this more standard + end +end +# TBD: Can these methods be merge by having a function as an arguement instead? +# TODO: Add inference test that show where things break and how this rewrite fixes it. +function apply_transpose_with_region(op, v, boundary_width::Integer, dim_size::Integer, i) + if 0 < i <= boundary_width + return LazyTensors.apply_transpose(op,v,Index(i,Lower)) + elseif boundary_width < i <= dim_size-boundary_width + return LazyTensors.apply_transpose(op,v,Index(i,Interior)) + elseif dim_size-boundary_width < i <= dim_size + return LazyTensors.apply_transpose(op,v,Index(i,Upper)) + else + error("Bounds error") # TODO: Make this more standard + end +end
--- a/src/LazyTensors/tensor_mapping.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/LazyTensors/tensor_mapping.jl Tue Mar 15 21:42:52 2022 +0100 @@ -1,10 +1,16 @@ +export TensorMapping +export apply +export apply_transpose +export range_dim, domain_dim +export range_size, domain_size + """ TensorMapping{T,R,D} Describes a mapping of a `D` dimension tensor to an `R` dimension tensor. The action of the mapping is implemented through the method ```julia - apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T} + apply(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg) where {R,D,T} ``` The size of the range and domain that the operator works with should be returned by @@ -21,23 +27,20 @@ ``` """ abstract type TensorMapping{T,R,D} end -export TensorMapping """ - apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T} + apply(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg) where {R,D,T} Return the result of the mapping for a given index. """ function apply end -export apply """ - apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg) where {R,D,T} + apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg) where {R,D,T} Return the result of the transposed mapping for a given index. """ function apply_transpose end -export apply_transpose """ range_dim(::TensorMapping) @@ -51,7 +54,6 @@ """ domain_dim(::TensorMapping{T,R,D}) where {T,R,D} = D -export range_dim, domain_dim """ range_size(M::TensorMapping) @@ -67,7 +69,6 @@ """ function domain_size end -export range_size, domain_size """ eltype(::TensorMapping{T})
--- a/src/RegionIndices/RegionIndices.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/RegionIndices/RegionIndices.jl Tue Mar 15 21:42:52 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 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/SbpOperators.jl Tue Mar 15 21:42:52 2022 +0100 @@ -19,12 +19,15 @@ even = 1 end +export closure_size + include("stencil.jl") include("readoperator.jl") include("volumeops/volume_operator.jl") include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/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 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Tue Mar 15 21:42:52 2022 +0100 @@ -62,28 +62,27 @@ LazyTensors.range_size(op::BoundaryOperator) = () LazyTensors.domain_size(op::BoundaryOperator) = (op.size,) -function LazyTensors.apply(op::BoundaryOperator{T,Lower}, v::AbstractVector{T}) where T +function LazyTensors.apply(op::BoundaryOperator{<:Any,Lower}, v::AbstractVector) apply_stencil(op.stencil,v,1) end -function LazyTensors.apply(op::BoundaryOperator{T,Upper}, v::AbstractVector{T}) where T +function LazyTensors.apply(op::BoundaryOperator{<:Any,Upper}, v::AbstractVector) apply_stencil_backwards(op.stencil,v,op.size) end -function LazyTensors.apply_transpose(op::BoundaryOperator{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T +function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Lower}, v::AbstractArray{<:Any,0}, i::Index{Lower}) return op.stencil[Int(i)-1]*v[] end -function LazyTensors.apply_transpose(op::BoundaryOperator{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T +function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Upper}, v::AbstractArray{<:Any,0}, i::Index{Upper}) return op.stencil[op.size[1] - Int(i)]*v[] end # Catch all combinations of Lower, Upper and Interior not caught by the two previous methods. -function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i::Index) where T - return zero(T) +function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i::Index) + return zero(eltype(v)) end -function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T - r = getregion(i, closure_size(op), op.size) - apply_transpose(op, v, Index(i,r)) +function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i) + 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 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/operators/standard_diagonal.toml Tue Mar 15 21:42:52 2022 +0100 @@ -20,6 +20,10 @@ H.inner = "1" H.closure = ["1/2"] +e.closure = ["1"] +d1.closure = {s = ["3/2", "-2", "1/2"], c = 1} + + D1.inner_stencil = ["-1/2", "0", "1/2"] D1.closure_stencils = [ {s = ["-1", "1"], c = 1}, @@ -30,8 +34,10 @@ {s = ["1", "-2", "1"], c = 1}, ] -e.closure = ["1"] -d1.closure = {s = ["3/2", "-2", "1/2"], c = 1} +D2variable.inner_stencil = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]] +D2variable.closure_stencils = [ + {s = [["2", "-1", "0"],["-3", "1", "0"],["1","0","0"]], c = 1}, +] [[stencil_set]] @@ -40,6 +46,9 @@ H.inner = "1" H.closure = ["17/48", "59/48", "43/48", "49/48"] +e.closure = ["1"] +d1.closure = {s = ["11/6", "-3", "3/2", "-1/3"], c = 1} + D1.inner_stencil = ["1/12","-2/3","0","2/3","-1/12"] D1.closure_stencils = [ {s = [ "-24/17", "59/34", "-4/17", "-3/34", "0", "0"], c = 1}, @@ -56,5 +65,79 @@ {s = [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"], c = 4}, ] +D2variable.inner_stencil = [ + ["-1/8", "1/6", "-1/8", "0", "0" ], + [ "1/6", "1/2", "1/2", "1/6", "0" ], + ["-1/24", "-5/6", "-3/4", "-5/6", "-1/24"], + [ "0", "1/6", "1/2", "1/2", "1/6" ], + [ "0", "0", "-1/8", "1/6", "-1/8" ], +] +D2variable.closure_stencils = [ + {c = 1, s = [ + [ "920/289", "-59/68", "-81031200387/366633756146", "-69462376031/733267512292", "0", "0", "0", "0" ], + ["-1740/289", "0", "6025413881/7482321554", "1612249989/7482321554", "0", "0", "0", "0" ], + [ "1128/289", "59/68", "-6251815797/8526366422", "-639954015/17052732844", "0", "0", "0", "0" ], + [ "-308/289", "0", "1244724001/7482321554", "-752806667/7482321554", "0", "0", "0", "0" ], + [ "0", "0", "-148737261/10783345769", "148737261/10783345769", "0", "0", "0", "0" ], + [ "0", "0", "-3/833", "3/833", "0", "0", "0", "0" ], + ]}, + {c = 2, s = [ + [ "12/17", "0", "102125659/440136562", "27326271/440136562", "0", "0", "0", "0" ], + [ "-59/68", "0", "-156920047993625/159775733917868", "-12001237118451/79887866958934", "0", "0", "0", "0" ], + [ "2/17", "0", "1489556735319/1857857371138", "149729180391/1857857371138", "0", "0", "0", "0" ], + [ "3/68", "0", "-13235456910147/159775733917868", "3093263736297/79887866958934", "0", "0", "0", "0" ], + [ "0", "0", "67535018271/2349643145851", "-67535018271/2349643145851", "0", "0", "0", "0" ], + [ "0", "0", "441/181507", "-441/181507", "0", "0", "0", "0" ], + ]}, + {c = 3, s = [ + [ "-96/731", "59/172", "-6251815797/21566691538", "-639954015/43133383076", "0", "0", "0", "0" ], + [ "118/731", "0", "87883847383821/79887866958934", "8834021643069/79887866958934", "0", "0", "0", "0" ], + [ "-16/731", "-59/172", "-1134866646907639536627/727679167377258785038", "-13777050223300597/23487032885926596", "-26254/557679", "0", "0", "0" ], + [ "-6/731", "0", "14509020271326561681/14850595252597118062", "17220493277981/79887866958934", "1500708/7993399", "0", "0", "0" ], + [ "0", "0", "-4841930283098652915/21402328452272317207", "31597236232005/115132514146699", "-26254/185893", "0", "0", "0" ], + [ "0", "0", "-2318724711/1653303156799", "960119/1147305747", "13564/23980197", "0", "0", "0" ], + ]}, + {c = 4, s = [ + [ "-36/833", "0", "1244724001/21566691538", "-752806667/21566691538", "0", "0", "0", "0" ], + [ "177/3332", "0", "-780891957698673/7829010961975532", "3724542049827/79887866958934", "0", "0", "0", "0" ], + [ "-6/833", "0", "14509020271326561681/16922771334354855466", "2460070468283/13005001597966", "1500708/9108757", "0", "0", "0" ], + [ "-9/3332", "0", "-217407431400324796377/207908333536359652868", "-1950062198436997/3914505480987766", "-7476412/9108757", "-2/49", "0", "0" ], + [ "0", "0", "4959271814984644613/21402328452272317207", "47996144728947/115132514146699", "4502124/9108757", "8/49", "0", "0" ], + [ "0", "0", "-2258420001/1653303156799", "-1063649/8893843", "1473580/9108757", "-6/49", "0", "0" ], + ]}, + {c = 5, s = [ + [ "0", "0", "-49579087/10149031312", "49579087/10149031312", "0", "0", "0", "0" ], + [ "0", "0", "1328188692663/37594290333616", "-1328188692663/37594290333616", "0", "0", "0", "0" ], + [ "0", "0", "-1613976761032884305/7963657098519931984", "10532412077335/42840005263888", "-564461/4461432", "0", "0", "0" ], + [ "0", "0", "4959271814984644613/20965546238960637264", "15998714909649/37594290333616", "375177/743572", "1/6", "0", "0" ], + [ "0", "0", "-8386761355510099813/128413970713633903242", "-2224717261773437/2763180339520776", "-280535/371786", "-5/6", "-1/24", "0" ], + [ "0", "0", "13091810925/13226425254392", "35039615/213452232", "1118749/2230716", "1/2", "1/6", "0" ], + [ "0", "0", "0", "0", "-1/8", "1/6", "-1/8", "0" ], + ]}, + {c = 6, s = [ + [ "0", "0", "-1/784", "1/784", "0", "0", "0", "0" ], + [ "0", "0", "8673/2904112", "-8673/2904112", "0", "0", "0", "0" ], + [ "0", "0", "-33235054191/26452850508784", "960119/1280713392", "3391/6692148", "0", "0", "0" ], + [ "0", "0", "-752806667/539854092016", "-1063649/8712336", "368395/2230716", "-1/8", "0", "0" ], + [ "0", "0", "13091810925/13226425254392", "35039615/213452232", "1118749/2230716", "1/2", "1/6", "0" ], + [ "0", "0", "-660204843/13226425254392", "-3290636/80044587", "-5580181/6692148", "-3/4", "-5/6", "-1/24"], + [ "0", "0", "0", "0", "1/6", "1/2", "1/2", "1/6" ], + [ "0", "0", "0", "0", "0", "-1/8", "1/6", "-1/8" ], + ]} +] + + + +[[stencil_set]] + +order = 6 + +H.inner = "1" +H.closure = ["13649/43200", "12013/8640", "2711/4320", "5359/4320", "7877/8640", "43801/43200"] + + + + + e.closure = ["1"] -d1.closure = {s = ["11/6", "-3", "3/2", "-1/3"], c = 1} +d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"]
--- a/src/SbpOperators/readoperator.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/readoperator.jl Tue Mar 15 21:42:52 2022 +0100 @@ -4,6 +4,7 @@ export get_stencil_set export parse_stencil +export parse_nested_stencil export parse_scalar export parse_tuple @@ -106,6 +107,33 @@ end end + +""" + parse_nested_stencil(parsed_toml) + +Accept parsed TOML and read it as a nested tuple. + +See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref). +""" +function parse_nested_stencil(parsed_toml) + if parsed_toml isa Array + weights = parse_stencil.(parsed_toml) + return CenteredNestedStencil(weights...) + end + + center = parsed_toml["c"] + weights = parse_tuple.(parsed_toml["s"]) + return NestedStencil(weights...; center) +end + +""" + parse_nested_stencil(T, parsed_toml) + +Parse the input as a nested stencil with element type `T`. +""" +parse_nested_stencil(T, parsed_toml) = NestedStencil{T}(parse_nested_stencil(parsed_toml)) + + """ parse_scalar(parsed_toml)
--- a/src/SbpOperators/stencil.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/stencil.jl Tue Mar 15 21:42:52 2022 +0100 @@ -1,11 +1,12 @@ export CenteredStencil +export CenteredNestedStencil struct Stencil{T,N} - range::Tuple{Int,Int} + range::UnitRange{Int64} weights::NTuple{N,T} - function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N} - @assert range[2]-range[1]+1 == N + function Stencil(range::UnitRange,weights::NTuple{N,T}) where {T, N} + @assert length(range) == N new{T,N}(range,weights) end end @@ -15,27 +16,30 @@ Create a stencil with the given weights with element `center` as the center of the stencil. """ -function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error +function Stencil(weights...; center::Int) + weights = promote(weights...) N = length(weights) - range = (1, N) .- center + range = (1:N) .- center return Stencil(range, weights) end -function Stencil{T}(s::Stencil) where T - return Stencil(s.range, T.(s.weights)) -end +Stencil{T,N}(s::Stencil{S,N}) where {T,S,N} = Stencil(s.range, T.(s.weights)) +Stencil{T}(s::Stencil) where T = Stencil{T,length(s)}(s) -Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil) +Base.convert(::Type{Stencil{T1,N}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s) +Base.convert(::Type{Stencil{T1}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s) -function CenteredStencil(weights::Vararg) +Base.promote_rule(::Type{Stencil{T1,N}}, ::Type{Stencil{T2,N}}) where {T1,T2,N} = Stencil{promote_type(T1,T2),N} + +function CenteredStencil(weights...) if iseven(length(weights)) throw(ArgumentError("a centered stencil must have an odd number of weights.")) end r = length(weights) ÷ 2 - return Stencil((-r, r), weights) + return Stencil(-r:r, weights) end @@ -48,7 +52,8 @@ return Stencil(s.range, a.*s.weights) end -Base.eltype(::Stencil{T}) where T = T +Base.eltype(::Stencil{T,N}) where {T,N} = T +Base.length(::Stencil{T,N}) where {T,N} = N function flip(s::Stencil) range = (-s.range[2], -s.range[1]) @@ -57,24 +62,103 @@ # Provides index into the Stencil based on offset for the root element @inline function Base.getindex(s::Stencil, i::Int) - @boundscheck if i < s.range[1] || s.range[2] < i + @boundscheck if i ∉ s.range return zero(eltype(s)) end return s.weights[1 + i - s.range[1]] end -Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N} - w = s.weights[1]*v[i + s.range[1]] - @simd for k ∈ 2:N - w += s.weights[k]*v[i + s.range[1] + k-1] +Base.@propagate_inbounds @inline function apply_stencil(s::Stencil, v::AbstractVector, i::Int) + w = zero(promote_type(eltype(s),eltype(v))) + @simd for k ∈ 1:length(s) + w += s.weights[k]*v[i + s.range[k]] + end + + return w +end + +Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil, v::AbstractVector, i::Int) + w = zero(promote_type(eltype(s),eltype(v))) + @simd for k ∈ length(s):-1:1 + w += s.weights[k]*v[i - s.range[k]] end return w end -Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N} - w = s.weights[N]*v[i - s.range[2]] - @simd for k ∈ N-1:-1:1 - w += s.weights[k]*v[i - s.range[1] - k + 1] - end - return w + +struct NestedStencil{T,N,M} + s::Stencil{Stencil{T,N},M} +end + +# Stencil input +NestedStencil(s::Vararg{Stencil}; center) = NestedStencil(Stencil(s... ; center)) +CenteredNestedStencil(s::Vararg{Stencil}) = NestedStencil(CenteredStencil(s...)) + +# Tuple input +function NestedStencil(weights::Vararg{NTuple{N,Any}}; center) where N + inner_stencils = map(w -> Stencil(w...; center), weights) + return NestedStencil(Stencil(inner_stencils... ; center)) +end +function CenteredNestedStencil(weights::Vararg{NTuple{N,Any}}) where N + inner_stencils = map(w->CenteredStencil(w...), weights) + return CenteredNestedStencil(inner_stencils...) +end + + +# Conversion +function NestedStencil{T,N,M}(ns::NestedStencil{S,N,M}) where {T,S,N,M} + return NestedStencil(Stencil{Stencil{T}}(ns.s)) +end + +function NestedStencil{T}(ns::NestedStencil{S,N,M}) where {T,S,N,M} + NestedStencil{T,N,M}(ns) +end + +function Base.convert(::Type{NestedStencil{T,N,M}}, s::NestedStencil{S,N,M}) where {T,S,N,M} + return NestedStencil{T,N,M}(s) +end +Base.convert(::Type{NestedStencil{T}}, stencil) where T = NestedStencil{T}(stencil) + +function Base.promote_rule(::Type{NestedStencil{T,N,M}}, ::Type{NestedStencil{S,N,M}}) where {T,S,N,M} + return NestedStencil{promote_type(T,S),N,M} end + +Base.eltype(::NestedStencil{T}) where T = T + +function scale(ns::NestedStencil, a) + range = ns.s.range + weights = ns.s.weights + + return NestedStencil(Stencil(range, scale.(weights,a))) +end + +function flip(ns::NestedStencil) + s_flip = flip(ns.s) + return NestedStencil(Stencil(s_flip.range, flip.(s_flip.weights))) +end + +Base.getindex(ns::NestedStencil, i::Int) = ns.s[i] + +"Apply inner stencils to `c` and get a concrete stencil" +Base.@propagate_inbounds function apply_inner_stencils(ns::NestedStencil, c::AbstractVector, i::Int) + weights = apply_stencil.(ns.s.weights, Ref(c), i) + return Stencil(ns.s.range, weights) +end + +"Apply the whole nested stencil" +Base.@propagate_inbounds function apply_stencil(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int) + s = apply_inner_stencils(ns,c,i) + return apply_stencil(s, v, i) +end + +"Apply inner stencils backwards to `c` and get a concrete stencil" +Base.@propagate_inbounds @inline function apply_inner_stencils_backwards(ns::NestedStencil, c::AbstractVector, i::Int) + weights = apply_stencil_backwards.(ns.s.weights, Ref(c), i) + return Stencil(ns.s.range, weights) +end + +"Apply the whole nested stencil backwards" +Base.@propagate_inbounds @inline function apply_stencil_backwards(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int) + s = apply_inner_stencils_backwards(ns,c,i) + return apply_stencil_backwards(s, v, i) +end
--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Tue Mar 15 21:42:52 2022 +0100 @@ -28,21 +28,20 @@ LazyTensors.domain_size(op::ConstantInteriorScalingOperator) = (op.size,) # TBD: @inbounds in apply methods? -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Lower}) return op.closure_weights[Int(i)]*v[Int(i)] end -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Interior}) return op.interior_weight*v[Int(i)] end -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Upper}) return op.closure_weights[op.size[1]-Int(i)+1]*v[Int(i)] end -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T - r = getregion(i, closure_size(op), op.size[1]) - return LazyTensors.apply(op, v, Index(i, r)) +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i) + 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 Tue Mar 15 21:42:52 2022 +0100 @@ -0,0 +1,195 @@ +export SecondDerivativeVariable + +# REVIEW: Fixa docs +""" + SecondDerivativeVariable{Dir,T,D,...} <: TensorMapping{T,D,D} + +A second derivative operator in direction `Dir` with a variable coefficient. +""" +struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: TensorMapping{T,D,D} + inner_stencil::IStencil + closure_stencils::NTuple{M,CStencil} + size::NTuple{D,Int} + coefficient::TArray + + function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, size::NTuple{D,Int}, coefficient::AbstractArray) where {Dir,T,D,M} + IStencil = typeof(inner_stencil) + CStencil = eltype(closure_stencils) + TArray = typeof(coefficient) + return new{Dir,T,D,M,IStencil,CStencil,TArray}(inner_stencil,closure_stencils,size, coefficient) + end +end + +function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir) + check_coefficient(grid, coeff) + + Δxᵢ = spacing(grid)[dir] + scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2) + scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2) + return SecondDerivativeVariable{dir, dimension(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff) +end + +function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils) + return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1) +end + +@doc raw""" + SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir) + +Create a `TensorMapping` for the second derivative with a variable coefficient +`coeff` on `grid` from the stencils in `stencil_set`. The direction is +determined by `dir`. + +`coeff` is a grid function on `grid`. + +# Example +With +``` +D = SecondDerivativeVariable(g, c, stencil_set, 2) +``` +then `D*u` approximates +```math +\frac{\partial}{\partial y} c(x,y) \frac{\partial u}{\partial y}, +``` +on ``(0,1)⨯(0,1)`` represented by `g`. +""" +function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir::Int) + inner_stencil = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"]) + closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"]) + + return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, dir) +end + +function check_coefficient(grid, coeff) + if dimension(grid) != ndims(coeff) + throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(dimension(grid))")) + end + + if size(grid) != size(coeff) + throw(DimensionMismatch("the size $(size(coeff)) of the coefficient does not match the size $(size(grid)) of the grid")) + end +end + +derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir + +closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils) + +LazyTensors.range_size(op::SecondDerivativeVariable) = op.size +LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size + + +function derivative_view(op, a, I) + d = derivative_direction(op) + + Iview = Base.setindex(I,:,d) + return @view a[Iview...] +end + +function apply_lower(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i) +end + +function apply_interior(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + return apply_stencil(op.inner_stencil, c̃, ṽ, i) +end + +function apply_upper(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I::Vararg{Index}) + if I[derivative_direction(op)] isa Index{Lower} + return apply_lower(op, v, Int.(I)...) + elseif I[derivative_direction(op)] isa Index{Upper} + return apply_upper(op, v, Int.(I)...) + elseif I[derivative_direction(op)] isa Index{Interior} + return apply_interior(op, v, Int.(I)...) + else + error("Invalid region") + end +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...) + dir = derivative_direction(op) + + i = I[dir] + + I = map(i->Index(i, Interior), I) + if 0 < i <= closure_size(op) + I = Base.setindex(I, Index(i, Lower), dir) + return LazyTensors.apply(op, v, I...) + elseif closure_size(op) < i <= op.size[dir]-closure_size(op) + I = Base.setindex(I, Index(i, Interior), dir) + return LazyTensors.apply(op, v, I...) + elseif op.size[dir]-closure_size(op) < i <= op.size[dir] + I = Base.setindex(I, Index(i, Upper), dir) + return LazyTensors.apply(op, v, I...) + else + error("Bounds error") # TODO: Make this more standard + end +end + + +## 2D Specific implementations to avoid instability +## TODO: Should really be solved by fixing the general methods instead + + +## x-direction +function apply_lower(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i) +end + +function apply_interior(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, i) +end + +function apply_upper(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) +end + + +## y-direction +function apply_lower(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + return @inbounds apply_stencil(op.closure_stencils[j], c̃, ṽ, j) +end + +function apply_interior(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, j) +end + +function apply_upper(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + stencil = op.closure_stencils[op.size[derivative_direction(op)]-j+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, j) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/inference_trouble.txt Tue Mar 15 21:42:52 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 15 21:30:56 2022 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Tue Mar 15 21:42:52 2022 +0100 @@ -43,19 +43,18 @@ LazyTensors.range_size(op::VolumeOperator) = op.size LazyTensors.domain_size(op::VolumeOperator) = op.size -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Lower}) return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i)) end -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Interior}) return apply_stencil(op.inner_stencil, v, Int(i)) end -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Upper}) return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i)) end -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T - r = getregion(i, closure_size(op), op.size[1]) - return LazyTensors.apply(op, v, Index(i, r)) +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i) + return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i) end
--- a/test/LazyTensors/lazy_tensor_operations_test.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Tue Mar 15 21:42:52 2022 +0100 @@ -36,10 +36,8 @@ m = SizeDoublingMapping{Int, 1, 1}((3,)) v = [0,1,2] - @test m*v isa AbstractVector{Int} @test size(m*v) == 2 .*size(v) @test (m*v)[0] == (:apply,v,(0,)) - @test m*m*v isa AbstractVector{Int} @test (m*m*v)[1] == (:apply,m*v,(1,)) @test (m*m*v)[3] == (:apply,m*v,(3,)) @test (m*m*v)[6] == (:apply,m*v,(6,)) @@ -80,6 +78,43 @@ v = [[1 2];[3 4]] @test m*v == [[2 4];[6 8]] @test (m*v)[2,1] == 6 + + @testset "Type calculation" begin + m = ScalingOperator{Int,1}(2,(3,)) + v = [1.,2.,3.] + @test m*v isa AbstractVector{Float64} + @test m*v == [2.,4.,6.] + @inferred m*v + @inferred (m*v)[1] + + m = ScalingOperator{Int,2}(2,(2,2)) + v = [[1. 2.];[3. 4.]] + @test m*v == [[2. 4.];[6. 8.]] + @test (m*v)[2,1] == 6. + @inferred m*v + @inferred (m*v)[1] + + m = ScalingOperator{ComplexF64,1}(2. +2. *im,(3,)) + v = [1.,2.,3.] + @test m*v isa AbstractVector{ComplexF64} + @test m*v == [2. + 2. *im, 4. + 4. *im, 6. + 6. *im] + @inferred m*v + @inferred (m*v)[1] + + m = ScalingOperator{ComplexF64,1}(1,(3,)) + v = [2. + 2. *im, 4. + 4. *im, 6. + 6. *im] + @test m*v isa AbstractVector{ComplexF64} + @test m*v == [2. + 2. *im, 4. + 4. *im, 6. + 6. *im] + @inferred m*v + @inferred (m*v)[1] + + m = ScalingOperator{Float64,1}(2., (3,)) + v = [[1,2,3], [3,2,1],[1,3,1]] + @test m*v isa AbstractVector{Vector{Float64}} + @test m*v == [[2.,4.,6.], [6.,4.,2.],[2.,6.,2.]] + @inferred m*v + @inferred (m*v)[1] + end end @testset "TensorMapping binary operations" begin @@ -107,6 +142,8 @@ @test range_size(A+B) == range_size(A) == range_size(B) @test domain_size(A+B) == domain_size(A) == domain_size(B) + + @test ((A+B)*ComplexF64[1.1,1.2,1.3])[3] isa ComplexF64 end @@ -129,6 +166,9 @@ v = rand(2) @test (Ã∘B̃)'*v ≈ B'*A'*v rtol=1e-14 + + @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64 + @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64 end @testset "LazyLinearMap" begin @@ -194,6 +234,10 @@ @test I*v == v @test I'*v == v + v = rand(ComplexF64,sz...) + @test I*v == v + @test I'*v == v + @test range_size(I) == sz @test domain_size(I) == sz end @@ -322,6 +366,16 @@ end end + @testset "application to other type" begin + tm = InflatedTensorMapping(I(3,2), A, I(4)) + + v = rand(ComplexF64, domain_size(tm)...) + @test (tm*v)[1,2,3,1] isa ComplexF64 + + v = rand(ComplexF64, domain_size(tm')...) + @test (tm'*v)[1,2,2,1] isa ComplexF64 + end + @testset "Inference of application" begin struct ScalingOperator{T,D} <: TensorMapping{T,D,D} λ::T
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Tue Mar 15 21:42:52 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)) @@ -66,6 +66,14 @@ @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] + + v = evalOn(g_1D, x->1. +x*im) + @test (op_l*v)[] isa ComplexF64 + + u = fill(1. +im) + @test (op_l'*u)[1] isa ComplexF64 + @test (op_l'*u)[5] isa ComplexF64 + @test (op_l'*u)[11] isa ComplexF64 end @testset "2D" begin
--- a/test/SbpOperators/readoperator_test.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/test/SbpOperators/readoperator_test.jl Tue Mar 15 21:42:52 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 15 21:30:56 2022 +0100 +++ b/test/SbpOperators/stencil_test.jl Tue Mar 15 21:42:52 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
--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Tue Mar 15 21:42:52 2022 +0100 @@ -24,6 +24,9 @@ @test a*v == [.1,.2,.5,.5,.5,.2,.1] @test a'*v == [.1,.2,.5,.5,.5,.2,.1] + @test (a*rand(ComplexF64, domain_size(a)... ))[1] isa ComplexF64 + @test (a'*rand(ComplexF64, domain_size(a')...))[1] isa ComplexF64 + @test range_size(a) == (7,) @test domain_size(a) == (7,)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl Tue Mar 15 21:42:52 2022 +0100 @@ -0,0 +1,187 @@ +using Test + +using Sbplib.Grids +using Sbplib.LazyTensors +using Sbplib.SbpOperators +using Sbplib.RegionIndices +using Sbplib.SbpOperators: NestedStencil, CenteredNestedStencil + +using LinearAlgebra + +@testset "SecondDerivativeVariable" begin + interior_stencil = CenteredNestedStencil((1/2, 1/2, 0.),(-1/2, -1., -1/2),( 0., 1/2, 1/2)) + closure_stencils = [ + NestedStencil(( 2., -1., 0.),(-3., 1., 0.), (1., 0., 0.), center = 1), + ] + + @testset "1D" begin + g = EquidistantGrid(11, 0., 1.) + c = [ 1., 3., 6., 10., 15., 21., 28., 36., 45., 55., 66.] + @testset "Constructors" begin + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) isa TensorMapping + + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) + @test range_dim(D₂ᶜ) == 1 + @test domain_dim(D₂ᶜ) == 1 + + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test SecondDerivativeVariable(g, c, stencil_set, 1) isa SecondDerivativeVariable + + @testset "checking c" begin + c_short = rand(5) + c_long = rand(16) + c_higher_dimension = rand(11,11) + + @test_throws DimensionMismatch("the size (5,) of the coefficient does not match the size (11,) of the grid") SecondDerivativeVariable(g, c_short, interior_stencil, closure_stencils) + @test_throws DimensionMismatch("the size (16,) of the coefficient does not match the size (11,) of the grid") SecondDerivativeVariable(g, c_long, interior_stencil, closure_stencils) + @test_throws ArgumentError("The coefficient has dimension 2 while the grid is dimension 1") SecondDerivativeVariable(g, c_higher_dimension, interior_stencil, closure_stencils,1) + end + end + + @testset "sizes" begin + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) + @test closure_size(D₂ᶜ) == 1 + @test range_size(D₂ᶜ) == (11,) + @test domain_size(D₂ᶜ) == (11,) + end + + @testset "application" begin + + function apply_to_functions(; v, c) + g = EquidistantGrid(11, 0., 10.) # h = 1 + c̄ = evalOn(g,c) + v̄ = evalOn(g,v) + + D₂ᶜ = SecondDerivativeVariable(g, c̄, interior_stencil, closure_stencils) + return D₂ᶜ*v̄ + end + + @test apply_to_functions(v=x->1., c=x-> -1.) == zeros(11) + @test apply_to_functions(v=x->1., c=x-> -x ) == zeros(11) + @test apply_to_functions(v=x->x, c=x-> 1.) == zeros(11) + @test apply_to_functions(v=x->x, c=x-> -x ) == -ones(11) + @test apply_to_functions(v=x->x^2, c=x-> 1.) == 2ones(11) + end + + @testset "type stability" begin + g = EquidistantGrid(11, 0., 10.) # h = 1 + c̄ = evalOn(g,x-> -1) + v̄ = evalOn(g,x->1.) + + D₂ᶜ = SecondDerivativeVariable(g, c̄, interior_stencil, closure_stencils) + + @inferred SbpOperators.apply_lower(D₂ᶜ, v̄, 1) + @inferred SbpOperators.apply_interior(D₂ᶜ, v̄, 5) + @inferred SbpOperators.apply_upper(D₂ᶜ, v̄, 11) + @inferred (D₂ᶜ*v̄)[Index(1,Lower)] + end + end + + @testset "2D" begin + g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1 + c = evalOn(g, (x,y)->x+y) + @testset "Constructors" begin + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) isa TensorMapping + @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) isa TensorMapping + + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) + @test range_dim(D₂ᶜ) == 2 + @test domain_dim(D₂ᶜ) == 2 + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test SecondDerivativeVariable(g, c, stencil_set, 1) isa SecondDerivativeVariable + end + + @testset "sizes" begin + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) + @test range_size(D₂ᶜ) == (11,9) + @test domain_size(D₂ᶜ) == (11,9) + @test closure_size(D₂ᶜ) == 1 + + D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) + @test range_size(D₂ᶜ) == (11,9) + @test domain_size(D₂ᶜ) == (11,9) + @test closure_size(D₂ᶜ) == 1 + end + + @testset "application" begin + function apply_to_functions(dir; v, c) + g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1 + c̄ = evalOn(g,c) + v̄ = evalOn(g,v) + + D₂ᶜ = SecondDerivativeVariable(g, c̄, interior_stencil, closure_stencils,dir) + return D₂ᶜ*v̄ + end + + # x-direction + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)-> -1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)->- x ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->x, c=(x,y)-> -x ) == -ones(11,9) + @test apply_to_functions(1,v=(x,y)->x^2, c=(x,y)-> 1.) == 2ones(11,9) + + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)->- y ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y, c=(x,y)-> -y ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y^2, c=(x,y)-> 1.) == zeros(11,9) + + # y-direction + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)-> -1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)->- y ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->y, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->y, c=(x,y)-> -y ) == -ones(11,9) + @test apply_to_functions(2,v=(x,y)->y^2, c=(x,y)-> 1.) == 2ones(11,9) + + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)->- x ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> -x ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x^2, c=(x,y)-> 1.) == zeros(11,9) + + + @testset "standard diagonal operators" begin + c(x,y) = exp(x) + exp(1.5(1-y)) + v(x,y) = sin(x) + cos(1.5(1-y)) + + Dxv(x,y) = cos(x)*exp(x) - (exp(x) + exp(1.5 - 1.5y))*sin(x) + Dyv(x,y) = -1.5(1.5exp(x) + 1.5exp(1.5 - 1.5y))*cos(1.5 - 1.5y) - 2.25exp(1.5 - 1.5y)*sin(1.5 - 1.5y) + + g₁ = EquidistantGrid((60,67), (0.,0.), (1.,2.)) + g₂ = refine(g₁,2) + + c̄₁ = evalOn(g₁, c) + c̄₂ = evalOn(g₂, c) + + v̄₁ = evalOn(g₁, v) + v̄₂ = evalOn(g₂, v) + + + function convergence_rate_estimate(stencil_set, dir, Dv_true) + D₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, dir) + D₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, dir) + + Dv̄₁ = D₁*v̄₁ + Dv̄₂ = D₂*v̄₂ + + Dv₁ = evalOn(g₁,Dv_true) + Dv₂ = evalOn(g₂,Dv_true) + + e₁ = norm(Dv̄₁ - Dv₁)/norm(Dv₁) + e₂ = norm(Dv̄₂ - Dv₂)/norm(Dv₂) + + return log2(e₁/e₂) + end + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 1.5 rtol = 1e-1 + @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 1.5 rtol = 1e-1 + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) + @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 2.5 rtol = 1e-1 + @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 2.5 rtol = 2e-1 + end + end + end +end +
--- a/test/SbpOperators/volumeops/volume_operator_test.jl Tue Mar 15 21:30:56 2022 +0100 +++ b/test/SbpOperators/volumeops/volume_operator_test.jl Tue Mar 15 21:42:52 2022 +0100 @@ -89,6 +89,8 @@ @testset "Application" begin @test op_x*v ≈ rx rtol = 1e-14 @test op_y*v ≈ ry rtol = 1e-14 + + @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64 end @testset "Regions" begin