Mercurial > repos > public > sbplib_julia

changeset 269:ccef055233a2 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move general methods for D2 to ConstantStencilOperator and increase verbosity of method names for clarity
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 05 Dec 2019 10:48:31 +0100
parents f67ce2eb6019
children 12b738f260a0
files DiffOps/src/laplace.jl SbpOperators/src/constantstenciloperator.jl SbpOperators/src/d2.jl SbpOperators/src/stencil.jl
diffstat 4 files changed, 109 insertions(+), 105 deletions(-) [+]
line wrap: on
line diff
--- a/DiffOps/src/laplace.jl	Thu Dec 05 09:44:34 2019 +0100
+++ b/DiffOps/src/laplace.jl	Thu Dec 05 10:48:31 2019 +0100
@@ -10,17 +10,17 @@
 
 # u = L*v
 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
-    uᵢ = L.a * SbpOperators.apply(L.op, inverse_spacing(L.grid)[1], v, i)
+    uᵢ = L.a * SbpOperators.apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], v, i)
     return uᵢ
 end
 
 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = L.a*SbpOperators.apply(L.op, inverse_spacing(L.grid)[1], vx , I[1])
+    @inbounds uᵢ = L.a*SbpOperators.apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], vx , I[1])
     # 2nd y-derivative
     @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += L.a*SbpOperators.apply(L.op, inverse_spacing(L.grid)[2], vy, I[2])
+    @inbounds uᵢ += L.a*SbpOperators.apply_2nd_derivative(L.op, inverse_spacing(L.grid)[2], vy, I[2])
     # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
     return uᵢ
 end
@@ -113,12 +113,12 @@
     i = I[dim(e.bId)]
     j = I[3-dim(e.bId)]
     N_i = size(e.grid)[dim(e.bId)]
-    return apply_e(e.op, v[j], N_i, i, region(e.bId))
+    return apply_boundary_value(e.op, v[j], N_i, i, region(e.bId))
 end
 
 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
     u = selectdim(v,3-dim(e.bId),I[1])
-    return apply_e_T(e.op, u, region(e.bId))
+    return apply_boundary_value_transpose(e.op, u, region(e.bId))
 end
 
 """
@@ -145,12 +145,12 @@
     j = I[3-dim(d.bId)]
     N_i = size(d.grid)[dim(d.bId)]
     h_inv = inverse_spacing(d.grid)[dim(d.bId)]
-    return apply_d(d.op, h_inv, v[j], N_i, i, region(d.bId))
+    return apply_normal_derivative(d.op, h_inv, v[j], N_i, i, region(d.bId))
 end
 
 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
     u = selectdim(v,3-dim(d.bId),I[1])
-    return apply_d_T(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
+    return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
 end
 
 """
--- a/SbpOperators/src/constantstenciloperator.jl	Thu Dec 05 09:44:34 2019 +0100
+++ b/SbpOperators/src/constantstenciloperator.jl	Thu Dec 05 10:48:31 2019 +0100
@@ -1,39 +1,124 @@
-export apply
-
 abstract type ConstantStencilOperator end
 
 # Apply for different regions Lower/Interior/Upper or Unknown region
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
-    return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Lower})
+    return @inbounds h_inv*h_inv*apply_stencil(op.closureStencils[Int(i)], v, Int(i))
 end
 
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
-    return @inbounds h*h*apply(op.innerStencil, v, Int(i))
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Interior})
+    return @inbounds h_inv*h_inv*apply_stencil(op.innerStencil, v, Int(i))
 end
 
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Upper})
     N = length(v)
-    return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
+    return @inbounds h_inv*h_inv*Int(op.parity)*apply_stencil_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
 end
 
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
     cSize = closuresize(op)
     N = length(v)
 
     i = Int(index)
 
     if 0 < i <= cSize
-        return apply(op, h, v, Index{Lower}(i))
+        return apply_2nd_derivative(op, h_inv, v, Index{Lower}(i))
     elseif cSize < i <= N-cSize
-        return apply(op, h, v, Index{Interior}(i))
+        return apply_2nd_derivative(op, h_inv, v, Index{Interior}(i))
     elseif N-cSize < i <= N
