Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/constantstenciloperator.jl @ 291:0f94dc29c4bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in branch boundary_conditions
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 22 Jun 2020 21:43:05 +0200
parents fe9e8737ddfa
children 36206445018b
comparison
equal deleted inserted replaced
231:fbabfd4e8f20 291:0f94dc29c4bf
1 abstract type ConstantStencilOperator end
2
3 # Apply for different regions Lower/Interior/Upper or Unknown region
4 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Lower})
5 return @inbounds h_inv*h_inv*apply_stencil(op.closureStencils[Int(i)], v, Int(i))
6 end
7
8 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Interior})
9 return @inbounds h_inv*h_inv*apply_stencil(op.innerStencil, v, Int(i))
10 end
11
12 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Upper})
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))
15 end
16
17 @inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
18 N = length(v)
19 r = getregion(Int(index), closuresize(op), N)
20 i = Index(Int(index), r)
21 return apply_2nd_derivative(op, h_inv, v, i)
22 end
23 export apply_2nd_derivative
24
25 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Lower}, N::Integer) where T = v*h*op.quadratureClosure[Int(i)]
26 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Upper}, N::Integer) where T = v*h*op.quadratureClosure[N-Int(i)+1]
27 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h
28
29 function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, index::Index{Unknown}, N::Integer) where T
30 r = getregion(Int(index), closuresize(op), N)
31 i = Index(Int(index), r)
32 return apply_quadrature(op, h, v, i, N)
33 end
34 export apply_quadrature
35
36 # TODO: Evaluate if divisions affect performance
37 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Lower}, N::Integer) where T = h_inv*v/op.quadratureClosure[Int(i)]
38 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]
39 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h_inv
40
41 function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, index::Index{Unknown}, N::Integer) where T
42 r = getregion(Int(index), closuresize(op), N)
43 i = Index(Int(index), r)
44 return apply_inverse_quadrature(op, h_inv, v, i, N)
45 end
46
47 export apply_inverse_quadrature
48
49 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
50 @boundscheck if length(v) < closuresize(op)
51 throw(BoundsError())
52 end
53 apply_stencil(op.eClosure,v,1)
54 end
55
56 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper})
57 @boundscheck if length(v) < closuresize(op)
58 throw(BoundsError())
59 end
60 apply_stencil_backwards(op.eClosure,v,length(v))
61 end
62 export apply_boundary_value_transpose
63
64 function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Lower})
65 @boundscheck if !(0<length(Int(i)) <= N)
66 throw(BoundsError())
67 end
68 op.eClosure[Int(i)-1]*v
69 end
70
71 function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Upper})
72 @boundscheck if !(0<length(Int(i)) <= N)
73 throw(BoundsError())
74 end
75 op.eClosure[N-Int(i)]*v
76 end
77 export apply_boundary_value
78
79 function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Lower})
80 @boundscheck if length(v) < closuresize(op)
81 throw(BoundsError())
82 end
83 h_inv*apply_stencil(op.dClosure,v,1)
84 end
85
86 function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Upper})
87 @boundscheck if length(v) < closuresize(op)
88 throw(BoundsError())
89 end
90 -h_inv*apply_stencil_backwards(op.dClosure,v,length(v))
91 end
92
93 export apply_normal_derivative_transpose
94
95 function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Lower})
96 @boundscheck if !(0<length(Int(i)) <= N)
97 throw(BoundsError())
98 end
99 h_inv*op.dClosure[Int(i)-1]*v
100 end
101
102 function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Upper})
103 @boundscheck if !(0<length(Int(i)) <= N)
104 throw(BoundsError())
105 end
106 -h_inv*op.dClosure[N-Int(i)]*v
107 end
108
109 export apply_normal_derivative