Mercurial > repos > public > sbplib_julia

annotate sbpD2.jl @ 91:c0f33eccd527 cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 29 Jan 2019 14:32:28 +0100
parents 8d505e9bc715
children 93df72e2b135
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
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
18 Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::InteriorIndex)
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
19 return apply(op.innerStencil, 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
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
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
22 Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::LowerClosureIndex)
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
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
26 Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::UpperClosureIndex)
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
27 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
28 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
29
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
30 @enum Parity begin
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
31 odd = -1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
32 even = 1
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
33 end
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
34
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
35 struct D2{T,N,M,K} <: ConstantStencilOperator
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
36 quadratureClosure::Vector{T}
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
37 innerStencil::Stencil{T,N}
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
38 closureStencils::NTuple{M, Stencil{T,K}}
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
39 eClosure::Vector{T}
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
40 dClosure::Vector{T}
47
8c4cfa680e5a Use enum for parity
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
41 parity::Parity
1
a286be5f8a32 Add type to sbpD2.jl
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
42 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
43
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
44 function closureSize(D::D2)::Int
43
ef060ab3b035 remove stride and remove some bugs
Ylva Rydin <ylva.rydin@telia.com>
parents: 40
diff changeset
45 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
46 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
47
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
48 function readOperator(D2fn, Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
49 d = readSectionedFile(D2fn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
50 h = readSectionedFile(Hfn)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
51
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
52 # Create inner stencil
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
53 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
54 width = length(innerStencilWeights)
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
55 r = (-div(width,2), div(width,2))
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
56
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
57 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
58
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
59 # Create boundary stencils
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
60 boundarySize = length(d["boundary_stencils"])
68
d485da6e3a77 Make D2 more type stable
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
61 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
62
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
63 for i ∈ 1:boundarySize
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
64 stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
65 width = length(stencilWeights)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
66 r = (1-i,width-i)
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 47
diff changeset
67 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
68 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
69
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
70 d2 = D2(
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
71 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
72 innerStencil,
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
73 closureStencils,
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
74 stringToVector(Float64, d["e"][1]),
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
75 stringToVector(Float64, d["d1"][1]),
34
bb841977d198 Move stencil operator application to its own function
Jonatan Werpers <jonatan@werpers.com>
parents: 24
diff changeset
76 even
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
77 )
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
78
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
79 return d2
24
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
80 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
81
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 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
84 f = open(filename)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
85 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
86 currentKey = ""
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
87
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
88 for ln ∈ eachline(f)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
89 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
90 continue
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
91 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
92
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
93 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
94 if ~haskey(sections, ln)
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
95 sections[ln] = Vector{String}()
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
96 end
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
97 currentKey = ln
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
98 continue
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
55fea1ceb6aa Start implementing reading 1D operator stencils from file into struct
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
101 push!(sections[currentKey], ln)
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 return sections
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
35
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
107 function stringToVector(T::DataType, s::String)
91e662512e9a Fix bugs in sbpD2
Ylva Rydin <ylva.rydin@telia.com>
parents: 24
diff changeset
108 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
109 end