Mercurial > repos > public > sbplib_julia
changeset 1101:1e8270c18edb feature/lazy_tensors/pretty_printing
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 12 May 2022 21:52:47 +0200 |
| parents | 67969bd7e642 (current diff) 74c54996de6a (diff) |
| children | 84820d4780fa |
| files | src/LazyTensors/lazy_tensor_operations.jl src/LazyTensors/tensor_types.jl test/LazyTensors/lazy_tensor_operations_test.jl test/LazyTensors/tensor_types_test.jl |
| diffstat | 31 files changed, 1001 insertions(+), 611 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Mon Mar 21 09:32:26 2022 +0100 +++ b/Manifest.toml Thu May 12 21:52:47 2022 +0200 @@ -1,13 +1,13 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.7.0" +julia_version = "1.7.1" manifest_format = "2.0" [[deps.Adapt]] deps = ["LinearAlgebra"] -git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7" +git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.3.1" +version = "3.3.3" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
--- a/Notes.md Mon Mar 21 09:32:26 2022 +0100 +++ b/Notes.md Thu May 12 21:52:47 2022 +0200 @@ -136,12 +136,9 @@ ## Reasearch and thinking - [ ] Use a trait to indicate that a LazyTensor har the same range and domain? - - [ ] Rename all the Tensor stuff to just LazyOperator, LazyApplication and so on? - [ ] Check how the native julia doc generator works - [ ] Check if Vidars design docs fit in there - [ ] Create a macro @lazy which replaces a binary op (+,-) by its lazy equivalent? Would be a neat way to indicate which evaluations are lazy without cluttering/confusing with special characters. - - [ ] Specificera operatorer i TOML eller något liknande? - H.. H_gamma etc.) - [ ] Dispatch on Lower() instead of the type Lower so `::Lower` instead of `::Type{Lower}` ??? Seems better unless there is some specific reason to use the type instead of the value. - [ ] How do we handle mixes of periodic and non-periodic grids? Seems it should be supported on the grid level and on the 1d operator level. Between there it should be transparent. @@ -149,12 +146,52 @@ - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input? - [ ] Figure out how to treat the borrowing parameters of operators. Include in into the struct? Expose via function dispatched on the operator type and grid? +## Identifiers for regions +The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probabily be included in the grid module. This allows new grid types to come with their own regions. + ## Regions and tensormappings - [ ] Use a trait to indicate if a LazyTensor uses indices with regions. The default should be that they do NOT. - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle? - [ ] Figure out repeated application of regioned LazyTensors. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work? +### Ideas for information sharing functions +```julia +using StaticArrays + +function regions(op::SecondDerivativeVariable) + t = ntuple(i->(Interior(),),range_dim(op)) + return Base.setindex(t, (Lower(), Interior(), Upper()), derivative_direction(op)) +end + +function regionsizes(op::SecondDerivativeVariable) + sz = tuple.(range_size(op)) + + cl = closuresize(op) + return Base.setindex(sz, (cl, n-2cl, cl), derivative_direction(op)) +end + + +g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1 +c = evalOn(g, (x,y)->x+y) + +D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) +@test regions(D₂ᶜ) == ( + (Lower(), Interior(), Upper()), + (Interior(),), +) +@test regionsizes(D₂ᶜ) == ((1,9,1),(9,)) + + +D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2) +@test regions(D₂ᶜ) == ( + (Interior(),), + (Lower(), Interior(), Upper()), +) +@test regionsizes(D₂ᶜ) == ((11,),(1,7,1)) +``` + + ## Boundschecking and dimension checking Does it make sense to have boundschecking only in getindex methods? This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway.
--- a/TODO.md Mon Mar 21 09:32:26 2022 +0100 +++ b/TODO.md Thu May 12 21:52:47 2022 +0200 @@ -1,19 +1,15 @@ # 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 - - [ ] Add boundschecking in LazyTensorApplication - - [ ] Start renaming things in LazyTensors - [ ] Clean up RegionIndices 1. [ ] Write tests for how things should work 2. [ ] Update RegionIndices accordingly @@ -22,13 +18,11 @@ - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`? - [ ] Add custom pretty printing to LazyTensors/SbpOperators to enhance readability of e.g error messages. See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing) + - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe - [ ] Move export statements to top of each module - - [ ] 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 Mon Mar 21 09:32:26 2022 +0100 +++ b/src/LazyTensors/LazyTensors.jl Thu May 12 21:52:47 2022 +0200 @@ -1,12 +1,13 @@ module LazyTensors -export LazyTensorApplication -export LazyTensorTranspose -export LazyTensorComposition -export LazyLinearMap +export TensorApplication +export TensorTranspose +export TensorComposition export IdentityTensor export ScalingTensor -export InflatedLazyTensor +export DiagonalTensor +export DenseTensor +export InflatedTensor export LazyOuterProduct export ⊗ export DomainSizeMismatch @@ -19,16 +20,16 @@ include("tuple_manipulation.jl") # Applying lazy tensors to vectors -Base.:*(a::LazyTensor, v::AbstractArray) = LazyTensorApplication(a,v) +Base.:*(a::LazyTensor, v::AbstractArray) = TensorApplication(a,v) Base.:*(a::LazyTensor, b::LazyTensor) = throw(MethodError(Base.:*,(a,b))) Base.:*(a::LazyTensor, args::Union{LazyTensor, AbstractArray}...) = foldr(*,(a,args...)) # Addition and subtraction of lazy tensors -Base.:+(s::LazyTensor, t::LazyTensor) = LazyTensorBinaryOperation{:+}(s,t) -Base.:-(s::LazyTensor, t::LazyTensor) = LazyTensorBinaryOperation{:-}(s,t) +Base.:+(s::LazyTensor, t::LazyTensor) = ElementwiseTensorOperation{:+}(s,t) +Base.:-(s::LazyTensor, t::LazyTensor) = ElementwiseTensorOperation{:-}(s,t) # Composing lazy tensors -Base.:∘(s::LazyTensor, t::LazyTensor) = LazyTensorComposition(s,t) +Base.:∘(s::LazyTensor, t::LazyTensor) = TensorComposition(s,t) # Outer products of tensors ⊗(a::LazyTensor, b::LazyTensor) = LazyOuterProduct(a,b)
--- a/src/LazyTensors/lazy_array.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/LazyTensors/lazy_array.jl Thu May 12 21:52:47 2022 +0200 @@ -36,7 +36,7 @@ function Base.getindex(lfa::LazyFunctionArray{F,T,D}, I::Vararg{Int,D}) where {F,T,D} @boundscheck checkbounds(lfa, I...) - return lfa.f(I...) + return @inbounds lfa.f(I...) end @@ -64,14 +64,14 @@ Base.size(v::LazyElementwiseOperation) = size(v.a) -evaluate(leo::LazyElementwiseOperation{T,D,:+}, I::Vararg{Int,D}) where {T,D} = @inbounds leo.a[I...] + leo.b[I...] -evaluate(leo::LazyElementwiseOperation{T,D,:-}, I::Vararg{Int,D}) where {T,D} = @inbounds leo.a[I...] - leo.b[I...] -evaluate(leo::LazyElementwiseOperation{T,D,:*}, I::Vararg{Int,D}) where {T,D} = @inbounds leo.a[I...] * leo.b[I...] -evaluate(leo::LazyElementwiseOperation{T,D,:/}, I::Vararg{Int,D}) where {T,D} = @inbounds leo.a[I...] / leo.b[I...] +Base.@propagate_inbounds evaluate(leo::LazyElementwiseOperation{T,D,:+}, I::Vararg{Int,D}) where {T,D} = leo.a[I...] + leo.b[I...] +Base.@propagate_inbounds evaluate(leo::LazyElementwiseOperation{T,D,:-}, I::Vararg{Int,D}) where {T,D} = leo.a[I...] - leo.b[I...] +Base.@propagate_inbounds evaluate(leo::LazyElementwiseOperation{T,D,:*}, I::Vararg{Int,D}) where {T,D} = leo.a[I...] * leo.b[I...] +Base.@propagate_inbounds evaluate(leo::LazyElementwiseOperation{T,D,:/}, I::Vararg{Int,D}) where {T,D} = leo.a[I...] / leo.b[I...] function Base.getindex(leo::LazyElementwiseOperation{T,D}, I::Vararg{Int,D}) where {T,D} @boundscheck checkbounds(leo.a, I...) - return evaluate(leo, I...) + return @inbounds evaluate(leo, I...) end # Define lazy operations for AbstractArrays. Operations constructs a LazyElementwiseOperation which
--- a/src/LazyTensors/lazy_tensor_operations.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Thu May 12 21:52:47 2022 +0200 @@ -1,17 +1,17 @@ """ - LazyTensorApplication{T,R,D} <: LazyArray{T,R} + TensorApplication{T,R,D} <: LazyArray{T,R} Struct for lazy application of a LazyTensor. Created using `*`. Allows the result of a `LazyTensor` applied to a vector to be treated as an `AbstractArray`. -With a mapping `m` and a vector `v` the LazyTensorApplication object can be created by `m*v`. +With a mapping `m` and a vector `v` the TensorApplication object can be created by `m*v`. The actual result will be calcualted when indexing into `m*v`. """ -struct LazyTensorApplication{T,R,D, TM<:LazyTensor{<:Any,R,D}, AA<:AbstractArray{<:Any,D}} <: LazyArray{T,R} +struct TensorApplication{T,R,D, TM<:LazyTensor{<:Any,R,D}, AA<:AbstractArray{<:Any,D}} <: LazyArray{T,R} t::TM o::AA - function LazyTensorApplication(t::LazyTensor{<:Any,R,D}, o::AbstractArray{<:Any,D}) where {R,D} + function TensorApplication(t::LazyTensor{<:Any,R,D}, o::AbstractArray{<:Any,D}) where {R,D} @boundscheck check_domain_size(t, size(o)) I = ntuple(i->1, range_dim(t)) T = typeof(apply(t,o,I...)) @@ -19,16 +19,16 @@ end end -function Base.getindex(ta::LazyTensorApplication{T,R}, I::Vararg{Any,R}) where {T,R} +function Base.getindex(ta::TensorApplication{T,R}, I::Vararg{Any,R}) where {T,R} @boundscheck checkbounds(ta, Int.(I)...) - return apply(ta.t, ta.o, I...) + return @inbounds apply(ta.t, ta.o, I...) end -Base.getindex(ta::LazyTensorApplication{T,1} where T, I::CartesianIndex{1}) = ta[Tuple(I)...] # Would otherwise be caught in the previous method. -Base.size(ta::LazyTensorApplication) = range_size(ta.t) +Base.@propagate_inbounds Base.getindex(ta::TensorApplication{T,1} where T, I::CartesianIndex{1}) = ta[Tuple(I)...] # Would otherwise be caught in the previous method. +Base.size(ta::TensorApplication) = range_size(ta.t) """ - LazyTensorTranspose{T,R,D} <: LazyTensor{T,D,R} + TensorTranspose{T,R,D} <: LazyTensor{T,D,R} Struct for lazy transpose of a LazyTensor. @@ -36,102 +36,104 @@ the transpose of mapping `m` by using `m'`. `m'` will work as a regular LazyTensor lazily calling the appropriate methods of `m`. """ -struct LazyTensorTranspose{T,R,D, TM<:LazyTensor{T,R,D}} <: LazyTensor{T,D,R} +struct TensorTranspose{T,R,D, TM<:LazyTensor{T,R,D}} <: LazyTensor{T,D,R} tm::TM end # # 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 LazyTensor even if it doesn't implement `apply_transpose`? -Base.adjoint(tm::LazyTensor) = LazyTensorTranspose(tm) -Base.adjoint(tmt::LazyTensorTranspose) = tmt.tm +Base.adjoint(tm::LazyTensor) = TensorTranspose(tm) +Base.adjoint(tmt::TensorTranspose) = tmt.tm -apply(tmt::LazyTensorTranspose{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) -apply_transpose(tmt::LazyTensorTranspose{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...) +apply(tmt::TensorTranspose{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) +apply_transpose(tmt::TensorTranspose{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...) -range_size(tmt::LazyTensorTranspose) = domain_size(tmt.tm) -domain_size(tmt::LazyTensorTranspose) = range_size(tmt.tm) +range_size(tmt::TensorTranspose) = domain_size(tmt.tm) +domain_size(tmt::TensorTranspose) = range_size(tmt.tm) -struct LazyTensorBinaryOperation{Op,T,R,D,T1<:LazyTensor{T,R,D},T2<:LazyTensor{T,R,D}} <: LazyTensor{T,D,R} +struct ElementwiseTensorOperation{Op,T,R,D,T1<:LazyTensor{T,R,D},T2<:LazyTensor{T,R,D}} <: LazyTensor{T,D,R} tm1::T1 tm2::T2 - function LazyTensorBinaryOperation{Op,T,R,D}(tm1::T1,tm2::T2) where {Op,T,R,D, T1<:LazyTensor{T,R,D},T2<:LazyTensor{T,R,D}} + function ElementwiseTensorOperation{Op,T,R,D}(tm1::T1,tm2::T2) where {Op,T,R,D, T1<:LazyTensor{T,R,D},T2<:LazyTensor{T,R,D}} @boundscheck check_domain_size(tm2, domain_size(tm1)) @boundscheck check_range_size(tm2, range_size(tm1)) return new{Op,T,R,D,T1,T2}(tm1,tm2) end end -LazyTensorBinaryOperation{Op}(s,t) where Op = LazyTensorBinaryOperation{Op,eltype(s), range_dim(s), domain_dim(s)}(s,t) +ElementwiseTensorOperation{Op}(s,t) where Op = ElementwiseTensorOperation{Op,eltype(s), range_dim(s), domain_dim(s)}(s,t) -apply(tmBinOp::LazyTensorBinaryOperation{:+,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::LazyTensorBinaryOperation{:-,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::ElementwiseTensorOperation{:+,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::ElementwiseTensorOperation{:-,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::LazyTensorBinaryOperation) = range_size(tmBinOp.tm1) -domain_size(tmBinOp::LazyTensorBinaryOperation) = domain_size(tmBinOp.tm1) +range_size(tmBinOp::ElementwiseTensorOperation) = range_size(tmBinOp.tm1) +domain_size(tmBinOp::ElementwiseTensorOperation) = domain_size(tmBinOp.tm1) """ - LazyTensorComposition{T,R,K,D} + TensorComposition{T,R,K,D} Lazily compose two `LazyTensor`s, so that they can be handled as a single `LazyTensor`. """ -struct LazyTensorComposition{T,R,K,D, TM1<:LazyTensor{T,R,K}, TM2<:LazyTensor{T,K,D}} <: LazyTensor{T,R,D} +struct TensorComposition{T,R,K,D, TM1<:LazyTensor{T,R,K}, TM2<:LazyTensor{T,K,D}} <: LazyTensor{T,R,D} t1::TM1 t2::TM2 - function LazyTensorComposition(t1::LazyTensor{T,R,K}, t2::LazyTensor{T,K,D}) where {T,R,K,D} + function TensorComposition(t1::LazyTensor{T,R,K}, t2::LazyTensor{T,K,D}) where {T,R,K,D} @boundscheck check_domain_size(t1, range_size(t2)) return new{T,R,K,D, typeof(t1), typeof(t2)}(t1,t2) end end -range_size(tm::LazyTensorComposition) = range_size(tm.t1) -domain_size(tm::LazyTensorComposition) = domain_size(tm.t2) +range_size(tm::TensorComposition) = range_size(tm.t1) +domain_size(tm::TensorComposition) = domain_size(tm.t2) -function apply(c::LazyTensorComposition{T,R,K,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,K,D} +function apply(c::TensorComposition{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::LazyTensorComposition{T,R,K,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,K,D} +function apply_transpose(c::TensorComposition{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 """ - LazyTensorComposition(tm, tmi::IdentityTensor) - LazyTensorComposition(tmi::IdentityTensor, tm) + TensorComposition(tm, tmi::IdentityTensor) + TensorComposition(tmi::IdentityTensor, tm) Composes a `Tensormapping` `tm` with an `IdentityTensor` `tmi`, by returning `tm` """ -function LazyTensorComposition(tm::LazyTensor{T,R,D}, tmi::IdentityTensor{T,D}) where {T,R,D} +function TensorComposition(tm::LazyTensor{T,R,D}, tmi::IdentityTensor{T,D}) where {T,R,D} @boundscheck check_domain_size(tm, range_size(tmi)) return tm end -function LazyTensorComposition(tmi::IdentityTensor{T,R}, tm::LazyTensor{T,R,D}) where {T,R,D} +function TensorComposition(tmi::IdentityTensor{T,R}, tm::LazyTensor{T,R,D}) where {T,R,D} @boundscheck check_domain_size(tmi, range_size(tm)) return tm end # Specialization for the case where tm is an IdentityTensor. Required to resolve ambiguity. -function LazyTensorComposition(tm::IdentityTensor{T,D}, tmi::IdentityTensor{T,D}) where {T,D} +function TensorComposition(tm::IdentityTensor{T,D}, tmi::IdentityTensor{T,D}) where {T,D} @boundscheck check_domain_size(tm, range_size(tmi)) return tmi end +Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm) +Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm """ - InflatedLazyTensor{T,R,D} <: LazyTensor{T,R,D} + InflatedTensor{T,R,D} <: LazyTensor{T,R,D} An inflated `LazyTensor` with dimensions added before and afer its actual dimensions. """ -struct InflatedLazyTensor{T,R,D,D_before,R_middle,D_middle,D_after, TM<:LazyTensor{T,R_middle,D_middle}} <: LazyTensor{T,R,D} +struct InflatedTensor{T,R,D,D_before,R_middle,D_middle,D_after, TM<:LazyTensor{T,R_middle,D_middle}} <: LazyTensor{T,R,D} before::IdentityTensor{T,D_before} tm::TM after::IdentityTensor{T,D_after} - function InflatedLazyTensor(before, tm::LazyTensor{T}, after) where T + function InflatedTensor(before, tm::LazyTensor{T}, after) where T R_before = range_dim(before) R_middle = range_dim(tm) R_after = range_dim(after) @@ -146,33 +148,35 @@ end """ - InflatedLazyTensor(before, tm, after) - InflatedLazyTensor(before,tm) - InflatedLazyTensor(tm,after) + InflatedTensor(before, tm, after) + InflatedTensor(before,tm) + InflatedTensor(tm,after) The outer product of `before`, `tm` and `after`, where `before` and `after` are `IdentityTensor`s. If one of `before` or `after` is left out, a 0-dimensional `IdentityTensor` is used as the default value. -If `tm` already is an `InflatedLazyTensor`, `before` and `after` will be extended instead of -creating a nested `InflatedLazyTensor`. +If `tm` already is an `InflatedTensor`, `before` and `after` will be extended instead of +creating a nested `InflatedTensor`. """ -InflatedLazyTensor(::IdentityTensor, ::LazyTensor, ::IdentityTensor) +InflatedTensor(::IdentityTensor, ::LazyTensor, ::IdentityTensor) -function InflatedLazyTensor(before, itm::InflatedLazyTensor, after) - return InflatedLazyTensor( +function InflatedTensor(before, itm::InflatedTensor, after) + return InflatedTensor( IdentityTensor(before.size..., itm.before.size...), itm.tm, IdentityTensor(itm.after.size..., after.size...), ) end -InflatedLazyTensor(before::IdentityTensor, tm::LazyTensor{T}) where T = InflatedLazyTensor(before,tm,IdentityTensor{T}()) -InflatedLazyTensor(tm::LazyTensor{T}, after::IdentityTensor) where T = InflatedLazyTensor(IdentityTensor{T}(),tm,after) +InflatedTensor(before::IdentityTensor, tm::LazyTensor{T}) where T = InflatedTensor(before,tm,IdentityTensor{T}()) +InflatedTensor(tm::LazyTensor{T}, after::IdentityTensor) where T = InflatedTensor(IdentityTensor{T}(),tm,after) # Resolve ambiguity between the two previous methods -InflatedLazyTensor(I1::IdentityTensor{T}, I2::IdentityTensor{T}) where T = InflatedLazyTensor(I1,I2,IdentityTensor{T}()) +InflatedTensor(I1::IdentityTensor{T}, I2::IdentityTensor{T}) where T = InflatedTensor(I1,I2,IdentityTensor{T}()) -function range_size(itm::InflatedLazyTensor) +# TODO: Implement some pretty printing in terms of ⊗. E.g InflatedTensor(I(3),B,I(2)) -> I(3)⊗B⊗I(2) + +function range_size(itm::InflatedTensor) return flatten_tuple( range_size(itm.before), range_size(itm.tm), @@ -180,7 +184,7 @@ ) end -function domain_size(itm::InflatedLazyTensor) +function domain_size(itm::InflatedTensor) return flatten_tuple( domain_size(itm.before), domain_size(itm.tm), @@ -188,7 +192,7 @@ ) end -function apply(itm::InflatedLazyTensor{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} +function apply(itm::InflatedTensor{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) @@ -200,7 +204,7 @@ return apply(itm.tm, v_inner, inner_index...) end -function apply_transpose(itm::InflatedLazyTensor{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} +function apply_transpose(itm::InflatedTensor{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) @@ -212,7 +216,7 @@ return apply_transpose(itm.tm, v_inner, inner_index...) end -function Base.show(io::IO, ::MIME"text/plain", tm::InflatedLazyTensor{T}) where T +function Base.show(io::IO, ::MIME"text/plain", tm::InflatedTensor{T}) where T show(IOContext(io, :compact=>true), MIME("text/plain"), tm.before) print(io, "⊗") # if get(io, :compact, false) @@ -225,7 +229,7 @@ @doc raw""" LazyOuterProduct(tms...) -Creates a `LazyTensorComposition` for the outerproduct of `tms...`. +Creates a `TensorComposition` for the outerproduct of `tms...`. This is done by separating the outer product into regular products of outer products involving only identity mappings and one non-identity mapping. First let @@ -262,15 +266,15 @@ function LazyOuterProduct end function LazyOuterProduct(tm1::LazyTensor{T}, tm2::LazyTensor{T}) where T - itm1 = InflatedLazyTensor(tm1, IdentityTensor{T}(range_size(tm2))) - itm2 = InflatedLazyTensor(IdentityTensor{T}(domain_size(tm1)),tm2) + itm1 = InflatedTensor(tm1, IdentityTensor{T}(range_size(tm2))) + itm2 = InflatedTensor(IdentityTensor{T}(domain_size(tm1)),tm2) return itm1∘itm2 end LazyOuterProduct(t1::IdentityTensor{T}, t2::IdentityTensor{T}) where T = IdentityTensor{T}(t1.size...,t2.size...) -LazyOuterProduct(t1::LazyTensor, t2::IdentityTensor) = InflatedLazyTensor(t1, t2) -LazyOuterProduct(t1::IdentityTensor, t2::LazyTensor) = InflatedLazyTensor(t1, t2) +LazyOuterProduct(t1::LazyTensor, t2::IdentityTensor) = InflatedTensor(t1, t2) +LazyOuterProduct(t1::IdentityTensor, t2::LazyTensor) = InflatedTensor(t1, t2) LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms)
--- a/src/LazyTensors/tensor_types.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/LazyTensors/tensor_types.jl Thu May 12 21:52:47 2022 +0200 @@ -2,7 +2,7 @@ IdentityTensor{T,D} <: LazyTensor{T,D,D} The lazy identity LazyTensor for a given size. Usefull for building up higher dimensional tensor mappings from lower -dimensional ones through outer products. Also used in the Implementation for InflatedLazyTensor. +dimensional ones through outer products. Also used in the Implementation for InflatedTensor. """ struct IdentityTensor{T,D} <: LazyTensor{T,D,D} size::NTuple{D,Int} @@ -60,20 +60,38 @@ end """ - LazyLinearMap{T,R,D,...}(A, range_indicies, domain_indicies) + DiagonalTensor{T,D,...} <: LazyTensor{T,D,D} + DiagonalTensor(a::AbstractArray) + +A lazy tensor with diagonal `a`. +""" +struct DiagonalTensor{T,D,AT<:AbstractArray{T,D}} <: LazyTensor{T,D,D} + diagonal::AT +end + +range_size(tm::DiagonalTensor) = size(tm.diagonal) +domain_size(tm::DiagonalTensor) = size(tm.diagonal) + + +LazyTensors.apply(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...] +LazyTensors.apply_transpose(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...] + + +""" + DenseTensor{T,R,D,...}(A, range_indicies, domain_indicies) LazyTensor defined by the AbstractArray A. `range_indicies` and `domain_indicies` define which indicies of A should be considerd the range and domain of the LazyTensor. Each set of indices must be ordered in ascending order. -For instance, if A is a m x n matrix, and range_size = (1,), domain_size = (2,), then the LazyLinearMap performs the +For instance, if A is a m x n matrix, and range_size = (1,), domain_size = (2,), then the DenseTensor performs the standard matrix-vector product on vectors of size n. """ -struct LazyLinearMap{T,R,D, RD, AA<:AbstractArray{T,RD}} <: LazyTensor{T,R,D} +struct DenseTensor{T,R,D, RD, AA<:AbstractArray{T,RD}} <: LazyTensor{T,R,D} A::AA range_indicies::NTuple{R,Int} domain_indicies::NTuple{D,Int} - function LazyLinearMap(A::AA, range_indicies::NTuple{R,Int}, domain_indicies::NTuple{D,Int}) where {T,R,D, RD, AA<:AbstractArray{T,RD}} + function DenseTensor(A::AA, range_indicies::NTuple{R,Int}, domain_indicies::NTuple{D,Int}) where {T,R,D, RD, AA<:AbstractArray{T,RD}} if !issorted(range_indicies) || !issorted(domain_indicies) throw(DomainError("range_indicies and domain_indicies must be sorted in ascending order")) end @@ -82,10 +100,10 @@ end end -range_size(llm::LazyLinearMap) = size(llm.A)[[llm.range_indicies...]] -domain_size(llm::LazyLinearMap) = size(llm.A)[[llm.domain_indicies...]] +range_size(llm::DenseTensor) = size(llm.A)[[llm.range_indicies...]] +domain_size(llm::DenseTensor) = size(llm.A)[[llm.domain_indicies...]] -function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} +function apply(llm::DenseTensor{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]) @@ -94,6 +112,6 @@ return sum(A_view.*v) end -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...) +function apply_transpose(llm::DenseTensor{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} + apply(DenseTensor(llm.A, llm.domain_indicies, llm.range_indicies), v, I...) end
--- a/src/SbpOperators/SbpOperators.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/SbpOperators.jl Thu May 12 21:52:47 2022 +0200 @@ -1,5 +1,15 @@ module SbpOperators +# Stencil set +export StencilSet +export read_stencil_set +export get_stencil_set +export parse_stencil +export parse_scalar +export parse_tuple +export sbp_operators_path + +# Operators export boundary_quadrature export boundary_restriction export inner_product @@ -20,7 +30,7 @@ end include("stencil.jl") -include("readoperator.jl") +include("stencil_set.jl") include("volumeops/volume_operator.jl") include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/first_derivative.jl")
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Thu May 12 21:52:47 2022 +0200 @@ -1,6 +1,3 @@ -# TODO: The type parameter closure_stencil::Stencil is required since there isnt any suitable type -# for stencil_set. We should consider adding type ::StencilSet and dispatch on that instead. -# The same goes for other operators """ boundary_restriction(grid, closure_stencil::Stencil, boundary) @@ -13,7 +10,7 @@ See also: [`boundary_operator`](@ref). """ -function boundary_restriction(grid, closure_stencil::Stencil, boundary) +function boundary_restriction(grid, closure_stencil, boundary) converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil) return SbpOperators.boundary_operator(grid, converted_stencil, boundary) end @@ -21,7 +18,6 @@ """ boundary_restriction(grid, stencil_set, boundary) -Creates a `boundary_restriction` operator on `grid` given a parsed TOML -`stencil_set`. +Creates a `boundary_restriction` operator on `grid` given a `stencil_set`. """ -boundary_restriction(grid, stencil_set, boundary) = boundary_restriction(grid, parse_stencil(stencil_set["e"]["closure"]), boundary) +boundary_restriction(grid, stencil_set::StencilSet, boundary) = boundary_restriction(grid, parse_stencil(stencil_set["e"]["closure"]), boundary)
--- a/src/SbpOperators/boundaryops/normal_derivative.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Thu May 12 21:52:47 2022 +0200 @@ -10,7 +10,7 @@ See also: [`boundary_operator`](@ref). """ -function normal_derivative(grid, closure_stencil::Stencil, boundary) +function normal_derivative(grid, closure_stencil, boundary) direction = dim(boundary) h_inv = inverse_spacing(grid)[direction] return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary) @@ -19,7 +19,6 @@ """ normal_derivative(grid, stencil_set, boundary) -Creates a `normal_derivative` operator on `grid` given a parsed TOML -`stencil_set`. +Creates a `normal_derivative` operator on `grid` given a `stencil_set`. """ -normal_derivative(grid, stencil_set, boundary) = normal_derivative(grid, parse_stencil(stencil_set["d1"]["closure"]), boundary) +normal_derivative(grid, stencil_set::StencilSet, boundary) = normal_derivative(grid, parse_stencil(stencil_set["d1"]["closure"]), boundary)
--- a/src/SbpOperators/readoperator.jl Mon Mar 21 09:32:26 2022 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,160 +0,0 @@ -using TOML - -export read_stencil_set -export get_stencil_set - -export parse_stencil -export parse_scalar -export parse_tuple - -export sbp_operators_path - - -""" - read_stencil_set(filename; filters) - -Picks out a stencil set from a TOML file based on some key-value -filters. If more than one set matches the filters an error is raised. The -returned stencil set contains parsed TOML intended for functions like -`parse_scalar` and `parse_stencil`. - -The stencil set is not parsed beyond the inital TOML parse. To get usable -stencils use the `parse_stencil` functions on the fields of the stencil set. - -The reason for this is that since stencil sets are intended to be very -general, and currently do not include any way to specify how to parse a given -section, the exact parsing is left to the user. - -For more information see [Operator file format](@ref) in the documentation. - -See also [`sbp_operators_path`](@ref), [`get_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref),. -""" -read_stencil_set(filename; filters...) = get_stencil_set(TOML.parsefile(filename); filters...) - -""" - get_stencil_set(parsed_toml; filters...) - -Picks out a stencil set from an already parsed TOML based on some key-value -filters. - -See also [`read_stencil_set`](@ref). -""" -function get_stencil_set(parsed_toml; filters...) - matches = findall(parsed_toml["stencil_set"]) do set - for (key, val) ∈ filters - if set[string(key)] != val - return false - end - end - - return true - end - - if length(matches) != 1 - throw(ArgumentError("filters must pick out a single stencil set")) - end - - i = matches[1] - return parsed_toml["stencil_set"][i] -end - -""" - parse_stencil(parsed_toml) - -Accepts parsed TOML and reads it as a stencil. - -See also [`read_stencil_set`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref). -""" -function parse_stencil(parsed_toml) - check_stencil_toml(parsed_toml) - - if parsed_toml isa Array - weights = parse_rational.(parsed_toml) - return CenteredStencil(weights...) - end - - weights = parse_rational.(parsed_toml["s"]) - return Stencil(weights..., center = parsed_toml["c"]) -end - -""" - parse_stencil(T, parsed_toml) - -Parses the input as a stencil with element type `T`. -""" -parse_stencil(T, parsed_toml) = Stencil{T}(parse_stencil(parsed_toml)) - -function check_stencil_toml(parsed_toml) - if !(parsed_toml isa Dict || parsed_toml isa Vector{String}) - throw(ArgumentError("the TOML for a stencil must be a vector of strings or a table.")) - end - - if parsed_toml isa Vector{String} - return - end - - if !(haskey(parsed_toml, "s") && haskey(parsed_toml, "c")) - throw(ArgumentError("the table form of a stencil must have fields `s` and `c`.")) - end - - if !(parsed_toml["s"] isa Vector{String}) - throw(ArgumentError("a stencil must be specified as a vector of strings.")) - end - - if !(parsed_toml["c"] isa Int) - throw(ArgumentError("the center of a stencil must be specified as an integer.")) - end -end - -""" - parse_scalar(parsed_toml) - -Parse a scalar, represented as a string or a number in the TOML, and return it as a `Rational` - -See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref) [`parse_tuple`](@ref). -""" -function parse_scalar(parsed_toml) - try - return parse_rational(parsed_toml) - catch e - throw(ArgumentError("must be a number or a string representing a number.")) - end -end - -""" - parse_tuple(parsed_toml) - -Parse an array as a tuple of scalars. - -See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref). -""" -function parse_tuple(parsed_toml) - if !(parsed_toml isa Array) - throw(ArgumentError("argument must be an array")) - end - return Tuple(parse_scalar.(parsed_toml)) -end - - -""" - parse_rational(parsed_toml) - -Parse a string or a number as a rational. -""" -function parse_rational(parsed_toml) - if parsed_toml isa String - expr = Meta.parse(replace(parsed_toml, "/"=>"//")) - return eval(:(Rational($expr))) - else - return Rational(parsed_toml) - end -end - -""" - sbp_operators_path() - -Calculate the path for the operators folder with included stencil sets. - -See also [`read_stencil_set`](@ref) -""" -sbp_operators_path() = (@__DIR__) * "/operators/"
--- a/src/SbpOperators/stencil.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/stencil.jl Thu May 12 21:52:47 2022 +0200 @@ -1,11 +1,12 @@ export CenteredStencil +export CenteredNestedStencil struct Stencil{T,N} - range::Tuple{Int,Int} + range::UnitRange{Int64} weights::NTuple{N,T} - function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N} - @assert range[2]-range[1]+1 == N + function Stencil(range::UnitRange,weights::NTuple{N,T}) where {T, N} + @assert length(range) == N new{T,N}(range,weights) end end @@ -15,27 +16,30 @@ Create a stencil with the given weights with element `center` as the center of the stencil. """ -function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error +function Stencil(weights...; center::Int) + weights = promote(weights...) N = length(weights) - range = (1, N) .- center + range = (1:N) .- center return Stencil(range, weights) end -function Stencil{T}(s::Stencil) where T - return Stencil(s.range, T.(s.weights)) -end +Stencil{T,N}(s::Stencil{S,N}) where {T,S,N} = Stencil(s.range, T.(s.weights)) +Stencil{T}(s::Stencil) where T = Stencil{T,length(s)}(s) -Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil) +Base.convert(::Type{Stencil{T1,N}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s) +Base.convert(::Type{Stencil{T1}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s) -function CenteredStencil(weights::Vararg) +Base.promote_rule(::Type{Stencil{T1,N}}, ::Type{Stencil{T2,N}}) where {T1,T2,N} = Stencil{promote_type(T1,T2),N} + +function CenteredStencil(weights...) if iseven(length(weights)) throw(ArgumentError("a centered stencil must have an odd number of weights.")) end r = length(weights) ÷ 2 - return Stencil((-r, r), weights) + return Stencil(-r:r, weights) end @@ -48,7 +52,8 @@ return Stencil(s.range, a.*s.weights) end -Base.eltype(::Stencil{T}) where T = T +Base.eltype(::Stencil{T,N}) where {T,N} = T +Base.length(::Stencil{T,N}) where {T,N} = N function flip(s::Stencil) range = (-s.range[2], -s.range[1]) @@ -57,24 +62,103 @@ # Provides index into the Stencil based on offset for the root element @inline function Base.getindex(s::Stencil, i::Int) - @boundscheck if i < s.range[1] || s.range[2] < i + @boundscheck if i ∉ s.range return zero(eltype(s)) end return s.weights[1 + i - s.range[1]] end -Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N} - w = s.weights[1]*v[i + s.range[1]] - @simd for k ∈ 2:N - w += s.weights[k]*v[i + s.range[1] + k-1] +Base.@propagate_inbounds @inline function apply_stencil(s::Stencil, v::AbstractVector, i::Int) + w = zero(promote_type(eltype(s),eltype(v))) + @simd for k ∈ 1:length(s) + w += s.weights[k]*v[i + s.range[k]] + end + + return w +end + +Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil, v::AbstractVector, i::Int) + w = zero(promote_type(eltype(s),eltype(v))) + @simd for k ∈ length(s):-1:1 + w += s.weights[k]*v[i - s.range[k]] end return w end -Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N} - w = s.weights[N]*v[i - s.range[2]] - @simd for k ∈ N-1:-1:1 - w += s.weights[k]*v[i - s.range[1] - k + 1] - end - return w + +struct NestedStencil{T,N,M} + s::Stencil{Stencil{T,N},M} +end + +# Stencil input +NestedStencil(s::Vararg{Stencil}; center) = NestedStencil(Stencil(s... ; center)) +CenteredNestedStencil(s::Vararg{Stencil}) = NestedStencil(CenteredStencil(s...)) + +# Tuple input +function NestedStencil(weights::Vararg{NTuple{N,Any}}; center) where N + inner_stencils = map(w -> Stencil(w...; center), weights) + return NestedStencil(Stencil(inner_stencils... ; center)) +end +function CenteredNestedStencil(weights::Vararg{NTuple{N,Any}}) where N + inner_stencils = map(w->CenteredStencil(w...), weights) + return CenteredNestedStencil(inner_stencils...) +end + + +# Conversion +function NestedStencil{T,N,M}(ns::NestedStencil{S,N,M}) where {T,S,N,M} + return NestedStencil(Stencil{Stencil{T}}(ns.s)) +end + +function NestedStencil{T}(ns::NestedStencil{S,N,M}) where {T,S,N,M} + NestedStencil{T,N,M}(ns) +end + +function Base.convert(::Type{NestedStencil{T,N,M}}, s::NestedStencil{S,N,M}) where {T,S,N,M} + return NestedStencil{T,N,M}(s) +end +Base.convert(::Type{NestedStencil{T}}, stencil) where T = NestedStencil{T}(stencil) + +function Base.promote_rule(::Type{NestedStencil{T,N,M}}, ::Type{NestedStencil{S,N,M}}) where {T,S,N,M} + return NestedStencil{promote_type(T,S),N,M} end + +Base.eltype(::NestedStencil{T}) where T = T + +function scale(ns::NestedStencil, a) + range = ns.s.range + weights = ns.s.weights + + return NestedStencil(Stencil(range, scale.(weights,a))) +end + +function flip(ns::NestedStencil) + s_flip = flip(ns.s) + return NestedStencil(Stencil(s_flip.range, flip.(s_flip.weights))) +end + +Base.getindex(ns::NestedStencil, i::Int) = ns.s[i] + +"Apply inner stencils to `c` and get a concrete stencil" +Base.@propagate_inbounds function apply_inner_stencils(ns::NestedStencil, c::AbstractVector, i::Int) + weights = apply_stencil.(ns.s.weights, Ref(c), i) + return Stencil(ns.s.range, weights) +end + +"Apply the whole nested stencil" +Base.@propagate_inbounds function apply_stencil(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int) + s = apply_inner_stencils(ns,c,i) + return apply_stencil(s, v, i) +end + +"Apply inner stencils backwards to `c` and get a concrete stencil" +Base.@propagate_inbounds @inline function apply_inner_stencils_backwards(ns::NestedStencil, c::AbstractVector, i::Int) + weights = apply_stencil_backwards.(ns.s.weights, Ref(c), i) + return Stencil(ns.s.range, weights) +end + +"Apply the whole nested stencil backwards" +Base.@propagate_inbounds @inline function apply_stencil_backwards(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int) + s = apply_inner_stencils_backwards(ns,c,i) + return apply_stencil_backwards(s, v, i) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/stencil_set.jl Thu May 12 21:52:47 2022 +0200 @@ -0,0 +1,164 @@ +using TOML + + +""" + StencilSet + +A `StencilSet` contains a set of associated stencils. The stencils +are are stored in a table, and can be accesed by indexing into the `StencilSet`. +""" +struct StencilSet + table +end +Base.getindex(set::StencilSet,I...) = set.table[I...] + + +""" +read_stencil_set(filename; filters) + +Creates a `StencilSet` from a TOML file based on some key-value +filters. If more than one set matches the filters an error is raised. The +table of the `StencilSet` is a parsed TOML intended for functions like +`parse_scalar` and `parse_stencil`. + +The `StencilSet` table is not parsed beyond the inital TOML parse. To get usable +stencils use the `parse_stencil` functions on the fields of the stencil set. + +The reason for this is that since stencil sets are intended to be very +general, and currently do not include any way to specify how to parse a given +section, the exact parsing is left to the user. + +For more information see [Operator file format](@ref) in the documentation. + +See also [`StencilSet`](@ref), [`sbp_operators_path`](@ref), [`get_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref). +""" +read_stencil_set(filename; filters...) = StencilSet(get_stencil_set(TOML.parsefile(filename); filters...)) + + +""" + get_stencil_set(parsed_toml; filters...) + +Picks out a stencil set from an already parsed TOML based on some key-value +filters. + +See also [`read_stencil_set`](@ref). +""" +function get_stencil_set(parsed_toml; filters...) + matches = findall(parsed_toml["stencil_set"]) do set + for (key, val) ∈ filters + if set[string(key)] != val + return false + end + end + + return true + end + + if length(matches) != 1 + throw(ArgumentError("filters must pick out a single stencil set")) + end + + i = matches[1] + return parsed_toml["stencil_set"][i] +end + +""" + parse_stencil(parsed_toml) + +Accepts parsed TOML and reads it as a stencil. + +See also [`read_stencil_set`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref). +""" +function parse_stencil(parsed_toml) + check_stencil_toml(parsed_toml) + + if parsed_toml isa Array + weights = parse_rational.(parsed_toml) + return CenteredStencil(weights...) + end + + weights = parse_rational.(parsed_toml["s"]) + return Stencil(weights..., center = parsed_toml["c"]) +end + +""" + parse_stencil(T, parsed_toml) + +Parses the input as a stencil with element type `T`. +""" +parse_stencil(T, parsed_toml) = Stencil{T}(parse_stencil(parsed_toml)) + +function check_stencil_toml(parsed_toml) + if !(parsed_toml isa Dict || parsed_toml isa Vector{String}) + throw(ArgumentError("the TOML for a stencil must be a vector of strings or a table.")) + end + + if parsed_toml isa Vector{String} + return + end + + if !(haskey(parsed_toml, "s") && haskey(parsed_toml, "c")) + throw(ArgumentError("the table form of a stencil must have fields `s` and `c`.")) + end + + if !(parsed_toml["s"] isa Vector{String}) + throw(ArgumentError("a stencil must be specified as a vector of strings.")) + end + + if !(parsed_toml["c"] isa Int) + throw(ArgumentError("the center of a stencil must be specified as an integer.")) + end +end + +""" + parse_scalar(parsed_toml) + +Parse a scalar, represented as a string or a number in the TOML, and return it as a `Rational` + +See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref) [`parse_tuple`](@ref). +""" +function parse_scalar(parsed_toml) + try + return parse_rational(parsed_toml) + catch e + throw(ArgumentError("must be a number or a string representing a number.")) + end +end + +""" + parse_tuple(parsed_toml) + +Parse an array as a tuple of scalars. + +See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref). +""" +function parse_tuple(parsed_toml) + if !(parsed_toml isa Array) + throw(ArgumentError("argument must be an array")) + end + return Tuple(parse_scalar.(parsed_toml)) +end + + +""" + parse_rational(parsed_toml) + +Parse a string or a number as a rational. +""" +function parse_rational(parsed_toml) + if parsed_toml isa String + expr = Meta.parse(replace(parsed_toml, "/"=>"//")) + return eval(:(Rational($expr))) + else + return Rational(parsed_toml) + end +end + +""" + sbp_operators_path() + +Calculate the path for the operators folder with included stencil sets. + +See also [`StencilSet`](@ref), [`read_stencil_set`](@ref). +""" +sbp_operators_path() = (@__DIR__) * "/operators/"
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl Thu May 12 21:52:47 2022 +0200 @@ -16,16 +16,31 @@ h_inv = inverse_spacing(grid)[direction] return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction) end -first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid,inner_stencil,closure_stencils,1) + """ - first_derivative(grid, stencil_set, direction) + first_derivative(grid, inner_stencil, closure_stencils) + +Creates a `first_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`. +""" +first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid, inner_stencil, closure_stencils, 1) + -Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a parsed TOML -`stencil_set`. """ -function first_derivative(grid::EquidistantGrid, stencil_set, direction) + first_derivative(grid, stencil_set::StencilSet, direction) + +Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`. +""" +function first_derivative(grid::EquidistantGrid, stencil_set::StencilSet, direction) inner_stencil = parse_stencil(stencil_set["D1"]["inner_stencil"]) closure_stencils = parse_stencil.(stencil_set["D1"]["closure_stencils"]) first_derivative(grid,inner_stencil,closure_stencils,direction); end + + +""" + first_derivative(grid, stencil_set) + +Creates a `first_derivative` operator on a 1D `grid` given a `stencil_set`. +""" +first_derivative(grid::EquidistantGrid{1}, stencil_set::StencilSet) = first_derivative(grid, stencil_set, 1)
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl Thu May 12 21:52:47 2022 +0200 @@ -16,16 +16,31 @@ h_inv = inverse_spacing(grid)[direction] return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction) end -second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1) + + +""" + second_derivative(grid, inner_stencil, closure_stencils) + +Creates a `second_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`. +""" +second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid, inner_stencil, closure_stencils,1) + """ second_derivative(grid, stencil_set, direction) -Creates a `second_derivative` operator on `grid` along coordinate dimension `direction` given a parsed TOML -`stencil_set`. +Creates a `second_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`. """ -function second_derivative(grid::EquidistantGrid, stencil_set, direction) +function second_derivative(grid::EquidistantGrid, stencil_set::StencilSet, direction) inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) second_derivative(grid,inner_stencil,closure_stencils,direction); -end +end + + +""" + second_derivative(grid, stencil_set) + +Creates a `second_derivative` operator on a 1D `grid` given a `stencil_set`. +""" +second_derivative(grid::EquidistantGrid{1}, stencil_set::StencilSet) = second_derivative(grid, stencil_set, 1) \ No newline at end of file
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl Thu May 12 21:52:47 2022 +0200 @@ -37,10 +37,9 @@ """ inner_product(grid, stencil_set) -Creates a `inner_product` operator on `grid` given a parsed TOML -`stencil_set`. +Creates a `inner_product` operator on `grid` given a `stencil_set`. """ -function inner_product(grid, stencil_set) +function inner_product(grid, stencil_set::StencilSet) inner_stencil = parse_scalar(stencil_set["H"]["inner"]) closure_stencils = parse_tuple(stencil_set["H"]["closure"]) return inner_product(grid, inner_stencil, closure_stencils)
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Thu May 12 21:52:47 2022 +0200 @@ -33,10 +33,9 @@ """ inverse_inner_product(grid, stencil_set) -Creates a `inverse_inner_product` operator on `grid` given a parsed TOML -`stencil_set`. +Creates a `inverse_inner_product` operator on `grid` given a `stencil_set`. """ -function inverse_inner_product(grid, stencil_set) +function inverse_inner_product(grid, stencil_set::StencilSet) inner_stencil = parse_scalar(stencil_set["H"]["inner"]) closure_stencils = parse_tuple(stencil_set["H"]["closure"]) return inverse_inner_product(grid, inner_stencil, closure_stencils)
--- a/src/SbpOperators/volumeops/laplace/laplace.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Thu May 12 21:52:47 2022 +0200 @@ -2,21 +2,22 @@ Laplace{T, Dim, TM} <: LazyTensor{T, Dim, Dim} Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a -`LazyTensor`. Additionally `Laplace` stores the stencil set (parsed from TOML) -used to construct the `LazyTensor`. +`LazyTensor`. Additionally `Laplace` stores the `StencilSet` +used to construct the `LazyTensor `. """ struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim} D::TM # Difference operator - stencil_set # Stencil set of the operator + stencil_set::StencilSet # Stencil set of the operator end """ Laplace(grid::Equidistant, stencil_set) -Creates the `Laplace` operator `Δ` on `grid` given a parsed TOML -`stencil_set`. See also [`laplace`](@ref). +Creates the `Laplace` operator `Δ` on `grid` given a `stencil_set`. + +See also [`laplace`](@ref). """ -function Laplace(grid::EquidistantGrid, stencil_set) +function Laplace(grid::EquidistantGrid, stencil_set::StencilSet) inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) Δ = laplace(grid, inner_stencil,closure_stencils)
--- a/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Thu May 12 21:52:47 2022 +0200 @@ -24,7 +24,7 @@ end """ - VolumeOperator{T,N,M,K} <: TensorOperator{T,1} + VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1} Implements a one-dimensional constant coefficients volume operator """ struct VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
--- a/test/LazyTensors/lazy_tensor_operations_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Thu May 12 21:52:47 2022 +0200 @@ -36,7 +36,7 @@ end -@testset "LazyTensorApplication" begin +@testset "TensorApplication" begin m = SizeDoublingMapping{Int, 1, 1}((3,)) mm = SizeDoublingMapping{Int, 1, 1}((6,)) v = [0,1,2] @@ -159,14 +159,14 @@ end -@testset "LazyTensorComposition" begin +@testset "TensorComposition" begin A = rand(2,3) B = rand(3,4) - Ã = LazyLinearMap(A, (1,), (2,)) - B̃ = LazyLinearMap(B, (1,), (2,)) + Ã = DenseTensor(A, (1,), (2,)) + B̃ = DenseTensor(B, (1,), (2,)) - @test Ã∘B̃ isa LazyTensorComposition + @test Ã∘B̃ isa TensorComposition @test range_size(Ã∘B̃) == (2,) @test domain_size(Ã∘B̃) == (4,) @test_throws DomainSizeMismatch B̃∘Ã @@ -181,114 +181,112 @@ @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64 @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64 + + a = 2. + v = rand(3) + @test a*Ã isa TensorComposition + @test a*Ã == Ã*a + @test range_size(a*Ã) == range_size(Ã) + @test domain_size(a*Ã) == domain_size(Ã) + @test a*Ã*v == a.*A*v end -@testset "InflatedLazyTensor" begin +@testset "InflatedTensor" begin I(sz...) = IdentityTensor(sz...) Ã = rand(4,2) B̃ = rand(4,2,3) C̃ = rand(4,2,3) - A = LazyLinearMap(Ã,(1,),(2,)) - B = LazyLinearMap(B̃,(1,2),(3,)) - C = LazyLinearMap(C̃,(1,),(2,3)) + A = DenseTensor(Ã,(1,),(2,)) + B = DenseTensor(B̃,(1,2),(3,)) + C = DenseTensor(C̃,(1,),(2,3)) @testset "Constructors" begin - @test InflatedLazyTensor(I(3,2), A, I(4)) isa LazyTensor{Float64, 4, 4} - @test InflatedLazyTensor(I(3,2), B, I(4)) isa LazyTensor{Float64, 5, 4} - @test InflatedLazyTensor(I(3), C, I(2,3)) isa LazyTensor{Float64, 4, 5} - @test InflatedLazyTensor(C, I(2,3)) isa LazyTensor{Float64, 3, 4} - @test InflatedLazyTensor(I(3), C) isa LazyTensor{Float64, 2, 3} - @test InflatedLazyTensor(I(3), I(2,3)) isa LazyTensor{Float64, 3, 3} + @test InflatedTensor(I(3,2), A, I(4)) isa LazyTensor{Float64, 4, 4} + @test InflatedTensor(I(3,2), B, I(4)) isa LazyTensor{Float64, 5, 4} + @test InflatedTensor(I(3), C, I(2,3)) isa LazyTensor{Float64, 4, 5} + @test InflatedTensor(C, I(2,3)) isa LazyTensor{Float64, 3, 4} + @test InflatedTensor(I(3), C) isa LazyTensor{Float64, 2, 3} + @test InflatedTensor(I(3), I(2,3)) isa LazyTensor{Float64, 3, 3} end @testset "Range and domain size" begin - @test range_size(InflatedLazyTensor(I(3,2), A, I(4))) == (3,2,4,4) - @test domain_size(InflatedLazyTensor(I(3,2), A, I(4))) == (3,2,2,4) + @test range_size(InflatedTensor(I(3,2), A, I(4))) == (3,2,4,4) + @test domain_size(InflatedTensor(I(3,2), A, I(4))) == (3,2,2,4) - @test range_size(InflatedLazyTensor(I(3,2), B, I(4))) == (3,2,4,2,4) - @test domain_size(InflatedLazyTensor(I(3,2), B, I(4))) == (3,2,3,4) + @test range_size(InflatedTensor(I(3,2), B, I(4))) == (3,2,4,2,4) + @test domain_size(InflatedTensor(I(3,2), B, I(4))) == (3,2,3,4) - @test range_size(InflatedLazyTensor(I(3), C, I(2,3))) == (3,4,2,3) - @test domain_size(InflatedLazyTensor(I(3), C, I(2,3))) == (3,2,3,2,3) + @test range_size(InflatedTensor(I(3), C, I(2,3))) == (3,4,2,3) + @test domain_size(InflatedTensor(I(3), C, I(2,3))) == (3,2,3,2,3) - @inferred range_size(InflatedLazyTensor(I(3,2), A, I(4))) == (3,2,4,4) - @inferred domain_size(InflatedLazyTensor(I(3,2), A, I(4))) == (3,2,2,4) + @inferred range_size(InflatedTensor(I(3,2), A, I(4))) == (3,2,4,4) + @inferred domain_size(InflatedTensor(I(3,2), A, I(4))) == (3,2,2,4) end @testset "Application" begin # Testing regular application and transposed application with inflation "before", "after" and "before and after". # The inflated tensor mappings are chosen to preserve, reduce and increase the dimension of the result compared to the input. - tests = [ + cases = [ ( - InflatedLazyTensor(I(3,2), A, I(4)), + InflatedTensor(I(3,2), A, I(4)), (v-> @tullio res[a,b,c,d] := Ã[c,i]*v[a,b,i,d]), # Expected result of apply (v-> @tullio res[a,b,c,d] := Ã[i,c]*v[a,b,i,d]), # Expected result of apply_transpose ), ( - InflatedLazyTensor(I(3,2), B, I(4)), + InflatedTensor(I(3,2), B, I(4)), (v-> @tullio res[a,b,c,d,e] := B̃[c,d,i]*v[a,b,i,e]), (v-> @tullio res[a,b,c,d] := B̃[i,j,c]*v[a,b,i,j,d]), ), ( - InflatedLazyTensor(I(3,2), C, I(4)), + InflatedTensor(I(3,2), C, I(4)), (v-> @tullio res[a,b,c,d] := C̃[c,i,j]*v[a,b,i,j,d]), (v-> @tullio res[a,b,c,d,e] := C̃[i,c,d]*v[a,b,i,e]), ), ( - InflatedLazyTensor(I(3,2), A), + InflatedTensor(I(3,2), A), (v-> @tullio res[a,b,c] := Ã[c,i]*v[a,b,i]), (v-> @tullio res[a,b,c] := Ã[i,c]*v[a,b,i]), ), ( - InflatedLazyTensor(I(3,2), B), + InflatedTensor(I(3,2), B), (v-> @tullio res[a,b,c,d] := B̃[c,d,i]*v[a,b,i]), (v-> @tullio res[a,b,c] := B̃[i,j,c]*v[a,b,i,j]), ), ( - InflatedLazyTensor(I(3,2), C), + InflatedTensor(I(3,2), C), (v-> @tullio res[a,b,c] := C̃[c,i,j]*v[a,b,i,j]), (v-> @tullio res[a,b,c,d] := C̃[i,c,d]*v[a,b,i]), ), ( - InflatedLazyTensor(A,I(4)), + InflatedTensor(A,I(4)), (v-> @tullio res[a,b] := Ã[a,i]*v[i,b]), (v-> @tullio res[a,b] := Ã[i,a]*v[i,b]), ), ( - InflatedLazyTensor(B,I(4)), + InflatedTensor(B,I(4)), (v-> @tullio res[a,b,c] := B̃[a,b,i]*v[i,c]), (v-> @tullio res[a,b] := B̃[i,j,a]*v[i,j,b]), ), ( - InflatedLazyTensor(C,I(4)), + InflatedTensor(C,I(4)), (v-> @tullio res[a,b] := C̃[a,i,j]*v[i,j,b]), (v-> @tullio res[a,b,c] := C̃[i,a,b]*v[i,c]), ), ] - @testset "apply" begin - for i ∈ 1:length(tests) - tm = tests[i][1] - v = rand(domain_size(tm)...) - true_value = tests[i][2](v) - @test tm*v ≈ true_value rtol=1e-14 - end - end + @testset "$tm" for (tm, true_apply, true_apply_transpose) ∈ cases + v = rand(domain_size(tm)...) + @test tm*v ≈ true_apply(v) rtol=1e-14 - @testset "apply_transpose" begin - for i ∈ 1:length(tests) - tm = tests[i][1] - v = rand(range_size(tm)...) - true_value = tests[i][3](v) - @test tm'*v ≈ true_value rtol=1e-14 - end + v = rand(range_size(tm)...) + @test tm'*v ≈ true_apply_transpose(v) rtol=1e-14 end @testset "application to other type" begin - tm = InflatedLazyTensor(I(3,2), A, I(4)) + tm = InflatedTensor(I(3,2), A, I(4)) v = rand(ComplexF64, domain_size(tm)...) @test (tm*v)[1,2,3,1] isa ComplexF64 @@ -298,7 +296,7 @@ end @testset "Inference of application" begin - tm = InflatedLazyTensor(I(2,3),ScalingTensor(2.0, (3,2)),I(3,4)) + tm = InflatedTensor(I(2,3),ScalingTensor(2.0, (3,2)),I(3,4)) v = rand(domain_size(tm)...) @inferred apply(tm,v,1,2,3,2,2,4) @@ -306,19 +304,19 @@ end end - @testset "InflatedLazyTensor of InflatedLazyTensor" begin + @testset "InflatedTensor of InflatedTensor" begin A = ScalingTensor(2.0,(2,3)) - itm = InflatedLazyTensor(I(3,2), A, I(4)) - @test InflatedLazyTensor(I(4), itm, I(2)) == InflatedLazyTensor(I(4,3,2), A, I(4,2)) - @test InflatedLazyTensor(itm, I(2)) == InflatedLazyTensor(I(3,2), A, I(4,2)) - @test InflatedLazyTensor(I(4), itm) == InflatedLazyTensor(I(4,3,2), A, I(4)) + itm = InflatedTensor(I(3,2), A, I(4)) + @test InflatedTensor(I(4), itm, I(2)) == InflatedTensor(I(4,3,2), A, I(4,2)) + @test InflatedTensor(itm, I(2)) == InflatedTensor(I(3,2), A, I(4,2)) + @test InflatedTensor(I(4), itm) == InflatedTensor(I(4,3,2), A, I(4)) - @test InflatedLazyTensor(I(2), I(2), I(2)) isa InflatedLazyTensor # The constructor should always return its type. + @test InflatedTensor(I(2), I(2), I(2)) isa InflatedTensor # The constructor should always return its type. end @testset "Pretty printing" begin cases = [ - InflatedLazyTensor(I(4), ScalingTensor(2., (3,2)), I(2)) => ( + InflatedTensor(I(4), ScalingTensor(2., (3,2)), I(2)) => ( regular="I(4)⊗ScalingTensor{Float64}(2.0, (3, 2))⊗I(2)", compact="I(4)⊗2.0*I(3,2)⊗I(2)" ) @@ -359,8 +357,8 @@ v₁ = rand(2,4,3) v₂ = rand(4,3,2) - Ã = LazyLinearMap(A,(1,),(2,)) - B̃ = LazyLinearMap(B,(1,),(2,3)) + Ã = DenseTensor(A,(1,),(2,)) + B̃ = DenseTensor(B,(1,),(2,3)) ÃB̃ = LazyOuterProduct(Ã,B̃) @tullio ABv[i,k] := A[i,j]*B[k,l,m]*v₁[j,l,m] @@ -373,12 +371,12 @@ @testset "Indentity mapping arguments" begin @test LazyOuterProduct(IdentityTensor(3,2), IdentityTensor(1,2)) == IdentityTensor(3,2,1,2) - Ã = LazyLinearMap(A,(1,),(2,)) - @test LazyOuterProduct(IdentityTensor(3,2), Ã) == InflatedLazyTensor(IdentityTensor(3,2),Ã) - @test LazyOuterProduct(Ã, IdentityTensor(3,2)) == InflatedLazyTensor(Ã,IdentityTensor(3,2)) + Ã = DenseTensor(A,(1,),(2,)) + @test LazyOuterProduct(IdentityTensor(3,2), Ã) == InflatedTensor(IdentityTensor(3,2),Ã) + @test LazyOuterProduct(Ã, IdentityTensor(3,2)) == InflatedTensor(Ã,IdentityTensor(3,2)) I1 = IdentityTensor(3,2) I2 = IdentityTensor(4) - @test I1⊗Ã⊗I2 == InflatedLazyTensor(I1, Ã, I2) + @test I1⊗Ã⊗I2 == InflatedTensor(I1, Ã, I2) end end
--- a/test/LazyTensors/tensor_types_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/LazyTensors/tensor_types_test.jl Thu May 12 21:52:47 2022 +0200 @@ -1,5 +1,6 @@ using Test using Sbplib.LazyTensors +using BenchmarkTools @testset "IdentityTensor" begin @test IdentityTensor{Float64}((4,5)) isa IdentityTensor{T,2} where T @@ -32,7 +33,7 @@ @inferred domain_dim(I) à = rand(4,2) - A = LazyLinearMap(Ã,(1,),(2,)) + A = DenseTensor(Ã,(1,),(2,)) I1 = IdentityTensor{Float64}(2) I2 = IdentityTensor{Float64}(4) @test A∘I1 == A @@ -74,16 +75,57 @@ end end +@testset "DiagonalTensor" begin + @test DiagonalTensor([1,2,3,4]) isa LazyTensor{Int,1,1} + @test DiagonalTensor([1 2 3; 4 5 6]) isa LazyTensor{Int,2,2} + @test DiagonalTensor([1. 2. 3.; 4. 5. 6.]) isa LazyTensor{Float64,2,2} -@testset "LazyLinearMap" begin + @test range_size(DiagonalTensor([1,2,3,4])) == (4,) + @test domain_size(DiagonalTensor([1,2,3,4])) == (4,) + + @test range_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3) + @test domain_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3) + + @testset "apply size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)] + diag = rand(sz...) + tm = DiagonalTensor(diag) + + v = rand(sz...) + + @test tm*v == diag.*v + @test tm'*v == diag.*v + end + + + @testset "allocations size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)] + diag = rand(sz...) + tm = DiagonalTensor(diag) + + v = rand(sz...) + + @test tm*v == diag.*v + @test tm'*v == diag.*v + end + + sz = (3,2) + diag = rand(sz...) + tm = DiagonalTensor(diag) + + v = rand(sz...) + LazyTensors.apply(tm,v, 2,1) + @test (@ballocated LazyTensors.apply($tm,$v, 2,1)) == 0 +end + + +@testset "DenseTensor" begin # Test a standard matrix-vector product # mapping vectors of size 4 to vectors of size 3. A = rand(3,4) - à = LazyLinearMap(A, (1,), (2,)) + à = DenseTensor(A, (1,), (2,)) v = rand(4) w = rand(3) - @test à isa LazyLinearMap{T,1,1} where T + @test à isa DenseTensor{T,1,1} where T @test à isa LazyTensor{T,1,1} where T @test range_size(Ã) == (3,) @test domain_size(Ã) == (4,) @@ -93,12 +135,12 @@ @test Ã'*w ≈ A'*w A = rand(2,3,4) - @test_throws DomainError LazyLinearMap(A, (3,1), (2,)) + @test_throws DomainError DenseTensor(A, (3,1), (2,)) # Test more exotic mappings B = rand(3,4,2) # Map vectors of size 2 to matrices of size (3,4) - B̃ = LazyLinearMap(B, (1,2), (3,)) + B̃ = DenseTensor(B, (1,2), (3,)) v = rand(2) @test range_size(B̃) == (3,4) @@ -108,7 +150,7 @@ @test B̃*v ≈ B[:,:,1]*v[1] + B[:,:,2]*v[2] atol=5e-13 # Map matrices of size (3,2) to vectors of size 4 - B̃ = LazyLinearMap(B, (2,), (1,3)) + B̃ = DenseTensor(B, (2,), (1,3)) v = rand(3,2) @test range_size(B̃) == (4,)
--- a/test/Manifest.toml Mon Mar 21 09:32:26 2022 +0100 +++ b/test/Manifest.toml Thu May 12 21:52:47 2022 +0200 @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.7.0" +julia_version = "1.7.1" manifest_format = "2.0" [[deps.ArgTools]] @@ -12,11 +12,17 @@ [[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +[[deps.BenchmarkTools]] +deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] +git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93" +uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +version = "1.3.1" + [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47" +git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.11.2" +version = "1.14.0" [[deps.ChangesOfVariables]] deps = ["ChainRulesCore", "LinearAlgebra", "Test"] @@ -26,9 +32,9 @@ [[deps.Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582" +git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.41.0" +version = "3.42.0" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -48,10 +54,10 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" [[deps.DiffRules]] -deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] -git-tree-sha1 = "9bc5dac3c8b6706b58ad5ce24cffd9861f07c94f" +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.9.0" +version = "1.10.0" [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] @@ -78,9 +84,9 @@ [[deps.InverseFunctions]] deps = ["Test"] -git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65" +git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d" uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.2" +version = "0.1.3" [[deps.IrrationalConstants]] git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" @@ -89,9 +95,15 @@ [[deps.JLLWrappers]] deps = ["Preferences"] -git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e" +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.3.0" +version = "1.4.1" + +[[deps.JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.3" [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] @@ -118,9 +130,9 @@ [[deps.LogExpFunctions]] deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1" +git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.6" +version = "0.3.10" [[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -140,9 +152,9 @@ uuid = "14a3606d-f60d-562e-9121-12d972cd8159" [[deps.NaNMath]] -git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283" +git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "0.3.6" +version = "1.0.0" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" @@ -161,20 +173,30 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" +[[deps.Parsers]] +deps = ["Dates"] +git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.2.3" + [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [[deps.Preferences]] deps = ["TOML"] -git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a" +git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.2.2" +version = "1.2.5" [[deps.Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[[deps.Profile]] +deps = ["Printf"] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" + [[deps.REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" @@ -185,9 +207,9 @@ [[deps.Requires]] deps = ["UUIDs"] -git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a" +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.2.0" +version = "1.3.0" [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -208,9 +230,9 @@ [[deps.SpecialFunctions]] deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e" +git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.0.0" +version = "2.1.4" [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -236,9 +258,9 @@ [[deps.Tullio]] deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"] -git-tree-sha1 = "0288b7a395fc412952baf756fac94e4f28bfec65" +git-tree-sha1 = "7830c974acc69437a3fee35dd7b510a74cbc862d" uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" -version = "0.3.2" +version = "0.3.3" [[deps.UUIDs]] deps = ["Random", "SHA"]
--- a/test/Project.toml Mon Mar 21 09:32:26 2022 +0100 +++ b/test/Project.toml Thu May 12 21:52:47 2022 +0200 @@ -1,4 +1,5 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Thu May 12 21:52:47 2022 +0200 @@ -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)) @@ -28,7 +28,7 @@ @testset "2D" begin e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}()) - @test e_w isa InflatedLazyTensor + @test e_w isa InflatedTensor @test e_w isa LazyTensor{T,1,2} where T end end
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl Thu May 12 21:52:47 2022 +0200 @@ -4,7 +4,7 @@ using Sbplib.Grids using Sbplib.LazyTensors using Sbplib.RegionIndices -import Sbplib.SbpOperators.BoundaryOperator +using Sbplib.SbpOperators: BoundaryOperator, Stencil @testset "boundary_restriction" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) @@ -30,7 +30,7 @@ @testset "2D" begin e_w = boundary_restriction(g_2D,e_closure,CartesianBoundary{1,Upper}()) @test e_w == boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}()) - @test e_w isa InflatedLazyTensor + @test e_w isa InflatedTensor @test e_w isa LazyTensor{T,1,2} where T end end
--- a/test/SbpOperators/readoperator_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,172 +0,0 @@ -using Test - -using TOML -using Sbplib.SbpOperators - -import Sbplib.SbpOperators.Stencil - -@testset "readoperator" begin - toml_str = """ - [meta] - authors = "Ken Mattson" - description = "Standard operators for equidistant grids" - type = "equidistant" - cite = "A paper a long time ago in a galaxy far far away." - - [[stencil_set]] - - order = 2 - test = 2 - - H.inner = ["1"] - H.closure = ["1/2"] - - D1.inner_stencil = ["-1/2", "0", "1/2"] - D1.closure_stencils = [ - {s = ["-1", "1"], c = 1}, - ] - - D2.inner_stencil = ["1", "-2", "1"] - D2.closure_stencils = [ - {s = ["1", "-2", "1"], c = 1}, - ] - - e.closure = ["1"] - d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1} - - [[stencil_set]] - - order = 4 - test = 1 - H.inner = ["1"] - H.closure = ["17/48", "59/48", "43/48", "49/48"] - - D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"] - D2.closure_stencils = [ - {s = [ "2", "-5", "4", "-1", "0", "0"], c = 1}, - {s = [ "1", "-2", "1", "0", "0", "0"], c = 2}, - {s = [ "-4/43", "59/43", "-110/43", "59/43", "-4/43", "0"], c = 3}, - {s = [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"], c = 4}, - ] - - e.closure = ["1"] - d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1} - - [[stencil_set]] - order = 4 - test = 2 - - H.closure = ["-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"] - """ - - parsed_toml = TOML.parse(toml_str) - - @testset "get_stencil_set" begin - @test get_stencil_set(parsed_toml; order = 2) isa Dict - @test get_stencil_set(parsed_toml; order = 2) == parsed_toml["stencil_set"][1] - @test get_stencil_set(parsed_toml; test = 1) == parsed_toml["stencil_set"][2] - @test get_stencil_set(parsed_toml; order = 4, test = 2) == parsed_toml["stencil_set"][3] - - @test_throws ArgumentError get_stencil_set(parsed_toml; test = 2) - @test_throws ArgumentError get_stencil_set(parsed_toml; order = 4) - end - - @testset "parse_stencil" begin - toml = """ - s1 = ["-1/12","4/3","-5/2","4/3","-1/12"] - s2 = {s = ["2", "-5", "4", "-1", "0", "0"], c = 1} - s3 = {s = ["1", "-2", "1", "0", "0", "0"], c = 2} - s4 = "not a stencil" - s5 = [-1, 4, 3] - s6 = {k = ["1", "-2", "1", "0", "0", "0"], c = 2} - s7 = {s = [-1, 4, 3], c = 2} - s8 = {s = ["1", "-2", "1", "0", "0", "0"], c = [2,2]} - """ - - @test parse_stencil(TOML.parse(toml)["s1"]) == CenteredStencil(-1//12, 4//3, -5//2, 4//3, -1//12) - @test parse_stencil(TOML.parse(toml)["s2"]) == Stencil(2//1, -5//1, 4//1, -1//1, 0//1, 0//1; center=1) - @test parse_stencil(TOML.parse(toml)["s3"]) == Stencil(1//1, -2//1, 1//1, 0//1, 0//1, 0//1; center=2) - - @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s4"]) - @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s5"]) - @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s6"]) - @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s7"]) - @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s8"]) - - stencil_set = get_stencil_set(parsed_toml; order = 4, test = 1) - - @test parse_stencil.(stencil_set["D2"]["closure_stencils"]) == [ - Stencil( 2//1, -5//1, 4//1, -1//1, 0//1, 0//1; center=1), - Stencil( 1//1, -2//1, 1//1, 0//1, 0//1, 0//1; center=2), - Stencil(-4//43, 59//43, -110//43, 59//43, -4//43, 0//1; center=3), - Stencil(-1//49, 0//1, 59//49, -118//49, 64//49, -4//49; center=4), - ] - - - @test parse_stencil(Float64, TOML.parse(toml)["s1"]) == CenteredStencil(-1/12, 4/3, -5/2, 4/3, -1/12) - @test parse_stencil(Float64, TOML.parse(toml)["s2"]) == Stencil(2/1, -5/1, 4/1, -1/1, 0/1, 0/1; center=1) - @test parse_stencil(Float64, TOML.parse(toml)["s3"]) == Stencil(1/1, -2/1, 1/1, 0/1, 0/1, 0/1; center=2) - end - - @testset "parse_scalar" begin - toml = TOML.parse(""" - a1 = 1 - a2 = 1.5 - a3 = 1.0 - a4 = 10 - a5 = "1/2" - a6 = "1.5" - - e1 = [1,2,3] - e2 = "a string value" - """) - - @test parse_scalar(toml["a1"]) == 1//1 - @test parse_scalar(toml["a2"]) == 3//2 - @test parse_scalar(toml["a3"]) == 1//1 - @test parse_scalar(toml["a4"]) == 10//1 - @test parse_scalar(toml["a5"]) == 1//2 - @test parse_scalar(toml["a6"]) == 3//2 - - @test_throws ArgumentError parse_scalar(toml["e1"]) - @test_throws ArgumentError parse_scalar(toml["e2"]) - end - - @testset "parse_tuple" begin - toml = TOML.parse(""" - t1 = [1,3,4] - t2 = ["1/2","3/4","2/1"] - - e1 = "not a tuple" - e2.a="1" - e3 = 1 - e4 = ["1/2","3/4","not a number"] - """) - - @test parse_tuple(toml["t1"]) == (1//1,3//1,4//1) - @test parse_tuple(toml["t2"]) == (1//2,3//4,2//1) - - @test_throws ArgumentError parse_tuple(toml["e1"]) - @test_throws ArgumentError parse_tuple(toml["e2"]) - @test_throws ArgumentError parse_tuple(toml["e3"]) - @test_throws ArgumentError parse_tuple(toml["e4"]) - end -end - -@testset "parse_rational" begin - @test SbpOperators.parse_rational("1") isa Rational - @test SbpOperators.parse_rational("1") == 1//1 - @test SbpOperators.parse_rational("1/2") isa Rational - @test SbpOperators.parse_rational("1/2") == 1//2 - @test SbpOperators.parse_rational("37/13") isa Rational - @test SbpOperators.parse_rational("37/13") == 37//13 - - @test SbpOperators.parse_rational(0.5) isa Rational - @test SbpOperators.parse_rational(0.5) == 1//2 - - @test SbpOperators.parse_rational("0.5") isa Rational - @test SbpOperators.parse_rational("0.5") == 1//2 - - @test SbpOperators.parse_rational(2) isa Rational - @test SbpOperators.parse_rational(2) == 2//1 -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/stencil_set_test.jl Thu May 12 21:52:47 2022 +0200 @@ -0,0 +1,172 @@ +using Test + +using TOML +using Sbplib.SbpOperators + +import Sbplib.SbpOperators.Stencil + +@testset "readoperator" begin + toml_str = """ + [meta] + authors = "Ken Mattson" + description = "Standard operators for equidistant grids" + type = "equidistant" + cite = "A paper a long time ago in a galaxy far far away." + + [[stencil_set]] + + order = 2 + test = 2 + + H.inner = ["1"] + H.closure = ["1/2"] + + D1.inner_stencil = ["-1/2", "0", "1/2"] + D1.closure_stencils = [ + {s = ["-1", "1"], c = 1}, + ] + + D2.inner_stencil = ["1", "-2", "1"] + D2.closure_stencils = [ + {s = ["1", "-2", "1"], c = 1}, + ] + + e.closure = ["1"] + d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1} + + [[stencil_set]] + + order = 4 + test = 1 + H.inner = ["1"] + H.closure = ["17/48", "59/48", "43/48", "49/48"] + + D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"] + D2.closure_stencils = [ + {s = [ "2", "-5", "4", "-1", "0", "0"], c = 1}, + {s = [ "1", "-2", "1", "0", "0", "0"], c = 2}, + {s = [ "-4/43", "59/43", "-110/43", "59/43", "-4/43", "0"], c = 3}, + {s = [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"], c = 4}, + ] + + e.closure = ["1"] + d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1} + + [[stencil_set]] + order = 4 + test = 2 + + H.closure = ["-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"] + """ + + parsed_toml = TOML.parse(toml_str) + + @testset "get_stencil_set" begin + @test get_stencil_set(parsed_toml; order = 2) isa Dict + @test get_stencil_set(parsed_toml; order = 2) == parsed_toml["stencil_set"][1] + @test get_stencil_set(parsed_toml; test = 1) == parsed_toml["stencil_set"][2] + @test get_stencil_set(parsed_toml; order = 4, test = 2) == parsed_toml["stencil_set"][3] + + @test_throws ArgumentError get_stencil_set(parsed_toml; test = 2) + @test_throws ArgumentError get_stencil_set(parsed_toml; order = 4) + end + + @testset "parse_stencil" begin + toml = """ + s1 = ["-1/12","4/3","-5/2","4/3","-1/12"] + s2 = {s = ["2", "-5", "4", "-1", "0", "0"], c = 1} + s3 = {s = ["1", "-2", "1", "0", "0", "0"], c = 2} + s4 = "not a stencil" + s5 = [-1, 4, 3] + s6 = {k = ["1", "-2", "1", "0", "0", "0"], c = 2} + s7 = {s = [-1, 4, 3], c = 2} + s8 = {s = ["1", "-2", "1", "0", "0", "0"], c = [2,2]} + """ + + @test parse_stencil(TOML.parse(toml)["s1"]) == CenteredStencil(-1//12, 4//3, -5//2, 4//3, -1//12) + @test parse_stencil(TOML.parse(toml)["s2"]) == Stencil(2//1, -5//1, 4//1, -1//1, 0//1, 0//1; center=1) + @test parse_stencil(TOML.parse(toml)["s3"]) == Stencil(1//1, -2//1, 1//1, 0//1, 0//1, 0//1; center=2) + + @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s4"]) + @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s5"]) + @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s6"]) + @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s7"]) + @test_throws ArgumentError parse_stencil(TOML.parse(toml)["s8"]) + + stencil_set = get_stencil_set(parsed_toml; order = 4, test = 1) + + @test parse_stencil.(stencil_set["D2"]["closure_stencils"]) == [ + Stencil( 2//1, -5//1, 4//1, -1//1, 0//1, 0//1; center=1), + Stencil( 1//1, -2//1, 1//1, 0//1, 0//1, 0//1; center=2), + Stencil(-4//43, 59//43, -110//43, 59//43, -4//43, 0//1; center=3), + Stencil(-1//49, 0//1, 59//49, -118//49, 64//49, -4//49; center=4), + ] + + + @test parse_stencil(Float64, TOML.parse(toml)["s1"]) == CenteredStencil(-1/12, 4/3, -5/2, 4/3, -1/12) + @test parse_stencil(Float64, TOML.parse(toml)["s2"]) == Stencil(2/1, -5/1, 4/1, -1/1, 0/1, 0/1; center=1) + @test parse_stencil(Float64, TOML.parse(toml)["s3"]) == Stencil(1/1, -2/1, 1/1, 0/1, 0/1, 0/1; center=2) + end + + @testset "parse_scalar" begin + toml = TOML.parse(""" + a1 = 1 + a2 = 1.5 + a3 = 1.0 + a4 = 10 + a5 = "1/2" + a6 = "1.5" + + e1 = [1,2,3] + e2 = "a string value" + """) + + @test parse_scalar(toml["a1"]) == 1//1 + @test parse_scalar(toml["a2"]) == 3//2 + @test parse_scalar(toml["a3"]) == 1//1 + @test parse_scalar(toml["a4"]) == 10//1 + @test parse_scalar(toml["a5"]) == 1//2 + @test parse_scalar(toml["a6"]) == 3//2 + + @test_throws ArgumentError parse_scalar(toml["e1"]) + @test_throws ArgumentError parse_scalar(toml["e2"]) + end + + @testset "parse_tuple" begin + toml = TOML.parse(""" + t1 = [1,3,4] + t2 = ["1/2","3/4","2/1"] + + e1 = "not a tuple" + e2.a="1" + e3 = 1 + e4 = ["1/2","3/4","not a number"] + """) + + @test parse_tuple(toml["t1"]) == (1//1,3//1,4//1) + @test parse_tuple(toml["t2"]) == (1//2,3//4,2//1) + + @test_throws ArgumentError parse_tuple(toml["e1"]) + @test_throws ArgumentError parse_tuple(toml["e2"]) + @test_throws ArgumentError parse_tuple(toml["e3"]) + @test_throws ArgumentError parse_tuple(toml["e4"]) + end +end + +@testset "parse_rational" begin + @test SbpOperators.parse_rational("1") isa Rational + @test SbpOperators.parse_rational("1") == 1//1 + @test SbpOperators.parse_rational("1/2") isa Rational + @test SbpOperators.parse_rational("1/2") == 1//2 + @test SbpOperators.parse_rational("37/13") isa Rational + @test SbpOperators.parse_rational("37/13") == 37//13 + + @test SbpOperators.parse_rational(0.5) isa Rational + @test SbpOperators.parse_rational(0.5) == 1//2 + + @test SbpOperators.parse_rational("0.5") isa Rational + @test SbpOperators.parse_rational("0.5") == 1//2 + + @test SbpOperators.parse_rational(2) isa Rational + @test SbpOperators.parse_rational(2) == 2//1 +end
--- a/test/SbpOperators/stencil_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/SbpOperators/stencil_test.jl Thu May 12 21:52:47 2022 +0200 @@ -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/derivatives/first_derivative_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Thu May 12 21:52:47 2022 +0200 @@ -29,12 +29,14 @@ @test first_derivative(g₁, stencil_set, 1) isa LazyTensor{Float64,1,1} @test first_derivative(g₂, stencil_set, 2) isa LazyTensor{Float64,2,2} + @test first_derivative(g₁, stencil_set, 1) == first_derivative(g₁, stencil_set) interior_stencil = CenteredStencil(-1,0,1) closure_stencils = [Stencil(-1,1, center=1)] @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa LazyTensor{Float64,1,1} @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa VolumeOperator + @test first_derivative(g₁, interior_stencil, closure_stencils, 1) == first_derivative(g₁, interior_stencil, closure_stencils) @test first_derivative(g₂, interior_stencil, closure_stencils, 2) isa LazyTensor{Float64,2,2} end @@ -69,6 +71,7 @@ end @testset "Accuracy on function" begin + # 1D g = EquidistantGrid(30, 0.,1.) v = evalOn(g, x->exp(x)) @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)] @@ -77,6 +80,18 @@ @test D₁*v ≈ v rtol=tol end + + # 2D + g = EquidistantGrid((30,60), (0.,0.),(1.,2.)) + v = evalOn(g, (x,y)->exp(0.8x+1.2*y)) + @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)] + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) + Dx = first_derivative(g, stencil_set, 1) + Dy = first_derivative(g, stencil_set, 2) + + @test Dx*v ≈ 0.8v rtol=tol + @test Dy*v ≈ 1.2v rtol=tol + end end end
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Thu May 12 21:52:47 2022 +0200 @@ -21,11 +21,12 @@ Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils,1) @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils) @test Dₓₓ == second_derivative(g_1D,stencil_set,1) + @test Dₓₓ == second_derivative(g_1D,stencil_set) @test Dₓₓ isa VolumeOperator end @testset "2D" begin Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1) - D2 = second_derivative(g_1D,inner_stencil,closure_stencils) + D2 = second_derivative(g_1D,inner_stencil,closure_stencils,1) I = IdentityTensor{Float64}(size(g_2D)[2]) @test Dₓₓ == D2⊗I @test Dₓₓ == second_derivative(g_2D,stencil_set,1) @@ -51,9 +52,7 @@ # implies that L*v should be exact for monomials up to order 2. @testset "2nd order" begin stencil_set = read_stencil_set(operator_path; order=2) - inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) - closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) - Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils) + Dₓₓ = second_derivative(g_1D,stencil_set) @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 @@ -64,9 +63,7 @@ # implies that L*v should be exact for monomials up to order 3. @testset "4th order" begin stencil_set = read_stencil_set(operator_path; order=4) - inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) - closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) - Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils) + Dₓₓ = second_derivative(g_1D,stencil_set) # NOTE: high tolerances for checking the "exact" differentiation # due to accumulation of round-off errors/cancellation errors? @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 @@ -92,9 +89,7 @@ # implies that L*v should be exact for binomials up to order 2. @testset "2nd order" begin stencil_set = read_stencil_set(operator_path; order=2) - inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) - closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) - Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2) + Dyy = second_derivative(g_2D,stencil_set,2) @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 @@ -105,9 +100,7 @@ # implies that L*v should be exact for binomials up to order 3. @testset "4th order" begin stencil_set = read_stencil_set(operator_path; order=4) - inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) - closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) - Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2) + Dyy = second_derivative(g_2D,stencil_set,2) # NOTE: high tolerances for checking the "exact" differentiation # due to accumulation of round-off errors/cancellation errors? @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl Mon Mar 21 09:32:26 2022 +0100 +++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl Thu May 12 21:52:47 2022 +0200 @@ -69,7 +69,7 @@ @testset "laplace" begin @testset "1D" begin Δ = laplace(g_1D, inner_stencil, closure_stencils) - @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils) + @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils, 1) @test Δ isa LazyTensor{T,1,1} where T end @testset "3D" begin