Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/constantstenciloperator.jl @ 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 396eadb652bd
children ecd49ffe0bc8
comparison
equal deleted inserted replaced
268:f67ce2eb6019 269:ccef055233a2
1 export apply
2
3 abstract type ConstantStencilOperator end 1 abstract type ConstantStencilOperator end
4 2
5 # Apply for different regions Lower/Interior/Upper or Unknown region 3 # Apply for different regions Lower/Interior/Upper or Unknown region
6 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower}) 4 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Lower})
7 return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i)) 5 return @inbounds h_inv*h_inv*apply_stencil(op.closureStencils[Int(i)], v, Int(i))
8 end 6 end
9 7
10 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior}) 8 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Interior})
11 return @inbounds h*h*apply(op.innerStencil, v, Int(i)) 9 return @inbounds h_inv*h_inv*apply_stencil(op.innerStencil, v, Int(i))
12 end 10 end
13 11
14 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper}) 12 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Upper})
15 N = length(v) 13 N = length(v)
16 return @inbounds h*h*Int(op.parity)*apply_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))
17 end 15 end
18 16
19 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown}) 17 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
20 cSize = closuresize(op) 18 cSize = closuresize(op)
21 N = length(v) 19 N = length(v)
22 20
23 i = Int(index) 21 i = Int(index)
24 22
25 if 0 < i <= cSize 23 if 0 < i <= cSize
26 return apply(op, h, v, Index{Lower}(i)) 24 return apply_2nd_derivative(op, h_inv, v, Index{Lower}(i))
27 elseif cSize < i <= N-cSize 25 elseif cSize < i <= N-cSize
28 return apply(op, h, v, Index{Interior}(i)) 26 return apply_2nd_derivative(op, h_inv, v, Index{Interior}(i))
29 elseif N-cSize < i <= N 27 elseif N-cSize < i <= N
30 return apply(op, h, v, Index{Upper}(i)) 28 return apply_2nd_derivative(op, h_inv, v, Index{Upper}(i))
31 else 29 else
32 error("Bounds error") # TODO: Make this more standard 30 error("Bounds error") # TODO: Make this more standard
33 end 31 end
34 end 32 end
35 33
36 # Wrapper functions for using regular indecies without specifying regions 34 # Wrapper functions for using regular indecies without specifying regions
37 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int) 35 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Int)
38 return apply(op, h, v, Index{Unknown}(i)) 36 return apply_2nd_derivative(op, h_inv, v, Index{Unknown}(i))
39 end 37 end
38 export apply_2nd_derivative
39
40 # TODO: Dispatch on Index{R}?
41 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Lower}) where T = v*h*op.quadratureClosure[i]
42 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Upper}) where T = v*h*op.quadratureClosure[N-i+1]
43 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h
44
45 # TODO: Avoid branching in inner loops
46 function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T
47 r = getregion(i, closuresize(op), N)
48 return apply_quadrature(op, h, v, i, N, r)
49 end
50 export apply_quadrature
51
52 # TODO: Dispatch on Index{R}?
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]
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]
55 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer, ::Type{Interior}) where T = v*h_inv
56
57 # TODO: Avoid branching in inner loops
58 function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Integer, N::Integer) where T
59 r = getregion(i, closuresize(op), N)
60 return apply_inverse_quadrature(op, h_inv, v, i, N, r)
61 end
62 export apply_inverse_quadrature
63
64 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
65 @boundscheck if length(v) < closuresize(op)
66 throw(BoundsError())
67 end
68 apply_stencil(op.eClosure,v,1)
69 end
70
71 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper})
72 @boundscheck if length(v) < closuresize(op)
73 throw(BoundsError())
74 end
75 apply_stencil_backwards(op.eClosure,v,length(v))
76 end
77 export apply_boundary_value_transpose
78
79 function apply_boundary_value(op::ConstantStencilOperator, v::Number, N::Integer, i::Integer, ::Type{Lower})
80 @boundscheck if !(0<length(i) <= N)
81 throw(BoundsError())
82 end
83 op.eClosure[i-1]*v
84 end
85
86 function apply_boundary_value(op::ConstantStencilOperator, v::Number, N::Integer, i::Integer, ::Type{Upper})
87 @boundscheck if !(0<length(i) <= N)
88 throw(BoundsError())
89 end
90 op.eClosure[N-i]*v
91 end
92 export apply_boundary_value
93
94 function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Lower})
95 @boundscheck if length(v) < closuresize(op)
96 throw(BoundsError())
97 end
98 h_inv*apply_stencil(op.dClosure,v,1)
99 end
100
101 function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Upper})
102 @boundscheck if length(v) < closuresize(op)
103 throw(BoundsError())
104 end
105 -h_inv*apply_stencil_backwards(op.dClosure,v,length(v))
106 end
107
108 export apply_normal_derivative_transpose
109
110 function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
111 @boundscheck if !(0<length(i) <= N)
112 throw(BoundsError())
113 end
114 h_inv*op.dClosure[i-1]*v
115 end
116
117 function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
118 @boundscheck if !(0<length(i) <= N)
119 throw(BoundsError())
120 end
121 -h_inv*op.dClosure[N-i]*v
122 end
123
124 export apply_normal_derivative