Line data Source code
1 : # TODO: migrate from FLoops to OhMyThreads?!
2 :
3 3584 : @inline function _diff(F, h, i::Tuple, ::Val{K}) where {K}
4 3584 : (i[K] == 1) && return _diffl(F, h, i, Val(K))
5 2944 : (i[K] == size(F, K)) && return _diffr(F, h, i, Val(K))
6 :
7 2304 : ir = MVector(i...)
8 2304 : ir[K] += 1
9 2304 : il = MVector(i...)
10 2304 : il[K] -= 1
11 :
12 : # NOTE: Need to convert back to SVector or won't optimize!
13 :
14 2304 : (F[SVector(ir)...] - F[SVector(il)...]) / 2h[K]
15 : end
16 :
17 : # left boundary
18 640 : function _diffl(F, h, i::Tuple, ::Val{K}) where {K}
19 640 : ia = MVector(i...)
20 640 : ia[K] += 1
21 640 : ib = MVector(i...)
22 640 : ib[K] += 1
23 :
24 640 : (-3F[SVector(i)...] + 4F[SVector(ia)...] - 1F[SVector(ib)...]) / 2h[K]
25 : end
26 :
27 : # right boundary
28 640 : function _diffr(F, h, i::Tuple, ::Val{K}) where {K}
29 640 : ia = MVector(i...)
30 640 : ia[K] -= 1
31 640 : ib = MVector(i...)
32 640 : ib[K] -= 1
33 :
34 640 : (+3F[SVector(i)...] - 4F[SVector(ia)...] + 1F[SVector(ib)...]) / 2h[K]
35 : end
36 :
37 : #-----------------------------------------------------------------------------
38 :
39 : # Helper that handles edge cases:
40 : #
41 : # eigvals(S) requires a symmetric matrix S and raises an error otherwise.
42 : # This includes the case that an element is NaN (which always implies that
43 : # S != S').
44 : #
45 : # _maximum_eigvals(S) returns NaN if S is not symmetric or any entry is not
46 : # finite (NaN or Inf).
47 : #
48 1536 : function _maximum_eigvals(S)
49 1536 : (S == S') || return eltype(S)(NaN)
50 1536 : maximum(eigvals(S))
51 : end
52 :
53 1024 : function _ftle(Φ, absτ, h::Tuple, i::Int, j::Int)
54 0 : @massert absτ > 0
55 : # Φx = (Φ[i+1, j] - Φ[i-1, j]) / 2h[1]
56 : # Φy = (Φ[i, j+1] - Φ[i, j-1]) / 2h[2]
57 :
58 1920 : Φx = _diff(Φ, h, (i, j), Val(1))
59 1920 : Φy = _diff(Φ, h, (i, j), Val(2))
60 :
61 1024 : ∇Φ = [ Φx Φy ] :: SMatrix
62 :
63 1024 : λmax = _maximum_eigvals(∇Φ'∇Φ)
64 :
65 1024 : log(λmax) / 2absτ
66 : end
67 :
68 512 : function _ftle(Φ, absτ, h::Tuple, i::Int, j::Int, k::Int)
69 0 : @massert absτ > 0
70 :
71 : # Φx = (Φ[i+1, j, k] - Φ[i-1, j, k]) / 2h[1]
72 : # Φy = (Φ[i, j+1, k] - Φ[i, j-1, k]) / 2h[2]
73 : # Φz = (Φ[i, j, k+1] - Φ[i, j, k+1]) / 2h[3]
74 :
75 896 : Φx = _diff(Φ, h, (i, j, k), Val(1))
76 896 : Φy = _diff(Φ, h, (i, j, k), Val(2))
77 896 : Φz = _diff(Φ, h, (i, j, k), Val(3))
78 :
79 512 : ∇Φ = [ Φx Φy Φz] :: SMatrix
80 :
81 1024 : λmax = _maximum_eigvals(∇Φ'∇Φ)
82 :
83 512 : log(λmax) / 2absτ
84 : end
85 :
86 : #-----------------------------------------------------------------------------
87 :
88 : """
89 : value = _ftle(Φ, absτ, h, i, j) # 2d
90 : value = _ftle(Φ, absτ, h, i, j, k) # 3d
91 :
92 : Helper that computes FTLE value for sampled flow map `Φ` and `absτ =
93 : abs(τ) > 0` and grid spacing `h = (hi, hi, ...)` at `i, j, ...`.
94 :
95 : See also [`ftle`](@ref), [`ftle!`](@ref)
96 : """ _ftle
97 :
98 12 : function ftle(p::Problem{N, T, S, Dy, Out},
99 : t0::T, τ::T, extent...) where {N, T, S, Dy, Out}
100 12 : @assert 2 <= N::Int <= 3 "require dimensions N in [2, 3]"
101 12 : sz = length.(extent)
102 12 : Φ = Array{SVector{N, T}, N}(undef, sz...)
103 12 : F = Array{T, N}(undef, sz...)
104 :
105 12 : ftle!(p, Φ, F, t0, τ, extent...)
106 : end
107 :
108 2 : function ftle(p::Problem{N, T, S, Dy, Out},
109 : t0::T, τs, extent...) where {N, T, S, Dy, Out}
110 2 : @assert 2 <= N::Int <= 3 "require dimensions N in [2, 3]"
111 2 : sz = length.(extent)
112 2 : n = length(τs)
113 2 : Φ = Array{SVector{N, T}, N + 1}(undef, n, sz...)
114 2 : F = Array{T, N + 1}(undef, n, sz...)
115 :
116 2 : ftle!(p, Φ, F, t0, τs, extent...)
117 : end
118 :
119 : # 2d single τ
120 16 : function ftle!(p::Problem{2, T, S, Dy, Out},
121 : Φ::AbstractArray{SVector{2, T}, 2}, F::Array{T, 2},
122 : t0::T, τ::T,
123 : x::StepRangeLen{T}, y::StepRangeLen{T};
124 : sample::Bool=true) where {T, S, Dy, Out}
125 8 : sz = (length(x), length(y))
126 8 : @assert size(Φ) == sz
127 8 : @assert size(F) == sz
128 :
129 8 : sample && sampleflowmap!(p, Φ, t0, t0 + τ, x, y)
130 :
131 8 : h = (x[2] - x[1], y[2] - y[1])
132 8 : τ = abs(τ)
133 :
134 8 : @floop for i in 1:sz[1], j in 1:sz[2] # TODO: [benchmark] parallel should be dim k?!
135 512 : F[i, j] = _ftle(Φ, τ, h, i, j)
136 0 : end
137 :
138 8 : F
139 : end
140 :
141 : # 2d multiple τs
142 2 : function ftle!(p::Problem{2, T, S, Dy, Out},
143 : Φ::AbstractArray{SVector{2, T}, 3}, F::AbstractArray{T, 3},
144 : t0::T, τs,
145 : x::StepRangeLen{T}, y::StepRangeLen{T};
146 : sample::Bool=true) where {T, S, Dy, Out}
147 1 : sz = (length(x), length(y))
148 1 : n = length(τs)
149 1 : @assert size(Φ) == (n, sz...)
150 1 : @assert size(F) == (n, sz...)
151 :
152 1 : t1 = t0 .+ τs
153 :
154 1 : sample && sampleflowmap!(p, Φ, t0, t1, x, y)
155 :
156 1 : h = (x[2] - x[1], y[2] - y[1])
157 :
158 1 : @floop for ti in 1:n
159 8 : τ = abs(τs[ti])
160 8 : Φh = @view Φ[ti, :, :]
161 :
162 8 : for i in 1:sz[1], j in 1:sz[2]
163 512 : F[ti, i, j] = _ftle(Φh, τ, h, i, j)
164 568 : end
165 0 : end
166 :
167 1 : F
168 : end
169 :
170 : # 3d single τ
171 8 : function ftle!(p::Problem{3, T, S, Dy, Out},
172 : Φ::Array{SVector{3, T}, 3}, F::Array{T, 3},
173 : t0::T, τ::T,
174 : x::StepRangeLen{T}, y::StepRangeLen{T}, z::StepRangeLen{T};
175 : sample::Bool=true) where {T, S, Dy, Out}
176 4 : sz = (length(x), length(y), length(z))
177 4 : @assert size(Φ) == sz
178 4 : @assert size(F) == sz
179 :
180 4 : sample && sampleflowmap!(p, Φ, t0, t0 + τ, x, y, z)
181 :
182 4 : h = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
183 4 : τ = abs(τ)
184 :
185 4 : @floop for i in 1:sz[1], j in 1:sz[2], k in 1:sz[3] # TODO: [benchmark] parallel should be dim k?!
186 256 : F[i, j, k] = _ftle(Φ, τ, h, i, j, k)
187 0 : end
188 :
189 4 : F
190 : end
191 :
192 : # 3d multiple τs
193 2 : function ftle!(p::Problem{3, T, S, Dy, Out},
194 : Φ::AbstractArray{SVector{3, T}, 4}, F::AbstractArray{T, 4},
195 : t0::T, τs,
196 : x::StepRangeLen{T}, y::StepRangeLen{T}, z::StepRangeLen{T};
197 : sample::Bool=true) where {T, S, Dy, Out}
198 1 : sz = (length(x), length(y),length(z))
199 1 : n = length(τs)
200 1 : @assert size(Φ) == (n, sz...)
201 1 : @assert size(F) == (n, sz...)
202 :
203 1 : t1 = t0 .+ τs
204 :
205 1 : sample && sampleflowmap!(p, Φ, t0, t1, x, y, z)
206 :
207 1 : h = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
208 :
209 1 : @floop for ti in 1:n
210 4 : τ = abs(τs[ti])
211 4 : Φh = @view Φ[ti, :, :, :]
212 :
213 4 : for i in 1:sz[1], j in 1:sz[2], k in 1:sz[3]
214 256 : F[ti, i, j, k] = _ftle(Φh, τ, h, i, j, k)
215 268 : end
216 0 : end
217 :
218 1 : F
219 : end
220 :
221 : """
222 : x = range(x0, x1, length = M)
223 : y = range(y0, y1, length = N)
224 : z = range(z0, z1, length = K)
225 :
226 : F = ftle(problem, t0, τ, x, y[, z])
227 : F = ftle(problem, t0, τs, x, y[, z])
228 :
229 : Sample 2d or 3d FTLE field in `F` for `problem` and integation from
230 : `t0` to `t0 + τ` (or for each `τ ∈ τs`) in domain `x × y[ ×z]`.
231 :
232 : !!! note
233 : Non-finite entries in `Φ` produce undefined FTLE values (`NaN`).
234 :
235 : !!! note
236 : FTLE is computed at the same resolution as the flow map with
237 : adpated rules for central differences at the boundary.
238 :
239 : !!! note
240 : `ftle` uses **multi-threading** ([`mtproblem`](@ref)).
241 :
242 : See also [`ftle!`](@ref)
243 : """ ftle
244 :
245 : """
246 : x = range(x0, x1, length = M)
247 : y = range(y0, y1, length = N)
248 : z = range(z0, z1, length = K)
249 :
250 : ftle!(problem, Φ, F, t0, τ, x, y[, z; sample=true])
251 : ftle!(problem, Φ, F, t0, τs, x, y[, z; sample=true])
252 :
253 : In-place variant of [`ftle`](@ref) that samples flow map in `Φ` and
254 : stores FTLE values in `F`.
255 :
256 : If `sample=false`, the flow map is not sampled and `Φ` is treated as
257 : input.
258 :
259 : !!! note
260 : Non-finite entries in `Φ` produce undefined FTLE values (`NaN`).
261 :
262 : !!! note
263 : Dimensions of `Φ` and `F` must match the lengths `M, N, K`.
264 :
265 : See also [`ftle`](@ref), [`sampleflowmap!`](@ref)
266 : """ ftle!
267 :
268 : # TODO: [test] test FTLE (against DifferentialEquations)
|
