comparison SbpOperators/src/SbpOperators.jl @ 249:7cb4492ccd60 boundary_conditions

Refactor package SbpOperators
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jun 2019 14:18:48 +0200
parents d9e262cb2e8d
children d1004b881da1
comparison
equal deleted inserted replaced
248:05e7bbe0af97 249:7cb4492ccd60
1 module SbpOperators 1 module SbpOperators
2 2
3 using RegionIndices 3 using RegionIndices
4 4
5 include("stencil.jl") 5 include("stencil.jl")
6 6 include("constantstenciloperator.jl")
7 export D2, closureSize, apply, readOperator, apply_e, apply_d, apply_e_T, apply_d_T 7 include("d2.jl")
8 8 include("readoperator.jl")
9 abstract type ConstantStencilOperator end
10
11 # Apply for different regions Lower/Interior/Upper or Unknown region
12 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
13 return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
14 end
15
16 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
17 return @inbounds h*h*apply(op.innerStencil, v, Int(i))
18 end
19
20 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
21 N = length(v)
22 return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
23 end
24
25 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
26 cSize = closureSize(op)
27 N = length(v)
28
29 i = Int(index)
30
31 if 0 < i <= cSize
32 return apply(op, h, v, Index{Lower}(i))
33 elseif cSize < i <= N-cSize
34 return apply(op, h, v, Index{Interior}(i))
35 elseif N-cSize < i <= N
36 return apply(op, h, v, Index{Upper}(i))
37 else
38 error("Bounds error") # TODO: Make this more standard
39 end
40 end
41
42
43 # Wrapper functions for using regular indecies without specifying regions
44 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
45 return apply(op, h, v, Index{Unknown}(i))
46 end
47
48 @enum Parity begin
49 odd = -1
50 even = 1
51 end
52
53 struct D2{T,N,M,K} <: ConstantStencilOperator
54 quadratureClosure::NTuple{M,T}
55 innerStencil::Stencil{T,N}
56 closureStencils::NTuple{M,Stencil{T,K}}
57 eClosure::Stencil{T,M}
58 dClosure::Stencil{T,M}
59 parity::Parity
60 end
61
62 function closureSize(D::D2)::Int
63 return length(D.quadratureClosure)
64 end
65
66 function readOperator(D2fn, Hfn)
67 d = readSectionedFile(D2fn)
68 h = readSectionedFile(Hfn)
69
70 # Create inner stencil
71 innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
72 width = length(innerStencilWeights)
73 r = (-div(width,2), div(width,2))
74
75 innerStencil = Stencil(r, innerStencilWeights)
76
77 # Create boundary stencils
78 boundarySize = length(d["boundary_stencils"])
79 closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
80
81 for i ∈ 1:boundarySize
82 stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
83 width = length(stencilWeights)
84 r = (1-i,width-i)
85 closureStencils = (closureStencils..., Stencil(r, stencilWeights))
86 end
87
88 quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
89 eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
90 dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
91
92 d2 = D2(
93 quadratureClosure,
94 innerStencil,
95 closureStencils,
96 eClosure,
97 dClosure,
98 even
99 )
100
101 return d2
102 end
103
104
105 function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower})
106 @boundscheck if length(v) < closureSize(op)
107 throw(BoundsError())
108 end
109 apply(op.eClosure,v,1)
110 end
111
112 function apply_e_T(op::D2, v::AbstractVector, ::Type{Upper})
113 @boundscheck if length(v) < closureSize(op)
114 throw(BoundsError())
115 end
116 apply(flip(op.eClosure),v,length(v))
117 end
118
119
120 function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Lower})
121 @boundscheck if !(0<length(i) <= N)
122 throw(BoundsError())
123 end
124 op.eClosure[i-1]*v
125 end
126
127 function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Upper})
128 @boundscheck if !(0<length(i) <= N)
129 throw(BoundsError())
130 end
131 op.eClosure[N-i]*v
132 end
133
134 function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
135 @boundscheck if length(v) < closureSize(op)
136 throw(BoundsError())
137 end
138 h_inv*apply(op.dClosure,v,1)
139 end
140
141 function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
142 @boundscheck if length(v) < closureSize(op)
143 throw(BoundsError())
144 end
145 -h_inv*apply(flip(op.dClosure),v,length(v))
146 end
147
148 function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
149 @boundscheck if !(0<length(i) <= N)
150 throw(BoundsError())
151 end
152 h_inv*op.dClosure[i-1]*v
153 end
154
155 function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
156 @boundscheck if !(0<length(i) <= N)
157 throw(BoundsError())
158 end
159 -h_inv*op.dClosure[N-i]*v
160 end
161
162
163 function readSectionedFile(filename)::Dict{String, Vector{String}}
164 f = open(filename)
165 sections = Dict{String, Vector{String}}()
166 currentKey = ""
167
168 for ln ∈ eachline(f)
169 if ln == "" || ln[1] == '#' # Skip comments and empty lines
170 continue
171 end
172
173 if isletter(ln[1]) # Found start of new section
174 if ~haskey(sections, ln)
175 sections[ln] = Vector{String}()
176 end
177 currentKey = ln
178 continue
179 end
180
181 push!(sections[currentKey], ln)
182 end
183
184 return sections
185 end
186
187 function stringToTuple(T::DataType, s::String)
188 return Tuple(stringToVector(T,s))
189 end
190
191 function stringToVector(T::DataType, s::String)
192 return T.(eval.(Meta.parse.(split(s))))
193 end
194
195
196 function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
197 if N >= n
198 return t
199 else
200 return pad_tuple((t..., zero(T)), n)
201 end
202 end
203
204 sbp_operators_path() = (@__DIR__) * "/../operators/"
205 export sbp_operators_path
206 9
207 end # module 10 end # module