Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/constantstenciloperator.jl @ 271:ecd49ffe0bc8 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Dispatch apply_quadrature and apply_inverse_quadrature on Index-type. Add convenience functions dispatching on integer-indices
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 05 Dec 2019 11:54:56 +0100
parents ccef055233a2
children fe9e8737ddfa
comparison
equal deleted inserted replaced
270:12b738f260a0 271:ecd49ffe0bc8
13 N = length(v) 13 N = length(v)
14 return @inbounds h_inv*h_inv*Int(op.parity)*apply_stencil_backwards(op.closureStencils[N-Int(i)+1], v, Int(i)) 14 return @inbounds h_inv*h_inv*Int(op.parity)*apply_stencil_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
15 end 15 end
16 16
17 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown}) 17 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
18 cSize = closuresize(op)
19 N = length(v) 18 N = length(v)
20 19 r = getregion(Int(index), closuresize(op), N)
21 i = Int(index) 20 i = Index(Int(index), r)
22 21 return apply_2nd_derivative(op, h_inv, v, i)
23 if 0 < i <= cSize
24 return apply_2nd_derivative(op, h_inv, v, Index{Lower}(i))
25 elseif cSize < i <= N-cSize
26 return apply_2nd_derivative(op, h_inv, v, Index{Interior}(i))
27 elseif N-cSize < i <= N
28 return apply_2nd_derivative(op, h_inv, v, Index{Upper}(i))
29 else
30 error("Bounds error") # TODO: Make this more standard
31 end
32 end 22 end
33 23
34 # Wrapper functions for using regular indecies without specifying regions 24 # Wrapper functions for using regular indecies without specifying regions
35 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Int) 25 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Int)
36 return apply_2nd_derivative(op, h_inv, v, Index{Unknown}(i)) 26 return apply_2nd_derivative(op, h_inv, v, Index{Unknown}(i))
37 end 27 end
38 export apply_2nd_derivative 28 export apply_2nd_derivative
39 29
40 # TODO: Dispatch on Index{R}? 30 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Lower}, N::Integer) where T = v*h*op.quadratureClosure[Int(i)]
41 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h*op.quadratureClosure[i] 31 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Upper}, N::Integer) where T = v*h*op.quadratureClosure[N-Int(i)+1]
42 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h*op.quadratureClosure[N-i+1] 32 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h
43 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h
44 33
45 # TODO: Avoid branching in inner loops 34 function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, index::Index{Unknown}, N::Integer) where T
35 r = getregion(Int(index), closuresize(op), N)
36 i = Index(Int(index), r)
37 return apply_quadrature(op, h, v, i, N)
38 end
39
40 # Wrapper functions for using regular indecies without specifying regions
46 function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T 41 function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T
47 r = getregion(i, closuresize(op), N) 42 return apply_quadrature(op, h, v, Index{Unknown}(i), N)
48 return apply_quadrature(op, h, v, i, N, r)
49 end 43 end
50 export apply_quadrature 44 export apply_quadrature
51 45
52 # TODO: Dispatch on Index{R}? 46 # TODO: Evaluate if divisions affect performance
53 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h_inv*op.inverseQuadratureClosure[i] 47 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Lower}, N::Integer) where T = h_inv*v/op.quadratureClosure[Int(i)]
54 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] 48 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Upper}, N::Integer) where T = h_inv*v/op.quadratureClosure[N-Int(i)+1]
55 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h_inv 49 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h_inv
56 50
57 # TODO: Avoid branching in inner loops 51 function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, index::Index{Unknown}, N::Integer) where T
58 function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer) where T 52 r = getregion(Int(index), closuresize(op), N)
59 r = getregion(i, closuresize(op), N) 53 i = Index(Int(index), r)
60 return apply_inverse_quadrature(op, h_inv, v, i, N, r) 54 return apply_inverse_quadrature(op, h_inv, v, i, N)
55 end
56
57 # Wrapper functions for using regular indecies without specifying regions
58 function apply_inverse_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T
59 return apply_inverse_quadrature(op, h, v, Index{Unknown}(i), N)
61 end 60 end
62 export apply_inverse_quadrature 61 export apply_inverse_quadrature
63 62
64 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower}) 63 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
65 @boundscheck if length(v) < closuresize(op) 64 @boundscheck if length(v) < closuresize(op)