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 (current diff) a78bda7084f6 (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 | 8 files changed, 321 insertions(+), 232 deletions(-) [+] |
line wrap: on
line diff
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/SbpOperators.jl --- a/src/SbpOperators/SbpOperators.jl Thu Dec 31 08:41:07 2020 +0100 +++ b/src/SbpOperators/SbpOperators.jl Fri Jan 01 16:39:57 2021 +0100 @@ -11,10 +11,8 @@ include("volumeops/volume_operator.jl") include("volumeops/derivatives/secondderivative.jl") include("volumeops/laplace/laplace.jl") -include("quadrature/diagonal_inner_product.jl") -include("quadrature/quadrature.jl") -include("quadrature/inverse_diagonal_inner_product.jl") -include("quadrature/inverse_quadrature.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")
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/quadrature/diagonal_inner_product.jl --- a/src/SbpOperators/quadrature/diagonal_inner_product.jl Thu Dec 31 08:41:07 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ -export DiagonalInnerProduct, closuresize -""" - DiagonalInnerProduct{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping -""" -struct DiagonalInnerProduct{T,M} <: TensorMapping{T,1,1} - h::T - quadratureClosure::NTuple{M,T} - size::Tuple{Int} -end - -function DiagonalInnerProduct(g::EquidistantGrid{1}, quadratureClosure) - return DiagonalInnerProduct(spacing(g)[1], quadratureClosure, size(g)) -end - -LazyTensors.range_size(H::DiagonalInnerProduct) = H.size -LazyTensors.domain_size(H::DiagonalInnerProduct) = H.size - -function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Lower}) where T - return @inbounds H.h*H.quadratureClosure[Int(i)]*v[Int(i)] -end - -function LazyTensors.apply(H::DiagonalInnerProduct{T},v::AbstractVector{T}, i::Index{Upper}) where T - N = length(v); - return @inbounds H.h*H.quadratureClosure[N-Int(i)+1]*v[Int(i)] -end - -function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Interior}) where T - return @inbounds H.h*v[Int(i)] -end - -function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T - N = length(v); - r = getregion(i, closuresize(H), N) - return LazyTensors.apply(H, v, Index(i, r)) -end - -LazyTensors.apply_transpose(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(H,v,i) - -closuresize(H::DiagonalInnerProduct{T,M}) where {T,M} = M
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/quadrature/diagonal_quadrature.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/quadrature/diagonal_quadrature.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,92 @@ +""" +diagonal_quadrature(g,quadrature_closure) + +Constructs the diagonal quadrature operator `H` on a grid of `Dim` dimensions as +a `TensorMapping`. The one-dimensional operator is a `DiagonalQuadrature`, while +the multi-dimensional operator is the outer-product of the +one-dimensional operators in each coordinate direction. +""" +function diagonal_quadrature(g::EquidistantGrid{Dim}, quadrature_closure) where Dim + H = DiagonalQuadrature(restrict(g,1), quadrature_closure) + for i ∈ 2:Dim + H = H⊗DiagonalQuadrature(restrict(g,i), quadrature_closure) + end + return H +end +export diagonal_quadrature + +""" + DiagonalQuadrature{T,M} <: TensorMapping{T,1,1} + +Implements the one-dimensional diagonal quadrature operator as a `TensorMapping` +The quadrature is defined by the quadrature interval length `h`, the quadrature +closure weights `closure` and the number of quadrature intervals `size`. The +interior stencil has the weight 1. +""" +struct DiagonalQuadrature{T,M} <: TensorMapping{T,1,1} + h::T + closure::NTuple{M,T} + size::Tuple{Int} +end +export DiagonalQuadrature + +""" + DiagonalQuadrature(g, quadrature_closure) + +Constructs the `DiagonalQuadrature` on the `EquidistantGrid` `g` with +closure given by `quadrature_closure`. +""" +function DiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure) + return DiagonalQuadrature(spacing(g)[1], quadrature_closure, size(g)) +end + +""" + range_size(H::DiagonalQuadrature) + +The size of an object in the range of `H` +""" +LazyTensors.range_size(H::DiagonalQuadrature) = H.size + +""" + domain_size(H::DiagonalQuadrature) + +The size of an object in the domain of `H` +""" +LazyTensors.domain_size(H::DiagonalQuadrature) = H.size + +""" + apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i) where T +Implements the application `(H*v)[i]` an `Index{R}` where `R` is one of the regions +`Lower`,`Interior`,`Upper`. If `i` is another type of index (e.g an `Int`) it will first +be converted to an `Index{R}`. +""" +function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Lower}) where T + return @inbounds H.h*H.closure[Int(i)]*v[Int(i)] +end + +function LazyTensors.apply(H::DiagonalQuadrature{T},v::AbstractVector{T}, i::Index{Upper}) where T + N = length(v); #TODO: Use dim_size here? + return @inbounds H.h*H.closure[N-Int(i)+1]*v[Int(i)] +end + +function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Interior}) where T + return @inbounds H.h*v[Int(i)] +end + +function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i) where T + N = length(v); #TODO: Use dim_size here? + r = getregion(i, closure_size(H), N) + return LazyTensors.apply(H, v, Index(i, r)) +end + +""" + apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, I::Index) where T +Implements the application (H'*v)[I]. The operator is self-adjoint. +""" +LazyTensors.apply_transpose(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(H,v,i) + +""" + closure_size(H) +Returns the size of the closure stencil of a DiagonalQuadrature `H`. +""" +closure_size(H::DiagonalQuadrature{T,M}) where {T,M} = M
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl --- a/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl Thu Dec 31 08:41:07 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,43 +0,0 @@ -export InverseDiagonalInnerProduct, closuresize -""" - InverseDiagonalInnerProduct{Dim,T<:Real,M} <: TensorMapping{T,1,1} - -Implements the inverse diagonal inner product operator `Hi` of as a 1D TensorOperator -""" -struct InverseDiagonalInnerProduct{T<:Real,M} <: TensorMapping{T,1,1} - h_inv::T - inverseQuadratureClosure::NTuple{M,T} - size::Tuple{Int} -end - -function InverseDiagonalInnerProduct(g::EquidistantGrid{1}, quadratureClosure) - return InverseDiagonalInnerProduct(inverse_spacing(g)[1], 1 ./ quadratureClosure, size(g)) -end - -LazyTensors.range_size(Hi::InverseDiagonalInnerProduct) = Hi.size -LazyTensors.domain_size(Hi::InverseDiagonalInnerProduct) = Hi.size - - -function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Lower}) where T - return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[Int(i)]*v[Int(i)] -end - -function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Upper}) where T - N = length(v); - return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[N-Int(i)+1]*v[Int(i)] -end - -function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Interior}) where T - return @inbounds Hi.h_inv*v[Int(i)] -end - -function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T - N = length(v); - r = getregion(i, closuresize(Hi), N) - return LazyTensors.apply(Hi, v, Index(i, r)) -end - -LazyTensors.apply_transpose(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(Hi,v,i) - - -closuresize(Hi::InverseDiagonalInnerProduct{T,M}) where {T,M} = M
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/quadrature/inverse_diagonal_quadrature.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/quadrature/inverse_diagonal_quadrature.jl Fri Jan 01 16:39:57 2021 +0100 @@ -0,0 +1,89 @@ +""" +inverse_diagonal_quadrature(g,quadrature_closure) + +Constructs the inverse `Hi` of a `DiagonalQuadrature` on a grid of `Dim` dimensions as +a `TensorMapping`. The one-dimensional operator is a `InverseDiagonalQuadrature`, while +the multi-dimensional operator is the outer-product of the one-dimensional operators +in each coordinate direction. +""" +function inverse_diagonal_quadrature(g::EquidistantGrid{Dim}, quadrature_closure) where Dim + Hi = InverseDiagonalQuadrature(restrict(g,1), quadrature_closure) + for i ∈ 2:Dim + Hi = Hi⊗InverseDiagonalQuadrature(restrict(g,i), quadrature_closure) + end + return Hi +end +export inverse_diagonal_quadrature + + +""" + InverseDiagonalQuadrature{T,M} <: TensorMapping{T,1,1} + +Implements the inverse of a one-dimensional `DiagonalQuadrature` as a `TensorMapping` +The operator is defined by the reciprocal of the quadrature interval length `h_inv`, the +reciprocal of the quadrature closure weights `closure` and the number of quadrature intervals `size`. The +interior stencil has the weight 1. +""" +struct InverseDiagonalQuadrature{T<:Real,M} <: TensorMapping{T,1,1} + h_inv::T + closure::NTuple{M,T} + size::Tuple{Int} +end +export InverseDiagonalQuadrature + +""" + InverseDiagonalQuadrature(g, quadrature_closure) + +Constructs the `InverseDiagonalQuadrature` on the `EquidistantGrid` `g` with +closure given by the reciprocal of `quadrature_closure`. +""" +function InverseDiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure) + return InverseDiagonalQuadrature(inverse_spacing(g)[1], 1 ./ quadrature_closure, size(g)) +end + +""" + domain_size(Hi::InverseDiagonalQuadrature) + +The size of an object in the range of `Hi` +""" +LazyTensors.range_size(Hi::InverseDiagonalQuadrature) = Hi.size + +""" + domain_size(Hi::InverseDiagonalQuadrature) + +The size of an object in the domain of `Hi` +""" +LazyTensors.domain_size(Hi::InverseDiagonalQuadrature) = Hi.size + +""" + apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i) where T +Implements the application `(Hi*v)[i]` an `Index{R}` where `R` is one of the regions +`Lower`,`Interior`,`Upper`. If `i` is another type of index (e.g an `Int`) it will first +be converted to an `Index{R}`. +""" +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Lower}) where T + return @inbounds Hi.h_inv*Hi.closure[Int(i)]*v[Int(i)] +end + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Upper}) where T + N = length(v); + return @inbounds Hi.h_inv*Hi.closure[N-Int(i)+1]*v[Int(i)] +end + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Interior}) where T + return @inbounds Hi.h_inv*v[Int(i)] +end + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i) where T + N = length(v); + r = getregion(i, closure_size(Hi), N) + return LazyTensors.apply(Hi, v, Index(i, r)) +end + +LazyTensors.apply_transpose(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(Hi,v,i) + +""" + closure_size(Hi) +Returns the size of the closure stencil of a InverseDiagonalQuadrature `Hi`. +""" +closure_size(Hi::InverseDiagonalQuadrature{T,M}) where {T,M} = M
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/quadrature/inverse_quadrature.jl --- a/src/SbpOperators/quadrature/inverse_quadrature.jl Thu Dec 31 08:41:07 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -export InverseQuadrature -""" - InverseQuadrature{Dim,T<:Real,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the inverse quadrature operator `Qi` of Dim dimension as a TensorMapping -The multi-dimensional tensor operator consists of a tuple of 1D InverseDiagonalInnerProduct -tensor operators. -""" -struct InverseQuadrature{Dim,T<:Real,M} <: TensorMapping{T,Dim,Dim} - Hi::NTuple{Dim,InverseDiagonalInnerProduct{T,M}} -end - - -function InverseQuadrature(g::EquidistantGrid{Dim}, quadratureClosure) where Dim - Hi = () - for i ∈ 1:Dim - Hi = (Hi..., InverseDiagonalInnerProduct(restrict(g,i), quadratureClosure)) - end - - return InverseQuadrature(Hi) -end - -LazyTensors.range_size(Hi::InverseQuadrature) = getindex.(range_size.(Hi.Hi),1) -LazyTensors.domain_size(Hi::InverseQuadrature) = getindex.(domain_size.(Hi.Hi),1) - -LazyTensors.domain_size(Qi::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size - -function LazyTensors.apply(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim} - error("not implemented") -end - -@inline function LazyTensors.apply(Qi::InverseQuadrature{1,T}, v::AbstractVector{T}, i) where T - @inbounds q = apply(Qi.Hi[1], v , i) - return q -end - -@inline function LazyTensors.apply(Qi::InverseQuadrature{2,T}, v::AbstractArray{T,2}, i,j) where T - # InverseQuadrature in x direction - @inbounds vx = view(v, :, Int(j)) - @inbounds qx_inv = apply(Qi.Hi[1], vx , i) - # InverseQuadrature in y-direction - @inbounds vy = view(v, Int(i), :) - @inbounds qy_inv = apply(Qi.Hi[2], vy, j) - return qx_inv*qy_inv -end - -LazyTensors.apply_transpose(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {Dim,T} = LazyTensors.apply(Qi,v,I...)
diff -r bf8b66c596f7 -r fb5ac62563aa src/SbpOperators/quadrature/quadrature.jl --- a/src/SbpOperators/quadrature/quadrature.jl Thu Dec 31 08:41:07 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,44 +0,0 @@ -export Quadrature -""" - Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the quadrature operator `Q` of Dim dimension as a TensorMapping -The multi-dimensional tensor operator consists of a tuple of 1D DiagonalInnerProduct H -tensor operators. -""" -struct Quadrature{Dim,T<:Real,M} <: TensorMapping{T,Dim,Dim} - H::NTuple{Dim,DiagonalInnerProduct{T,M}} -end - -function Quadrature(g::EquidistantGrid{Dim}, quadratureClosure) where Dim - H = () - for i ∈ 1:Dim - H = (H..., DiagonalInnerProduct(restrict(g,i), quadratureClosure)) - end - - return Quadrature(H) -end - -LazyTensors.range_size(H::Quadrature) = getindex.(range_size.(H.H),1) -LazyTensors.domain_size(H::Quadrature) = getindex.(domain_size.(H.H),1) - -function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim} - error("not implemented") -end - -function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, i) where T - @inbounds q = apply(Q.H[1], v , i) - return q -end - -function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, i, j) where T - # Quadrature in x direction - @inbounds vx = view(v, :, Int(j)) - @inbounds qx = apply(Q.H[1], vx , i) - # Quadrature in y-direction - @inbounds vy = view(v, Int(i), :) - @inbounds qy = apply(Q.H[2], vy, j) - return qx*qy -end - -LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)
diff -r bf8b66c596f7 -r fb5ac62563aa test/testSbpOperators.jl --- a/test/testSbpOperators.jl Thu Dec 31 08:41:07 2020 +0100 +++ b/test/testSbpOperators.jl Fri Jan 01 16:39:57 2021 +0100 @@ -405,67 +405,152 @@ end end -@testset "DiagonalInnerProduct" begin +@testset "DiagonalQuadrature" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) - L = 2.3 - g = EquidistantGrid(77, 0.0, L) - H = DiagonalInnerProduct(g,op.quadratureClosure) - v = ones(Float64, size(g)) + Lx = π/2. + Ly = Float64(π) + g_1D = EquidistantGrid(77, 0.0, Lx) + g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) + integral(H,v) = sum(H*v) + @testset "Constructors" begin + # 1D + H_x = DiagonalQuadrature(spacing(g_1D)[1],op.quadratureClosure,size(g_1D)); + @test H_x == DiagonalQuadrature(g_1D,op.quadratureClosure) + @test H_x == diagonal_quadrature(g_1D,op.quadratureClosure) + @test H_x isa TensorMapping{T,1,1} where T + @test H_x' isa TensorMapping{T,1,1} where T + # 2D + H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) + @test H_xy isa TensorMappingComposition + @test H_xy isa TensorMapping{T,2,2} where T + @test H_xy' isa TensorMapping{T,2,2} where T + end + + @testset "Sizes" begin + # 1D + H_x = diagonal_quadrature(g_1D,op.quadratureClosure) + @test domain_size(H_x) == size(g_1D) + @test range_size(H_x) == size(g_1D) + # 2D + H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) + @test domain_size(H_xy) == size(g_2D) + @test range_size(H_xy) == size(g_2D) + end - @test H isa TensorMapping{T,1,1} where T - @test H' isa TensorMapping{T,1,1} where T - @test sum(H*v) ≈ L - @test H*v == H'*v + @testset "Application" begin + # 1D + H_x = diagonal_quadrature(g_1D,op.quadratureClosure) + a = 3.2 + v_1D = a*ones(Float64, size(g_1D)) + u_1D = evalOn(g_1D,x->sin(x)) + @test integral(H_x,v_1D) ≈ a*Lx rtol = 1e-13 + @test integral(H_x,u_1D) ≈ 1. rtol = 1e-8 + @test H_x*v_1D == H_x'*v_1D + # 2D + H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) + b = 2.1 + v_2D = b*ones(Float64, size(g_2D)) + u_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y)) + @test integral(H_xy,v_2D) ≈ b*Lx*Ly rtol = 1e-13 + @test integral(H_xy,u_2D) ≈ π rtol = 1e-8 + @test H_xy*v_2D ≈ H_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason? + end + + @testset "Accuracy" begin + v = () + for i = 0:4 + f_i(x) = 1/factorial(i)*x^i + v = (v...,evalOn(g_1D,f_i)) + end + # TODO: Bug in readOperator for 2nd order + # # 2nd order + # op2 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) + # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure) + # for i = 1:3 + # @test integral(H2,v[i]) ≈ v[i+1] rtol = 1e-14 + # end + + # 4th order + op4 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + H4 = diagonal_quadrature(g_1D,op4.quadratureClosure) + for i = 1:4 + @test integral(H4,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 + end + end + + @testset "Inferred" begin + H_x = diagonal_quadrature(g_1D,op.quadratureClosure) + H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) + v_1D = ones(Float64, size(g_1D)) + v_2D = ones(Float64, size(g_2D)) + @inferred H_x*v_1D + @inferred H_x'*v_1D + @inferred H_xy*v_2D + @inferred H_xy'*v_2D + end end -@testset "Quadrature" begin - op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) - Lx = 2.3 - Ly = 5.2 - g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) - - Q = Quadrature(g, op.quadratureClosure) - - @test Q isa TensorMapping{T,2,2} where T - @test Q' isa TensorMapping{T,2,2} where T - - v = ones(Float64, size(g)) - @test sum(Q*v) ≈ Lx*Ly - - v = 2*ones(Float64, size(g)) - @test_broken sum(Q*v) ≈ 2*Lx*Ly - - @test Q*v == Q'*v -end - -@testset "InverseDiagonalInnerProduct" begin +@testset "InverseDiagonalQuadrature" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) - L = 2.3 - g = EquidistantGrid(77, 0.0, L) - H = DiagonalInnerProduct(g, op.quadratureClosure) - Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure) - v = evalOn(g, x->sin(x)) + Lx = π/2. + Ly = Float64(π) + g_1D = EquidistantGrid(77, 0.0, Lx) + g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) + @testset "Constructors" begin + # 1D + Hi_x = InverseDiagonalQuadrature(inverse_spacing(g_1D)[1], 1. ./ op.quadratureClosure, size(g_1D)); + @test Hi_x == InverseDiagonalQuadrature(g_1D,op.quadratureClosure) + @test Hi_x == inverse_diagonal_quadrature(g_1D,op.quadratureClosure) + @test Hi_x isa TensorMapping{T,1,1} where T + @test Hi_x' isa TensorMapping{T,1,1} where T - @test Hi isa TensorMapping{T,1,1} where T - @test Hi' isa TensorMapping{T,1,1} where T - @test Hi*H*v ≈ v - @test Hi*v == Hi'*v -end + # 2D + Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) + @test Hi_xy isa TensorMappingComposition + @test Hi_xy isa TensorMapping{T,2,2} where T + @test Hi_xy' isa TensorMapping{T,2,2} where T + end + + @testset "Sizes" begin + # 1D + Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) + @test domain_size(Hi_x) == size(g_1D) + @test range_size(Hi_x) == size(g_1D) + # 2D + Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) + @test domain_size(Hi_xy) == size(g_2D) + @test range_size(Hi_xy) == size(g_2D) + end -@testset "InverseQuadrature" begin - op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) - Lx = 7.3 - Ly = 8.2 - g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) + @testset "Application" begin + # 1D + H_x = diagonal_quadrature(g_1D,op.quadratureClosure) + Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) + v_1D = evalOn(g_1D,x->sin(x)) + u_1D = evalOn(g_1D,x->x^3-x^2+1) + @test Hi_x*H_x*v_1D ≈ v_1D rtol = 1e-15 + @test Hi_x*H_x*u_1D ≈ u_1D rtol = 1e-15 + @test Hi_x*v_1D == Hi_x'*v_1D + # 2D + H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) + Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) + v_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y)) + u_2D = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) + @test Hi_xy*H_xy*v_2D ≈ v_2D rtol = 1e-15 + @test Hi_xy*H_xy*u_2D ≈ u_2D rtol = 1e-15 + @test Hi_xy*v_2D ≈ Hi_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason? + end - Q = Quadrature(g, op.quadratureClosure) - Qinv = InverseQuadrature(g, op.quadratureClosure) - v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) - - @test Qinv isa TensorMapping{T,2,2} where T - @test Qinv' isa TensorMapping{T,2,2} where T - @test_broken Qinv*(Q*v) ≈ v - @test Qinv*v == Qinv'*v + @testset "Inferred" begin + Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) + Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) + v_1D = ones(Float64, size(g_1D)) + v_2D = ones(Float64, size(g_2D)) + @inferred Hi_x*v_1D + @inferred Hi_x'*v_1D + @inferred Hi_xy*v_2D + @inferred Hi_xy'*v_2D + end end @testset "BoundaryOperator" begin