Mercurial > repos > public > sbplib_julia
changeset 331:c8bbb4a83518
Merge
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Fri, 25 Sep 2020 09:36:06 +0200 |
| parents | 8d1a830b0c22 (diff) 41c3c25e4e3b (current diff) |
| children | 535f1bff4bcc 3c2238c681b5 |
| files | |
| diffstat | 9 files changed, 195 insertions(+), 199 deletions(-) [+] |
line wrap: on
line diff
--- a/SbpOperators/src/InverseQuadrature.jl Thu Sep 24 22:31:04 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 - -""" - InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the quadrature 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/Quadrature.jl Thu Sep 24 22:31:04 2020 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -# At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings -""" - 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 DiagonalNorm H -tensor operators. -""" -export Quadrature -struct Quadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim} - H::NTuple{Dim,DiagonalNorm{T,N,M}} -end - -LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size - -function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} - error("not implemented") -end - -LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I) - -@inline 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 - -@inline 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 - -""" - DiagonalNorm{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} - -Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping -""" -export DiagonalNorm, closuresize, LazyTensors.apply -struct DiagonalNorm{T<:Real,N,M} <: TensorOperator{T,1} - h::T # The grid spacing could be included in the stencil already. Preferable? - closure::NTuple{M,T} - #TODO: Write a nice constructor -end - -@inline function LazyTensors.apply(H::DiagonalNorm{T}, v::AbstractVector{T}, I::Index) where T - return @inbounds apply(H, v, I) -end - -LazyTensors.apply_transpose(H::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::Index) where T = LazyTensors.apply(H,v,I) - -@inline LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, I::Index{Lower}) where T - return @inbounds H.h*H.closure[Int(I)]*v[Int(I)] -end -@inline LazyTensors.apply(H::DiagonalNorm,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 - -@inline LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, I::Index{Interior}) where T - return @inbounds H.h*v[Int(I)] -end - -function LazyTensors.apply(H::DiagonalNorm, 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 - -function closuresize(H::DiagonalNorm{T<:Real,N,M}) where {T,N,M} - return M -end
--- a/SbpOperators/src/SbpOperators.jl Thu Sep 24 22:31:04 2020 +0200 +++ b/SbpOperators/src/SbpOperators.jl Fri Sep 25 09:36:06 2020 +0200 @@ -9,4 +9,8 @@ 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") end # module
--- a/SbpOperators/src/laplace/laplace.jl Thu Sep 24 22:31:04 2020 +0200 +++ b/SbpOperators/src/laplace/laplace.jl Fri Sep 25 09:36:06 2020 +0200 @@ -12,7 +12,7 @@ #TODO: Write a good constructor end -LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size +LazyTensors.domain_size(L::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} error("not implemented") @@ -24,8 +24,7 @@ return u end -# TODO: Fix dispatch on tuples! -@inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T +function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T # 2nd x-derivative @inbounds vx = view(v, :, Int(J)) @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , I) @@ -37,9 +36,7 @@ return uᵢ end -@inline function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} - return LazyTensors.apply(L, v, I) -end +LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} = LazyTensors.apply(L, v, I...) # quadrature(L::Laplace) = Quadrature(L.op, L.grid) # inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
--- a/SbpOperators/src/laplace/secondderivative.jl Thu Sep 24 22:31:04 2020 +0200 +++ b/SbpOperators/src/laplace/secondderivative.jl Fri Sep 25 09:36:06 2020 +0200 @@ -20,32 +20,26 @@ # I thought I::Vararg{Index,R} fell back to just Index for R = 1 # Apply for different regions Lower/Interior/Upper or Unknown region -@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Lower}) where T +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 -@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Interior}) where T +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 -@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T +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*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I)) end -@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T +function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T N = length(v) # TODO: Use domain_size here instead? r = getregion(Int(index), closuresize(D2), N) I = Index(Int(index), r) return LazyTensors.apply(D2, v, I) end - -@inline function LazyTensors.apply_transpose(D2::SecondDerivative, v::AbstractVector, I::Index) - return LazyTensors.apply(D2, v, I) -end +LazyTensors.apply_transpose(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index) where {T} = LazyTensors.apply(D2, v, I) - -function closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} - return M -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/SbpOperators/src/quadrature/diagonal_inner_product.jl Fri Sep 25 09:36:06 2020 +0200 @@ -0,0 +1,41 @@ +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} <: TensorOperator{T,1} + h::T # The grid spacing could be included in the stencil already. Preferable? + closure::NTuple{M,T} + #TODO: Write a nice constructor +end + +LazyTensors.domain_size(H::DiagonalInnerProduct, range_size::NTuple{1,Integer}) = range_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.closure[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.closure[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/SbpOperators/src/quadrature/inverse_quadrature.jl Fri Sep 25 09:36:06 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...)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SbpOperators/src/quadrature/quadrature.jl Fri Sep 25 09:36:06 2020 +0200 @@ -0,0 +1,34 @@ +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} <: TensorOperator{T,Dim} + H::NTuple{Dim,DiagonalInnerProduct{T,M}} +end + +LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size + +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/SbpOperators/test/runtests.jl Thu Sep 24 22:31:04 2020 +0200 +++ b/SbpOperators/test/runtests.jl Fri Sep 25 09:36:06 2020 +0200 @@ -119,35 +119,79 @@ @test sqrt(prod(h)*sum(collect(e4.^2))) <= accuracytol @test sqrt(prod(h)*sum(collect(e5.^2))) <= accuracytol end -# -# @testset "Quadrature" begin -# op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") -# Lx = 2.3 -# Ly = 5.2 -# g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) -# H = Quadrature(op,g) -# v = ones(Float64, size(g)) -# -# @test H isa TensorOperator{T,2} where T -# @test H' isa TensorMapping{T,2,2} where T -# @test sum(collect(H*v)) ≈ (Lx*Ly) -# @test collect(H*v) == collect(H'*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 "DiagonalInnerProduct" 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) + v = ones(Float64, size(g)) + + @test H isa TensorOperator{T,1} where T + @test H' isa TensorMapping{T,1,1} where T + @test sum(collect(H*v)) ≈ L + @test collect(H*v) == collect(H'*v) +end + +@testset "Quadrature" begin + op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") + Lx = 2.3 + Ly = 5.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)) + + v = ones(Float64, size(g)) + + @test Q isa TensorOperator{T,2} where T + @test Q' isa TensorMapping{T,2,2} where T + @test sum(collect(Q*v)) ≈ (Lx*Ly) + @test collect(Q*v) == collect(Q'*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")