-        return apply(op, h, v, Index{Upper}(i))
+        return apply_2nd_derivative(op, h_inv, v, Index{Upper}(i))
     else
         error("Bounds error") # TODO: Make this more standard
     end
 end
 
 # Wrapper functions for using regular indecies without specifying regions
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
-    return apply(op, h, v, Index{Unknown}(i))
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Int)
+    return apply_2nd_derivative(op, h_inv, v, Index{Unknown}(i))
+end
+export apply_2nd_derivative
+
+# TODO: Dispatch on Index{R}?
+apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h*op.quadratureClosure[i]
+apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h*op.quadratureClosure[N-i+1]
+apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h
+
+# TODO: Avoid branching in inner loops
+function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T
+    r = getregion(i, closuresize(op), N)
+    return apply_quadrature(op, h, v, i, N, r)
+end
+export apply_quadrature
+
+# TODO: Dispatch on Index{R}?
+apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h_inv*op.inverseQuadratureClosure[i]
+apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h_inv*op.inverseQuadratureClosure[N-i+1]
+apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h_inv
+
+# TODO: Avoid branching in inner loops
+function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer) where T
+    r = getregion(i, closuresize(op), N)
+    return apply_inverse_quadrature(op, h_inv, v, i, N, r)
+end
+export apply_inverse_quadrature
+
+function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    apply_stencil(op.eClosure,v,1)
+end
+
+function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    apply_stencil_backwards(op.eClosure,v,length(v))
 end
+export apply_boundary_value_transpose
+
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, N::Integer, i::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[i-1]*v
+end
+
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, N::Integer, i::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[N-i]*v
+end
+export apply_boundary_value
+
+function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    h_inv*apply_stencil(op.dClosure,v,1)
+end
+
+function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    -h_inv*apply_stencil_backwards(op.dClosure,v,length(v))
+end
+
+export apply_normal_derivative_transpose
+
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    h_inv*op.dClosure[i-1]*v
+end
+
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    -h_inv*op.dClosure[N-i]*v
+end
+
+export apply_normal_derivative
--- a/SbpOperators/src/d2.jl	Thu Dec 05 09:44:34 2019 +0100
+++ b/SbpOperators/src/d2.jl	Thu Dec 05 10:48:31 2019 +0100
@@ -1,4 +1,4 @@
-export D2, closuresize, readOperator, apply_e, apply_d, apply_e_T, apply_d_T
+export D2, closuresize, readOperator
 
 @enum Parity begin
     odd = -1
@@ -18,84 +18,3 @@
 function closuresize(D::D2)::Int
     return length(D.quadratureClosure)
 end
-
-# TODO: Dispatch on Index{R}?
-apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h*op.quadratureClosure[i]
-apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h*op.quadratureClosure[N-i+1]
-apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h
-
-# TODO: Avoid branching in inner loops
-function apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer) where T
-    r = getregion(i, closuresize(op), N)
-    return apply_quadrature(op, h, v, i, N, r)
-end
-export apply_quadrature
-
-# TODO: Dispatch on Index{R}?
-apply_inverse_quadrature(op::D2{T}, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h_inv*op.inverseQuadratureClosure[i]
-apply_inverse_quadrature(op::D2{T}, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h_inv*op.inverseQuadratureClosure[N-i+1]
-apply_inverse_quadrature(op::D2{T}, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h_inv
-
-# TODO: Avoid branching in inner loops
-function apply_inverse_quadrature(op::D2{T}, h_inv::Real, v::T, i::Integer, N::Integer) where T
-    r = getregion(i, closuresize(op), N)
-    return apply_inverse_quadrature(op, h_inv, v, i, N, r)
-end
-export apply_inverse_quadrature
-
-function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    apply(op.eClosure,v,1)
-end
-
-function apply_e_T(op::D2, v::AbstractVector, ::Type{Upper})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    apply(flip(op.eClosure),v,length(v))
-end
-
-
-function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    op.eClosure[i-1]*v
-end
-
-function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    op.eClosure[N-i]*v
-end
-
-function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    h_inv*apply(op.dClosure,v,1)
-end
-
-function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    -h_inv*apply(flip(op.dClosure),v,length(v))
-end
-
-function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    h_inv*op.dClosure[i-1]*v
-end
-
-function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    -h_inv*op.dClosure[N-i]*v
-end
--- a/SbpOperators/src/stencil.jl	Thu Dec 05 09:44:34 2019 +0100
+++ b/SbpOperators/src/stencil.jl	Thu Dec 05 10:48:31 2019 +0100
@@ -21,7 +21,7 @@
     return s.weights[1 + i - s.range[1]]
 end
 
-Base.@propagate_inbounds @inline function apply(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
     w = s.weights[1]*v[i + s.range[1]]
     @simd for k ∈ 2:N
         w += s.weights[k]*v[i + s.range[1] + k-1]
@@ -29,7 +29,7 @@
     return w
 end
 
-Base.@propagate_inbounds @inline function apply_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
     w = s.weights[N]*v[i - s.range[2]]
     @simd for k ∈ N-1:-1:1
         w += s.weights[k]*v[i - s.range[1] - k + 1]