Line data Source code
1 : """
2 : _patch!(solver, tend, y, dy)
3 :
4 : Helper that sets partial solver state returns modified `solver`. The
5 : patched solver will *output* constant `y`, `dy` from `solver.t` until
6 : `tend`.
7 :
8 : !!! note
9 : The sole purpose of `_patch!` is changing the solver state *after
10 : the final step* and before a potential final output.
11 :
12 : Example: If the solver stopped early, i.e., `!isdone(solver)` (e.g.,
13 : due to `OutOfDomain`), you may still want to have output between
14 : `solver.t` and `solver.tfinal`, e.g., undefined values or the last
15 : state reached. You can use `_patch!` before calling `ooutput(solver)`
16 : to achieve this.
17 :
18 : !!! warning
19 : If any step is taken after `_patch!`, the solver state is
20 : *undefined*.
21 : """
22 8 : function _patch!(solver::Solver{N, T, S}, t, y, dy) where {N, T, S}
23 8 : solver.tnew = t
24 8 : solver.y = solver.ynew = y
25 8 : solver.dy = solver.dynew = dy
26 :
27 8 : solver
28 : end
29 :
30 : # TODO: check need for type annotations (field)
31 :
32 1732 : function _solvefill!(p::Problem{N, T, S, Dy, Out0},
33 : t0::T, t1::T, y0::SVector{N, T},
34 : out::Out) where {N, T, S, Dy, Out0, Out}
35 1732 : rv = solve!(p, t0, t1, y0, out)
36 :
37 1732 : s = getsolver(p)
38 1732 : if !isdone(s)
39 8 : yundef = filled(s.y, NaN, first(pos(s)); fill=NaN)
40 8 : dyundef = filled(s.y; fill=0) # NOTE: Require 0 or _interpolate fails
41 :
42 8 : s1 = _patch!(copy(s), s.tfinal, yundef, dyundef) # NOTE: strictly don't need copy
43 8 : out(s1)
44 : end
45 :
46 1732 : rv
47 : end
48 :
49 1728 : function _solvefill!(ps::TaskLocalValue{Problem{N, T, S, Dy, Out0}},
50 : t0::T, t1::T, y0::SVector{N, T},
51 : out::Out) where {N, T, S, Dy, Out0, Out}
52 1728 : _solvefill!(ps[], t0, t1, y0, out)
53 : end
54 :
55 : """
56 : Helper for [`sampleflowmap!`](@ref).
57 :
58 : Same as [`solve!`](@ref) but fill remaining output slots if solver
59 : stopped early (but `!isdone(solver)`).
60 :
61 : Uses [`_patch!`](@ref) to fill with undefined (`NaN`) vectors, where
62 : the last valid time (`first(pos(solver))`) is put into the second
63 : vector component.
64 :
65 : See also [`_patch!`](@ref).
66 : """ _solvefill!
67 :
68 : #
69 : # -----------------------------------------------------------------------------
70 : #
71 :
72 3 : function sampleflowmap(p::Problem{N, T, S, Dy, Out},
73 : t0::T, t1, extent...) where {N, T, S, Dy, Out}
74 3 : sz = (length(t1), length.(extent)...)
75 6 : Φ = Array{SVector{N, T}, N::Int+1}(undef, sz)
76 1287 : fill!(Φ, SVector{N, T}(ntuple(_->T(NaN), N)...))
77 :
78 3 : @assert issorted(t1)
79 :
80 3 : sampleflowmap!(p,Φ, t0, t1, extent...)
81 : end
82 :
83 20 : function sampleflowmap(p::Problem{N, T, S, Dy, Out},
84 : t0::T, t1::T, extent...) where {N, T, S, Dy, Out}
85 20 : sz = length.(extent)
86 40 : Φ = Array{SVector{N, T}, N::Int}(undef, sz)
87 1324 : fill!(Φ, SVector{N, T}(ntuple(_->T(NaN), N)...))
88 :
89 20 : sampleflowmap!(p,Φ, t0, t1, extent...)
90 : end
91 :
92 4 : function sampleflowmap(p::Problem{N, T, S, Dy, Out},
93 : t0, t1, extent...) where {N, T, S, Dy, Out}
94 4 : sz = (length(t0), length(t1), length.(extent)...)
95 8 : Φ = Array{SVector{N, T}, N::Int+2}(undef, sz)
96 4104 : fill!(Φ, SVector{N, T}(ntuple(_->T(NaN), N)...))
97 :
98 4 : @assert issorted(t0)
99 :
100 4 : sampleflowmap!(p, Φ, t0, t1, extent...)
101 : end
102 :
103 : # map coordinate from CartesianIndex to grid coordinate in sampleflowmap!
104 15616 : function _yfromindex(ysi)
105 15616 : ys, i = ysi
106 15616 : ys[i]
107 : end
108 :
109 : # get bounds in order and order as Bool
110 21 : function _temporalbounds(t1)
111 21 : if first(t1) <= last(t1)
112 21 : first(t1), last(t1), true
113 : else
114 0 : last(t1), first(t1), false
115 : end
116 : end
117 :
118 : # Is t0 is within times spanned by t1? Return index for splitting t1 or 0.
119 : #
120 : # Splitting at j means integration in one direction (forward or
121 : # backward) for t1[1:j] and into the other direction for t1[j+1:end]
122 21 : function _checkoverlap(t0::Real, t1)
123 21 : t1min, t1max, ascending = _temporalbounds(t1)
124 :
125 21 : if (t1min < t0 <= t1max)
126 7 : searchsortedlast(t1, t0; rev=!ascending)
127 14 : elseif t0 <= t1min
128 0 : 0
129 : else
130 2 : length(t1)
131 : end
132 : end
133 :
134 : # Partition indices of t1 into parts left and right such that t0 is
135 : # not inside t1[left] and not in t1[right]. t0 may by be equal to
136 : # the bound t1[left[end]].
137 : #
138 21 : function _partition(t0::Real, t1)
139 21 : n = length(t1)
140 21 : jsplit = _checkoverlap(t0, t1)
141 :
142 30 : @assert jsplit == 0 || (1 <= jsplit <= length(t1))
143 :
144 21 : left= 1:jsplit
145 24 : right = jsplit+1:n
146 :
147 21 : @assert length(left) + length(right) == n
148 :
149 21 : left, right
150 : end
151 :
152 21 : function sampleflowmap!(p::Problem{N, T, S, Dy, Out},
153 : Φ::AbstractArray{SVector{N, T}, M},
154 : t0::T, t1, extent...) where {N, T, S, Dy, Out, M}
155 21 : n = length(t1)
156 21 : dims = length.(extent)
157 :
158 21 : @assert length(extent) == N::Int
159 21 : @assert size(Φ) == (n, dims...)
160 :
161 21 : tend = t1[end]
162 :
163 21 : @assert first(t1) <= last(t1) "times t1 must be increasing"
164 : # NOTE: We already required/checked issorted(t0)!
165 : # NOTE: Allowing decreasing t1 requires extra care!
166 :
167 21 : left, right = _partition(t0, t1)
168 :
169 21 : ps = mtproblem(p)
170 :
171 21 : cart = CartesianIndices(dims)
172 :
173 21 : rv_t1 = reverse(t1)
174 41 : rv_left = reverse(left)
175 21 : rv_t1left = t1[rv_left]
176 42 : rv_idx = n:-1:1
177 :
178 : # for i in eachindex(cart)
179 21 : OhMyThreads.@tasks for i in eachindex(cart)
180 0 : @set begin
181 : scheduler=:dynamic
182 : # chunksize=4
183 : end
184 :
185 1344 : I = cart[i]
186 1344 : y0 = SVector(map(_yfromindex, zip(extent, Tuple(I)))...)
187 :
188 1344 : if isempty(left)
189 : # This is the standard case "forward" w/o overlap, i.e., t0 <= first(t1)
190 768 : _solvefill!(ps, t0, tend, y0, sample(view(Φ, :, Tuple(I)...), t1))
191 576 : elseif isempty(right)
192 : # This is the standard case "backward" w/o overlap, i.e., last(t1) < t0
193 192 : _solvefill!(ps, t0, t1[1], y0, sample(view(Φ, rv_idx, Tuple(I)...), rv_t1))
194 : else
195 384 : _solvefill!(ps, t0, t1[1], y0, sample(view(Φ, rv_left, Tuple(I)...), rv_t1left))
196 :
197 384 : _solvefill!(ps, t0, tend, y0, sample(view(Φ, right, Tuple(I)...), t1[right]))
198 : end
199 0 : end
200 :
201 : # NOTE: OhMyThreads seems to schedule better if the number of
202 : # threads is high (number of ht cores).
203 :
204 21 : Φ
205 : end
206 :
207 :
208 32 : function sampleflowmap!(p::Problem{N, T, S, Dy, Out},
209 : Φ::AbstractArray{SVector{N, T}, M},
210 : t0::T, tend::T, extent...) where {N, T, S, Dy, Out, M}
211 32 : dims = length.(extent)
212 32 : @assert length(extent) == N::Int
213 32 : @assert size(Φ) == dims
214 :
215 32 : ps = mtproblem(p)
216 :
217 32 : cart = CartesianIndices(dims)
218 :
219 32 : OhMyThreads.tmap!(Φ, eachindex(Φ); scheduler=:dynamic) do i
220 2048 : I = cart[i]
221 2048 : y0 = SVector(map(_yfromindex, zip(extent, Tuple(I)))...)
222 2048 : solve!(ps, t0, tend, y0)
223 2048 : isdone(ps[]) && return pos(ps[])[2]
224 :
225 32 : filled(SVector{N, T}, NaN, first(pos(ps[])); fill=NaN)
226 : # NOTE: tmap! cannot leave entry "untouched"
227 : end
228 :
229 32 : Φ
230 : end
231 :
232 4 : function sampleflowmap!(p::Problem{N, T, S, Dy, Out},
233 : Φ::AbstractArray{SVector{N, T}, M},
234 : t0, t1, extent...) where {N, T, S, Dy, Out, M}
235 4 : n0 = length(t0)
236 4 : n1 = length(t1)
237 4 : dims = length.(extent)
238 :
239 4 : @assert length(extent) == N::Int
240 4 : @assert size(Φ) == (n0, n1, dims...)
241 :
242 : # NOTE: This loop is not executed in multiple threads. MT is used by
243 : # the callee sampleflowmap!. This gives better load balancing
244 : # as cost to the calls to solve! differ depending on wether
245 : # times "overlap".
246 :
247 4 : for i in eachindex(t0)
248 64 : sampleflowmap!(p, view(Φ, i, ntuple(_ -> Colon(), N + 1)...),
249 : t0[i], t1, extent...)
250 28 : end
251 :
252 4 : Φ
253 : end
254 :
255 : """
256 : Φ = simple_sampleflowmap!(problem, t0, t1, X[, Y, ,...])
257 :
258 : Sample flow map by calling [`advect`](@ref) for every entry in `Φ`,
259 : i.e., starting a "full" integration for every sample.
260 :
261 : !!! note
262 : This method is *inefficient*. It serves as ground truth for
263 : testing [`sampleflowmap!`](@ref). Expect small deviations due
264 : to interpolation from different polynomial pieces.
265 :
266 : See also [`sampleflowmap!`](@ref), [`advect`](@ref)
267 : """
268 4 : function simple_sampleflowmap!(p::Problem{N, T, S, Dy, Out},
269 : Φ::AbstractArray{SVector{N, T}, M},
270 : t0, t1, extent...) where {N, T, S, Dy, Out, M}
271 4 : n0 = length(t0)
272 4 : n1 = length(t1)
273 4 : dims = length.(extent)
274 :
275 4 : @assert length(extent) == N::Int
276 4 : @assert size(Φ) == (n0, n1, dims...)
277 :
278 8 : for i in CartesianIndices(Φ)
279 4096 : j, k = i[1], i[2]
280 4096 : y0 = SVector(map(_yfromindex, zip(extent, Tuple(i)[3:end]))...)
281 :
282 4096 : Φ[i] = advect(p, t0[j], t1[k], y0)
283 4100 : end
284 :
285 4 : Φ
286 : end
287 :
288 : # NOTE: entries on diagonal are wrong for sampleflowmap! in case(t0=t1), there is more wrong...
289 :
290 : """
291 : Φ = sampleflowmap(problem, t0::T, T1, X[, Y, ,...])
292 :
293 : Sample flow map for [`problem`](@ref) in spatial grid `X × Y × ...`
294 : for times `T1`. Returns `Array{SVector{T, N}, N+1}` where `N` is
295 : the number of spatial dimensions and time steps are stored in the
296 : first dimension `Φ[:, ix, iy, ...]`.
297 :
298 : The sample times `T1` and dimensions `X, Y,...` can be given as a
299 : `Vector`, range or any generator.
300 :
301 : Φ = sampleflowmap(problem, T0, T1, X[, Y, ...])
302 :
303 : If sample times `T0` are also given as a `Vector`, etc. The result is
304 : an `N + 2`-dimensional array where the first and second dimensions
305 : refer to times `T0` and `T1` followed by the `N` spatial dimensions,
306 : i.e., the array is indexed as `Φ[it0, it1, ix, iy, ...]` with
307 : `it0 in eachindex(T0)` and `it1 in eachindex(T1)`.
308 :
309 : Φ = sampleflowmap(problem, t0::T, t1::T, X[, Y, ...])
310 :
311 : If start time `t0` and end time `t1` are given as *numbers*, Φ is
312 : returned as `Array{SVector{T, N}, N}`. This variant does not require a
313 : [`sample`](@ref) and is more efficient for single start and end times.
314 :
315 : !!! note "undefined entries"
316 : `sampleflowmap` yields undefined vectors if the solver terminated
317 : early (e.g., because `dy` returned `OutOfDomain`). In this case,
318 : the first component of the vector is `NaN`, and the *second*
319 : component stores the *time* when the solver terminated, i.e., the
320 : last time that the solver could yield a well-defined result (e.g.,
321 : the time when the trajectory intersected the domain boundary).
322 :
323 : !!! warning "precondition"
324 : Times `T0` and `T1` must be **monotonically increasing**.
325 :
326 : Times `T0` and `T1` may be represented as any iterable collection with
327 : fixed length, e.g., ranges or vectors. The time ranges *may overlap*.
328 :
329 : !!! note
330 : `sampleflowmap` uses **multi-threading** ([`mtproblem`](@ref)).
331 :
332 : [`sampleflowmap!`](@ref) provides in-place version.
333 :
334 : See also [`sampleflowmap!`](@ref), [`solve!`](@ref)
335 : """
336 : sampleflowmap
337 :
338 : """
339 : sampleflowmap!(problem, Φ, t0::T, T1, X[, Y, ,...])
340 : sampleflowmap!(problem, Φ, T0, T1, X[, Y, ,...])
341 : sampleflowmap!(problem, Φ, t0::T, t1::T, X[, Y, ,...])
342 :
343 : In-place version of [`sampleflowmap`](@ref) that modifies `Φ`. The
344 : dimensions of the array `Φ` must match times `T0`, `T1` and the
345 : spatial extent `X, Y, ...`.
346 :
347 : !!! note
348 : `sampleflowmap!` uses **multi-threading** ([`mtproblem`](@ref)).
349 :
350 : See also [`sampleflowmap`](@ref)
351 : """ sampleflowmap!
352 :
353 : # TODO: [feature] interpolate flow map Φ given dy (Hermite interpolation in time)
|
