Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/d2.jl @ 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 4ca3794fffef
children ccef055233a2
comparison
equal deleted inserted replaced
259:5571d2c5bf0f 260:f89718833620
5 even = 1 5 even = 1
6 end 6 end
7 7
8 struct D2{T,N,M,K} <: ConstantStencilOperator 8 struct D2{T,N,M,K} <: ConstantStencilOperator
9 quadratureClosure::NTuple{M,T} 9 quadratureClosure::NTuple{M,T}
10 inverseQuadratureClosure::NTuple{M,T}
10 innerStencil::Stencil{T,N} 11 innerStencil::Stencil{T,N}
11 closureStencils::NTuple{M,Stencil{T,K}} 12 closureStencils::NTuple{M,Stencil{T,K}}
12 eClosure::Stencil{T,M} 13 eClosure::Stencil{T,M}
13 dClosure::Stencil{T,M} 14 dClosure::Stencil{T,M}
14 parity::Parity 15 parity::Parity
16 17
17 function closuresize(D::D2)::Int 18 function closuresize(D::D2)::Int
18 return length(D.quadratureClosure) 19 return length(D.quadratureClosure)
19 end 20 end
20 21
22 # TODO: Dispatch on Index{R}?
21 apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h*op.quadratureClosure[i] 23 apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h*op.quadratureClosure[i]
22 apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h*op.quadratureClosure[N-i+1] 24 apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h*op.quadratureClosure[N-i+1]
23 apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h 25 apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h
24 26
27 # TODO: Avoid branching in inner loops
25 function apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer) where T 28 function apply_quadrature(op::D2{T}, h::Real, v::T, i::Integer, N::Integer) where T
26 r = getregion(i, closuresize(op), N) 29 r = getregion(i, closuresize(op), N)
27 return apply_quadrature(op, h, v, i, N, r) 30 return apply_quadrature(op, h, v, i, N, r)
28 end 31 end
29 export apply_quadrature 32 export apply_quadrature
33
34 # TODO: Dispatch on Index{R}?
35 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]
36 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]
37 apply_inverse_quadrature(op::D2{T}, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h_inv
38
39 # TODO: Avoid branching in inner loops
40 function apply_inverse_quadrature(op::D2{T}, h_inv::Real, v::T, i::Integer, N::Integer) where T
41 r = getregion(i, closuresize(op), N)
42 return apply_inverse_quadrature(op, h_inv, v, i, N, r)
43 end
44 export apply_inverse_quadrature
30 45
31 function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower}) 46 function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower})
32 @boundscheck if length(v) < closuresize(op) 47 @boundscheck if length(v) < closuresize(op)
33 throw(BoundsError()) 48 throw(BoundsError())
34 end 49 end