Mercurial > repos > public > sbplib_julia
changeset 329:408c37b295c2
Refactor 1D tensor mapping in inverse quadrature to separate file, InverseDiagonalNorm. Add tests
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Fri, 25 Sep 2020 09:34:37 +0200 |
| parents | 9cc5d1498b2d |
| children | 8d1a830b0c22 |
| files | SbpOperators/src/InverseQuadrature.jl SbpOperators/src/SbpOperators.jl SbpOperators/src/quadrature/inversequadrature.jl SbpOperators/test/runtests.jl |
| diffstat | 4 files changed, 76 insertions(+), 91 deletions(-) [+] |
line wrap: on
line diff
--- a/SbpOperators/src/InverseQuadrature.jl Thu Sep 24 22:31:48 2020 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -""" - InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the inverse quadrature operator `Qi` of Dim dimension as a TensorOperator -The multi-dimensional tensor operator consists of a tuple of 1D InverseDiagonalNorm -tensor operators. -""" -export InverseQuadrature -struct InverseQuadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim} - Hi::NTuple{Dim,InverseDiagonalNorm{T,N,M}} -end - -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 - -LazyTensors.apply_transpose(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I) - -@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 - -""" - DiagonalNorm{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the inverse diagnoal norm operator `Hi` of Dim dimension as a TensorMapping -""" -export InverseDiagonalNorm, closuresize -struct InverseDiagonalNorm{T<:Real,N,M} <: TensorOperator{T,1} - h_inv::T # The reciprocl grid spacing could be included in the stencil already. Preferable? - closure::NTuple{M,T} - #TODO: Write a nice constructor -end - -@inline function LazyTensors.apply(Hi::InverseDiagonalNorm{T}, v::AbstractVector{T}, I:Index) where T - return @inbounds apply(Hi, v, I) -end - -LazyTensors.apply_transpose(Hi::InverseQuadrature{Dim,T}, v::AbstractArray{T,2}, I::Index) where T = LazyTensors.apply(Hi,v,I) - -@inline LazyTensors.apply(Hi::InverseDiagonalNorm, v::AbstractVector{T}, I::Index{Lower}) where T - return @inbounds Hi.h_inv*Hi.closure[Int(i)]*v[Int(I)] -end -@inline LazyTensors.apply(Hi::InverseDiagonalNorm,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 - -@inline LazyTensors.apply(Hi::InverseDiagonalNorm, v::AbstractVector{T}, I::Index{Interior}) where T - return @inbounds Hi.h_inv*v[Int(I)] -end - -function LazyTensors.apply(Hi::InverseDiagonalNorm, 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 -export LazyTensors.apply - -function closuresize(Hi::InverseDiagonalNorm{T<:Real,N,M}) where {T,N,M} - return M -end
--- a/SbpOperators/src/SbpOperators.jl Thu Sep 24 22:31:48 2020 +0200 +++ b/SbpOperators/src/SbpOperators.jl Fri Sep 25 09:34:37 2020 +0200 @@ -11,4 +11,6 @@ include("laplace/laplace.jl") include("quadrature/diagonal_inner_product.jl") include("quadrature/quadrature.jl") +include("quadrature/inverse_diagonal_inner_product.jl") +include("quadrature/inversequadrature.jl") end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SbpOperators/src/quadrature/inversequadrature.jl Fri Sep 25 09:34:37 2020 +0200 @@ -0,0 +1,34 @@ +export InverseQuadrature +""" + InverseQuadrature{Dim,T<:Real,M,K} <: TensorMapping{T,Dim,Dim} + +Implements the inverse quadrature operator `Qi` of Dim dimension as a TensorOperator +The multi-dimensional tensor operator consists of a tuple of 1D InverseDiagonalInnerProduct +tensor operators. +""" +struct InverseQuadrature{Dim,T<:Real,M} <: TensorOperator{T,Dim} + Hi::NTuple{Dim,InverseDiagonalInnerProduct{T,M}} +end + +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/SbpOperators/test/runtests.jl Thu Sep 24 22:31:48 2020 +0200 +++ b/SbpOperators/test/runtests.jl Fri Sep 25 09:34:37 2020 +0200 @@ -152,21 +152,46 @@ @test sum(collect(Q*v)) ≈ (Lx*Ly) @test collect(Q*v) == collect(Q'*v) end -# -# @testset "InverseQuadrature" begin -# op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") -# Lx = 7.3 -# Ly = 8.2 -# g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) -# H = Quadrature(op,g) -# Hinv = InverseQuadrature(op,g) -# v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) -# -# @test Hinv isa TensorOperator{T,2} where T -# @test Hinv' isa TensorMapping{T,2,2} where T -# @test collect(Hinv*H*v) ≈ v -# @test collect(Hinv*v) == collect(Hinv'*v) -# end + +@testset "InverseDiagonalInnerProduct" 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 = spacing(g) + H = DiagonalInnerProduct(h[1],op.quadratureClosure) + + h_i = inverse_spacing(g) + Hi = InverseDiagonalInnerProduct(h_i[1],1 ./ op.quadratureClosure) + v = evalOn(g, x->sin(x)) + + @test Hi isa TensorOperator{T,1} where T + @test Hi' isa TensorMapping{T,1,1} where T + @test collect(Hi*H*v) ≈ v + @test collect(Hi*v) == collect(Hi'*v) +end + +@testset "InverseQuadrature" begin + op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") + Lx = 7.3 + Ly = 8.2 + g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) + + h = spacing(g) + Hx = DiagonalInnerProduct(h[1], op.quadratureClosure); + Hy = DiagonalInnerProduct(h[2], op.quadratureClosure); + Q = Quadrature((Hx,Hy)) + + hi = inverse_spacing(g) + Hix = InverseDiagonalInnerProduct(hi[1], 1 ./ op.quadratureClosure); + Hiy = InverseDiagonalInnerProduct(hi[2], 1 ./ op.quadratureClosure); + Qinv = InverseQuadrature((Hix,Hiy)) + v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) + + @test Qinv isa TensorOperator{T,2} where T + @test Qinv' isa TensorMapping{T,2,2} where T + @test collect(Qinv*Q*v) ≈ v + @test collect(Qinv*v) == collect(Qinv'*v) +end # # @testset "BoundaryValue" begin # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")