Mercurial > repos > public > sbplib_julia

comparison sbpD2.jl @ 134:79699dda29be

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