Mercurial > repos > public > sbplib_julia
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 |