Mercurial > repos > public > sbplib_julia

changeset 304:5645021683d3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Sep 2020 20:41:31 +0200
parents d5475ad78b28 (current diff) 6fa2ba769ae3 (diff)
children 7a7d9daa9eb7 d244f2e5f822
files
diffstat 4 files changed, 164 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/InverseQuadrature.jl	Wed Sep 09 20:41:31 2020 +0200
@@ -0,0 +1,75 @@
+"""
+    Quadrature{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.
+"""
+struct Quadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim}
+    Hi::NTuple{Dim,InverseDiagonalNorm{T,N,M}}
+end
+export Quadrature
+
+LazyTensors.domain_size(Qi::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
+
+function LazyTensors.apply(Qi::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
+    error("not implemented")
+end
+
+LazyTensors.apply_transpose(Qi::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where {Dim,T} = LazyTensors.apply(Q,v,I)
+
+@inline function LazyTensors.apply(Qi::Quadrature{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
+    @inbounds q = apply(Qi.Hi[1], v , I[1])
+    return q
+end
+
+@inline function LazyTensors.apply(Qi::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
+    # Quadrature in x direction
+    @inbounds vx = view(v, :, Int(I[2]))
+    @inbounds qx_inv = apply(Qi.Hi[1], vx , I[1])
+    # Quadrature in y-direction
+    @inbounds vy = view(v, Int(I[1]), :)
+    @inbounds qy_inv = apply(Qi.Hi[2], vy, I[2])
+    return qx_inv*qy_inv
+end
+
+"""
+    Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+
+Implements the quadrature operator `Hi` of Dim dimension as a TensorMapping
+"""
+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::NTuple{1,Index}) where T
+    return @inbounds apply(Hi, v, I[1])
+end
+
+LazyTensors.apply_transpose(Hi::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::NTuple{2,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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/Quadrature.jl	Wed Sep 09 20:41:31 2020 +0200
@@ -0,0 +1,76 @@
+# 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.
+"""
+struct Quadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim}
+    H::NTuple{Dim,DiagonalNorm{T,N,M}}
+end
+export Quadrature
+
+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::NTuple{Dim,Index}) where {T,Dim}
+    error("not implemented")
+end
+
+LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where {Dim,T} = LazyTensors.apply(Q,v,I)
+
+@inline function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
+    @inbounds q = apply(Q.H[1], v , I[1])
+    return q
+end
+
+@inline function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
+    # Quadrature in x direction
+    @inbounds vx = view(v, :, Int(I[2]))
+    @inbounds qx = apply(Q.H[1], vx , I[1])
+    # Quadrature in y-direction
+    @inbounds vy = view(v, Int(I[1]), :)
+    @inbounds qy = apply(Q.H[2], vy, I[2])
+    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
+"""
+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::NTuple{1,Index}) where T
+    return @inbounds apply(H, v, I[1])
+end
+
+LazyTensors.apply_transpose(H::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::NTuple{2,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
+export LazyTensors.apply
+
+function closuresize(H::DiagonalNorm{T<:Real,N,M}) where {T,N,M}
+    return M
+end
--- a/SbpOperators/src/laplace/laplace.jl	Wed Sep 09 20:41:12 2020 +0200
+++ b/SbpOperators/src/laplace/laplace.jl	Wed Sep 09 20:41:31 2020 +0200
@@ -2,8 +2,8 @@
     Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
 
 Implements the Laplace operator `L` in Dim dimensions as a tensor operator
-The multi-dimensional tensor operator simply consists of a tuple of the 1D
-Laplace tensor operator as defined by ConstantLaplaceOperator.
+The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
+tensor operators.
 """
 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
     D2::NTuple(Dim,SecondDerivative{T,N,M,K})
@@ -17,20 +17,23 @@
     error("not implemented")
 end
 
+function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} = LazyTensors.apply(L, v, I)
+
 # u = L*v
 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
-    return apply(L.D2[1],v,I)
+    @inbounds u = apply(L.D2[1],v,I)
+    return u
 end
 
 
 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = apply(L.D2[1], vx , (I[1],)) #Tuple conversion here is ugly. How to do it? Should we use indexing here?
+    @inbounds uᵢ = apply(L.D2[1], vx , I[1])
 
     # 2nd y-derivative
     @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += apply(L.D2[2], vy , (I[2],)) #Tuple conversion here is ugly. How to do it?
+    @inbounds uᵢ += apply(L.D2[2], vy , I[2])
 
     return uᵢ
 end
@@ -42,56 +45,6 @@
 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
 export quadrature
 
-# 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 `H` of Dim dimension as a TensorMapping
-"""
-struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
-    op::D2{T,N,M,K}
-    grid::EquidistantGrid{Dim,T}
-end
-export Quadrature
-
-LazyTensors.domain_size(H::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
-
-@inline function LazyTensors.apply(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
-    N = size(H.grid)
-    # Quadrature in x direction
-    @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1])
-    # Quadrature in y-direction
-    @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2])
-    return q
-end
-
-LazyTensors.apply_transpose(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H,v,I)
-
-
-"""
-    InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
-"""
-struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
-    op::D2{T,N,M,K}
-    grid::EquidistantGrid{Dim,T}
-end
-export InverseQuadrature
-
-LazyTensors.domain_size(H_inv::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
-
-@inline function LazyTensors.apply(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
-    N = size(H_inv.grid)
-    # Inverse quadrature in x direction
-    @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1])
-    # Inverse quadrature in y-direction
-    @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[2], q_inv, I[2], N[2])
-    return q_inv
-end
-
-LazyTensors.apply_transpose(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H_inv,v,I)
-
 """
     BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
 
--- a/SbpOperators/src/laplace/secondderivative.jl	Wed Sep 09 20:41:12 2020 +0200
+++ b/SbpOperators/src/laplace/secondderivative.jl	Wed Sep 09 20:41:31 2020 +0200
@@ -18,10 +18,12 @@
 
 LazyTensors.domain_size(D2::SecondDerivative, range_size::NTuple{1,Integer}) = range_size
 
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
-    return apply(D2, v, I[1])
+@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
+    return @inbounds apply(D2, v, I[1])
 end
 
+function LazyTensors.apply_transpose(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T = LazyTensors.apply(D2, v, I)
+
 # Apply for different regions Lower/Interior/Upper or Unknown region
 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Lower})
     return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i))
@@ -43,6 +45,6 @@
     return apply(D2, v, i)
 end
 
-function closuresize(D2::SecondDerivative{T<:Real,N,M,K}) where T,N,M,K
+function closuresize(D2::SecondDerivative{T<:Real,N,M,K}) where {T,N,M,K}
     return M
 end