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, 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));

    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 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: */
    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;

    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);

    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: }

    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: }