Mercurial > repos > public > sbplib_julia

annotate sbpD2.jl @ 133:c4cc2a1a68b9 stencil_index

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Close branch
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 21 Feb 2019 16:20:55 +0100
parents 84b1ad5a3755
children
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
78
fbf7398f8154 Inline and inbounds everything
Ylva Rydin <ylva.rydin@telia.com>
parents: 68
diff changeset
3 @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
4 cSize = closureSize(op)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
5 N = length(v)
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
6
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
7 if i ∈ range(1; length=cSize)
72
4640839b1616 Make apply operator more efficient using @inbounds and @inline
Jonatan Werpers <jonatan@werpers.com>
parents: 68
diff changeset
8 @inbounds uᵢ = apply(op.closureStencils[i], v, i)/h^2
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
9 elseif i ∈ range(N - cSize+1, length=cSize)
72
4640839b1616 Make apply operator more efficient using @inbounds and @inline
Jonatan Werpers <jonatan@werpers.com>
parents: 68
diff changeset
10 @inbounds uᵢ = Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
11 else
72
4640839b1616 Make apply operator more efficient using @inbounds and @inline
Jonatan Werpers <jonatan@werpers.com>
parents: 68
diff changeset
12 @inbounds uᵢ = apply(op.innerStencil, v, i)/h^2
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
13 end
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
14
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 47
diff changeset
15 return uᵢ
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
16 end
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
17
94
84b1ad5a3755 Made everything work(?) but also go really slow. And also not type-stable.
Ylva Rydin <ylva.rydin@telia.com>
parents: 93
diff changeset
18 Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::StencilIndex{Interior})
93
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 91
diff changeset
19 return apply(op.innerStencil, v, i)/h^2
91
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
20 end
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
21
94
84b1ad5a3755 Made everything work(?) but also go really slow. And also not type-stable.
Ylva Rydin <ylva.rydin@telia.com>
parents: 93
diff changeset
22 Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::StencilIndex{Lower})
91
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
23 return apply(op.closureStencils[i], v, i)/h^2
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
24 end
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
25
94
84b1ad5a3755 Made everything work(?) but also go really slow. And also not type-stable.
Ylva Rydin <ylva.rydin@telia.com>
parents: 93
diff changeset
26 Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::StencilIndex{Upper})
93
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 91
diff changeset
27 N = length(v)
91
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
28 return Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2 #TODO: Write an applybackwards instead?
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
29 end
c0f33eccd527 Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
30
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
31 @enum Parity begin
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
32 odd = -1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
33 even = 1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
34 end
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
35
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
36 struct D2{T,N,M,K} <: ConstantStencilOperator
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
37 quadratureClosure::Vector{T}
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
38 innerStencil::Stencil{T,N}
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
39 closureStencils::NTuple{M, Stencil{T,K}}
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
40 eClosure::Vector{T}
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
41 dClosure::Vector{T}
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
42 parity::Parity
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
43 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
44
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
45 function closureSize(D::D2)::Int
43
ef060ab3b035 remove stride and remove some bugs
Ylva Rydin <ylva.rydin@telia.com>
parents: 40
diff changeset
46 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
47 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
48
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
49 function readOperator(D2fn, Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
50 d = readSectionedFile(D2fn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
51 h = readSectionedFile(Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
52
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
53 # Create inner stencil
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
54 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
55 width = length(innerStencilWeights)
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
56 r = (-div(width,2), div(width,2))
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
57
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
58 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
59
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
60 # Create boundary stencils
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
61 boundarySize = length(d["boundary_stencils"])
68
d485da6e3a77 Make D2 more type stable
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
62 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
63
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
64 for i ∈ 1:boundarySize
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
65 stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
66 width = length(stencilWeights)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
67 r = (1-i,width-i)
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
68 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
69 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
70
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
71 d2 = D2(
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
72 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
73 innerStencil,
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
74 closureStencils,
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
75 stringToVector(Float64, d["e"][1]),
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
76 stringToVector(Float64, d["d1"][1]),
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
77 even
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
78 )
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
79
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
80 return d2
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
81 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
82
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
83
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
84 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
85 f = open(filename)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
86 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
87 currentKey = ""
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
88
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
89 for ln ∈ eachline(f)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
90 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
91 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
92 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
93
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
94 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
95 if ~haskey(sections, ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
96 sections[ln] = Vector{String}()
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
97 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
98 currentKey = ln
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
99 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
100 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
101
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
102 push!(sections[currentKey], ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
103 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
104
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
105 return sections
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
106 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
107
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
108 function stringToVector(T::DataType, s::String)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
109 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
110 end