Mercurial > repos > public > sbplib_julia

changeset 260:f89718833620 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Store the inverse quadrature closure for D2. Implement the stencil application for the inverse quadrature
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 26 Nov 2019 08:18:23 -0800
parents 5571d2c5bf0f
children 01017d2b46b0
files SbpOperators/src/d2.jl SbpOperators/src/readoperator.jl
diffstat 2 files changed, 17 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/SbpOperators/src/d2.jl	Fri Jun 28 14:17:13 2019 +0200
+++ b/SbpOperators/src/d2.jl	Tue Nov 26 08:18:23 2019 -0800
@@ -7,6 +7,7 @@
 
 struct D2{T,N,M,K} <: ConstantStencilOperator
     quadratureClosure::NTuple{M,T}
+    inverseQuadratureClosure::NTuple{M,T}
     innerStencil::Stencil{T,N}
     closureStencils::NTuple{M,Stencil{T,K}}
     eClosure::Stencil{T,M}
@@ -18,16 +19,30 @@
     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())
--- a/SbpOperators/src/readoperator.jl	Fri Jun 28 14:17:13 2019 +0200
+++ b/SbpOperators/src/readoperator.jl	Tue Nov 26 08:18:23 2019 -0800
@@ -21,11 +21,13 @@
     end
 
     quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
+    inverseQuadratureClosure = 1.0 ./ quadratureClosure
     eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
     dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
 
     d2 = D2(
         quadratureClosure,
+        inverseQuadratureClosure,
         innerStencil,
         closureStencils,
         eClosure,