Mercurial > repos > public > sbplib_julia

annotate sbpD2.jl @ 97:8324c82c2dfb cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make D2 support sbp.Index and specification of region
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 05 Feb 2019 22:30:46 +0100
parents 8d505e9bc715
children 50273f745f05
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
97
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
3 # Apply for different regions Lower/Interior/Upper
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
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
17 # Wrapper functions for using regular indecies without specifying regions
78
fbf7398f8154 Inline and inbounds everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 68
diff changeset
18 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
19 cSize = closureSize(op)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
20 N = length(v)
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
21
97
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
22 if 0 < i <= cSize
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
23 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
24 elseif cSize < i <= N-cSize
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
25 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
26 elseif N-cSize < i <= N
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
27 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
28 else
97
8324c82c2dfb Make D2 support sbp.Index and specification of region
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
29 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
30 end
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
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
33 @enum Parity begin
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
34 odd = -1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
35 even = 1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
36 end
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
37
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
38 struct D2{T,N,M,K} <: ConstantStencilOperator
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
39 quadratureClosure::Vector{T}
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
40 innerStencil::Stencil{T,N}
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
41 closureStencils::NTuple{M, Stencil{T,K}}
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
42 eClosure::Vector{T}
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
43 dClosure::Vector{T}
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
44 parity::Parity
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
45 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
46
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
47 function closureSize(D::D2)::Int
43
ef060ab3b035 remove stride and remove some bugs
Ylva Rydin <ylva.rydin@telia.com>
parents: 40
diff changeset
48 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
49 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
50
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
51 function readOperator(D2fn, Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
52 d = readSectionedFile(D2fn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
53 h = readSectionedFile(Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
54
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
55 # Create inner stencil
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
56 innerStencilWeights = stringToVector(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
57 width = length(innerStencilWeights)
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
58 r = (-div(width,2), div(width,2))
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
59
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
60 innerStencil = Stencil(r, Tuple(innerStencilWeights))
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
61
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
62 # Create boundary stencils
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
63 boundarySize = length(d["boundary_stencils"])
68
d485da6e3a77 Make D2 more type stable
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
64 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
65
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
66 for i ∈ 1:boundarySize
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
67 stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
68 width = length(stencilWeights)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
69 r = (1-i,width-i)
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
70 closureStencils = (closureStencils..., Stencil(r, Tuple(stencilWeights)))
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
71 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
72
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
73 d2 = D2(
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
74 stringToVector(Float64, h["closure"][1]),
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
75 innerStencil,
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
76 closureStencils,
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
77 stringToVector(Float64, d["e"][1]),
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
78 stringToVector(Float64, d["d1"][1]),
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
79 even
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
80 )
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
81
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
82 return d2
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
83 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
84
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
85
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
86 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
87 f = open(filename)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
88 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
89 currentKey = ""
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 for ln ∈ eachline(f)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
92 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
93 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
94 end
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 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
97 if ~haskey(sections, ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
98 sections[ln] = Vector{String}()
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
99 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
100 currentKey = ln
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
101 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
102 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
103
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
104 push!(sections[currentKey], ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
105 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
106
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
107 return sections
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
108 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
109
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
110 function stringToVector(T::DataType, s::String)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
111 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
112 end