annotate sbpD2.jl @ 109:b3fbef345810 cell_based_test

Introduce apply! for given indices
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 08 Feb 2019 13:08:07 +0100
parents 19031733bbbf
children 6c6979ff17f4 f01b70b81e95
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
1 abstract type ConstantStencilOperator end
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
2
98
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
3 # Apply for different regions Lower/Interior/Upper or Unknown region
97
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
4 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
5 return @inbounds apply(op.closureStencils[Int(i)], v, Int(i))/h^2
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
6 end
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
7
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
8 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
9 return @inbounds apply(op.innerStencil, v, Int(i))/h^2
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
10 end
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
11
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
12 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
13 N = length(v)
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
14 return @inbounds Int(op.parity)*apply(flip(op.closureStencils[N-Int(i)+1]), v, Int(i))/h^2
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
15 end
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
16
98
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
17 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
18 cSize = closureSize(op)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
19 N = length(v)
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
20
98
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
21 i = Int(index)
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
22
97
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
23 if 0 < i <= cSize
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
24 return apply(op, h, v, Index{Lower}(i))
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
25 elseif cSize < i <= N-cSize
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
26 return apply(op, h, v, Index{Interior}(i))
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
27 elseif N-cSize < i <= N
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
28 return apply(op, h, v, Index{Upper}(i))
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
29 else
97
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
30 error("Bounds error") # TODO: Make this more standard
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
31 end
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
32 end
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
33
98
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
34 # Wrapper functions for using regular indecies without specifying regions
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
35 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
36 return apply(op, h, v, Index{Unknown}(i))
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
37 end
50273f745f05 Add "Unknown" region and implement D2 for it
Jonatan Werpers <jonatan@werpers.com>
parents: 97
diff changeset
38
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
39 @enum Parity begin
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
40 odd = -1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
41 even = 1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
42 end
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
43
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
44 struct D2{T,N,M,K} <: ConstantStencilOperator
101
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
45 quadratureClosure::NTuple{M,T}
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
46 innerStencil::Stencil{T,N}
101
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
47 closureStencils::NTuple{M,Stencil{T,K}}
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
48 eClosure::NTuple{M,T}
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
49 dClosure::NTuple{M,T}
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
50 parity::Parity
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
51 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
52
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
53 function closureSize(D::D2)::Int
43
ef060ab3b035 remove stride and remove some bugs
Ylva Rydin <ylva.rydin@telia.com>
parents: 40
diff changeset
54 return length(D.quadratureClosure)
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
55 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
56
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
57 function readOperator(D2fn, Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
58 d = readSectionedFile(D2fn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
59 h = readSectionedFile(Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
60
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
61 # Create inner stencil
100
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
62 innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
63 width = length(innerStencilWeights)
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
64 r = (-div(width,2), div(width,2))
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
65
100
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
66 innerStencil = Stencil(r, innerStencilWeights)
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
67
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
68 # Create boundary stencils
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
69 boundarySize = length(d["boundary_stencils"])
68
d485da6e3a77 Make D2 more type stable
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
70 closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
71
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
72 for i ∈ 1:boundarySize
100
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
73 stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
74 width = length(stencilWeights)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
75 r = (1-i,width-i)
100
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
76 closureStencils = (closureStencils..., Stencil(r, stencilWeights))
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
77 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
78
101
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
79 quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
80 eClosure = pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize)
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
81 dClosure = pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize)
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
82
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
83 d2 = D2(
101
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
84 quadratureClosure,
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
85 innerStencil,
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
86 closureStencils,
101
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
87 eClosure,
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
88 dClosure,
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
89 even
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
90 )
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
91
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
92 return d2
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
93 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
94
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
95
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
96 function readSectionedFile(filename)::Dict{String, Vector{String}}
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
97 f = open(filename)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
98 sections = Dict{String, Vector{String}}()
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
99 currentKey = ""
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
100
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
101 for ln ∈ eachline(f)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
102 if ln == "" || ln[1] == '#' # Skip comments and empty lines
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
103 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
104 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
105
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
106 if isletter(ln[1]) # Found start of new section
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
107 if ~haskey(sections, ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
108 sections[ln] = Vector{String}()
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
109 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
110 currentKey = ln
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
111 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
112 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
113
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
114 push!(sections[currentKey], ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
115 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
116
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
117 return sections
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
118 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
119
100
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
120 function stringToTuple(T::DataType, s::String)
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
121 return Tuple(stringToVector(T,s))
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
122 end
49796ca2dfa0 Introduce and use the function stringToTuple in sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 98
diff changeset
123
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
124 function stringToVector(T::DataType, s::String)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
125 return T.(eval.(Meta.parse.(split(s))))
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
126 end
101
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
127
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
128
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
129 function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
130 if N >= n
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
131 return t
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
132 else
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
133 return pad_tuple((t..., zero(T)), n)
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
134 end
19031733bbbf Convert D2 to only use tuples
Jonatan Werpers <jonatan@werpers.com>
parents: 100
diff changeset
135 end