Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/SbpOperators.jl @ 219:69a6049e14d9 package_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create package SbpOperators
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2019 13:12:38 +0200
parents sbpD2.jl@6ba2238a9687
children 235f0a771c8f
comparison
equal deleted inserted replaced
218:03375aa30edd 219:69a6049e14d9
1 module SbpOperators
2
3 using RegionIndices
4
5 include("stencil.jl")
6
7 abstract type ConstantStencilOperator end
8
9 # Apply for different regions Lower/Interior/Upper or Unknown region
10 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
11 return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
12 end
13
14 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
15 return @inbounds h*h*apply(op.innerStencil, v, Int(i))
16 end
17
18 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
19 N = length(v)
20 return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
21 end
22
23 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
24 cSize = closureSize(op)
25 N = length(v)
26
27 i = Int(index)
28
29 if 0 < i <= cSize
30 return apply(op, h, v, Index{Lower}(i))
31 elseif cSize < i <= N-cSize
32 return apply(op, h, v, Index{Interior}(i))
33 elseif N-cSize < i <= N
34 return apply(op, h, v, Index{Upper}(i))
35 else
36 error("Bounds error") # TODO: Make this more standard
37 end
38 end
39
40
41 # Wrapper functions for using regular indecies without specifying regions
42 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
43 return apply(op, h, v, Index{Unknown}(i))
44 end
45
46 @enum Parity begin
47 odd = -1
48 even = 1
49 end
50
51 struct D2{T,N,M,K} <: ConstantStencilOperator
52 quadratureClosure::NTuple{M,T}
53 innerStencil::Stencil{T,N}
54 closureStencils::NTuple{M,Stencil{T,K}}
55 eClosure::Stencil{T,M}
56 dClosure::Stencil{T,M}
57 parity::Parity
58 end
59
60 function closureSize(D::D2)::Int
61 return length(D.quadratureClosure)
62 end
63
64 function readOperator(D2fn, Hfn)
65 d = readSectionedFile(D2fn)
66 h = readSectionedFile(Hfn)
67
68 # Create inner stencil
69 innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
70 width = length(innerStencilWeights)
71 r = (-div(width,2), div(width,2))
72
73 innerStencil = Stencil(r, innerStencilWeights)
74
75 # Create boundary stencils
76 boundarySize = length(d["boundary_stencils"])
77 closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
78
79 for i ∈ 1:boundarySize
80 stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
81 width = length(stencilWeights)
82 r = (1-i,width-i)
83 closureStencils = (closureStencils..., Stencil(r, stencilWeights))
84 end
85
86 quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
87 eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
88 dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
89
90 d2 = D2(
91 quadratureClosure,
92 innerStencil,
93 closureStencils,
94 eClosure,
95 dClosure,
96 even
97 )
98
99 return d2
100 end
101
102
103 function apply_e(op::D2, v::AbstractVector, ::Type{Lower})
104 apply(op.eClosure,v,1)
105 end
106
107 function apply_e(op::D2, v::AbstractVector, ::Type{Upper})
108 apply(flip(op.eClosure),v,length(v))
109 end
110
111
112 function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
113 -h_inv*apply(op.dClosure,v,1)
114 end
115
116 function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
117 -h_inv*apply(flip(op.dClosure),v,length(v))
118 end
119
120 function readSectionedFile(filename)::Dict{String, Vector{String}}
121 f = open(filename)
122 sections = Dict{String, Vector{String}}()
123 currentKey = ""
124
125 for ln ∈ eachline(f)
126 if ln == "" || ln[1] == '#' # Skip comments and empty lines
127 continue
128 end
129
130 if isletter(ln[1]) # Found start of new section
131 if ~haskey(sections, ln)
132 sections[ln] = Vector{String}()
133 end
134 currentKey = ln
135 continue
136 end
137
138 push!(sections[currentKey], ln)
139 end
140
141 return sections
142 end
143
144 function stringToTuple(T::DataType, s::String)
145 return Tuple(stringToVector(T,s))
146 end
147
148 function stringToVector(T::DataType, s::String)
149 return T.(eval.(Meta.parse.(split(s))))
150 end
151
152
153 function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
154 if N >= n
155 return t
156 else
157 return pad_tuple((t..., zero(T)), n)
158 end
159 end
160
161 end # module