1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2
3 /*
4 Fills array sign[] with signs of array f[]. If abs(f[i]) < vtol[i], the zero sign is taken.
5 All arrays should have length 'nev'
6 */
TSEventCalcSigns(PetscInt nev,const PetscReal * f,const PetscReal * vtol,PetscInt * sign)7 static inline void TSEventCalcSigns(PetscInt nev, const PetscReal *f, const PetscReal *vtol, PetscInt *sign)
8 {
9 for (PetscInt i = 0; i < nev; i++) {
10 if (PetscAbsReal(f[i]) < vtol[i]) sign[i] = 0;
11 else sign[i] = PetscSign(f[i]);
12 }
13 }
14
15 /*
16 TSEventInitialize - Initializes TSEvent for TSSolve
17 */
TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U)18 PetscErrorCode TSEventInitialize(TSEvent event, TS ts, PetscReal t, Vec U)
19 {
20 PetscFunctionBegin;
21 if (!event) PetscFunctionReturn(PETSC_SUCCESS);
22 PetscAssertPointer(event, 1);
23 PetscValidHeaderSpecific(ts, TS_CLASSID, 2);
24 PetscValidHeaderSpecific(U, VEC_CLASSID, 4);
25 event->ptime_prev = t;
26 event->iterctr = 0;
27 event->processing = PETSC_FALSE;
28 event->revisit_right = PETSC_FALSE;
29 PetscCallBack("TSEvent indicator", (*event->indicator)(ts, t, U, event->fvalue_prev, event->ctx));
30 TSEventCalcSigns(event->nevents, event->fvalue_prev, event->vtol, event->fsign_prev); // by this time event->vtol should have been defined
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
TSEventDestroy(TSEvent * event)34 PetscErrorCode TSEventDestroy(TSEvent *event)
35 {
36 PetscFunctionBegin;
37 PetscAssertPointer(event, 1);
38 if (!*event) PetscFunctionReturn(PETSC_SUCCESS);
39 if (--(*event)->refct > 0) {
40 *event = NULL;
41 PetscFunctionReturn(PETSC_SUCCESS);
42 }
43
44 PetscCall(PetscFree((*event)->fvalue_prev));
45 PetscCall(PetscFree((*event)->fvalue));
46 PetscCall(PetscFree((*event)->fvalue_right));
47 PetscCall(PetscFree((*event)->fsign_prev));
48 PetscCall(PetscFree((*event)->fsign));
49 PetscCall(PetscFree((*event)->fsign_right));
50 PetscCall(PetscFree((*event)->side));
51 PetscCall(PetscFree((*event)->side_prev));
52 PetscCall(PetscFree((*event)->justrefined_AB));
53 PetscCall(PetscFree((*event)->gamma_AB));
54 PetscCall(PetscFree((*event)->direction));
55 PetscCall(PetscFree((*event)->terminate));
56 PetscCall(PetscFree((*event)->events_zero));
57 PetscCall(PetscFree((*event)->vtol));
58
59 for (PetscInt i = 0; i < (*event)->recsize; i++) PetscCall(PetscFree((*event)->recorder.eventidx[i]));
60 PetscCall(PetscFree((*event)->recorder.eventidx));
61 PetscCall(PetscFree((*event)->recorder.nevents));
62 PetscCall(PetscFree((*event)->recorder.stepnum));
63 PetscCall(PetscFree((*event)->recorder.time));
64
65 PetscCall(PetscViewerDestroy(&(*event)->monitor));
66 PetscCall(PetscFree(*event));
67 PetscFunctionReturn(PETSC_SUCCESS);
68 }
69
70 // PetscClangLinter pragma disable: -fdoc-section-header-unknown
71 /*@
72 TSSetPostEventStep - Set the first time step to use after the event
73
74 Logically Collective
75
76 Input Parameters:
77 + ts - time integration context
78 - dt1 - first post event step
79
80 Options Database Key:
81 . -ts_event_post_event_step <dt1> - first time step after the event
82
83 Level: advanced
84
85 Notes:
86 `TSSetPostEventStep()` allows one to set a time step to use immediately following an event.
87 Note, if `TSAdapt` is allowed to interfere and reject steps, a large 'dt1' set by `TSSetPostEventStep()` may get truncated,
88 resulting in a smaller actual post-event step. See also the warning below regarding the `TSAdapt`.
89
90 The post-event time steps should be selected based on the post-event dynamics.
91 If the dynamics are stiff, or a significant jump in the equations or the state vector has taken place at the event,
92 conservative (small) steps should be employed. If not, then larger time steps may be appropriate.
93
94 This function accepts either a numerical value for `dt1`, or `PETSC_DECIDE`. The special value `PETSC_DECIDE` signals the event handler to follow
95 the originally planned trajectory, and is assumed by default.
96
97 To describe the way `PETSC_DECIDE` affects the post-event steps, consider a trajectory of time points t1 -> t2 -> t3 -> t4.
98 Suppose the `TS` has reached and calculated the solution at point t3, and has planned the next move: t3 -> t4.
99 At this moment, an event between t2 and t3 is detected, and after a few iterations it is resolved at point `te`, t2 < te < t3.
100 After event `te`, two post-event steps can be specified: the first one dt1 (`TSSetPostEventStep()`),
101 and the second one dt2 (`TSSetPostEventSecondStep()`). Both post-event steps can be either `PETSC_DECIDE`, or a number.
102 Four different combinations are possible\:
103
104 1. dt1 = `PETSC_DECIDE`, dt2 = `PETSC_DECIDE`. Then, after `te` `TS` goes to t3, and then to t4. This is the all-default behaviour.
105
106 2. dt1 = `PETSC_DECIDE`, dt2 = x2 (numerical). Then, after `te` `TS` goes to t3, and then to t3+x2.
107
108 3. dt1 = x1 (numerical), dt2 = x2 (numerical). Then, after `te` `TS` goes to te+x1, and then to te+x1+x2.
109
110 4. dt1 = x1 (numerical), dt2 = `PETSC_DECIDE`. Then, after `te` `TS` goes to te+x1, and event handler does not interfere to the subsequent steps.
111
112 In the special case when `te` == t3 with a good precision, the post-event step te -> t3 is not performed, so behaviour of (1) and (2) becomes\:
113
114 1a. After `te` `TS` goes to t4, and event handler does not interfere to the subsequent steps.
115
116 2a. After `te` `TS` goes to t4, and then to t4+x2.
117
118 Warning! When the second post-event step (either `PETSC_DECIDE` or a numerical value) is managed by the event handler, i.e. in cases 1, 2, 3 and 2a,
119 `TSAdapt` will never analyse (and never do a reasonable rejection of) the first post-event step. The first post-event step will always be accepted.
120 In this situation, it is the user's responsibility to make sure the step size is appropriate!
121 In cases 4 and 1a, however, `TSAdapt` will analyse the first post-event step, and is allowed to reject it.
122
123 This function can be called not only in the initial setup, but also inside the `postevent()` callback set with `TSSetEventHandler()`,
124 affecting the post-event steps for the current event, and the subsequent ones.
125 Thus, the strategy of the post-event time step definition can be adjusted on the fly.
126 In case several events are triggered in the given time point, only a single postevent handler is invoked,
127 and the user is to determine what post-event time step is more appropriate in this situation.
128
129 The default value is `PETSC_DECIDE`.
130
131 Developer Notes:
132 Event processing starts after visiting point t3, which means ts->adapt->dt_span_cached has been set to whatever value is required
133 when planning the step t3 -> t4.
134
135 .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`, `TSSetPostEventSecondStep()`
136 @*/
TSSetPostEventStep(TS ts,PetscReal dt1)137 PetscErrorCode TSSetPostEventStep(TS ts, PetscReal dt1)
138 {
139 PetscFunctionBegin;
140 ts->event->timestep_postevent = dt1;
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
144 // PetscClangLinter pragma disable: -fdoc-section-header-unknown
145 /*@
146 TSSetPostEventSecondStep - Set the second time step to use after the event
147
148 Logically Collective
149
150 Input Parameters:
151 + ts - time integration context
152 - dt2 - second post event step
153
154 Options Database Key:
155 . -ts_event_post_event_second_step <dt2> - second time step after the event
156
157 Level: advanced
158
159 Notes:
160 `TSSetPostEventSecondStep()` allows one to set the second time step after the event.
161
162 The post-event time steps should be selected based on the post-event dynamics.
163 If the dynamics are stiff, or a significant jump in the equations or the state vector has taken place at the event,
164 conservative (small) steps should be employed. If not, then larger time steps may be appropriate.
165
166 This function accepts either a numerical value for `dt2`, or `PETSC_DECIDE` (default).
167
168 To describe the way `PETSC_DECIDE` affects the post-event steps, consider a trajectory of time points t1 -> t2 -> t3 -> t4.
169 Suppose the `TS` has reached and calculated the solution at point t3, and has planned the next move: t3 -> t4.
170 At this moment, an event between t2 and t3 is detected, and after a few iterations it is resolved at point `te`, t2 < te < t3.
171 After event `te`, two post-event steps can be specified: the first one dt1 (`TSSetPostEventStep()`),
172 and the second one dt2 (`TSSetPostEventSecondStep()`). Both post-event steps can be either `PETSC_DECIDE`, or a number.
173 Four different combinations are possible\:
174
175 1. dt1 = `PETSC_DECIDE`, dt2 = `PETSC_DECIDE`. Then, after `te` `TS` goes to t3, and then to t4. This is the all-default behaviour.
176
177 2. dt1 = `PETSC_DECIDE`, dt2 = x2 (numerical). Then, after `te` `TS` goes to t3, and then to t3+x2.
178
179 3. dt1 = x1 (numerical), dt2 = x2 (numerical). Then, after `te` `TS` goes to te+x1, and then to te+x1+x2.
180
181 4. dt1 = x1 (numerical), dt2 = `PETSC_DECIDE`. Then, after `te` `TS` goes to te+x1, and event handler does not interfere to the subsequent steps.
182
183 In the special case when `te` == t3 with a good precision, the post-event step te -> t3 is not performed, so behaviour of (1) and (2) becomes\:
184
185 1a. After `te` `TS` goes to t4, and event handler does not interfere to the subsequent steps.
186
187 2a. After `te` `TS` goes to t4, and then to t4+x2.
188
189 Warning! When the second post-event step (either `PETSC_DECIDE` or a numerical value) is managed by the event handler, i.e. in cases 1, 2, 3 and 2a,
190 `TSAdapt` will never analyse (and never do a reasonable rejection of) the first post-event step. The first post-event step will always be accepted.
191 In this situation, it is the user's responsibility to make sure the step size is appropriate!
192 In cases 4 and 1a, however, `TSAdapt` will analyse the first post-event step, and is allowed to reject it.
193
194 This function can be called not only in the initial setup, but also inside the `postevent()` callback set with `TSSetEventHandler()`,
195 affecting the post-event steps for the current event, and the subsequent ones.
196
197 The default value is `PETSC_DECIDE`.
198
199 .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`, `TSSetPostEventStep()`
200 @*/
TSSetPostEventSecondStep(TS ts,PetscReal dt2)201 PetscErrorCode TSSetPostEventSecondStep(TS ts, PetscReal dt2)
202 {
203 PetscFunctionBegin;
204 ts->event->timestep_2nd_postevent = dt2;
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
208 /*@
209 TSSetEventTolerances - Set tolerances for event (indicator function) zero crossings
210
211 Logically Collective
212
213 Input Parameters:
214 + ts - time integration context
215 . tol - tolerance, `PETSC_CURRENT` to leave the current value
216 - vtol - array of tolerances or `NULL`, used in preference to `tol` if present
217
218 Options Database Key:
219 . -ts_event_tol <tol> - tolerance for event (indicator function) zero crossing
220
221 Level: beginner
222
223 Notes:
224 One must call `TSSetEventHandler()` before setting the tolerances.
225
226 The size of `vtol` should be equal to the number of events on the given process.
227
228 This function can be also called from the `postevent()` callback set with `TSSetEventHandler()`,
229 to adjust the tolerances on the fly.
230
231 .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`
232 @*/
TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[])233 PetscErrorCode TSSetEventTolerances(TS ts, PetscReal tol, PetscReal vtol[])
234 {
235 TSEvent event;
236
237 PetscFunctionBegin;
238 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
239 if (vtol) PetscAssertPointer(vtol, 3);
240 PetscCheck(ts->event, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set the events first by calling TSSetEventHandler()");
241
242 event = ts->event;
243 if (vtol) {
244 for (PetscInt i = 0; i < event->nevents; i++) event->vtol[i] = vtol[i];
245 } else {
246 if (tol != (PetscReal)PETSC_CURRENT) {
247 for (PetscInt i = 0; i < event->nevents; i++) event->vtol[i] = tol;
248 }
249 }
250 PetscFunctionReturn(PETSC_SUCCESS);
251 }
252
253 /*@C
254 TSSetEventHandler - Sets functions and parameters used for indicating events and handling them
255
256 Logically Collective
257
258 Input Parameters:
259 + ts - the `TS` context obtained from `TSCreate()`
260 . nevents - number of local events (i.e. managed by the given MPI process)
261 . direction - direction of zero crossing to be detected (one for each local event).
262 `-1` => zero crossing in negative direction,
263 `+1` => zero crossing in positive direction, `0` => both ways
264 . terminate - flag to indicate whether time stepping should be terminated after
265 an event is detected (one for each local event)
266 . indicator - callback defininig the user indicator functions whose sign changes (see `direction`) mark presence of the events
267 . postevent - [optional] user post-event callback; it can change the solution, ODE etc at the time of the event
268 - ctx - [optional] user-defined context for private data for the
269 `indicator()` and `postevent()` routines (use `NULL` if no
270 context is desired)
271
272 Calling sequence of `indicator`:
273 + ts - the `TS` context
274 . t - current time
275 . U - current solution
276 . fvalue - output array with values of local indicator functions (length == `nevents`) for time t and state-vector U
277 - ctx - the context passed as the final argument to `TSSetEventHandler()`
278
279 Calling sequence of `postevent`:
280 + ts - the `TS` context
281 . nevents_zero - number of triggered local events (whose indicator function is marked as crossing zero, and direction is appropriate)
282 . events_zero - indices of the triggered local events
283 . t - current time
284 . U - current solution
285 . forwardsolve - flag to indicate whether `TS` is doing a forward solve (`PETSC_TRUE`) or adjoint solve (`PETSC_FALSE`)
286 - ctx - the context passed as the final argument to `TSSetEventHandler()`
287
288 Options Database Keys:
289 + -ts_event_tol <tol> - tolerance for zero crossing check of indicator functions
290 . -ts_event_monitor - print choices made by event handler
291 . -ts_event_recorder_initial_size <recsize> - initial size of event recorder
292 . -ts_event_post_event_step <dt1> - first time step after event
293 . -ts_event_post_event_second_step <dt2> - second time step after event
294 - -ts_event_dt_min <dt> - minimum time step considered for TSEvent
295
296 Level: intermediate
297
298 Notes:
299 The indicator functions should be defined in the `indicator` callback using the components of solution `U` and/or time `t`.
300 Note that `U` is `PetscScalar`-valued, and the indicator functions are `PetscReal`-valued. It is the user's responsibility to
301 properly handle this difference, e.g. by applying `PetscRealPart()` or other appropriate conversion means.
302
303 The full set of events is distributed (by the user design) across MPI processes, with each process defining its own local sub-set of events.
304 However, the `postevent()` callback invocation is performed synchronously on all processes, including
305 those processes which have not currently triggered any events.
306
307 .seealso: [](ch_ts), `TSEvent`, `TSCreate()`, `TSSetTimeStep()`, `TSSetConvergedReason()`
308 @*/
TSSetEventHandler(TS ts,PetscInt nevents,PetscInt direction[],PetscBool terminate[],PetscErrorCode (* indicator)(TS ts,PetscReal t,Vec U,PetscReal fvalue[],PetscCtx ctx),PetscErrorCode (* postevent)(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,PetscCtx ctx),PetscCtx ctx)309 PetscErrorCode TSSetEventHandler(TS ts, PetscInt nevents, PetscInt direction[], PetscBool terminate[], PetscErrorCode (*indicator)(TS ts, PetscReal t, Vec U, PetscReal fvalue[], PetscCtx ctx), PetscErrorCode (*postevent)(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx), PetscCtx ctx)
310 {
311 TSAdapt adapt;
312 PetscReal hmin;
313 TSEvent event;
314 PetscBool flg;
315 #if defined(PETSC_USE_REAL_SINGLE)
316 PetscReal tol = 1e-4;
317 #else
318 PetscReal tol = 1e-6;
319 #endif
320
321 PetscFunctionBegin;
322 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
323 if (nevents) {
324 PetscAssertPointer(direction, 3);
325 PetscAssertPointer(terminate, 4);
326 }
327 PetscCall(PetscNew(&event));
328 PetscCall(PetscMalloc1(nevents, &event->fvalue_prev));
329 PetscCall(PetscMalloc1(nevents, &event->fvalue));
330 PetscCall(PetscMalloc1(nevents, &event->fvalue_right));
331 PetscCall(PetscMalloc1(nevents, &event->fsign_prev));
332 PetscCall(PetscMalloc1(nevents, &event->fsign));
333 PetscCall(PetscMalloc1(nevents, &event->fsign_right));
334 PetscCall(PetscMalloc1(nevents, &event->side));
335 PetscCall(PetscMalloc1(nevents, &event->side_prev));
336 PetscCall(PetscMalloc1(nevents, &event->justrefined_AB));
337 PetscCall(PetscMalloc1(nevents, &event->gamma_AB));
338 PetscCall(PetscMalloc1(nevents, &event->direction));
339 PetscCall(PetscMalloc1(nevents, &event->terminate));
340 PetscCall(PetscMalloc1(nevents, &event->events_zero));
341 PetscCall(PetscMalloc1(nevents, &event->vtol));
342 for (PetscInt i = 0; i < nevents; i++) {
343 event->direction[i] = direction[i];
344 event->terminate[i] = terminate[i];
345 event->justrefined_AB[i] = PETSC_FALSE;
346 event->gamma_AB[i] = 1;
347 event->side[i] = 2;
348 event->side_prev[i] = 0;
349 }
350 event->iterctr = 0;
351 event->processing = PETSC_FALSE;
352 event->revisit_right = PETSC_FALSE;
353 event->nevents = nevents;
354 event->indicator = indicator;
355 event->postevent = postevent;
356 event->ctx = ctx;
357 event->timestep_postevent = PETSC_DECIDE;
358 event->timestep_2nd_postevent = PETSC_DECIDE;
359 PetscCall(TSGetAdapt(ts, &adapt));
360 PetscCall(TSAdaptGetStepLimits(adapt, &hmin, NULL));
361 event->timestep_min = hmin;
362
363 event->recsize = 8; /* Initial size of the recorder */
364 PetscOptionsBegin(((PetscObject)ts)->comm, ((PetscObject)ts)->prefix, "TS Event options", "TS");
365 {
366 PetscCall(PetscOptionsReal("-ts_event_tol", "Tolerance for zero crossing check of indicator functions", "TSSetEventTolerances", tol, &tol, NULL));
367 PetscCall(PetscOptionsName("-ts_event_monitor", "Print choices made by event handler", "", &flg));
368 PetscCall(PetscOptionsInt("-ts_event_recorder_initial_size", "Initial size of event recorder", "", event->recsize, &event->recsize, NULL));
369 PetscCall(PetscOptionsDeprecated("-ts_event_post_eventinterval_step", "-ts_event_post_event_second_step", "3.21", NULL));
370 PetscCall(PetscOptionsReal("-ts_event_post_event_step", "First time step after event", "", event->timestep_postevent, &event->timestep_postevent, NULL));
371 PetscCall(PetscOptionsReal("-ts_event_post_event_second_step", "Second time step after event", "", event->timestep_2nd_postevent, &event->timestep_2nd_postevent, NULL));
372 PetscCall(PetscOptionsReal("-ts_event_dt_min", "Minimum time step considered for TSEvent", "", event->timestep_min, &event->timestep_min, NULL));
373 }
374 PetscOptionsEnd();
375
376 PetscCall(PetscMalloc1(event->recsize, &event->recorder.time));
377 PetscCall(PetscMalloc1(event->recsize, &event->recorder.stepnum));
378 PetscCall(PetscMalloc1(event->recsize, &event->recorder.nevents));
379 PetscCall(PetscMalloc1(event->recsize, &event->recorder.eventidx));
380 for (PetscInt i = 0; i < event->recsize; i++) PetscCall(PetscMalloc1(event->nevents, &event->recorder.eventidx[i]));
381 /* Initialize the event recorder */
382 event->recorder.ctr = 0;
383
384 for (PetscInt i = 0; i < event->nevents; i++) event->vtol[i] = tol;
385 if (flg) PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &event->monitor));
386
387 PetscCall(TSEventDestroy(&ts->event));
388 ts->event = event;
389 ts->event->refct = 1;
390 PetscFunctionReturn(PETSC_SUCCESS);
391 }
392
393 /*
394 TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
395 is reached.
396 */
TSEventRecorderResize(TSEvent event)397 static PetscErrorCode TSEventRecorderResize(TSEvent event)
398 {
399 PetscReal *time;
400 PetscInt *stepnum, *nevents;
401 PetscInt **eventidx;
402 PetscInt fact = 2;
403
404 PetscFunctionBegin;
405 /* Create larger arrays */
406 PetscCall(PetscMalloc1(fact * event->recsize, &time));
407 PetscCall(PetscMalloc1(fact * event->recsize, &stepnum));
408 PetscCall(PetscMalloc1(fact * event->recsize, &nevents));
409 PetscCall(PetscMalloc1(fact * event->recsize, &eventidx));
410 for (PetscInt i = 0; i < fact * event->recsize; i++) PetscCall(PetscMalloc1(event->nevents, &eventidx[i]));
411
412 /* Copy over data */
413 PetscCall(PetscArraycpy(time, event->recorder.time, event->recsize));
414 PetscCall(PetscArraycpy(stepnum, event->recorder.stepnum, event->recsize));
415 PetscCall(PetscArraycpy(nevents, event->recorder.nevents, event->recsize));
416 for (PetscInt i = 0; i < event->recsize; i++) PetscCall(PetscArraycpy(eventidx[i], event->recorder.eventidx[i], event->recorder.nevents[i]));
417
418 /* Destroy old arrays */
419 for (PetscInt i = 0; i < event->recsize; i++) PetscCall(PetscFree(event->recorder.eventidx[i]));
420 PetscCall(PetscFree(event->recorder.eventidx));
421 PetscCall(PetscFree(event->recorder.nevents));
422 PetscCall(PetscFree(event->recorder.stepnum));
423 PetscCall(PetscFree(event->recorder.time));
424
425 /* Set pointers */
426 event->recorder.time = time;
427 event->recorder.stepnum = stepnum;
428 event->recorder.nevents = nevents;
429 event->recorder.eventidx = eventidx;
430
431 /* Update the size */
432 event->recsize *= fact;
433 PetscFunctionReturn(PETSC_SUCCESS);
434 }
435
436 /*
437 Helper routine to handle user postevents and recording
438 */
TSPostEvent(TS ts,PetscReal t,Vec U)439 static PetscErrorCode TSPostEvent(TS ts, PetscReal t, Vec U)
440 {
441 TSEvent event = ts->event;
442 PetscBool restart = PETSC_FALSE;
443 PetscBool terminate = PETSC_FALSE;
444 PetscBool statechanged = PETSC_FALSE;
445 PetscInt ctr, stepnum;
446 PetscBool inflag[3], outflag[3];
447 PetscBool forwardsolve = PETSC_TRUE; // Flag indicating that TS is doing a forward solve
448
449 PetscFunctionBegin;
450 if (event->postevent) {
451 PetscObjectState state_prev, state_post;
452 PetscCall(PetscObjectStateGet((PetscObject)U, &state_prev));
453 PetscCallBack("TSEvent post-event processing", (*event->postevent)(ts, event->nevents_zero, event->events_zero, t, U, forwardsolve, event->ctx)); // TODO update 'restart' here?
454 PetscCall(PetscObjectStateGet((PetscObject)U, &state_post));
455 if (state_prev != state_post) {
456 restart = PETSC_TRUE;
457 statechanged = PETSC_TRUE;
458 }
459 }
460
461 // Handle termination events and step restart
462 for (PetscInt i = 0; i < event->nevents_zero; i++)
463 if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
464 inflag[0] = restart;
465 inflag[1] = terminate;
466 inflag[2] = statechanged;
467 PetscCallMPI(MPIU_Allreduce(inflag, outflag, 3, MPI_C_BOOL, MPI_LOR, ((PetscObject)ts)->comm));
468 restart = outflag[0];
469 terminate = outflag[1];
470 statechanged = outflag[2];
471 if (restart) PetscCall(TSRestartStep(ts));
472 if (terminate) PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_EVENT));
473
474 /*
475 Recalculate the indicator functions and signs if the state has been changed by the user postevent callback.
476 Note! If the state HAS NOT changed, the existing event->fsign (equal to zero) is kept, which:
477 - might have been defined using the previous (now-possibly-overridden) event->vtol,
478 - might have been set to zero on reaching a small time step rather than using the vtol criterion.
479 This will enforce keeping event->fsign = 0 where the zero-crossings were actually marked,
480 resulting in a more consistent behaviour of fsign's.
481 */
482 if (statechanged) {
483 if (event->monitor) PetscCall(PetscPrintf(((PetscObject)ts)->comm, "TSEvent: at time %g the vector state has been changed by PostEvent, recalculating fvalues and signs\n", (double)t));
484 PetscCall(VecLockReadPush(U));
485 PetscCallBack("TSEvent indicator", (*event->indicator)(ts, t, U, event->fvalue, event->ctx));
486 PetscCall(VecLockReadPop(U));
487 TSEventCalcSigns(event->nevents, event->fvalue, event->vtol, event->fsign); // note, event->vtol might have been changed by the postevent()
488 }
489
490 // Record the event in the event recorder
491 PetscCall(TSGetStepNumber(ts, &stepnum));
492 ctr = event->recorder.ctr;
493 if (ctr == event->recsize) PetscCall(TSEventRecorderResize(event));
494 event->recorder.time[ctr] = t;
495 event->recorder.stepnum[ctr] = stepnum;
496 event->recorder.nevents[ctr] = event->nevents_zero;
497 for (PetscInt i = 0; i < event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
498 event->recorder.ctr++;
499 PetscFunctionReturn(PETSC_SUCCESS);
500 }
501
502 // PetscClangLinter pragma disable: -fdoc-sowing-chars
503 /*
504 (modified) Anderson-Bjorck variant of regula falsi method, refines [tleft, t] or [t, tright] based on 'side' (-1 or +1).
505 The scaling parameter is defined based on the 'justrefined' flag, the history of repeats of 'side', and the threshold.
506 To escape certain failure modes, the algorithm may drift towards the bisection rule.
507 The value pointed to by 'side_prev' gets updated.
508 This function returns the new time step.
509
510 The underlying pure Anderson-Bjorck algorithm was taken as described in
511 J.M. Fernandez-Diaz, C.O. Menendez-Perez "A common framework for modified Regula Falsi methods and new methods of this kind", 2023.
512 The modifications subsequently introduced have little effect on the behaviour for simple cases requiring only a few iterations
513 (some minor convergence slowdown may take place though), but the effect becomes more pronounced for tough cases requiring many iterations.
514 For the latter situation the speed-up may be order(s) of magnitude compared to the classical Anderson-Bjorck.
515 The modifications (the threshold trick, and the drift towards bisection) were tested and tweaked
516 based on a number of test functions from the mentioned paper.
517 */
RefineAndersonBjorck(PetscReal tleft,PetscReal t,PetscReal tright,PetscReal fleft,PetscReal f,PetscReal fright,PetscInt side,PetscInt * side_prev,PetscBool justrefined,PetscReal * gamma)518 static inline PetscReal RefineAndersonBjorck(PetscReal tleft, PetscReal t, PetscReal tright, PetscReal fleft, PetscReal f, PetscReal fright, PetscInt side, PetscInt *side_prev, PetscBool justrefined, PetscReal *gamma)
519 {
520 PetscReal new_dt, scal = 1.0, scalB = 1.0, threshold = 0.0, power;
521 PetscInt reps = 0;
522 const PetscInt REPS_CAP = 8; // an upper bound to be imposed on 'reps' (set to 8, somewhat arbitrary number, found after some tweaking)
523
524 // Preparations
525 if (justrefined) {
526 if (*side_prev * side > 0) *side_prev += side; // the side keeps repeating -> increment the side counter (-ve or +ve)
527 else *side_prev = side; // reset the counter
528 reps = PetscMin(*side_prev * side, REPS_CAP); // the length of the recent side-repeat series, including the current 'side'
529 threshold = PetscPowReal(0.5, reps) * 0.1; // ad-hoc strategy for threshold calculation (involved some tweaking)
530 } else *side_prev = side; // initial reset of the counter
531
532 // Time step calculation
533 if (side == -1) {
534 if (justrefined && fright != 0.0 && fleft != 0.0) {
535 scal = (fright - f) / fright;
536 scalB = -f / fleft;
537 }
538 } else { // must be side == +1
539 if (justrefined && fleft != 0.0 && fright != 0.0) {
540 scal = (fleft - f) / fleft;
541 scalB = -f / fright;
542 }
543 }
544
545 if (scal < threshold) scal = 0.5;
546 if (reps > 1) *gamma *= scal; // side did not switch since the last time, accumulate gamma
547 else *gamma = 1.0; // side switched -> reset gamma
548 power = PetscMax(0.0, (reps - 2.0) / (REPS_CAP - 2.0));
549 scal = PetscPowReal(scalB / *gamma, power) * (*gamma); // mix the Anderson-Bjorck scaling and Bisection scaling
550
551 if (side == -1) new_dt = scal * fleft / (scal * fleft - f) * (t - tleft);
552 else new_dt = f / (f - scal * fright) * (tright - t);
553 /*
554 In tough cases (e.g. a polynomial of high order), there is a failure mode for the standard Anderson-Bjorck,
555 when the new proposed point jumps from one end-point of the bracket to the other, however the bracket
556 is contracting very slowly. A larger threshold for 'scal' prevents entering this mode.
557 On the other hand, if the iteration gets stuck near one end-point of the bracket, and the 'side' does not switch for a while,
558 the 'scal' drifts towards the bisection approach (via scalB), ensuring stable convergence.
559 */
560 return new_dt;
561 }
562
563 /*
564 Checks if the current point (t) is the zero-crossing location, based on the indicator function signs and direction[]:
565 - using the dt_min criterion,
566 - using the vtol criterion.
567 The situation (fsign_prev, fsign) = (0, 0) is treated as staying in the near-zero-zone of the previous zero-crossing,
568 and is not marked as a new zero-crossing.
569 This function may update event->side[].
570 */
TSEventTestZero(TS ts,PetscReal t)571 static PetscErrorCode TSEventTestZero(TS ts, PetscReal t)
572 {
573 TSEvent event = ts->event;
574
575 PetscFunctionBegin;
576 for (PetscInt i = 0; i < event->nevents; i++) {
577 const PetscBool bracket_is_left = (event->fsign_prev[i] * event->fsign[i] < 0 && event->fsign[i] * event->direction[i] >= 0) ? PETSC_TRUE : PETSC_FALSE;
578
579 if (bracket_is_left && ((t - event->ptime_prev <= event->timestep_min) || event->revisit_right)) event->side[i] = 0; // mark zero-crossing from dt_min; 'bracket_is_left' accounts for direction
580 if (event->fsign[i] == 0 && event->fsign_prev[i] != 0 && event->fsign_prev[i] * event->direction[i] <= 0) event->side[i] = 0; // mark zero-crossing from vtol
581 }
582 PetscFunctionReturn(PETSC_SUCCESS);
583 }
584
585 /*
586 Checks if [fleft, f] or [f, fright] are 'brackets', i.e. intervals with the sign change, satisfying the 'direction'.
587 The right interval is only checked if iterctr > 0 (i.e. Anderson-Bjorck refinement has started).
588 The intervals like [0, x] and [x, 0] are not counted as brackets, i.e. intervals with the sign change.
589 The function returns the 'side' value: -1 (left, or both are brackets), +1 (only right one), +2 (neither).
590 */
TSEventTestBracket(PetscInt fsign_left,PetscInt fsign,PetscInt fsign_right,PetscInt direction,PetscInt iterctr)591 static inline PetscInt TSEventTestBracket(PetscInt fsign_left, PetscInt fsign, PetscInt fsign_right, PetscInt direction, PetscInt iterctr)
592 {
593 PetscInt side = 2;
594 if (fsign_left * fsign < 0 && fsign * direction >= 0) side = -1;
595 if (side != -1 && iterctr > 0 && fsign * fsign_right < 0 && fsign_right * direction >= 0) side = 1;
596 return side;
597 }
598
599 /*
600 Caps the time steps, accounting for evaluation time points.
601 It uses 'event->timestep_cache' as a time step to calculate the tolerance for eval_times points detection. This
602 is done since the event resolution may result in significant time step refinement, and we don't use these small steps for tolerances.
603 To enhance the consistency of eval_times points detection, tolerance 'eval_times->worktol' is reused later in the TSSolve iteration.
604 If a user-defined step is cut by this function, the input uncut step is saved to adapt->dt_span_cached.
605 Flag 'user_dt' indicates if the step was defined by user.
606 */
TSEvent_dt_cap(TS ts,PetscReal t,PetscReal dt,PetscBool user_dt)607 static inline PetscReal TSEvent_dt_cap(TS ts, PetscReal t, PetscReal dt, PetscBool user_dt)
608 {
609 PetscReal res = dt;
610 if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
611 PetscReal maxdt = ts->max_time - t; // this may be overridden by eval_times
612 PetscBool cut_made = PETSC_FALSE;
613 PetscReal eps = 10 * PETSC_MACHINE_EPSILON;
614 if (ts->eval_times) {
615 PetscInt idx = ts->eval_times->time_point_idx;
616 PetscInt Ns = ts->eval_times->num_time_points;
617 PetscReal *st = ts->eval_times->time_points;
618
619 if (ts->eval_times->worktol == 0) ts->eval_times->worktol = ts->eval_times->reltol * ts->event->timestep_cache + ts->eval_times->abstol; // in case TSAdaptChoose() has not defined it
620 if (idx < Ns && PetscIsCloseAtTol(t, st[idx], ts->eval_times->worktol, 0)) { // just hit a evaluation time point
621 if (idx + 1 < Ns) maxdt = st[idx + 1] - t; // ok to use the next evaluation time point
622 else maxdt = ts->max_time - t; // can't use the next evaluation time point: they have finished
623 } else if (idx < Ns) maxdt = st[idx] - t; // haven't hit a evaluation time point, use the nearest one
624 }
625 maxdt = PetscMin(maxdt, ts->max_time - t);
626 PetscCheck((maxdt > eps) || (PetscAbsReal(maxdt) <= eps && PetscIsCloseAtTol(t, ts->max_time, eps, 0)), PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state: bad maxdt (%g) in TSEvent_dt_cap()", (double)maxdt);
627
628 if (PetscIsCloseAtTol(dt, maxdt, eps, 0)) res = maxdt; // no cut
629 else {
630 if (dt > maxdt) {
631 res = maxdt; // yes cut
632 cut_made = PETSC_TRUE;
633 } else res = dt; // no cut
634 }
635 if (ts->adapt && user_dt) { // only update dt_span_cached for the user-defined step
636 if (cut_made) ts->adapt->dt_eval_times_cached = dt;
637 else ts->adapt->dt_eval_times_cached = 0;
638 }
639 }
640 return res;
641 }
642
643 /*
644 Updates the left-end values
645 */
TSEvent_update_left(TSEvent event,PetscReal t)646 static inline void TSEvent_update_left(TSEvent event, PetscReal t)
647 {
648 for (PetscInt i = 0; i < event->nevents; i++) {
649 event->fvalue_prev[i] = event->fvalue[i];
650 event->fsign_prev[i] = event->fsign[i];
651 }
652 event->ptime_prev = t;
653 }
654
655 /*
656 Updates the right-end values
657 */
TSEvent_update_right(TSEvent event,PetscReal t)658 static inline void TSEvent_update_right(TSEvent event, PetscReal t)
659 {
660 for (PetscInt i = 0; i < event->nevents; i++) {
661 event->fvalue_right[i] = event->fvalue[i];
662 event->fsign_right[i] = event->fsign[i];
663 }
664 event->ptime_right = t;
665 }
666
667 /*
668 Updates the current values from the right-end values
669 */
TSEvent_update_from_right(TSEvent event)670 static inline PetscReal TSEvent_update_from_right(TSEvent event)
671 {
672 for (PetscInt i = 0; i < event->nevents; i++) {
673 event->fvalue[i] = event->fvalue_right[i];
674 event->fsign[i] = event->fsign_right[i];
675 }
676 return event->ptime_right;
677 }
678
Not_PETSC_DECIDE(PetscReal dt)679 static inline PetscBool Not_PETSC_DECIDE(PetscReal dt)
680 {
681 return dt == PETSC_DECIDE ? PETSC_FALSE : PETSC_TRUE;
682 }
683
684 // PetscClangLinter pragma disable: -fdoc-section-spacing
685 // PetscClangLinter pragma disable: -fdoc-section-header-unknown
686 // PetscClangLinter pragma disable: -fdoc-section-header-spelling
687 // PetscClangLinter pragma disable: -fdoc-section-header-missing
688 // PetscClangLinter pragma disable: -fdoc-section-header-fishy-header
689 // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
690 // PetscClangLinter pragma disable: -fdoc-synopsis-missing-description
691 // PetscClangLinter pragma disable: -fdoc-sowing-chars
692 /*
693 TSEventHandler - the main function to perform a single iteration of event detection.
694
695 Developer notes:
696 A) The 'event->iterctr > 0' is used as an indicator that Anderson-Bjorck refinement has started.
697 B) If event->iterctr == 0, then justrefined_AB[i] is always false.
698 C) The right-end quantities: ptime_right, fvalue_right[i] and fsign_right[i] are only guaranteed to be valid
699 for event->iterctr > 0.
700 D) If event->iterctr > 0, then event->processing is PETSC_TRUE; the opposite may not hold.
701 E) When event->processing == PETSC_TRUE and event->iterctr == 0, the event handler iterations are complete, but
702 the event handler continues managing the 1st and 2nd post-event steps. In this case the 1st post-event step
703 proposed by the event handler is not checked by TSAdapt, and is always accepted (beware!).
704 However, if the 2nd post-event step is not managed by the event handler (e.g. 1st = numerical, 2nd = PETSC_DECIDE),
705 condition "E" does not hold, and TSAdapt may reject/adjust the 1st post-event step.
706 F) event->side[i] may take values: 0 <=> point t is a zero-crossing for indicator function i (via vtol/dt_min criterion);
707 -1/+1 <=> detected a bracket to the left/right of t for indicator function i; +2 <=> no brackets/zero-crossings.
708 G) The signs event->fsign[i] (with values 0/-1/+1) are calculated for each new point. Zero sign is set if the function value is
709 smaller than the tolerance. Besides, zero sign is enforced after marking a zero-crossing due to small bracket size criterion.
710
711 The intervals with the indicator function sign change (i.e. containing the potential zero-crossings) are called 'brackets'.
712 To find a zero-crossing, the algorithm first locates a bracket, and then sequentially subdivides it, generating a sequence
713 of brackets whose length tends to zero. The bracket subdivision involves the (modified) Anderson-Bjorck method.
714
715 Apart from the comments scattered throughout the code to clarify different lines and blocks,
716 a few tricky aspects of the algorithm (and the underlying reasoning) are discussed in detail below:
717
718 =Sign tracking=
719 When a zero-crossing is found, the sign variable (event->fsign[i]) is set to zero for the current point t.
720 This happens both for zero-crossings triggered via the vtol criterion, and those triggered via the dt_min
721 criterion. After the event, as the TS steps forward, the current sign values are handed over to event->fsign_prev[i].
722 The recalculation of signs is avoided if possible: e.g. if a 'vtol' criterion resulted in a zero-crossing at point t,
723 but the subsequent call to postevent() handler decreased 'vtol', making the indicator function no longer "close to zero"
724 at point t, the fsign[i] will still consistently keep the zero value. This allows avoiding the erroneous duplication
725 of events:
726
727 E.g. consider a bracket [t0, t2], where f0 < 0, f2 > 0, which resulted in a zero-crossing t1 with f1 < 0, abs(f1) < vtol.
728 Suppose the postevent() handler changes vtol to vtol*, such that abs(f1) > vtol*. The TS makes a step t1 -> t3, where
729 again f1 < 0, f3 > 0, and the event handler will find a new event near t1, which is actually a duplication of the
730 original event at t1. The duplications are avoided by NOT counting the sign progressions 0 -> +1, or 0 -> -1
731 as brackets. Tracking (instead of recalculating) the sign values makes this procedure work more consistently.
732
733 The sign values are however recalculated if the postevent() callback has changed the current solution vector U
734 (such a change resets everything).
735 The sign value is also set to zero if the dt_min criterion has triggered the event. This allows the algorithm to
736 work consistently, irrespective of the type of criterion involved (vtol/dt_min).
737
738 =Event from min bracket=
739 When the event handler ends up with a bracket [t0, t1] with size <= dt_min, a zero crossing is reported at t1,
740 and never at t0. If such a bracket is discovered when TS is staying at t0, one more step forward (to t1) is necessary
741 to mark the found event. This is the situation of revisiting t1, which is described below (see Revisiting).
742
743 Why t0 is not reported as event location? Suppose it is, and let f0 < 0, f1 > 0. Also suppose that the
744 postevent() handler has slightly changed the solution U, so the sign at t0 is recalculated: it equals -1. As the TS steps
745 further: t0 -> t2, with sign0 == -1, and sign2 == +1, the event handler will locate the bracket [t0, t2], eventually
746 resolving a new event near t1, i.e. finding a duplicate event.
747 This situation is avoided by reporting the event at t1 in the first place.
748
749 =Revisiting=
750 When handling the situation with small bracket size, the TS solver may happen to visit the same point twice,
751 but with different results.
752
753 E.g. originally it discovered a bracket with sign change [t0, t10], and started resolving the zero-crossing,
754 visiting the points t1,...,t9 : t0 < t1 < ... < t9 < t10. Suppose that at t9 the algorithm discovers
755 that [t9, t10] is a bracket with the sign change it was looking for, and that |t10 - t9| is too small.
756 So point t10 should be revisited and marked as the zero crossing (by the minimum bracket size criterion).
757 On re-visiting t10, via the refined sequence of steps t0,...,t10, the TS solver may arrive at a solution U*
758 different from the solution U it found at t10 originally. Hence, the indicator functions at t10 may become different,
759 and the condition of the sign change, which existed originally, may disappear, breaking the logic of the algorithm.
760
761 To handle such (-=unlikely=-, but possible) situations, two strategies can be considered:
762 1) [not used here] Allow the brackets with sign change to disappear during iterations. The algorithm should be able
763 to cleanly exit the iteration and leave all the objects/variables/caches involved in a valid state.
764 2) [ADOPTED HERE!] On revisiting t10, the event handler reuses the indicator functions previously calculated for the
765 original solution U. This U may be less precise than U*, but this trick does not allow the algorithm logic to break down.
766 HOWEVER, the original U is not stored anywhere, it is essentially lost since the TS performed the rollback from it.
767 On revisiting t10, the updated solution U* will inevitably be found and used everywhere EXCEPT the current
768 indicator functions calculation, e.g. U* will be used in the postevent() handler call. Since t10 is the event location,
769 the appropriate indicator-function-signs will be enforced to be 0 (regardless if the solution was U or U*).
770 If the solution is then changed by the postevent(), the indicator-function-signs will be recalculated.
771
772 Whether the algorithm is revisiting a point in the current TSEventHandler() call is flagged by 'event->revisit_right'.
773 */
TSEventHandler(TS ts)774 PetscErrorCode TSEventHandler(TS ts)
775 {
776 TSEvent event;
777 PetscReal t, dt_next = 0.0;
778 Vec U;
779 PetscInt minsidein = 2, minsideout = 2; // minsideout is sync on all ranks
780 PetscBool finished = PETSC_FALSE; // should stay sync on all ranks
781 PetscBool revisit_right_cache; // [sync] flag for inner consistency checks
782
783 PetscFunctionBegin;
784 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
785
786 if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
787 event = ts->event;
788 event->nevents_zero = 0;
789 revisit_right_cache = event->revisit_right;
790 for (PetscInt i = 0; i < event->nevents; i++) event->side[i] = 2; // side's are reset on each new iteration
791 if (event->iterctr == 0)
792 for (PetscInt i = 0; i < event->nevents; i++) event->justrefined_AB[i] = PETSC_FALSE;
793
794 PetscCall(TSGetTime(ts, &t));
795 if (!event->processing) { // update the caches
796 PetscReal dt;
797 PetscCall(TSGetTimeStep(ts, &dt));
798 event->ptime_cache = t;
799 event->timestep_cache = dt; // the next TS move is planned to be: t -> t+dt
800 }
801 if (event->processing && event->iterctr == 0 && Not_PETSC_DECIDE(event->timestep_2nd_postevent)) { // update the caches while processing the post-event steps
802 event->ptime_cache = t;
803 event->timestep_cache = event->timestep_2nd_postevent;
804 }
805
806 PetscCall(TSGetSolution(ts, &U)); // if revisiting, this will be the updated U* (see discussion on "Revisiting" in the Developer notes above)
807 if (event->revisit_right) {
808 PetscReal tr = TSEvent_update_from_right(event);
809 PetscCheck(PetscAbsReal(tr - t) < PETSC_SMALL, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Inconsistent time value when performing 'revisiting' in TSEventHandler()");
810 } else {
811 PetscCall(VecLockReadPush(U));
812 PetscCallBack("TSEvent indicator", (*event->indicator)(ts, t, U, event->fvalue, event->ctx)); // fill fvalue's at point 't'
813 PetscCall(VecLockReadPop(U));
814 TSEventCalcSigns(event->nevents, event->fvalue, event->vtol, event->fsign); // fill fvalue signs
815 }
816 PetscCall(TSEventTestZero(ts, t)); // check if the current point 't' is the event location; event->side[] may get updated
817
818 for (PetscInt i = 0; i < event->nevents; i++) { // check for brackets on the left/right of 't'
819 if (event->side[i] != 0) event->side[i] = TSEventTestBracket(event->fsign_prev[i], event->fsign[i], event->fsign_right[i], event->direction[i], event->iterctr);
820 minsidein = PetscMin(minsidein, event->side[i]);
821 }
822 PetscCallMPI(MPIU_Allreduce(&minsidein, &minsideout, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)ts)));
823 /*
824 minsideout (sync on all ranks) indicates the minimum of the following states:
825 -1 : [ptime_prev, t] is a bracket for some indicator-function-i
826 +1 : [t, ptime_right] is a bracket for some indicator-function-i
827 0 : t is a zero-crossing for some indicator-function-i
828 2 : none of the above
829 */
830 PetscCheck(!event->revisit_right || minsideout == 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "minsideout != 0 when performing 'revisiting' in TSEventHandler()");
831
832 if (minsideout == -1 || minsideout == +1) { // this if-branch will refine the left/right bracket
833 const PetscReal bracket_size = (minsideout == -1) ? t - event->ptime_prev : event->ptime_right - t; // sync on all ranks
834
835 if (minsideout == +1 && bracket_size <= event->timestep_min) { // check if the bracket (right) is small
836 // [--------------------|-]
837 dt_next = bracket_size; // need one more step to get to event->ptime_right
838 event->revisit_right = PETSC_TRUE;
839 TSEvent_update_left(event, t);
840 if (event->monitor)
841 PetscCall(PetscViewerASCIIPrintf(event->monitor, "[%d] TSEvent: iter %" PetscInt_FMT " - reached too small bracket [%g - %g], next stepping to its right end %g (revisiting)\n", PetscGlobalRank, event->iterctr, (double)event->ptime_prev,
842 (double)event->ptime_right, (double)(event->ptime_prev + dt_next)));
843 } else { // the bracket is not very small -> refine it
844 // [--------|-------------]
845 if (bracket_size <= 2 * event->timestep_min) dt_next = bracket_size / 2; // the bracket is almost small -> bisect it
846 else { // the bracket is not small -> use Anderson-Bjorck
847 PetscReal dti_min = PETSC_MAX_REAL;
848 for (PetscInt i = 0; i < event->nevents; i++) {
849 if (event->side[i] == minsideout) { // only refine the appropriate brackets
850 PetscReal dti = RefineAndersonBjorck(event->ptime_prev, t, event->ptime_right, event->fvalue_prev[i], event->fvalue[i], event->fvalue_right[i], event->side[i], &event->side_prev[i], event->justrefined_AB[i], &event->gamma_AB[i]);
851 dti_min = PetscMin(dti_min, dti);
852 }
853 }
854 PetscCallMPI(MPIU_Allreduce(&dti_min, &dt_next, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
855 if (dt_next < event->timestep_min) dt_next = event->timestep_min;
856 if (bracket_size - dt_next < event->timestep_min) dt_next = bracket_size - event->timestep_min;
857 }
858
859 if (minsideout == -1) { // minsideout == -1, update the right-end values, retain the left-end values
860 TSEvent_update_right(event, t);
861 PetscCall(TSRollBack(ts));
862 PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING)); // e.g. to override TS_CONVERGED_TIME on reaching ts->max_time
863 } else TSEvent_update_left(event, t); // minsideout == +1, update the left-end values, retain the right-end values
864
865 for (PetscInt i = 0; i < event->nevents; i++) { // update the "Anderson-Bjorck" flags
866 if (event->side[i] == minsideout) {
867 event->justrefined_AB[i] = PETSC_TRUE; // only for these i's Anderson-Bjorck was invoked
868 if (event->monitor)
869 PetscCall(PetscViewerASCIIPrintf(event->monitor, "[%d] TSEvent: iter %" PetscInt_FMT " - Event %" PetscInt_FMT " refining the bracket with sign change [%g - %g], next stepping to %g\n", PetscGlobalRank, event->iterctr, i, (double)event->ptime_prev,
870 (double)event->ptime_right, (double)(event->ptime_prev + dt_next)));
871 } else event->justrefined_AB[i] = PETSC_FALSE; // for these i's Anderson-Bjorck was not invoked
872 }
873 }
874 event->iterctr++;
875 event->processing = PETSC_TRUE;
876 } else if (minsideout == 0) { // found the appropriate zero-crossing (and no brackets to the left), finishing!
877 // [--------0-------------]
878 finished = PETSC_TRUE;
879 event->revisit_right = PETSC_FALSE;
880 for (PetscInt i = 0; i < event->nevents; i++)
881 if (event->side[i] == minsideout) {
882 event->events_zero[event->nevents_zero++] = i;
883 if (event->fsign[i] == 0) { // vtol was engaged
884 if (event->monitor)
885 PetscCall(PetscViewerASCIIPrintf(event->monitor, "[%d] TSEvent: iter %" PetscInt_FMT " - Event %" PetscInt_FMT " zero crossing located at time %g (tol=%g)\n", PetscGlobalRank, event->iterctr, i, (double)t, (double)event->vtol[i]));
886 } else { // dt_min was engaged
887 event->fsign[i] = 0; // sign = 0 is enforced further
888 if (event->monitor)
889 PetscCall(PetscViewerASCIIPrintf(event->monitor, "[%d] TSEvent: iter %" PetscInt_FMT " - Event %" PetscInt_FMT " accepting time %g as event location, due to reaching too small bracket [%g - %g]\n", PetscGlobalRank, event->iterctr, i, (double)t,
890 (double)event->ptime_prev, (double)t));
891 }
892 }
893 event->iterctr++;
894 event->processing = PETSC_TRUE;
895 } else { // minsideout == 2: no brackets, no zero-crossings
896 // [----------------------]
897 PetscCheck(event->iterctr == 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state (event->iterctr != 0) in TSEventHandler()");
898 if (event->processing) {
899 PetscReal dt2;
900 if (event->timestep_2nd_postevent == PETSC_DECIDE) dt2 = event->timestep_cache; // (1)
901 else dt2 = event->timestep_2nd_postevent; // (2), (2a), (3)
902 PetscCall(TSSetTimeStep(ts, TSEvent_dt_cap(ts, t, dt2, Not_PETSC_DECIDE(event->timestep_2nd_postevent)))); // set the second post-event step
903 }
904 event->processing = PETSC_FALSE;
905 }
906
907 // if 'revisit_right' was flagged before the current iteration started, the iteration is expected to finish
908 PetscCheck(!revisit_right_cache || finished, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state of 'revisit_right_cache' in TSEventHandler()");
909
910 if (finished) { // finished handling the current event
911 PetscCall(TSPostEvent(ts, t, U));
912
913 PetscReal dt1;
914 if (event->timestep_postevent == PETSC_DECIDE) { // (1), (2)
915 dt1 = event->ptime_cache - t;
916 event->processing = PETSC_TRUE;
917 if (PetscAbsReal(dt1) < PETSC_SMALL) { // (1a), (2a): the cached post-event point == event point
918 dt1 = event->timestep_cache;
919 event->processing = Not_PETSC_DECIDE(event->timestep_2nd_postevent);
920 }
921 } else { // (3), (4)
922 dt1 = event->timestep_postevent; // 1st post-event dt = user-provided value
923 event->processing = Not_PETSC_DECIDE(event->timestep_2nd_postevent);
924 }
925
926 PetscCall(TSSetTimeStep(ts, TSEvent_dt_cap(ts, t, dt1, Not_PETSC_DECIDE(event->timestep_postevent)))); // set the first post-event step
927 event->iterctr = 0;
928 } // if-finished
929
930 if (event->iterctr == 0) TSEvent_update_left(event, t); // not found an event, or finished the event
931 else {
932 PetscCall(TSGetTime(ts, &t)); // update 't' to account for potential rollback
933 PetscCall(TSSetTimeStep(ts, TSEvent_dt_cap(ts, t, dt_next, PETSC_FALSE))); // continue resolving the event
934 }
935 PetscFunctionReturn(PETSC_SUCCESS);
936 }
937
TSAdjointEventHandler(TS ts)938 PetscErrorCode TSAdjointEventHandler(TS ts)
939 {
940 TSEvent event;
941 PetscReal t;
942 Vec U;
943 PetscInt ctr;
944 PetscBool forwardsolve = PETSC_FALSE; // Flag indicating that TS is doing an adjoint solve
945
946 PetscFunctionBegin;
947 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
948 if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
949 event = ts->event;
950
951 PetscCall(TSGetTime(ts, &t));
952 PetscCall(TSGetSolution(ts, &U));
953
954 ctr = event->recorder.ctr - 1;
955 if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
956 // Call the user post-event function
957 if (event->postevent) {
958 PetscCallBack("TSEvent post-event processing", (*event->postevent)(ts, event->recorder.nevents[ctr], event->recorder.eventidx[ctr], t, U, forwardsolve, event->ctx));
959 event->recorder.ctr--;
960 }
961 }
962 PetscCall(PetscBarrier((PetscObject)ts));
963 PetscFunctionReturn(PETSC_SUCCESS);
964 }
965
966 /*@
967 TSGetNumEvents - Get the number of events defined on the given MPI process
968
969 Logically Collective
970
971 Input Parameter:
972 . ts - the `TS` context
973
974 Output Parameter:
975 . nevents - the number of local events on each MPI process
976
977 Level: intermediate
978
979 .seealso: [](ch_ts), `TSEvent`, `TSSetEventHandler()`
980 @*/
TSGetNumEvents(TS ts,PetscInt * nevents)981 PetscErrorCode TSGetNumEvents(TS ts, PetscInt *nevents)
982 {
983 PetscFunctionBegin;
984 *nevents = ts->event->nevents;
985 PetscFunctionReturn(PETSC_SUCCESS);
986 }
987