Mercurial > repos > public > sbplib_julia
changeset 634:fb5ac62563aa feature/volume_and_boundary_operators
Integrate feature/quadrature_as_outer_product into branch, before closing feature/quadrature_as_outer_product. (It is now obsolete apart from tests)
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Fri, 01 Jan 2021 16:39:57 +0100 |
| parents | bf8b66c596f7 (diff) a78bda7084f6 (current diff) |
| children | a1dfaf305f41 |
| files | src/SbpOperators/SbpOperators.jl src/SbpOperators/quadrature/diagonal_inner_product.jl src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl src/SbpOperators/quadrature/inverse_quadrature.jl src/SbpOperators/quadrature/quadrature.jl test/testSbpOperators.jl |
| diffstat | 13 files changed, 780 insertions(+), 504 deletions(-) [+] |
line wrap: on
line diff
--- a/TODO.md Fri Jan 01 16:34:55 2021 +0100 +++ b/TODO.md Fri Jan 01 16:39:57 2021 +0100 @@ -5,22 +5,23 @@ - [ ] Skriv tester ## Coding - - [ ] Add new Laplace opertor to DiffOps, probably named WaveEqOp(?!!?) + - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?) - [ ] Add 1D operators (D1, D2, e, d ... ) as TensorOperators - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation. - [ ] Add a quick and simple way of running all tests for all subpackages. - - [ ] Replace getindex hack for flatteing tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`) + - [ ] Replace getindex hack for flattening tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`) - [ ] Use `@inferred` in a lot of tests. - [ ] 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 convetions. For example i,j,... for indecies and I for multi-index + - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index - [ ] Add boundschecking in TensorMappingApplication - [ ] Start renaming things in LazyTensors - [ ] Clean up RegionIndices 1. [ ] Write tests for how things should work 2. [ ] Update RegionIndices accordingly 3. [ ] Fix the rest of the library - - [ ] Add posibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`? + Should getregion also work for getregion(::Colon,...) + - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`? ## Repo - [ ] Add Vidar to the authors list @@ -31,14 +32,3 @@ - [ ] Kolla att vi gör boundschecks överallt och att de är markerade med @boundscheck - [ ] Kolla att vi har @inline på rätt ställen - [ ] Profilera - - -# Old stuff todos (Are these still relevant?) -Borde det finns motsvarande apply_stencil för apply_quadrature, -apply_boundary_value och apply_normal_derivative? - -Borde man alltid skicka in N som parameter i apply_2nd_derivative, t.ex som i -apply_quadrature? - -Just nu agerar apply_normal_derivative, apply_boundary_value på inte på v som -en vektor, utan randvärdet plockas ut utanför. Känns inte konsistent med övrig design.
--- a/src/LazyTensors/lazy_tensor_operations.jl Fri Jan 01 16:34:55 2021 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Fri Jan 01 16:39:57 2021 +0100 @@ -331,7 +331,7 @@ split_tuple((1,2,3,4),Val(3)) -> (1,2,3), (4,) ``` """ -function split_tuple(t::NTuple{N},::Val{M}) where {N,M} +function split_tuple(t::NTuple{N,Any},::Val{M}) where {N,M} return slice_tuple(t,Val(1), Val(M)), slice_tuple(t,Val(M+1), Val(N)) end @@ -341,7 +341,7 @@ Same as `split_tuple(t::NTuple{N},::Val{M})` but splits the tuple in three parts. With the first two parts having lenght `M` and `K`. """ -function split_tuple(t::NTuple{N},::Val{M},::Val{K}) where {N,M,K} +function split_tuple(t::NTuple{N,Any},::Val{M},::Val{K}) where {N,M,K} p1, tail = split_tuple(t, Val(M)) p2, p3 = split_tuple(tail, Val(K)) return p1,p2,p3
--- a/src/SbpOperators/SbpOperators.jl Fri Jan 01 16:34:55 2021 +0100 +++ b/src/SbpOperators/SbpOperators.jl Fri Jan 01 16:39:57 2021 +0100 @@ -8,10 +8,13 @@ include("constantstenciloperator.jl") include("d2.jl") include("readoperator.jl") -include("laplace/secondderivative.jl") -include("laplace/laplace.jl") +include("volumeops/volume_operator.jl") +include("volumeops/derivatives/secondderivative.jl") +include("volumeops/laplace/laplace.jl") include("quadrature/diagonal_quadrature.jl") include("quadrature/inverse_diagonal_quadrature.jl") +include("boundaryops/boundary_operator.jl") include("boundaryops/boundary_restriction.jl") +include("boundaryops/normal_derivative.jl") end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,89 @@ +""" + boundary_operator(grid,closure_stencil,boundary) + +Creates a boundary operator on a `Dim`-dimensional grid for the +specified `boundary`. The action of the operator is determined by `closure_stencil`. + +When `Dim=1`, the corresponding `BoundaryOperator` tensor mapping is returned. +When `Dim>1`, the `BoundaryOperator` `op` is inflated by the outer product +of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`, +the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`. +""" +function boundary_operator(grid::EquidistantGrid{Dim,T}, closure_stencil::Stencil{T}, boundary::CartesianBoundary) where {Dim,T} + #TODO:Check that dim(boundary) <= Dim? + + # Create 1D boundary operator + r = region(boundary) + d = dim(boundary) + op = BoundaryOperator(restrict(grid, d), closure_stencil, r) + + # Create 1D IdentityMappings for each coordinate direction + one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) + Is = IdentityMapping{T}.(size.(one_d_grids)) + + # Formulate the correct outer product sequence of the identity mappings and + # the boundary operator + parts = Base.setindex(Is, op, d) + return foldl(⊗, parts) +end + +""" + BoundaryOperator{T,R,N} <: TensorMapping{T,0,1} + +Implements the boundary operator `op` for 1D as a `TensorMapping` + +`op` is the restriction of a grid function to the boundary using some closure `Stencil{T,N}`. +The boundary to restrict to is determined by `R`. +`op'` is the prolongation of a zero dimensional array to the whole grid using the same closure stencil. +""" +struct BoundaryOperator{T,R<:Region,N} <: TensorMapping{T,0,1} + stencil::Stencil{T,N} + size::Int +end + +BoundaryOperator{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryOperator{T,R,N}(stencil, size) + +""" + BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region) + +Constructs the BoundaryOperator with stencil `closure_stencil` for a one-dimensional `grid`, restricting to +to the boundary specified by `region`. +""" +function BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil::Stencil{T,N}, region::Region) where {T,N} + return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1]) +end + +""" + closure_size(::BoundaryOperator) +The size of the closure stencil. +""" +closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N + +LazyTensors.range_size(op::BoundaryOperator) = () +LazyTensors.domain_size(op::BoundaryOperator) = (op.size,) + +function LazyTensors.apply(op::BoundaryOperator{T,Lower}, v::AbstractVector{T}) where T + apply_stencil(op.stencil,v,1) +end + +function LazyTensors.apply(op::BoundaryOperator{T,Upper}, v::AbstractVector{T}) where T + apply_stencil_backwards(op.stencil,v,op.size) +end + +function LazyTensors.apply_transpose(op::BoundaryOperator{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T + return op.stencil[Int(i)-1]*v[] +end + +function LazyTensors.apply_transpose(op::BoundaryOperator{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T + return op.stencil[op.size[1] - Int(i)]*v[] +end + +# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods. +function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i::Index) where T + return zero(T) +end + +function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T + r = getregion(i, closure_size(op), op.size) + apply_transpose(op, v, Index(i,r)) +end
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Fri Jan 01 16:34:55 2021 +0100 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Fri Jan 01 16:39:57 2021 +0100 @@ -1,81 +1,15 @@ -""" - boundary_restriction(grid,closureStencil,boundary) - -Creates a boundary restriction operator on a `Dim`-dimensional grid for the -specified `boundary`. - -When `Dim=1`, the corresponding `BoundaryRestriction` tensor mapping is returned. -When `Dim>1`, the `BoundaryRestriction` `e` is inflated by the outer product -of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`, -the boundary restriction operator in the y-direction direction is `Ix⊗e⊗Iz`. -""" -function boundary_restriction(grid::EquidistantGrid{Dim,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary) where {Dim,T,M} - # Create 1D boundary restriction operator - r = region(boundary) - d = dim(boundary) - e = BoundaryRestriction(restrict(grid, d), closureStencil, r) - - # Create 1D IdentityMappings for each coordinate direction - one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) - Is = IdentityMapping{T}.(size.(one_d_grids)) - - # Formulate the correct outer product sequence of the identity mappings and - # the boundary restriction operator - parts = Base.setindex(Is, e, d) - return foldl(⊗, parts) -end - -export boundary_restriction - -""" - BoundaryRestriction{T,R,N} <: TensorMapping{T,0,1} - -Implements the boundary operator `e` for 1D as a `TensorMapping` - -`e` is the restriction of a grid function to the boundary using some `closureStencil`. -The boundary to restrict to is determined by `R`. - -`e'` is the prolongation of a zero dimensional array to the whole grid using the same `closureStencil`. """ -struct BoundaryRestriction{T,R<:Region,N} <: TensorMapping{T,0,1} - stencil::Stencil{T,N} - size::Int -end -export BoundaryRestriction - -BoundaryRestriction{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryRestriction{T,R,N}(stencil, size) + BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) + BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) -function BoundaryRestriction(grid::EquidistantGrid{1}, closureStencil::Stencil{T,N}, region::Region) where {T,N} - return BoundaryRestriction{T,typeof(region),N}(closureStencil,size(grid)[1]) -end - -closure_size(::BoundaryRestriction{T,R,N}) where {T,R,N} = N - -LazyTensors.range_size(e::BoundaryRestriction) = () -LazyTensors.domain_size(e::BoundaryRestriction) = (e.size,) - -function LazyTensors.apply(e::BoundaryRestriction{T,Lower}, v::AbstractVector{T}) where T - apply_stencil(e.stencil,v,1) -end +Creates the boundary restriction operator `e` as a `TensorMapping` -function LazyTensors.apply(e::BoundaryRestriction{T,Upper}, v::AbstractVector{T}) where T - apply_stencil_backwards(e.stencil,v,e.size) -end - -function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T - return e.stencil[Int(i)-1]*v[] -end +`e` is the restriction of a grid function to the boundary specified by `boundary` or `region` using some `closure_stencil`. +`e'` is the prolongation of a grid function on the boundary to the whole grid using the same `closure_stencil`. +On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of +a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details. +""" +BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary) +BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = BoundaryRestriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}()) -function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T - return e.stencil[e.size[1] - Int(i)]*v[] -end - -# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods. -function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i::Index) where T - return zero(T) -end - -function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i) where T - r = getregion(i, closure_size(e), e.size) - apply_transpose(e, v, Index(i,r)) -end +export BoundaryRestriction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,18 @@ +""" + NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) + NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) + +Creates the normal derivative boundary operator `d` as a `TensorMapping` + +`d` is the normal derivative of a grid function at the boundary specified by `boundary` or `region` using some `closure_stencil`. +`d'` is the prolongation of the normal derivative of a grid function to the whole grid using the same `closure_stencil`. +On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of +a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details. +""" +function NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) + direction = dim(boundary) + h_inv = inverse_spacing(grid)[direction] + return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary) +end +NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = NormalDerivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}()) +export NormalDerivative
--- a/src/SbpOperators/laplace/laplace.jl Fri Jan 01 16:34:55 2021 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,145 +0,0 @@ -export Laplace -""" - Laplace{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the Laplace operator `L` in Dim dimensions as a tensor operator -The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative -tensor operators. -""" -#export quadrature, inverse_quadrature, boundary_quadrature, boundary_value, normal_derivative -struct Laplace{Dim,T,N,M,K} <: TensorMapping{T,Dim,Dim} - D2::NTuple{Dim,SecondDerivative{T,N,M,K}} -end - -function Laplace(g::EquidistantGrid{Dim}, innerStencil, closureStencils) where Dim - D2 = () - for i ∈ 1:Dim - D2 = (D2..., SecondDerivative(restrict(g,i), innerStencil, closureStencils)) - end - - return Laplace(D2) -end - -LazyTensors.range_size(L::Laplace) = getindex.(range_size.(L.D2),1) -LazyTensors.domain_size(L::Laplace) = getindex.(domain_size.(L.D2),1) - -function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim} - error("not implemented") -end - -# u = L*v -function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, i) where T - @inbounds u = LazyTensors.apply(L.D2[1],v,i) - return u -end - -function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, i, j) where T - # 2nd x-derivative - @inbounds vx = view(v, :, Int(j)) - @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , i) - - # 2nd y-derivative - @inbounds vy = view(v, Int(i), :) - @inbounds uᵢ += LazyTensors.apply(L.D2[2], vy , j) - - return uᵢ -end - -# quadrature(L::Laplace) = Quadrature(L.op, L.grid) -# inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid) -# boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId) -# normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId) -# boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId) -# export NormalDerivative -# """ -# NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} -# -# Implements the boundary operator `d` as a TensorMapping -# """ -# struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} -# op::D2{T,N,M,K} -# grid::EquidistantGrid{2} -# bId::CartesianBoundary -# end -# -# # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue? -# # Can we give special treatment to TensorMappings that go to a higher dim? -# function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer}) -# if dim(e.bId) == 1 -# return (UnknownDim, domain_size[1]) -# elseif dim(e.bId) == 2 -# return (domain_size[1], UnknownDim) -# end -# end -# LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],) -# -# # TODO: Not type stable D:< -# # TODO: Make this independent of dimension -# function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T -# i = I[dim(d.bId)] -# j = I[3-dim(d.bId)] -# N_i = size(d.grid)[dim(d.bId)] -# h_inv = inverse_spacing(d.grid)[dim(d.bId)] -# return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId)) -# end -# -# function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T -# u = selectdim(v,3-dim(d.bId),Int(I[1])) -# return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId)) -# end -# -# """ -# BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} -# -# Implements the boundary operator `q` as a TensorOperator -# """ -# export BoundaryQuadrature -# struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} -# op::D2{T,N,M,K} -# grid::EquidistantGrid{2} -# bId::CartesianBoundary -# end -# -# -# # TODO: Make this independent of dimension -# function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T -# h = spacing(q.grid)[3-dim(q.bId)] -# N = size(v) -# return apply_quadrature(q.op, h, v[I[1]], I[1], N[1]) -# end -# -# LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I) -# -# -# -# -# struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end -# -# function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} -# e = boundary_value(L, Bid()) -# d = normal_derivative(L, Bid()) -# Hᵧ = boundary_quadrature(L, Bid()) -# H⁻¹ = inverse_quadrature(L) -# return (-H⁻¹*e*Hᵧ*(d'*v - g))[I] -# end -# -# struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition -# tau::Float64 -# end -# -# function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} -# e = boundary_value(L, Bid()) -# d = normal_derivative(L, Bid()) -# Hᵧ = boundary_quadrature(L, Bid()) -# H⁻¹ = inverse_quadrature(L) -# return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I] -# # Need to handle scalar multiplication and addition of TensorMapping -# end - -# function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D - # return apply(s.L, v, i) + -# sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + -# sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + -# sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + -# sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) -# end
--- a/src/SbpOperators/laplace/secondderivative.jl Fri Jan 01 16:34:55 2021 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,43 +0,0 @@ -""" - SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1} -Implements the Laplace tensor operator `L` with constant grid spacing and coefficients -in 1D dimension -""" - -struct SecondDerivative{T,N,M,K} <: TensorMapping{T,1,1} - h_inv::T # The grid spacing could be included in the stencil already. Preferable? - innerStencil::Stencil{T,N} - closureStencils::NTuple{M,Stencil{T,K}} - size::NTuple{1,Int} -end -export SecondDerivative - -function SecondDerivative(grid::EquidistantGrid{1}, innerStencil, closureStencils) - h_inv = inverse_spacing(grid)[1] - return SecondDerivative(h_inv, innerStencil, closureStencils, size(grid)) -end - -LazyTensors.range_size(D2::SecondDerivative) = D2.size -LazyTensors.domain_size(D2::SecondDerivative) = D2.size - -# Apply for different regions Lower/Interior/Upper or Unknown region -function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Lower}) where T - return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i)) -end - -function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Interior}) where T - return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(i)) -end - -function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Upper}) where T - N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v)) - return @inbounds D2.h_inv*D2.h_inv*apply_stencil_backwards(D2.closureStencils[N-Int(i)+1], v, Int(i)) -end - -function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i) where T - N = length(v) # TODO: Use domain_size here instead? - r = getregion(i, closuresize(D2), N) - return LazyTensors.apply(D2, v, Index(i, r)) -end - -closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} = M
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/derivatives/secondderivative.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,20 @@ +""" + SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) + SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) + +Creates the second-derivative operator `D2` as a `TensorMapping` + +`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by +`direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` +for the points in the closure regions. + +On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the outer product of the +one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections. +Also see the documentation of `SbpOperators.volume_operator(...)` for more details. +""" +function SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim + 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 +SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = SecondDerivative(grid,inner_stencil,closure_stencils,1) +export SecondDerivative
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,82 @@ +""" + Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) + +Creates the Laplace ooperator operator `Δ` as a `TensorMapping` + +`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using +the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` +for the points in the closure regions. + +On a one-dimensional `grid`, `Δ` is a `SecondDerivative`. On a multi-dimensional `grid`, `Δ` is the sum of +multi-dimensional `SecondDerivative`s where the sum is carried out lazily. +""" +function Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim + Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1) + for d = 2:Dim + Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d) + end + return Δ +end +export Laplace + +# quadrature(L::Laplace) = Quadrature(L.op, L.grid) +# inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid) +# boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId) +# normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId) +# boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId) + +# """ +# BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} +# +# Implements the boundary operator `q` as a TensorOperator +# """ +# export BoundaryQuadrature +# struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} +# op::D2{T,N,M,K} +# grid::EquidistantGrid{2} +# bId::CartesianBoundary +# end +# +# +# # TODO: Make this independent of dimension +# function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T +# h = spacing(q.grid)[3-dim(q.bId)] +# N = size(v) +# return apply_quadrature(q.op, h, v[I[1]], I[1], N[1]) +# end +# +# LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I) +# +# +# +# +# struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end +# +# function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} +# e = boundary_value(L, Bid()) +# d = normal_derivative(L, Bid()) +# Hᵧ = boundary_quadrature(L, Bid()) +# H⁻¹ = inverse_quadrature(L) +# return (-H⁻¹*e*Hᵧ*(d'*v - g))[I] +# end +# +# struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition +# tau::Float64 +# end +# +# function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} +# e = boundary_value(L, Bid()) +# d = normal_derivative(L, Bid()) +# Hᵧ = boundary_quadrature(L, Bid()) +# H⁻¹ = inverse_quadrature(L) +# return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I] +# # Need to handle scalar multiplication and addition of TensorMapping +# end + +# function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D + # return apply(s.L, v, i) + +# sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + +# sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + +# sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + +# sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) +# end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/volume_operator.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,60 @@ +""" + volume_operator(grid,inner_stencil,closure_stencils,parity,direction) +Creates a volume operator on a `Dim`-dimensional grid acting along the +specified coordinate `direction`. The action of the operator is determined by the +stencils `inner_stencil` and `closure_stencils`. +When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned. +When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product +of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`, +the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`. +""" +function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M} + #TODO: Check that direction <= Dim? + + # Create 1D volume operator in along coordinate direction + op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) + # Create 1D IdentityMappings for each coordinate direction + one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) + Is = IdentityMapping{T}.(size.(one_d_grids)) + # Formulate the correct outer product sequence of the identity mappings and + # the volume operator + parts = Base.setindex(Is, op, direction) + return foldl(⊗, parts) +end + +""" + VolumeOperator{T,N,M,K} <: TensorOperator{T,1} +Implements a one-dimensional constant coefficients volume operator +""" +struct VolumeOperator{T,N,M,K} <: TensorMapping{T,1,1} + inner_stencil::Stencil{T,N} + closure_stencils::NTuple{M,Stencil{T,K}} + size::NTuple{1,Int} + parity::Parity +end + +function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity) + return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity) +end + +closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M + +LazyTensors.range_size(op::VolumeOperator) = op.size +LazyTensors.domain_size(op::VolumeOperator) = op.size + +function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T + return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i)) +end + +function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T + return apply_stencil(op.inner_stencil, v, Int(i)) +end + +function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T + return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i)) +end + +function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T + r = getregion(i, closure_size(op), op.size[1]) + return LazyTensors.apply(op, v, Index(i, r)) +end
--- a/test/testLazyTensors.jl Fri Jan 01 16:34:55 2021 +0100 +++ b/test/testLazyTensors.jl Fri Jan 01 16:39:57 2021 +0100 @@ -487,17 +487,22 @@ @test LazyTensors.split_tuple((1,2,3,4),Val(3)) == ((1,2,3),(4,)) @test LazyTensors.split_tuple((1,2,3,4),Val(4)) == ((1,2,3,4),()) + @test LazyTensors.split_tuple((1,2,true,4),Val(3)) == ((1,2,true),(4,)) + @inferred LazyTensors.split_tuple((1,2,3,4),Val(3)) + @inferred LazyTensors.split_tuple((1,2,true,4),Val(3)) end @testset "3 parts" begin @test LazyTensors.split_tuple((),Val(0),Val(0)) == ((),(),()) @test LazyTensors.split_tuple((1,2,3),Val(1), Val(1)) == ((1,),(2,),(3,)) + @test LazyTensors.split_tuple((1,true,3),Val(1), Val(1)) == ((1,),(true,),(3,)) @test LazyTensors.split_tuple((1,2,3,4,5,6),Val(1),Val(2)) == ((1,),(2,3),(4,5,6)) @test LazyTensors.split_tuple((1,2,3,4,5,6),Val(3),Val(2)) == ((1,2,3),(4,5),(6,)) @inferred LazyTensors.split_tuple((1,2,3,4,5,6),Val(3),Val(2)) + @inferred LazyTensors.split_tuple((1,true,3),Val(1), Val(1)) end end
--- a/test/testSbpOperators.jl Fri Jan 01 16:34:55 2021 +0100 +++ b/test/testSbpOperators.jl Fri Jan 01 16:39:57 2021 +0100 @@ -7,6 +7,12 @@ using TOML import Sbplib.SbpOperators.Stencil +import Sbplib.SbpOperators.VolumeOperator +import Sbplib.SbpOperators.volume_operator +import Sbplib.SbpOperators.BoundaryOperator +import Sbplib.SbpOperators.boundary_operator +import Sbplib.SbpOperators.Parity + @testset "SbpOperators" begin @@ -129,86 +135,274 @@ # end # end -@testset "SecondDerivative" begin - op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) - L = 3.5 - g = EquidistantGrid(101, 0.0, L) - Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils) +@testset "VolumeOperator" begin + inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2) + closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2)) + g_1D = EquidistantGrid(11,0.,1.) + g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.)) + g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.)) + @testset "Constructors" begin + #TODO: How are even and odd in SbpOperators.Parity exposed? Currently constructing even as Parity(1) + @testset "1D" begin + op = VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1)) + @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,Parity(1)) + @test op == volume_operator(g_1D,inner_stencil,closure_stencils,Parity(1),1) + @test op isa TensorMapping{T,1,1} where T + end + @testset "2D" begin + op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1) + op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2) + Ix = IdentityMapping{Float64}((11,)) + Iy = IdentityMapping{Float64}((12,)) + @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))⊗Iy + @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),Parity(1)) + @test op_x isa TensorMapping{T,2,2} where T + @test op_y isa TensorMapping{T,2,2} where T + end + @testset "3D" begin + op_x = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),1) + op_y = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),2) + op_z = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),3) + Ix = IdentityMapping{Float64}((11,)) + Iy = IdentityMapping{Float64}((12,)) + Iz = IdentityMapping{Float64}((10,)) + @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))⊗Iy⊗Iz + @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),Parity(1))⊗Iz + @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),Parity(1)) + @test op_x isa TensorMapping{T,3,3} where T + @test op_y isa TensorMapping{T,3,3} where T + @test op_z isa TensorMapping{T,3,3} where T + end + end - f0(x) = 1. - f1(x) = x - f2(x) = 1/2*x^2 - f3(x) = 1/6*x^3 - f4(x) = 1/24*x^4 - f5(x) = sin(x) - f5ₓₓ(x) = -f5(x) + @testset "Sizes" begin + @testset "1D" begin + op = volume_operator(g_1D,inner_stencil,closure_stencils,Parity(1),1) + @test range_size(op) == domain_size(op) == size(g_1D) + end - v0 = evalOn(g,f0) - v1 = evalOn(g,f1) - v2 = evalOn(g,f2) - v3 = evalOn(g,f3) - v4 = evalOn(g,f4) - v5 = evalOn(g,f5) + @testset "2D" begin + op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1) + op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2) + @test range_size(op_y) == domain_size(op_y) == + range_size(op_x) == domain_size(op_x) == size(g_2D) + end + @testset "3D" begin + op_x = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),1) + op_y = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),2) + op_z = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),3) + @test range_size(op_z) == domain_size(op_z) == + range_size(op_y) == domain_size(op_y) == + range_size(op_x) == domain_size(op_x) == size(g_3D) + end + end - @test Dₓₓ isa TensorMapping{T,1,1} where T - @test Dₓₓ' isa TensorMapping{T,1,1} where T + # TODO: Test for other dimensions? + op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1) + op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(-1),2) + v = zeros(size(g_2D)) + Nx = size(g_2D)[1] + Ny = size(g_2D)[2] + for i = 1:Nx + v[i,:] .= i + end + rx = copy(v) + rx[1,:] .= 1.5 + rx[Nx,:] .= (2*Nx-1)/2 + ry = copy(v) + ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny] - # 4th order interior stencil, 2nd order boundary stencil, - # implies that L*v should be exact for v - monomial up to order 3. - # Exact differentiation is measured point-wise. For other grid functions - # the error is measured in the l2-norm. - @test norm(Dₓₓ*v0) ≈ 0.0 atol=5e-10 - @test norm(Dₓₓ*v1) ≈ 0.0 atol=5e-10 - @test Dₓₓ*v2 ≈ v0 atol=5e-11 - @test Dₓₓ*v3 ≈ v1 atol=5e-11 + @testset "Application" begin + @test op_x*v ≈ rx rtol = 1e-14 + @test op_y*v ≈ ry rtol = 1e-14 + end + + # TODO: Test for other dimensions? + @testset "Regions" begin + @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14 + @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14 + @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14 + @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14 + @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14 + + @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)] + @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)] - h = spacing(g)[1]; - l2(v) = sqrt(h*sum(v.^2)) - @test Dₓₓ*v4 ≈ v2 atol=5e-4 norm=l2 - @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2 + @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14 + @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14 + @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14 + @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14 + @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14 + + @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)] + @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)] + end + + # TODO: Test for other dimensions? + @testset "Inferred" begin + @inferred apply(op_x, v,1,1) + @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower)) + @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower)) + @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower)) + + @inferred apply(op_y, v,1,1) + @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower)) + @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior)) + @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper)) + end + end - -@testset "Laplace2D" begin +@testset "SecondDerivative" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) - Lx = 1.5 - Ly = 3.2 - g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly)) - L = Laplace(g, op.innerStencil, op.closureStencils) + Lx = 3.5 + Ly = 3. + g_1D = EquidistantGrid(121, 0.0, Lx) + g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly)) + + # TODO: These areant really constructors. Better name? + @testset "Constructors" begin + @testset "1D" begin + Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) + @test Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1) + @test Dₓₓ isa VolumeOperator + end + @testset "2D" begin + Dₓₓ = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,1) + D2 = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) + I = IdentityMapping{Float64}(size(g_2D)[2]) + @test Dₓₓ == D2⊗I + @test Dₓₓ isa TensorMapping{T,2,2} where T + end + end + + @testset "Accuracy" begin + @testset "1D" begin + monomials = () + maxOrder = 4; + for i = 0:maxOrder-1 + f_i(x) = 1/factorial(i)*x^i + monomials = (monomials...,evalOn(g_1D,f_i)) + end + l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2)); + #TODO: Error when reading second order stencil! + # @testset "2nd order" begin + # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) + # Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) + # @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-13 + # @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-13 + # end - f0(x,y) = 2. - f1(x,y) = x+y - f2(x,y) = 1/2*x^2 + 1/2*y^2 - f3(x,y) = 1/6*x^3 + 1/6*y^3 - f4(x,y) = 1/24*x^4 + 1/24*y^4 - f5(x,y) = sin(x) + cos(y) - f5ₓₓ(x,y) = -f5(x,y) + # 4th order interior stencil, 2nd order boundary stencil, + # implies that L*v should be exact for monomials up to order 3. + # Exact differentiation is measured point-wise. For other grid functions + # the error is measured in the l2-norm. + @testset "4th order" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) + # TODO: 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 + @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 + @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 + @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10 + @test Dₓₓ*evalOn(g_1D,x -> sin(x)) ≈ evalOn(g_1D,x -> -sin(x)) rtol = 5e-4 norm = l2 + end + end + + @testset "2D" begin + binomials = () + maxOrder = 4; + for i = 0:maxOrder-1 + f_i(x,y) = 1/factorial(i)*y^i + x^i + binomials = (binomials...,evalOn(g_2D,f_i)) + end + l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2)); + + #TODO: Error when reading second order stencil! + # @testset "2nd order" begin + # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) + # Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) + # @test Dyy*binomials[1] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 + # @test Dyy*binomials[2] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 + # end - v0 = evalOn(g,f0) - v1 = evalOn(g,f1) - v2 = evalOn(g,f2) - v3 = evalOn(g,f3) - v4 = evalOn(g,f4) - v5 = evalOn(g,f5) - v5ₓₓ = evalOn(g,f5ₓₓ) - - @test L isa TensorMapping{T,2,2} where T - @test L' isa TensorMapping{T,2,2} where T + # 4th order interior stencil, 2nd order boundary stencil, + # implies that L*v should be exact for binomials up to order 3. + # Exact differentiation is measured point-wise. For other grid functions + # the error is measured in the l2-norm. + @testset "4th order" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) + # TODO: 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 + @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 + @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9 + @test Dyy*evalOn(g_2D, (x,y) -> sin(x)+cos(y)) ≈ evalOn(g_2D,(x,y) -> -cos(y)) rtol = 5e-4 norm = l2 + end + end + end +end - # 4th order interior stencil, 2nd order boundary stencil, - # implies that L*v should be exact for v - monomial up to order 3. - # Exact differentiation is measured point-wise. For other grid functions - # the error is measured in the H-norm. - @test norm(L*v0) ≈ 0 atol=5e-10 - @test norm(L*v1) ≈ 0 atol=5e-10 - @test L*v2 ≈ v0 # Seems to be more accurate - @test L*v3 ≈ v1 atol=5e-10 +@testset "Laplace" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + g_1D = EquidistantGrid(101, 0.0, 1.) + #TODO: It's nice to verify that 3D works somewhere at least, but perhaps should keep 3D tests to a minimum to avoid + # long run time for test? + g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) + # TODO: These areant really constructors. Better name? + @testset "Constructors" begin + @testset "1D" begin + L = Laplace(g_1D, op.innerStencil, op.closureStencils) + @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils) + @test L isa TensorMapping{T,1,1} where T + end + @testset "3D" begin + L = Laplace(g_3D, op.innerStencil, op.closureStencils) + @test L isa TensorMapping{T,3,3} where T + Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1) + Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2) + Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3) + @test L == Dxx + Dyy + Dzz + end + end - h = spacing(g) - l2(v) = sqrt(prod(h)*sum(v.^2)) - @test L*v4 ≈ v2 atol=5e-4 norm=l2 - @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2 + @testset "Accuracy" begin + polynomials = () + maxOrder = 4; + for i = 0:maxOrder-1 + f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i) + polynomials = (polynomials...,evalOn(g_3D,f_i)) + end + l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2)); + + #TODO: Error when reading second order stencil! + # @testset "2nd order" begin + # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) + # Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) + # @test Dyy*binomials[1] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 + # @test Dyy*binomials[2] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 + # end + + # 4th order interior stencil, 2nd order boundary stencil, + # implies that L*v should be exact for binomials up to order 3. + # Exact differentiation is measured point-wise. For other grid functions + # the error is measured in the l2-norm. + @testset "4th order" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + L = Laplace(g_3D,op.innerStencil,op.closureStencils) + # TODO: high tolerances for checking the "exact" differentiation + # due to accumulation of round-off errors/cancellation errors? + @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 + @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 + @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 + @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9 + @test L*evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z)) ≈ evalOn(g_3D,(x,y,z) -> -sin(x)-cos(y) + exp(z)) rtol = 5e-4 norm = l2 + end + end end @testset "DiagonalQuadrature" begin @@ -359,66 +553,186 @@ end end -@testset "BoundaryRestrictrion" begin - op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) +@testset "BoundaryOperator" begin + closure_stencil = Stencil((0,2), (2.,1.,3.)) g_1D = EquidistantGrid(11, 0.0, 1.0) g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) @testset "Constructors" begin @testset "1D" begin - e_l = BoundaryRestriction{Lower}(op.eClosure,size(g_1D)[1]) - @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower()) - @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}()) - @test e_l isa TensorMapping{T,0,1} where T + op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1]) + @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower()) + @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}()) + @test op_l isa TensorMapping{T,0,1} where T - e_r = BoundaryRestriction{Upper}(op.eClosure,size(g_1D)[1]) - @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper()) - @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}()) - @test e_r isa TensorMapping{T,0,1} where T + op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1]) + @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper()) + @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}()) + @test op_r isa TensorMapping{T,0,1} where T end @testset "2D" begin - e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}()) + e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}()) @test e_w isa InflatedTensorMapping @test e_w isa TensorMapping{T,1,2} where T end end - e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}()) - e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}()) + op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}()) + op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}()) - e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}()) - e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}()) - e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}()) - e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}()) + op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}()) + op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}()) + op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}()) + op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}()) @testset "Sizes" begin @testset "1D" begin - @test domain_size(e_l) == (11,) - @test domain_size(e_r) == (11,) + @test domain_size(op_l) == (11,) + @test domain_size(op_r) == (11,) - @test range_size(e_l) == () - @test range_size(e_r) == () + @test range_size(op_l) == () + @test range_size(op_r) == () end @testset "2D" begin - @test domain_size(e_w) == (11,15) - @test domain_size(e_e) == (11,15) - @test domain_size(e_s) == (11,15) - @test domain_size(e_n) == (11,15) + @test domain_size(op_w) == (11,15) + @test domain_size(op_e) == (11,15) + @test domain_size(op_s) == (11,15) + @test domain_size(op_n) == (11,15) - @test range_size(e_w) == (15,) - @test range_size(e_e) == (15,) - @test range_size(e_s) == (11,) - @test range_size(e_n) == (11,) + @test range_size(op_w) == (15,) + @test range_size(op_e) == (15,) + @test range_size(op_s) == (11,) + @test range_size(op_n) == (11,) end end - @testset "Application" begin @testset "1D" begin v = evalOn(g_1D,x->1+x^2) u = fill(3.124) + @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3] + @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2] + @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] + @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] + @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] + end + + @testset "2D" begin + v = rand(size(g_2D)...) + u = fill(3.124) + @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14 + @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14 + @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14 + @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14 + + + g_x = rand(size(g_2D)[1]) + g_y = rand(size(g_2D)[2]) + + G_w = zeros(Float64, size(g_2D)...) + G_w[1,:] = 2*g_y + G_w[2,:] = g_y + G_w[3,:] = 3*g_y + + G_e = zeros(Float64, size(g_2D)...) + G_e[end,:] = 2*g_y + G_e[end-1,:] = g_y + G_e[end-2,:] = 3*g_y + + G_s = zeros(Float64, size(g_2D)...) + G_s[:,1] = 2*g_x + G_s[:,2] = g_x + G_s[:,3] = 3*g_x + + G_n = zeros(Float64, size(g_2D)...) + G_n[:,end] = 2*g_x + G_n[:,end-1] = g_x + G_n[:,end-2] = 3*g_x + + @test op_w'*g_y == G_w + @test op_e'*g_y == G_e + @test op_s'*g_x == G_s + @test op_n'*g_x == G_n + end + + @testset "Regions" begin + u = fill(3.124) + @test (op_l'*u)[Index(1,Lower)] == 2*u[] + @test (op_l'*u)[Index(2,Lower)] == u[] + @test (op_l'*u)[Index(6,Interior)] == 0 + @test (op_l'*u)[Index(10,Upper)] == 0 + @test (op_l'*u)[Index(11,Upper)] == 0 + + @test (op_r'*u)[Index(1,Lower)] == 0 + @test (op_r'*u)[Index(2,Lower)] == 0 + @test (op_r'*u)[Index(6,Interior)] == 0 + @test (op_r'*u)[Index(10,Upper)] == u[] + @test (op_r'*u)[Index(11,Upper)] == 2*u[] + end + end + + @testset "Inferred" begin + v = ones(Float64, 11) + u = fill(1.) + + @inferred apply(op_l, v) + @inferred apply(op_r, v) + + @inferred apply_transpose(op_l, u, 4) + @inferred apply_transpose(op_l, u, Index(1,Lower)) + @inferred apply_transpose(op_l, u, Index(2,Lower)) + @inferred apply_transpose(op_l, u, Index(6,Interior)) + @inferred apply_transpose(op_l, u, Index(10,Upper)) + @inferred apply_transpose(op_l, u, Index(11,Upper)) + + @inferred apply_transpose(op_r, u, 4) + @inferred apply_transpose(op_r, u, Index(1,Lower)) + @inferred apply_transpose(op_r, u, Index(2,Lower)) + @inferred apply_transpose(op_r, u, Index(6,Interior)) + @inferred apply_transpose(op_r, u, Index(10,Upper)) + @inferred apply_transpose(op_r, u, Index(11,Upper)) + end + +end + +@testset "BoundaryRestriction" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + g_1D = EquidistantGrid(11, 0.0, 1.0) + g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) + + # TODO: These areant really constructors. Better name? + @testset "Constructors" begin + @testset "1D" begin + e_l = BoundaryRestriction(g_1D,op.eClosure,Lower()) + @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}()) + @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower()) + @test e_l isa BoundaryOperator{T,Lower} where T + @test e_l isa TensorMapping{T,0,1} where T + + e_r = BoundaryRestriction(g_1D,op.eClosure,Upper()) + @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}()) + @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper()) + @test e_r isa BoundaryOperator{T,Upper} where T + @test e_r isa TensorMapping{T,0,1} where T + end + + @testset "2D" begin + e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}()) + @test e_w isa InflatedTensorMapping + @test e_w isa TensorMapping{T,1,2} where T + end + end + + @testset "Application" begin + @testset "1D" begin + e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}()) + e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}()) + + v = evalOn(g_1D,x->1+x^2) + u = fill(3.124) + @test (e_l*v)[] == v[1] @test (e_r*v)[] == v[end] @test (e_r*v)[1] == v[end] @@ -427,6 +741,11 @@ end @testset "2D" begin + e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}()) + e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}()) + e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}()) + e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}()) + v = rand(11, 15) u = fill(3.124) @@ -436,144 +755,88 @@ @test e_n*v == v[:,end] - g_x = rand(11) - g_y = rand(15) - - G_w = zeros(Float64, (11,15)) - G_w[1,:] = g_y + g_x = rand(11) + g_y = rand(15) - G_e = zeros(Float64, (11,15)) - G_e[end,:] = g_y + G_w = zeros(Float64, (11,15)) + G_w[1,:] = g_y - G_s = zeros(Float64, (11,15)) - G_s[:,1] = g_x - - G_n = zeros(Float64, (11,15)) - G_n[:,end] = g_x + G_e = zeros(Float64, (11,15)) + G_e[end,:] = g_y - @test e_w'*g_y == G_w - @test e_e'*g_y == G_e - @test e_s'*g_x == G_s - @test e_n'*g_x == G_n - end + G_s = zeros(Float64, (11,15)) + G_s[:,1] = g_x - @testset "Regions" begin - u = fill(3.124) - @test (e_l'*u)[Index(1,Lower)] == 3.124 - @test (e_l'*u)[Index(2,Lower)] == 0 - @test (e_l'*u)[Index(6,Interior)] == 0 - @test (e_l'*u)[Index(10,Upper)] == 0 - @test (e_l'*u)[Index(11,Upper)] == 0 + G_n = zeros(Float64, (11,15)) + G_n[:,end] = g_x - @test (e_r'*u)[Index(1,Lower)] == 0 - @test (e_r'*u)[Index(2,Lower)] == 0 - @test (e_r'*u)[Index(6,Interior)] == 0 - @test (e_r'*u)[Index(10,Upper)] == 0 - @test (e_r'*u)[Index(11,Upper)] == 3.124 + @test e_w'*g_y == G_w + @test e_e'*g_y == G_e + @test e_s'*g_x == G_s + @test e_n'*g_x == G_n end end +end - @testset "Inferred" begin - v = ones(Float64, 11) - u = fill(1.) +@testset "NormalDerivative" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0)) - @inferred apply(e_l, v) - @inferred apply(e_r, v) + d_w = NormalDerivative(g, op.dClosure, CartesianBoundary{1,Lower}()) + d_e = NormalDerivative(g, op.dClosure, CartesianBoundary{1,Upper}()) + d_s = NormalDerivative(g, op.dClosure, CartesianBoundary{2,Lower}()) + d_n = NormalDerivative(g, op.dClosure, CartesianBoundary{2,Upper}()) + - @inferred apply_transpose(e_l, u, 4) - @inferred apply_transpose(e_l, u, Index(1,Lower)) - @inferred apply_transpose(e_l, u, Index(2,Lower)) - @inferred apply_transpose(e_l, u, Index(6,Interior)) - @inferred apply_transpose(e_l, u, Index(10,Upper)) - @inferred apply_transpose(e_l, u, Index(11,Upper)) + v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) + v∂x = evalOn(g, (x,y)-> 2*x + y) + v∂y = evalOn(g, (x,y)-> 2*(y-1) + x) + + @test d_w isa TensorMapping{T,1,2} where T - @inferred apply_transpose(e_r, u, 4) - @inferred apply_transpose(e_r, u, Index(1,Lower)) - @inferred apply_transpose(e_r, u, Index(2,Lower)) - @inferred apply_transpose(e_r, u, Index(6,Interior)) - @inferred apply_transpose(e_r, u, Index(10,Upper)) - @inferred apply_transpose(e_r, u, Index(11,Upper)) + @test d_w*v ≈ v∂x[1,:] + @test d_e*v ≈ -v∂x[end,:] + @test d_s*v ≈ v∂y[:,1] + @test d_n*v ≈ -v∂y[:,end] + + + d_x_l = zeros(Float64, size(g)[1]) + d_x_u = zeros(Float64, size(g)[1]) + for i ∈ eachindex(d_x_l) + d_x_l[i] = op.dClosure[i-1] + d_x_u[i] = op.dClosure[length(d_x_u)-i] end + d_y_l = zeros(Float64, size(g)[2]) + d_y_u = zeros(Float64, size(g)[2]) + for i ∈ eachindex(d_y_l) + d_y_l[i] = op.dClosure[i-1] + d_y_u[i] = op.dClosure[length(d_y_u)-i] + end + + function prod_matrix(x,y) + G = zeros(Float64, length(x), length(y)) + for I ∈ CartesianIndices(G) + G[I] = x[I[1]]*y[I[2]] + end + + return G + end + + g_x = [1,2,3,4.0,5] + g_y = [5,4,3,2,1.0,11] + + G_w = prod_matrix(d_x_l, g_y) + G_e = prod_matrix(d_x_u, g_y) + G_s = prod_matrix(g_x, d_y_l) + G_n = prod_matrix(g_x, d_y_u) + + @test d_w'*g_y ≈ G_w + @test_broken d_e'*g_y ≈ G_e + @test d_s'*g_x ≈ G_s + @test_broken d_n'*g_x ≈ G_n end # -# @testset "NormalDerivative" begin -# op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) -# g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0)) -# -# d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}()) -# d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}()) -# d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}()) -# d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}()) -# -# -# v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) -# v∂x = evalOn(g, (x,y)-> 2*x + y) -# v∂y = evalOn(g, (x,y)-> 2*(y-1) + x) -# -# @test d_w isa TensorMapping{T,2,1} where T -# @test d_w' isa TensorMapping{T,1,2} where T -# -# @test domain_size(d_w, (3,2)) == (2,) -# @test domain_size(d_e, (3,2)) == (2,) -# @test domain_size(d_s, (3,2)) == (3,) -# @test domain_size(d_n, (3,2)) == (3,) -# -# @test size(d_w'*v) == (6,) -# @test size(d_e'*v) == (6,) -# @test size(d_s'*v) == (5,) -# @test size(d_n'*v) == (5,) -# -# @test d_w'*v .≈ v∂x[1,:] -# @test d_e'*v .≈ v∂x[5,:] -# @test d_s'*v .≈ v∂y[:,1] -# @test d_n'*v .≈ v∂y[:,6] -# -# -# d_x_l = zeros(Float64, 5) -# d_x_u = zeros(Float64, 5) -# for i ∈ eachindex(d_x_l) -# d_x_l[i] = op.dClosure[i-1] -# d_x_u[i] = -op.dClosure[length(d_x_u)-i] -# end -# -# d_y_l = zeros(Float64, 6) -# d_y_u = zeros(Float64, 6) -# for i ∈ eachindex(d_y_l) -# d_y_l[i] = op.dClosure[i-1] -# d_y_u[i] = -op.dClosure[length(d_y_u)-i] -# end -# -# function prod_matrix(x,y) -# G = zeros(Float64, length(x), length(y)) -# for I ∈ CartesianIndices(G) -# G[I] = x[I[1]]*y[I[2]] -# end -# -# return G -# end -# -# g_x = [1,2,3,4.0,5] -# g_y = [5,4,3,2,1.0,11] -# -# G_w = prod_matrix(d_x_l, g_y) -# G_e = prod_matrix(d_x_u, g_y) -# G_s = prod_matrix(g_x, d_y_l) -# G_n = prod_matrix(g_x, d_y_u) -# -# -# @test size(d_w*g_y) == (UnknownDim,6) -# @test size(d_e*g_y) == (UnknownDim,6) -# @test size(d_s*g_x) == (5,UnknownDim) -# @test size(d_n*g_x) == (5,UnknownDim) -# -# # These tests should be moved to where they are possible (i.e we know what the grid should be) -# @test_broken d_w*g_y .≈ G_w -# @test_broken d_e*g_y .≈ G_e -# @test_broken d_s*g_x .≈ G_s -# @test_broken d_n*g_x .≈ G_n -# end -# # @testset "BoundaryQuadrature" begin # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))