Actual source code: tsevent.c

  1: #include <petsc/private/tsimpl.h>

  3: /*
  4:   TSEventInitialize - Initializes TSEvent for TSSolve
  5: */
  6: PetscErrorCode TSEventInitialize(TSEvent event, TS ts, PetscReal t, Vec U)
  7: {
  8:   PetscFunctionBegin;
  9:   if (!event) PetscFunctionReturn(PETSC_SUCCESS);
 13:   event->ptime_prev = t;
 14:   event->iterctr    = 0;
 15:   PetscCall((*event->eventhandler)(ts, t, U, event->fvalue_prev, event->ctx));
 16:   PetscFunctionReturn(PETSC_SUCCESS);
 17: }

 19: PetscErrorCode TSEventDestroy(TSEvent *event)
 20: {
 21:   PetscInt i;

 23:   PetscFunctionBegin;
 25:   if (!*event) PetscFunctionReturn(PETSC_SUCCESS);
 26:   if (--(*event)->refct > 0) {
 27:     *event = NULL;
 28:     PetscFunctionReturn(PETSC_SUCCESS);
 29:   }

 31:   PetscCall(PetscFree((*event)->fvalue));
 32:   PetscCall(PetscFree((*event)->fvalue_prev));
 33:   PetscCall(PetscFree((*event)->fvalue_right));
 34:   PetscCall(PetscFree((*event)->zerocrossing));
 35:   PetscCall(PetscFree((*event)->side));
 36:   PetscCall(PetscFree((*event)->direction));
 37:   PetscCall(PetscFree((*event)->terminate));
 38:   PetscCall(PetscFree((*event)->events_zero));
 39:   PetscCall(PetscFree((*event)->vtol));

 41:   for (i = 0; i < (*event)->recsize; i++) PetscCall(PetscFree((*event)->recorder.eventidx[i]));
 42:   PetscCall(PetscFree((*event)->recorder.eventidx));
 43:   PetscCall(PetscFree((*event)->recorder.nevents));
 44:   PetscCall(PetscFree((*event)->recorder.stepnum));
 45:   PetscCall(PetscFree((*event)->recorder.time));

 47:   PetscCall(PetscViewerDestroy(&(*event)->monitor));
 48:   PetscCall(PetscFree(*event));
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: /*@
 53:   TSSetPostEventIntervalStep - Set the time-step used immediately following an event interval

 55:   Logically Collective

 57:   Input Parameters:
 58: + ts - time integration context
 59: - dt - post event interval step

 61:   Options Database Key:
 62: . -ts_event_post_eventinterval_step <dt> time-step after event interval

 64:   Level: advanced

 66:   Notes:
 67:  `TSSetPostEventIntervalStep()` allows one to set a time-step that is used immediately following an event interval.

 69:   This function should be called from the postevent function set with `TSSetEventHandler()`.

 71:   The post event interval time-step should be selected based on the dynamics following the event.
 72:   If the dynamics are stiff, a conservative (small) step should be used.
 73:   If not, then a larger time-step can be used.

 75: .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`
 76: @*/
 77: PetscErrorCode TSSetPostEventIntervalStep(TS ts, PetscReal dt)
 78: {
 79:   PetscFunctionBegin;
 80:   ts->event->timestep_posteventinterval = dt;
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

 84: /*@
 85:    TSSetEventTolerances - Set tolerances for event zero crossings

 87:    Logically Collective

 89:    Input Parameters:
 90: +  ts - time integration context
 91: .  tol - scalar tolerance, `PETSC_DECIDE` to leave current value
 92: -  vtol - array of tolerances or `NULL`, used in preference to tol if present

 94:    Options Database Key:
 95: .  -ts_event_tol <tol> - tolerance for event zero crossing

 97:    Level: beginner

 99:    Notes:
100:    Must call `TSSetEventHandler(`) before setting the tolerances.

102:    The size of `vtol` is equal to the number of events.

104:    The tolerance is some measure of how close the event function is to zero for the event detector to stop
105:    and declare the time of the event has been detected.

107: .seealso: [](ch_ts), `TS`, `TSEvent`, `TSSetEventHandler()`
108: @*/
109: PetscErrorCode TSSetEventTolerances(TS ts, PetscReal tol, PetscReal vtol[])
110: {
111:   TSEvent  event;
112:   PetscInt i;

114:   PetscFunctionBegin;
117:   PetscCheck(ts->event, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set the events first by calling TSSetEventHandler()");

119:   event = ts->event;
120:   if (vtol) {
121:     for (i = 0; i < event->nevents; i++) event->vtol[i] = vtol[i];
122:   } else {
123:     if ((tol != (PetscReal)PETSC_DECIDE) || (tol != (PetscReal)PETSC_DEFAULT)) {
124:       for (i = 0; i < event->nevents; i++) event->vtol[i] = tol;
125:     }
126:   }
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*@C
131:    TSSetEventHandler - Sets a function used for detecting events

133:    Logically Collective

135:    Input Parameters:
136: +  ts - the `TS` context obtained from `TSCreate()`
137: .  nevents - number of local events
138: .  direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction,
139:                +1 => Zero crossing in positive direction, 0 => both ways (one for each event)
140: .  terminate - flag to indicate whether time stepping should be terminated after
141:                event is detected (one for each event)
142: .  eventhandler - a change in sign of this function (see `direction`) is used to determine an even has occurred
143: .  postevent - [optional] post-event function, this function can change properties of the solution, ODE etc at the time of the event
144: -  ctx       - [optional] user-defined context for private data for the
145:                event detector and post event routine (use `NULL` if no
146:                context is desired)

148:    Calling sequence of `eventhandler`:
149: $   PetscErrorCode eventhandler(TS ts, PetscReal t, Vec U, PetscScalar fvalue[], void* ctx)
150: +  ts  - the `TS` context
151: .  t   - current time
152: .  U   - current iterate
153: .  ctx - [optional] context passed with eventhandler
154: -  fvalue    - function value of events at time t

156:    Calling sequence of `postevent`:
157: $   PetscErrorCode postevent(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
158: +  ts - the `TS` context
159: .  nevents_zero - number of local events whose event function is zero
160: .  events_zero  - indices of local events which have reached zero
161: .  t            - current time
162: .  U            - current solution
163: .  forwardsolve - Flag to indicate whether `TS` is doing a forward solve (1) or adjoint solve (0)
164: -  ctx          - the context passed with eventhandler

166:    Level: intermediate

168:    Note:
169:    The `eventhandler` is actually the event detector function and the `postevent` function actually handles the desired changes that
170:    should take place at the time of the event

172: .seealso: [](ch_ts), `TSEvent`, `TSCreate()`, `TSSetTimeStep()`, `TSSetConvergedReason()`
173: @*/
174: PetscErrorCode TSSetEventHandler(TS ts, PetscInt nevents, PetscInt direction[], PetscBool terminate[], PetscErrorCode (*eventhandler)(TS, PetscReal, Vec, PetscScalar[], void *), PetscErrorCode (*postevent)(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *), void *ctx)
175: {
176:   TSAdapt   adapt;
177:   PetscReal hmin;
178:   TSEvent   event;
179:   PetscInt  i;
180:   PetscBool flg;
181: #if defined PETSC_USE_REAL_SINGLE
182:   PetscReal tol = 1e-4;
183: #else
184:   PetscReal tol = 1e-6;
185: #endif

187:   PetscFunctionBegin;
189:   if (nevents) {
192:   }

194:   PetscCall(PetscNew(&event));
195:   PetscCall(PetscMalloc1(nevents, &event->fvalue));
196:   PetscCall(PetscMalloc1(nevents, &event->fvalue_prev));
197:   PetscCall(PetscMalloc1(nevents, &event->fvalue_right));
198:   PetscCall(PetscMalloc1(nevents, &event->zerocrossing));
199:   PetscCall(PetscMalloc1(nevents, &event->side));
200:   PetscCall(PetscMalloc1(nevents, &event->direction));
201:   PetscCall(PetscMalloc1(nevents, &event->terminate));
202:   PetscCall(PetscMalloc1(nevents, &event->vtol));
203:   for (i = 0; i < nevents; i++) {
204:     event->direction[i]    = direction[i];
205:     event->terminate[i]    = terminate[i];
206:     event->zerocrossing[i] = PETSC_FALSE;
207:     event->side[i]         = 0;
208:   }
209:   PetscCall(PetscMalloc1(nevents, &event->events_zero));
210:   event->nevents                    = nevents;
211:   event->eventhandler               = eventhandler;
212:   event->postevent                  = postevent;
213:   event->ctx                        = ctx;
214:   event->timestep_posteventinterval = ts->time_step;
215:   PetscCall(TSGetAdapt(ts, &adapt));
216:   PetscCall(TSAdaptGetStepLimits(adapt, &hmin, NULL));
217:   event->timestep_min = hmin;

219:   event->recsize = 8; /* Initial size of the recorder */
220:   PetscOptionsBegin(((PetscObject)ts)->comm, ((PetscObject)ts)->prefix, "TS Event options", "TS");
221:   {
222:     PetscCall(PetscOptionsReal("-ts_event_tol", "Scalar event tolerance for zero crossing check", "TSSetEventTolerances", tol, &tol, NULL));
223:     PetscCall(PetscOptionsName("-ts_event_monitor", "Print choices made by event handler", "", &flg));
224:     PetscCall(PetscOptionsInt("-ts_event_recorder_initial_size", "Initial size of event recorder", "", event->recsize, &event->recsize, NULL));
225:     PetscCall(PetscOptionsReal("-ts_event_post_eventinterval_step", "Time step after event interval", "", event->timestep_posteventinterval, &event->timestep_posteventinterval, NULL));
226:     PetscCall(PetscOptionsReal("-ts_event_post_event_step", "Time step after event", "", event->timestep_postevent, &event->timestep_postevent, NULL));
227:     PetscCall(PetscOptionsReal("-ts_event_dt_min", "Minimum time step considered for TSEvent", "", event->timestep_min, &event->timestep_min, NULL));
228:   }
229:   PetscOptionsEnd();

231:   PetscCall(PetscMalloc1(event->recsize, &event->recorder.time));
232:   PetscCall(PetscMalloc1(event->recsize, &event->recorder.stepnum));
233:   PetscCall(PetscMalloc1(event->recsize, &event->recorder.nevents));
234:   PetscCall(PetscMalloc1(event->recsize, &event->recorder.eventidx));
235:   for (i = 0; i < event->recsize; i++) PetscCall(PetscMalloc1(event->nevents, &event->recorder.eventidx[i]));
236:   /* Initialize the event recorder */
237:   event->recorder.ctr = 0;

239:   for (i = 0; i < event->nevents; i++) event->vtol[i] = tol;
240:   if (flg) PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF, "stdout", &event->monitor));

242:   PetscCall(TSEventDestroy(&ts->event));
243:   ts->event        = event;
244:   ts->event->refct = 1;
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: /*
249:   TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
250:                           is reached.
251: */
252: static PetscErrorCode TSEventRecorderResize(TSEvent event)
253: {
254:   PetscReal *time;
255:   PetscInt  *stepnum;
256:   PetscInt  *nevents;
257:   PetscInt **eventidx;
258:   PetscInt   i, fact = 2;

260:   PetscFunctionBegin;

262:   /* Create large arrays */
263:   PetscCall(PetscMalloc1(fact * event->recsize, &time));
264:   PetscCall(PetscMalloc1(fact * event->recsize, &stepnum));
265:   PetscCall(PetscMalloc1(fact * event->recsize, &nevents));
266:   PetscCall(PetscMalloc1(fact * event->recsize, &eventidx));
267:   for (i = 0; i < fact * event->recsize; i++) PetscCall(PetscMalloc1(event->nevents, &eventidx[i]));

269:   /* Copy over data */
270:   PetscCall(PetscArraycpy(time, event->recorder.time, event->recsize));
271:   PetscCall(PetscArraycpy(stepnum, event->recorder.stepnum, event->recsize));
272:   PetscCall(PetscArraycpy(nevents, event->recorder.nevents, event->recsize));
273:   for (i = 0; i < event->recsize; i++) PetscCall(PetscArraycpy(eventidx[i], event->recorder.eventidx[i], event->recorder.nevents[i]));

275:   /* Destroy old arrays */
276:   for (i = 0; i < event->recsize; i++) PetscCall(PetscFree(event->recorder.eventidx[i]));
277:   PetscCall(PetscFree(event->recorder.eventidx));
278:   PetscCall(PetscFree(event->recorder.nevents));
279:   PetscCall(PetscFree(event->recorder.stepnum));
280:   PetscCall(PetscFree(event->recorder.time));

282:   /* Set pointers */
283:   event->recorder.time     = time;
284:   event->recorder.stepnum  = stepnum;
285:   event->recorder.nevents  = nevents;
286:   event->recorder.eventidx = eventidx;

288:   /* Double size */
289:   event->recsize *= fact;

291:   PetscFunctionReturn(PETSC_SUCCESS);
292: }

294: /*
295:    Helper routine to handle user postevents and recording
296: */
297: static PetscErrorCode TSPostEvent(TS ts, PetscReal t, Vec U)
298: {
299:   TSEvent   event     = ts->event;
300:   PetscBool terminate = PETSC_FALSE;
301:   PetscBool restart   = PETSC_FALSE;
302:   PetscInt  i, ctr, stepnum;
303:   PetscBool inflag[2], outflag[2];
304:   PetscBool forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */

306:   PetscFunctionBegin;
307:   if (event->postevent) {
308:     PetscObjectState state_prev, state_post;
309:     PetscCall(PetscObjectStateGet((PetscObject)U, &state_prev));
310:     PetscCall((*event->postevent)(ts, event->nevents_zero, event->events_zero, t, U, forwardsolve, event->ctx));
311:     PetscCall(PetscObjectStateGet((PetscObject)U, &state_post));
312:     if (state_prev != state_post) restart = PETSC_TRUE;
313:   }

315:   /* Handle termination events and step restart */
316:   for (i = 0; i < event->nevents_zero; i++)
317:     if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
318:   inflag[0] = restart;
319:   inflag[1] = terminate;
320:   PetscCall(MPIU_Allreduce(inflag, outflag, 2, MPIU_BOOL, MPI_LOR, ((PetscObject)ts)->comm));
321:   restart   = outflag[0];
322:   terminate = outflag[1];
323:   if (restart) PetscCall(TSRestartStep(ts));
324:   if (terminate) PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_EVENT));
325:   event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;

327:   /* Reset event residual functions as states might get changed by the postevent callback */
328:   if (event->postevent) {
329:     PetscCall(VecLockReadPush(U));
330:     PetscCall((*event->eventhandler)(ts, t, U, event->fvalue, event->ctx));
331:     PetscCall(VecLockReadPop(U));
332:   }

334:   /* Cache current time and event residual functions */
335:   event->ptime_prev = t;
336:   for (i = 0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];

338:   /* Record the event in the event recorder */
339:   PetscCall(TSGetStepNumber(ts, &stepnum));
340:   ctr = event->recorder.ctr;
341:   if (ctr == event->recsize) PetscCall(TSEventRecorderResize(event));
342:   event->recorder.time[ctr]    = t;
343:   event->recorder.stepnum[ctr] = stepnum;
344:   event->recorder.nevents[ctr] = event->nevents_zero;
345:   for (i = 0; i < event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
346:   event->recorder.ctr++;
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: /* Uses Anderson-Bjorck variant of regula falsi method */
351: static inline PetscReal TSEventComputeStepSize(PetscReal tleft, PetscReal t, PetscReal tright, PetscScalar fleft, PetscScalar f, PetscScalar fright, PetscInt side, PetscReal dt)
352: {
353:   PetscReal new_dt, scal = 1.0;
354:   if (PetscRealPart(fleft) * PetscRealPart(f) < 0) {
355:     if (side == 1) {
356:       scal = (PetscRealPart(fright) - PetscRealPart(f)) / PetscRealPart(fright);
357:       if (scal < PETSC_SMALL) scal = 0.5;
358:     }
359:     new_dt = (scal * PetscRealPart(fleft) * t - PetscRealPart(f) * tleft) / (scal * PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
360:   } else {
361:     if (side == -1) {
362:       scal = (PetscRealPart(fleft) - PetscRealPart(f)) / PetscRealPart(fleft);
363:       if (scal < PETSC_SMALL) scal = 0.5;
364:     }
365:     new_dt = (PetscRealPart(f) * tright - scal * PetscRealPart(fright) * t) / (PetscRealPart(f) - scal * PetscRealPart(fright)) - t;
366:   }
367:   return PetscMin(dt, new_dt);
368: }

370: static PetscErrorCode TSEventDetection(TS ts)
371: {
372:   TSEvent   event = ts->event;
373:   PetscReal t;
374:   PetscInt  i;
375:   PetscInt  fvalue_sign, fvalueprev_sign;
376:   PetscInt  in, out;

378:   PetscFunctionBegin;
379:   PetscCall(TSGetTime(ts, &t));
380:   for (i = 0; i < event->nevents; i++) {
381:     if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
382:       if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
383:       event->status = TSEVENT_LOCATED_INTERVAL;
384:       if (event->monitor) {
385:         PetscCall(PetscViewerASCIIPrintf(event->monitor, "TSEvent: iter %" PetscInt_FMT " - Event %" PetscInt_FMT " interval detected due to zero value (tol=%g) [%g - %g]\n", event->iterctr, i, (double)event->vtol[i], (double)event->ptime_prev, (double)t));
386:       }
387:       continue;
388:     }
389:     if (PetscAbsScalar(event->fvalue_prev[i]) < event->vtol[i]) continue; /* avoid duplicative detection if the previous endpoint is an event location */
390:     fvalue_sign     = PetscSign(PetscRealPart(event->fvalue[i]));
391:     fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
392:     if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign)) {
393:       if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
394:       event->status = TSEVENT_LOCATED_INTERVAL;
395:       if (event->monitor) PetscCall(PetscViewerASCIIPrintf(event->monitor, "TSEvent: iter %" PetscInt_FMT " - Event %" PetscInt_FMT " interval detected due to sign change [%g - %g]\n", event->iterctr, i, (double)event->ptime_prev, (double)t));
396:     }
397:   }
398:   in = (PetscInt)event->status;
399:   PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ts)));
400:   event->status = (TSEventStatus)out;
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: static PetscErrorCode TSEventLocation(TS ts, PetscReal *dt)
405: {
406:   TSEvent   event = ts->event;
407:   PetscReal diff  = PetscAbsReal((event->ptime_right - event->ptime_prev) / 2);
408:   PetscInt  i;
409:   PetscReal t;
410:   PetscInt  fvalue_sign, fvalueprev_sign;
411:   PetscInt  rollback = 0, in[2], out[2];

413:   PetscFunctionBegin;
414:   PetscCall(TSGetTime(ts, &t));
415:   event->nevents_zero = 0;
416:   for (i = 0; i < event->nevents; i++) {
417:     if (event->zerocrossing[i]) {
418:       if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i] || *dt < event->timestep_min || PetscAbsReal(*dt) < diff * event->vtol[i]) { /* stopping criteria */
419:         event->status          = TSEVENT_ZERO;
420:         event->fvalue_right[i] = event->fvalue[i];
421:         continue;
422:       }
423:       /* Compute new time step */
424:       *dt             = TSEventComputeStepSize(event->ptime_prev, t, event->ptime_right, event->fvalue_prev[i], event->fvalue[i], event->fvalue_right[i], event->side[i], *dt);
425:       fvalue_sign     = PetscSign(PetscRealPart(event->fvalue[i]));
426:       fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
427:       switch (event->direction[i]) {
428:       case -1:
429:         if (fvalue_sign < 0) {
430:           rollback               = 1;
431:           event->fvalue_right[i] = event->fvalue[i];
432:           event->side[i]         = 1;
433:         }
434:         break;
435:       case 1:
436:         if (fvalue_sign > 0) {
437:           rollback               = 1;
438:           event->fvalue_right[i] = event->fvalue[i];
439:           event->side[i]         = 1;
440:         }
441:         break;
442:       case 0:
443:         if (fvalue_sign != fvalueprev_sign) { /* trigger rollback only when there is a sign change */
444:           rollback               = 1;
445:           event->fvalue_right[i] = event->fvalue[i];
446:           event->side[i]         = 1;
447:         }
448:         break;
449:       }
450:       if (event->status == TSEVENT_PROCESSING) event->side[i] = -1;
451:     }
452:   }
453:   in[0] = (PetscInt)event->status;
454:   in[1] = rollback;
455:   PetscCall(MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ts)));
456:   event->status = (TSEventStatus)out[0];
457:   rollback      = out[1];
458:   /* If rollback is true, the status will be overwritten so that an event at the endtime of current time step will be postponed to guarantee correct order */
459:   if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
460:   if (event->status == TSEVENT_ZERO) {
461:     for (i = 0; i < event->nevents; i++) {
462:       if (event->zerocrossing[i]) {
463:         if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i] || *dt < event->timestep_min || PetscAbsReal(*dt) < diff * event->vtol[i]) { /* stopping criteria */
464:           event->events_zero[event->nevents_zero++] = i;
465:           if (event->monitor) PetscCall(PetscViewerASCIIPrintf(event->monitor, "TSEvent: iter %" PetscInt_FMT " - Event %" PetscInt_FMT " zero crossing located at time %g\n", event->iterctr, i, (double)t));
466:           event->zerocrossing[i] = PETSC_FALSE;
467:         }
468:       }
469:       event->side[i] = 0;
470:     }
471:   }
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: PetscErrorCode TSEventHandler(TS ts)
476: {
477:   TSEvent   event;
478:   PetscReal t;
479:   Vec       U;
480:   PetscInt  i;
481:   PetscReal dt, dt_min, dt_reset = 0.0;

483:   PetscFunctionBegin;
485:   if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
486:   event = ts->event;

488:   PetscCall(TSGetTime(ts, &t));
489:   PetscCall(TSGetTimeStep(ts, &dt));
490:   PetscCall(TSGetSolution(ts, &U));

492:   if (event->status == TSEVENT_NONE) {
493:     event->timestep_prev = dt;
494:     event->ptime_end     = t;
495:   }
496:   if (event->status == TSEVENT_RESET_NEXTSTEP) {
497:     /* user has specified a PostEventInterval dt */
498:     dt = event->timestep_posteventinterval;
499:     if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
500:       PetscReal maxdt = ts->max_time - t;
501:       dt              = dt > maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
502:     }
503:     PetscCall(TSSetTimeStep(ts, dt));
504:     event->status = TSEVENT_NONE;
505:   }

