Mercurial > repos > public > sbplib_julia
changeset 504:21fba50cb5b0 feature/quadrature_as_outer_product
Use LazyOuterProduct to construct multi-dimensional quadratures. This change allwed to:
- Replace the types Quadrature and InverseQuadrature by functions returning outer products of the 1D operators.
- Avoid convoluted naming of types. DiagonalInnerProduct is now renamed to DiagonalQuadrature, similarly for InverseDiagonalInnerProduct.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Sat, 07 Nov 2020 13:07:31 +0100 |
| parents | fbbb3733650c |
| children | 26485066394a |
| files | src/SbpOperators/SbpOperators.jl src/SbpOperators/quadrature/diagonal_inner_product.jl src/SbpOperators/quadrature/diagonal_quadrature.jl src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl src/SbpOperators/quadrature/inverse_diagonal_quadrature.jl src/SbpOperators/quadrature/inverse_quadrature.jl src/SbpOperators/quadrature/quadrature.jl test/testSbpOperators.jl |
| diffstat | 8 files changed, 161 insertions(+), 216 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl Fri Nov 06 21:37:10 2020 +0100 +++ b/src/SbpOperators/SbpOperators.jl Sat Nov 07 13:07:31 2020 +0100 @@ -10,9 +10,7 @@ include("readoperator.jl") include("laplace/secondderivative.jl") include("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") end # module
--- a/src/SbpOperators/quadrature/diagonal_inner_product.jl Fri Nov 06 21:37:10 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,46 +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) where T - return @inbounds apply(H, v, I) -end - -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}, index::Index{Unknown}) where T - N = length(v); - r = getregion(Int(index), closuresize(H), N) - i = Index(Int(index), r) - return LazyTensors.apply(H, v, i) -end - -LazyTensors.apply_transpose(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(H,v,I) - -closuresize(H::DiagonalInnerProduct{T,M}) where {T,M} = M
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/quadrature/diagonal_quadrature.jl Sat Nov 07 13:07:31 2020 +0100 @@ -0,0 +1,64 @@ +""" +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{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} + +Implements the diagonal quadrature operator `H` of Dim dimension as a TensorMapping +""" +struct DiagonalQuadrature{T,M} <: TensorMapping{T,1,1} + h::T + closure::NTuple{M,T} + size::Tuple{Int} +end +export DiagonalQuadrature + +function DiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure) + return DiagonalQuadrature(spacing(g)[1], quadrature_closure, size(g)) +end + +LazyTensors.range_size(H::DiagonalQuadrature) = H.size +LazyTensors.domain_size(H::DiagonalQuadrature) = H.size + +function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, I::Index) where T + return @inbounds apply(H, v, I) +end + +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); + 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}, index::Index{Unknown}) where T + N = length(v); + r = getregion(Int(index), closuresize(H), N) + i = Index(Int(index), r) + return LazyTensors.apply(H, v, i) +end + +LazyTensors.apply_transpose(H::DiagonalQuadrature{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(H,v,I) + +closuresize(H::DiagonalQuadrature{T,M}) where {T,M} = M +export closuresize
--- a/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl Fri Nov 06 21:37:10 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,44 +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, v::AbstractVector{T}, index::Index{Unknown}) where T - N = length(v); - r = getregion(Int(index), closuresize(Hi), N) - i = Index(Int(index), r) - return LazyTensors.apply(Hi, v, i) -end - -LazyTensors.apply_transpose(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(Hi,v,I) - - -closuresize(Hi::InverseDiagonalInnerProduct{T,M}) where {T,M} = M
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/quadrature/inverse_diagonal_quadrature.jl Sat Nov 07 13:07:31 2020 +0100 @@ -0,0 +1,64 @@ +""" +inverse_diagonal_quadrature(g,quadrature_closure) + +Constructs the diagonal quadrature inverse operator `Qi` 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{Dim,T<:Real,M} <: TensorMapping{T,1,1} + +Implements the inverse diagonal inner product operator `Hi` of as a 1D TensorOperator +""" +struct InverseDiagonalQuadrature{T<:Real,M} <: TensorMapping{T,1,1} + h_inv::T + closure::NTuple{M,T} + size::Tuple{Int} +end +export InverseDiagonalQuadrature + +function InverseDiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure) + return InverseDiagonalQuadrature(inverse_spacing(g)[1], 1 ./ quadrature_closure, size(g)) +end + + +LazyTensors.range_size(Hi::InverseDiagonalQuadrature) = Hi.size +LazyTensors.domain_size(Hi::InverseDiagonalQuadrature) = Hi.size + + +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, v::AbstractVector{T}, index::Index{Unknown}) where T + N = length(v); + r = getregion(Int(index), closuresize(Hi), N) + i = Index(Int(index), r) + return LazyTensors.apply(Hi, v, i) +end + +LazyTensors.apply_transpose(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(Hi,v,I) + + +closuresize(Hi::InverseDiagonalQuadrature{T,M}) where {T,M} = M +export closuresize
--- a/src/SbpOperators/quadrature/inverse_quadrature.jl Fri Nov 06 21:37:10 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{Index,Dim}) where {T,Dim} - error("not implemented") -end - -@inline function LazyTensors.apply(Qi::InverseQuadrature{1,T}, v::AbstractVector{T}, I::Index) 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::Index, J::Index) 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{Index,Dim}) where {Dim,T} = LazyTensors.apply(Qi,v,I...)
--- a/src/SbpOperators/quadrature/quadrature.jl Fri Nov 06 21:37:10 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{Index,Dim}) where {T,Dim} - error("not implemented") -end - -function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::Index) 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::Index, J::Index) 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{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)
--- a/test/testSbpOperators.jl Fri Nov 06 21:37:10 2020 +0100 +++ b/test/testSbpOperators.jl Sat Nov 07 13:07:31 2020 +0100 @@ -110,67 +110,67 @@ @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2 end -@testset "DiagonalInnerProduct" begin +@testset "DiagonalQuadrature" begin op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") - L = 2.3 - g = EquidistantGrid(77, 0.0, L) - H = DiagonalInnerProduct(g,op.quadratureClosure) + Lx = 2.3 + g = EquidistantGrid(77, 0.0, Lx) + H = DiagonalQuadrature(g,op.quadratureClosure) v = ones(Float64, size(g)) @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 == diagonal_quadrature(g,op.quadratureClosure) + @test sum(H*v) ≈ Lx @test H*v == H'*v -end + @inferred H*v -@testset "Quadrature" begin - op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") - Lx = 2.3 + # Test multiple dimension Ly = 5.2 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) - Q = Quadrature(g, op.quadratureClosure) + H = diagonal_quadrature(g, op.quadratureClosure) - @test Q isa TensorMapping{T,2,2} where T - @test Q' isa TensorMapping{T,2,2} where T + @test H isa TensorMapping{T,2,2} where T + @test H' isa TensorMapping{T,2,2} where T v = ones(Float64, size(g)) - @test sum(Q*v) ≈ Lx*Ly + @test sum(H*v) ≈ Lx*Ly v = 2*ones(Float64, size(g)) - @test_broken sum(Q*v) ≈ 2*Lx*Ly + @test sum(H*v) ≈ 2*Lx*Ly - @test Q*v == Q'*v + @test_broken H*v == H'*v + + @inferred H*v end -@testset "InverseDiagonalInnerProduct" begin +@testset "InverseDiagonalQuadrature" begin op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") - L = 2.3 - g = EquidistantGrid(77, 0.0, L) - H = DiagonalInnerProduct(g, op.quadratureClosure) - Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure) + Lx = 2.3 + g = EquidistantGrid(77, 0.0, Lx) + H = DiagonalQuadrature(g, op.quadratureClosure) + Hi = InverseDiagonalQuadrature(g,op.quadratureClosure) v = evalOn(g, x->sin(x)) @test Hi isa TensorMapping{T,1,1} where T @test Hi' isa TensorMapping{T,1,1} where T + @test Hi == inverse_diagonal_quadrature(g,op.quadratureClosure) @test Hi*H*v ≈ v @test Hi*v == Hi'*v -end + @inferred Hi*v -@testset "InverseQuadrature" begin - op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") - Lx = 7.3 - Ly = 8.2 + # Test multiple dimension + Ly = 5.2 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) - Q = Quadrature(g, op.quadratureClosure) - Qinv = InverseQuadrature(g, op.quadratureClosure) + H = diagonal_quadrature(g, op.quadratureClosure) + Hi = inverse_diagonal_quadrature(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 + @test Hi isa TensorMapping{T,2,2} where T + @test Hi' isa TensorMapping{T,2,2} where T + @test Hi*H*v ≈ v + @test_broken Hi*v == Hi'*v end # # @testset "BoundaryValue" begin