Mercurial > repos > public > sbplib_julia
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 |
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 | 39 @enum Parity begin |
40 odd = -1 | |
41 even = 1 | |
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 | 50 parity::Parity |
1 | 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 | 64 r = (-div(width,2), div(width,2)) |
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 | 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 | 74 width = length(stencilWeights) |
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 | 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 | 124 function stringToVector(T::DataType, s::String) |
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 |