507:   PetscCall(VecLockReadPush(U));
508:   PetscCall((*event->eventhandler)(ts, t, U, event->fvalue, event->ctx));
509:   PetscCall(VecLockReadPop(U));

511:   /* Detect the events */
512:   PetscCall(TSEventDetection(ts));

514:   /* Locate the events */
515:   if (event->status == TSEVENT_LOCATED_INTERVAL || event->status == TSEVENT_PROCESSING) {
516:     /* Approach the zero crosing by setting a new step size */
517:     PetscCall(TSEventLocation(ts, &dt));
518:     /* Roll back when new events are detected */
519:     if (event->status == TSEVENT_LOCATED_INTERVAL) {
520:       PetscCall(TSRollBack(ts));
521:       PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
522:       event->iterctr++;
523:     }
524:     PetscCall(MPIU_Allreduce(&dt, &dt_min, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
525:     if (dt_reset > 0.0 && dt_reset < dt_min) dt_min = dt_reset;
526:     PetscCall(TSSetTimeStep(ts, dt_min));
527:     /* Found the zero crossing */
528:     if (event->status == TSEVENT_ZERO) {
529:       PetscCall(TSPostEvent(ts, t, U));

531:       dt = event->ptime_end - t;
532:       if (PetscAbsReal(dt) < PETSC_SMALL) { /* we hit the event, continue with the candidate time step */
533:         dt            = event->timestep_prev;
534:         event->status = TSEVENT_NONE;
535:       }
536:       if (event->timestep_postevent) { /* user has specified a PostEvent dt*/
537:         dt = event->timestep_postevent;
538:       }
539:       if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
540:         PetscReal maxdt = ts->max_time - t;
541:         dt              = dt > maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
542:       }
543:       PetscCall(TSSetTimeStep(ts, dt));
544:       event->iterctr = 0;
545:     }
546:     /* Have not found the zero crosing yet */
547:     if (event->status == TSEVENT_PROCESSING) {
548:       if (event->monitor) PetscCall(PetscViewerASCIIPrintf(event->monitor, "TSEvent: iter %" PetscInt_FMT " - Stepping forward as no event detected in interval [%g - %g]\n", event->iterctr, (double)event->ptime_prev, (double)t));
549:       event->iterctr++;
550:     }
551:   }
552:   if (event->status == TSEVENT_LOCATED_INTERVAL) { /* The step has been rolled back */
553:     event->status      = TSEVENT_PROCESSING;
554:     event->ptime_right = t;
555:   } else {
556:     for (i = 0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];
557:     event->ptime_prev = t;
558:   }
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: PetscErrorCode TSAdjointEventHandler(TS ts)
563: {
564:   TSEvent   event;
565:   PetscReal t;
566:   Vec       U;
567:   PetscInt  ctr;
568:   PetscBool forwardsolve = PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */

570:   PetscFunctionBegin;
572:   if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
573:   event = ts->event;

575:   PetscCall(TSGetTime(ts, &t));
576:   PetscCall(TSGetSolution(ts, &U));

578:   ctr = event->recorder.ctr - 1;
579:   if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
580:     /* Call the user postevent function */
581:     if (event->postevent) {
582:       PetscCall((*event->postevent)(ts, event->recorder.nevents[ctr], event->recorder.eventidx[ctr], t, U, forwardsolve, event->ctx));
583:       event->recorder.ctr--;
584:     }
585:   }

587:   PetscCall(PetscBarrier((PetscObject)ts));
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

591: /*@
592:   TSGetNumEvents - Get the numbers of events currently set to be detected

594:   Logically Collective

596:   Input Parameter:
597: . ts - the `TS` context

599:   Output Parameter:
600: . nevents - the number of events

602:   Level: intermediate

604: .seealso: [](ch_ts), `TSEvent`, `TSSetEventHandler()`
605: @*/
606: PetscErrorCode TSGetNumEvents(TS ts, PetscInt *nevents)
607: {
608:   PetscFunctionBegin;
609:   *nevents = ts->event->nevents;
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }