Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/stencil.jl @ 1032:11767fbb29f4 feature/dissipation_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add padding functions for stencils
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 22 Mar 2022 11:18:57 +0100
parents 4433be383840
children 14cb97284373
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
671
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
1 export CenteredStencil
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
2
806
9fc6d38da03f Remove type requirement for stencil weights
Jonatan Werpers <jonatan@werpers.com>
parents: 672
diff changeset
3 struct Stencil{T,N}
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 69
diff changeset
4 range::Tuple{Int,Int}
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 69
diff changeset
5 weights::NTuple{N,T}
126
66c239678a21 Add Assertion in stencil constructor
Ylva Rydin <ylva.rydin@telia.com>
parents: 85
diff changeset
6
806
9fc6d38da03f Remove type requirement for stencil weights
Jonatan Werpers <jonatan@werpers.com>
parents: 672
diff changeset
7 function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N}
126
66c239678a21 Add Assertion in stencil constructor
Ylva Rydin <ylva.rydin@telia.com>
parents: 85
diff changeset
8 @assert range[2]-range[1]+1 == N
66c239678a21 Add Assertion in stencil constructor
Ylva Rydin <ylva.rydin@telia.com>
parents: 85
diff changeset
9 new{T,N}(range,weights)
66c239678a21 Add Assertion in stencil constructor
Ylva Rydin <ylva.rydin@telia.com>
parents: 85
diff changeset
10 end
8
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 end
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12
584
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
13 """
595
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
14 Stencil(weights::NTuple; center::Int)
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
15
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
16 Create a stencil with the given weights with element `center` as the center of the stencil.
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
17 """
807
1fc6c4daec44 Add type parameter to Vararg to make better error messages when Stencil is called with incompatible weights
Jonatan Werpers <jonatan@werpers.com>
parents: 806
diff changeset
18 function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error
595
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
19 N = length(weights)
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
20 range = (1, N) .- center
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
21
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
22 return Stencil(range, weights)
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
23 end
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
24
826
4433be383840 Add stencil constructor to change the type of the weights along with a convert method
Jonatan Werpers <jonatan@werpers.com>
parents: 807
diff changeset
25 function Stencil{T}(s::Stencil) where T
4433be383840 Add stencil constructor to change the type of the weights along with a convert method
Jonatan Werpers <jonatan@werpers.com>
parents: 807
diff changeset
26 return Stencil(s.range, T.(s.weights))
4433be383840 Add stencil constructor to change the type of the weights along with a convert method
Jonatan Werpers <jonatan@werpers.com>
parents: 807
diff changeset
27 end
4433be383840 Add stencil constructor to change the type of the weights along with a convert method
Jonatan Werpers <jonatan@werpers.com>
parents: 807
diff changeset
28
4433be383840 Add stencil constructor to change the type of the weights along with a convert method
Jonatan Werpers <jonatan@werpers.com>
parents: 807
diff changeset
29 Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil)
4433be383840 Add stencil constructor to change the type of the weights along with a convert method
Jonatan Werpers <jonatan@werpers.com>
parents: 807
diff changeset
30
671
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
31 function CenteredStencil(weights::Vararg)
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
32 if iseven(length(weights))
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
33 throw(ArgumentError("a centered stencil must have an odd number of weights."))
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
34 end
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
35
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
36 r = length(weights) ÷ 2
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
37
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
38 return Stencil((-r, r), weights)
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
39 end
e14627e79a54 Add stencil constructor for centered stencils and change from tuple to vararg in stencil constructor taking cneter
Jonatan Werpers <jonatan@werpers.com>
parents: 595
diff changeset
40
672
59a81254fefc Formatting
Jonatan Werpers <jonatan@werpers.com>
parents: 671
diff changeset
41
595
03ef4d4740ab Add a constructor for Stencil where you can specify the center of the stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 584
diff changeset
42 """
584
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
43 scale(s::Stencil, a)
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
44
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
45 Scale the weights of the stencil `s` with `a` and return a new stencil.
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
46 """
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
47 function scale(s::Stencil, a)
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
48 return Stencil(s.range, a.*s.weights)
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
49 end
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
50
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
51 Base.eltype(::Stencil{T}) where T = T
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
52
8
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 function flip(s::Stencil)
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 range = (-s.range[2], -s.range[1])
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 69
diff changeset
55 return Stencil(range, reverse(s.weights))
8
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 end
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 # Provides index into the Stencil based on offset for the root element
128
7c0b9bb7ab4d Improve stencil application code to make it more friendly to compiler optimizations
Jonatan Werpers <jonatan@werpers.com>
parents: 122
diff changeset
59 @inline function Base.getindex(s::Stencil, i::Int)
71
18d0d794d3bb Make stencils respond to @ inbounds
Jonatan Werpers <jonatan@werpers.com>
parents: 67
diff changeset
60 @boundscheck if i < s.range[1] || s.range[2] < i
584
4aa7fe13a984 Add scale() and eltype() methods for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
61 return zero(eltype(s))
8
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 end
71
18d0d794d3bb Make stencils respond to @ inbounds
Jonatan Werpers <jonatan@werpers.com>
parents: 67
diff changeset
63 return s.weights[1 + i - s.range[1]]
8
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 end
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65
269
ccef055233a2 Move general methods for D2 to ConstantStencilOperator and increase verbosity of method names for clarity
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 219
diff changeset
66 Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
131
8569c637d923 Enable compiler loop-unrolling in apply_backwards
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
67 w = s.weights[1]*v[i + s.range[1]]
129
1aaeb46ba5f4 Improve efficiency of apply by the following:
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 128
diff changeset
68 @simd for k ∈ 2:N
131
8569c637d923 Enable compiler loop-unrolling in apply_backwards
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
69 w += s.weights[k]*v[i + s.range[1] + k-1]
8
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 end
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 return w
433008d3b7d3 Move stencil to its own file
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 end
122
6c6979ff17f4 Introduce and use apply_backwards for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
73
269
ccef055233a2 Move general methods for D2 to ConstantStencilOperator and increase verbosity of method names for clarity
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 219
diff changeset
74 Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
131
8569c637d923 Enable compiler loop-unrolling in apply_backwards
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
75 w = s.weights[N]*v[i - s.range[2]]
8569c637d923 Enable compiler loop-unrolling in apply_backwards
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
76 @simd for k ∈ N-1:-1:1
8569c637d923 Enable compiler loop-unrolling in apply_backwards
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
77 w += s.weights[k]*v[i - s.range[1] - k + 1]
122
6c6979ff17f4 Introduce and use apply_backwards for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
78 end
6c6979ff17f4 Introduce and use apply_backwards for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
79 return w
6c6979ff17f4 Introduce and use apply_backwards for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
80 end
1032
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
81
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
82
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
83 function left_pad(s::Stencil, N)
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
84 weights = LazyTensors.left_pad_tuple(s.weights, zero(eltype(s)), N)
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
85 range = (s.range[1] - (N - length(s.weights)) ,s.range[2])
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
86
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
87 return Stencil(range, weights)
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
88 end
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
89
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
90 function right_pad(s::Stencil, N)
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
91 weights = LazyTensors.right_pad_tuple(s.weights, zero(eltype(s)), N)
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
92 range = (s.range[1], s.range[2] + (N - length(s.weights)))
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
93
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
94 return Stencil(range, weights)
11767fbb29f4 Add padding functions for stencils
Jonatan Werpers <jonatan@werpers.com>
parents: 826
diff changeset
95 end