Actual source code: tsevent.c
1: #include <petsc/private/tsimpl.h>
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: */
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: }
15: /*
16: TSEventInitialize - Initializes TSEvent for TSSolve
17: */
18: PetscErrorCode TSEventInitialize(TSEvent event, TS ts, PetscReal t, Vec U)
19: {
20: PetscFunctionBegin;
21: if (!event) PetscFunctionReturn(PETSC_SUCCESS);
22: PetscAssertPointer(event, 1);
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: }
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: }
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));
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));
65: PetscCall(PetscViewerDestroy(&(*event)->monitor));
66: PetscCall(PetscFree(*event));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
71: /*@
72: TSSetPostEventStep - Set the first time step to use after the event
74: Logically Collective
76: Input Parameters:
77: + ts - time integration context
78: - dt1 - first post event step
80: Options Database Key:
81: . -ts_event_post_event_step <dt1> - first time step after the event
83: Level: advanced
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`.
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.
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.
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\:
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.
106: 2. dt1 = `PETSC_DECIDE`, dt2 = x2 (numerical). Then, after `te` TS goes to t3, and then to t3+x2.
108: 3. dt1 = x1 (numerical), dt2 = x2 (numerical). Then, after `te` TS goes to te+x1, and then to te+x1+x2.
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.
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\:
114: 1a. After `te` TS goes to t4, and event handler does not interfere to the subsequent steps.
116: 2a. After `te` TS goes to t4, and then to t4+x2.
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.
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.
129: The default value is `PETSC_DECIDE`.
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.
135: .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`, `TSSetPostEventSecondStep()`
136: @*/
137: PetscErrorCode TSSetPostEventStep(TS ts, PetscReal dt1)
138: {
139: PetscFunctionBegin;
140: ts->event->timestep_postevent = dt1;
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
145: /*@
146: TSSetPostEventSecondStep - Set the second time step to use after the event
148: Logically Collective
150: Input Parameters:
151: + ts - time integration context
152: - dt2 - second post event step
154: Options Database Key:
155: . -ts_event_post_event_second_step <dt2> - second time step after the event
157: Level: advanced
159: Notes:
160: `TSSetPostEventSecondStep()` allows one to set the second time step after the event.
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.
166: This function accepts either a numerical value for `dt2`, or `PETSC_DECIDE` (default).
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\:
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.
177: 2. dt1 = `PETSC_DECIDE`, dt2 = x2 (numerical). Then, after `te` TS goes to t3, and then to t3+x2.
179: 3. dt1 = x1 (numerical), dt2 = x2 (numerical). Then, after `te` TS goes to te+x1, and then to te+x1+x2.
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.
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\:
185: 1a. After `te` TS goes to t4, and event handler does not interfere to the subsequent steps.
187: 2a. After `te` TS goes to t4, and then to t4+x2.
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.
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.
197: The default value is `PETSC_DECIDE`.
199: .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`, `TSSetPostEventStep()`
200: @*/
201: PetscErrorCode TSSetPostEventSecondStep(TS ts, PetscReal dt2)
202: {
203: PetscFunctionBegin;
204: ts->event->timestep_2nd_postevent = dt2;
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: TSSetEventTolerances - Set tolerances for event (indicator function) zero crossings
211: Logically Collective
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
218: Options Database Key:
219: . -ts_event_tol <tol> - tolerance for event (indicator function) zero crossing
221: Level: beginner
223: Notes:
224: One must call `TSSetEventHandler()` before setting the tolerances.
226: The size of `vtol` should be equal to the number of events on the given process.
228: This function can be also called from the `postevent()` callback set with `TSSetEventHandler()`,
229: to adjust the tolerances on the fly.
231: .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`
232: @*/
233: PetscErrorCode TSSetEventTolerances(TS ts, PetscReal tol, PetscReal vtol[])
234: {
235: TSEvent event;
237: PetscFunctionBegin;
239: if (vtol) PetscAssertPointer(vtol, 3);
240: PetscCheck(ts->event, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set the events first by calling TSSetEventHandler()");
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: }
253: /*@C
254: TSSetEventHandler - Sets functions and parameters used for indicating events and handling them
256: Logically Collective
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)
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()`
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()`
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
296: Level: intermediate
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.
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.
307: .seealso: [](ch_ts), `TSEvent`, `TSCreate()`, `TSSetTimeStep()`, `TSSetConvergedReason()`
308: @*/
309: PetscErrorCode TSSetEventHandler(TS ts, PetscInt nevents, PetscInt direction[], PetscBool terminate[], PetscErrorCode (*indicator)(TS ts, PetscReal t, Vec U, PetscReal fvalue[], void *ctx), PetscErrorCode (*postevent)(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx), void *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
321: PetscFunctionBegin;
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;
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();
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;
384: for (PetscInt i = 0; i < event->nevents; i++) event->vtol[i] = tol;
385: if (flg) PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &event->monitor));
387: PetscCall(TSEventDestroy(&ts->event));
388: ts->event = event;
389: ts->event->refct = 1;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /*
394: TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
395: is reached.
396: */
397: static PetscErrorCode TSEventRecorderResize(TSEvent event)
398: {
399: PetscReal *time;
400: PetscInt *stepnum, *nevents;
401: PetscInt **eventidx;
402: PetscInt fact = 2;
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]));
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]));
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));
425: /* Set pointers */
426: event->recorder.time = time;
427: event->recorder.stepnum = stepnum;
428: event->recorder.nevents = nevents;
429: event->recorder.eventidx = eventidx;
431: /* Update the size */
432: event->recsize *= fact;
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /*
437: Helper routine to handle user postevents and recording
438: */
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
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: }
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, MPIU_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));
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: }
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: }
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.
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: */
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)
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
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: }
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
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: }
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: */
571: static PetscErrorCode TSEventTestZero(TS ts, PetscReal t)
572: {
573: TSEvent event = ts->event;
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;
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: }
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: */
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: }
599: /*
600: Caps the time steps, accounting for time span points.
601: It uses 'event->timestep_cache' as a time step to calculate the tolerance for tspan 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 tspan points detection, tolerance 'tspan->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: */
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 tspan
612: PetscBool cut_made = PETSC_FALSE;
613: PetscReal eps = 10 * PETSC_MACHINE_EPSILON;
614: if (ts->tspan) {
615: PetscInt ctr = ts->tspan->spanctr;
616: PetscInt Ns = ts->tspan->num_span_times;
617: PetscReal *st = ts->tspan->span_times;
619: if (ts->tspan->worktol == 0) ts->tspan->worktol = ts->tspan->reltol * ts->event->timestep_cache + ts->tspan->abstol; // in case TSAdaptChoose() has not defined it
620: if (ctr < Ns && PetscIsCloseAtTol(t, st[ctr], ts->tspan->worktol, 0)) { // just hit a time span point
621: if (ctr + 1 < Ns) maxdt = st[ctr + 1] - t; // ok to use the next time span point
622: else maxdt = ts->max_time - t; // can't use the next time span point: they have finished
623: } else if (ctr < Ns) maxdt = st[ctr] - t; // haven't hit a time span 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 in TSEvent_dt_cap()");
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_span_cached = dt;
637: else ts->adapt->dt_span_cached = 0;
638: }
639: }
640: return res;
641: }
643: /*
644: Updates the left-end values
645: */
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: }
655: /*
656: Updates the right-end values
657: */
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: }
667: /*
668: Updates the current values from the right-end values
669: */
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: }
679: static inline PetscBool Not_PETSC_DECIDE(PetscReal dt)
680: {
681: return dt == PETSC_DECIDE ? PETSC_FALSE : PETSC_TRUE;
682: }
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.
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.
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.
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:
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:
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.
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).
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).
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.
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.
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.
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.
772: Whether the algorithm is revisiting a point in the current TSEventHandler() call is flagged by 'event->revisit_right'.
773: */
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
783: PetscFunctionBegin;
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;
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: }
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
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()");
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
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: }
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
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: }
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()");
910: if (finished) { // finished handling the current event
911: PetscCall(TSPostEvent(ts, t, U));
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: }
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
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: }
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
946: PetscFunctionBegin;
948: if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
949: event = ts->event;
951: PetscCall(TSGetTime(ts, &t));
952: PetscCall(TSGetSolution(ts, &U));
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: }
966: /*@
967: TSGetNumEvents - Get the number of events defined on the given MPI process
969: Logically Collective
971: Input Parameter:
972: . ts - the `TS` context
974: Output Parameter:
975: . nevents - the number of local events on each MPI process
977: Level: intermediate
979: .seealso: [](ch_ts), `TSEvent`, `TSSetEventHandler()`
980: @*/
981: PetscErrorCode TSGetNumEvents(TS ts, PetscInt *nevents)
982: {
983: PetscFunctionBegin;
984: *nevents = ts->event->nevents;
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }