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);
10: PetscAssertPointer(event, 1);
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;
24: PetscAssertPointer(event, 1);
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;
116: if (vtol) PetscAssertPointer(vtol, 3);
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: 0.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
143: even has occurred
144: . postevent - [optional] post-event function, this function can change properties of the solution, ODE etc at the time of the event
145: - ctx - [optional] user-defined context for private data for the
146: event detector and post event routine (use `NULL` if no
147: context is desired)
149: Calling sequence of `eventhandler`:
150: + ts - the `TS` context
151: . t - current time
152: . U - current iterate
153: . fvalue - function value of events at time t
154: - ctx - [optional] context passed with eventhandler
156: Calling sequence of `postevent`:
157: + ts - the `TS` context
158: . nevents_zero - number of local events whose event function has been marked as crossing 0
159: . events_zero - indices of local events which have been marked as crossing 0
160: . t - current time
161: . U - current solution
162: . forwardsolve - Flag to indicate whether `TS` is doing a forward solve (1) or adjoint solve (0)
163: - ctx - the context passed with eventhandler
165: Level: intermediate
167: Note:
168: The `eventhandler` is actually the event detector function and the `postevent` function actually handles the desired changes that
169: should take place at the time of the event
171: .seealso: [](ch_ts), `TSEvent`, `TSCreate()`, `TSSetTimeStep()`, `TSSetConvergedReason()`
172: @*/
173: PetscErrorCode TSSetEventHandler(TS ts, PetscInt nevents, PetscInt direction[], PetscBool terminate[], PetscErrorCode (*eventhandler)(TS ts, PetscReal t, Vec U, PetscScalar fvalue[], void *ctx), PetscErrorCode (*postevent)(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx), void *ctx)
174: {
175: TSAdapt adapt;
176: PetscReal hmin;
177: TSEvent event;
178: PetscInt i;
179: PetscBool flg;
180: #if defined PETSC_USE_REAL_SINGLE
181: PetscReal tol = 1e-4;
182: #else
183: PetscReal tol = 1e-6;
184: #endif
186: PetscFunctionBegin;
188: if (nevents) {
189: PetscAssertPointer(direction, 3);
190: PetscAssertPointer(terminate, 4);
191: }
193: PetscCall(PetscNew(&event));
194: PetscCall(PetscMalloc1(nevents, &event->fvalue));
195: PetscCall(PetscMalloc1(nevents, &event->fvalue_prev));
196: PetscCall(PetscMalloc1(nevents, &event->fvalue_right));
197: PetscCall(PetscMalloc1(nevents, &event->zerocrossing));
198: PetscCall(PetscMalloc1(nevents, &event->side));
199: PetscCall(PetscMalloc1(nevents, &event->direction));
200: PetscCall(PetscMalloc1(nevents, &event->terminate));
201: PetscCall(PetscMalloc1(nevents, &event->vtol));
202: for (i = 0; i < nevents; i++) {
203: event->direction[i] = direction[i];
204: event->terminate[i] = terminate[i];
205: event->zerocrossing[i] = PETSC_FALSE;
206: event->side[i] = 0;
207: }
208: PetscCall(PetscMalloc1(nevents, &event->events_zero));
209: event->nevents = nevents;
210: event->eventhandler = eventhandler;
211: event->postevent = postevent;
212: event->ctx = ctx;
213: event->timestep_posteventinterval = ts->time_step;
214: PetscCall(TSGetAdapt(ts, &adapt));
215: PetscCall(TSAdaptGetStepLimits(adapt, &hmin, NULL));
216: event->timestep_min = hmin;
218: event->recsize = 8; /* Initial size of the recorder */
219: PetscOptionsBegin(((PetscObject)ts)->comm, ((PetscObject)ts)->prefix, "TS Event options", "TS");
220: {
221: PetscCall(PetscOptionsReal("-ts_event_tol", "Scalar event tolerance for zero crossing check", "TSSetEventTolerances", tol, &tol, NULL));
222: PetscCall(PetscOptionsName("-ts_event_monitor", "Print choices made by event handler", "", &flg));
223: PetscCall(PetscOptionsInt("-ts_event_recorder_initial_size", "Initial size of event recorder", "", event->recsize, &event->recsize, NULL));
224: PetscCall(PetscOptionsReal("-ts_event_post_eventinterval_step", "Time step after event interval", "", event->timestep_posteventinterval, &event->timestep_posteventinterval, NULL));
225: PetscCall(PetscOptionsReal("-ts_event_post_event_step", "Time step after event", "", event->timestep_postevent, &event->timestep_postevent, NULL));
226: PetscCall(PetscOptionsReal("-ts_event_dt_min", "Minimum time step considered for TSEvent", "", event->timestep_min, &event->timestep_min, NULL));
227: }
228: PetscOptionsEnd();
230: PetscCall(PetscMalloc1(event->recsize, &event->recorder.time));
231: PetscCall(PetscMalloc1(event->recsize, &event->recorder.stepnum));
232: PetscCall(PetscMalloc1(event->recsize, &event->recorder.nevents));
233: PetscCall(PetscMalloc1(event->recsize, &event->recorder.eventidx));
234: for (i = 0; i < event->recsize; i++) PetscCall(PetscMalloc1(event->nevents, &event->recorder.eventidx[i]));
235: /* Initialize the event recorder */
236: event->recorder.ctr = 0;
238: for (i = 0; i < event->nevents; i++) event->vtol[i] = tol;
239: if (flg) PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF, "stdout", &event->monitor));
241: PetscCall(TSEventDestroy(&ts->event));
242: ts->event = event;
243: ts->event->refct = 1;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*
248: TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
249: is reached.
250: */
251: static PetscErrorCode TSEventRecorderResize(TSEvent event)
252: {
253: PetscReal *time;
254: PetscInt *stepnum;
255: PetscInt *nevents;
256: PetscInt **eventidx;
257: PetscInt i, fact = 2;
259: PetscFunctionBegin;
261: /* Create large arrays */
262: PetscCall(PetscMalloc1(fact * event->recsize, &time));
263: PetscCall(PetscMalloc1(fact * event->recsize, &stepnum));
264: PetscCall(PetscMalloc1(fact * event->recsize, &nevents));
265: PetscCall(PetscMalloc1(fact * event->recsize, &eventidx));
266: for (i = 0; i < fact * event->recsize; i++) PetscCall(PetscMalloc1(event->nevents, &eventidx[i]));
268: /* Copy over data */
269: PetscCall(PetscArraycpy(time, event->recorder.time, event->recsize));
270: PetscCall(PetscArraycpy(stepnum, event->recorder.stepnum, event->recsize));
271: PetscCall(PetscArraycpy(nevents, event->recorder.nevents, event->recsize));
272: for (i = 0; i < event->recsize; i++) PetscCall(PetscArraycpy(eventidx[i], event->recorder.eventidx[i], event->recorder.nevents[i]));
274: /* Destroy old arrays */
275: for (i = 0; i < event->recsize; i++) PetscCall(PetscFree(event->recorder.eventidx[i]));
276: PetscCall(PetscFree(event->recorder.eventidx));
277: PetscCall(PetscFree(event->recorder.nevents));
278: PetscCall(PetscFree(event->recorder.stepnum));
279: PetscCall(PetscFree(event->recorder.time));
281: /* Set pointers */
282: event->recorder.time = time;
283: event->recorder.stepnum = stepnum;
284: event->recorder.nevents = nevents;
285: event->recorder.eventidx = eventidx;
287: /* Double size */
288: event->recsize *= fact;
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*
294: Helper routine to handle user postevents and recording
295: */
296: static PetscErrorCode TSPostEvent(TS ts, PetscReal t, Vec U)
297: {
298: TSEvent event = ts->event;
299: PetscBool terminate = PETSC_FALSE;
300: PetscBool restart = PETSC_FALSE;
301: PetscInt i, ctr, stepnum;
302: PetscBool inflag[2], outflag[2];
303: PetscBool forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */
305: PetscFunctionBegin;
306: if (event->postevent) {
307: PetscObjectState state_prev, state_post;
308: PetscCall(PetscObjectStateGet((PetscObject)U, &state_prev));
309: PetscCall((*event->postevent)(ts, event->nevents_zero, event->events_zero, t, U, forwardsolve, event->ctx));
310: PetscCall(PetscObjectStateGet((PetscObject)U, &state_post));
311: if (state_prev != state_post) restart = PETSC_TRUE;
312: }
314: /* Handle termination events and step restart */
315: for (i = 0; i < event->nevents_zero; i++)
316: if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
317: inflag[0] = restart;
318: inflag[1] = terminate;
319: PetscCall(MPIU_Allreduce(inflag, outflag, 2, MPIU_BOOL, MPI_LOR, ((PetscObject)ts)->comm));
320: restart = outflag[0];
321: terminate = outflag[1];
322: if (restart) PetscCall(TSRestartStep(ts));
323: if (terminate) PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_EVENT));
324: event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;
326: /* Reset event residual functions as states might get changed by the postevent callback */
327: if (event->postevent) {
328: PetscCall(VecLockReadPush(U));
329: PetscCall((*event->eventhandler)(ts, t, U, event->fvalue, event->ctx));
330: PetscCall(VecLockReadPop(U));
331: }
333: /* Cache current time and event residual functions */
334: event->ptime_prev = t;
335: for (i = 0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];
337: /* Record the event in the event recorder */
338: PetscCall(TSGetStepNumber(ts, &stepnum));
339: ctr = event->recorder.ctr;
340: if (ctr == event->recsize) PetscCall(TSEventRecorderResize(event));
341: event->recorder.time[ctr] = t;
342: event->recorder.stepnum[ctr] = stepnum;
343: event->recorder.nevents[ctr] = event->nevents_zero;
344: for (i = 0; i < event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
345: event->recorder.ctr++;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /* Uses Anderson-Bjorck variant of regula falsi method */
350: static inline PetscReal TSEventComputeStepSize(PetscReal tleft, PetscReal t, PetscReal tright, PetscScalar fleft, PetscScalar f, PetscScalar fright, PetscInt side, PetscReal dt)
351: {
352: PetscReal new_dt, scal = 1.0;
353: if (PetscRealPart(fleft) * PetscRealPart(f) < 0) {
354: if (side == 1) {
355: scal = (PetscRealPart(fright) - PetscRealPart(f)) / PetscRealPart(fright);
356: if (scal < PETSC_SMALL) scal = 0.5;
357: }
358: new_dt = (scal * PetscRealPart(fleft) * t - PetscRealPart(f) * tleft) / (scal * PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
359: } else {
360: if (side == -1) {
361: scal = (PetscRealPart(fleft) - PetscRealPart(f)) / PetscRealPart(fleft);
362: if (scal < PETSC_SMALL) scal = 0.5;
363: }
364: new_dt = (PetscRealPart(f) * tright - scal * PetscRealPart(fright) * t) / (PetscRealPart(f) - scal * PetscRealPart(fright)) - t;
365: }
366: return PetscMin(dt, new_dt);
367: }
369: static PetscErrorCode TSEventDetection(TS ts)
370: {
371: TSEvent event = ts->event;
372: PetscReal t;
373: PetscInt i;
374: PetscInt fvalue_sign, fvalueprev_sign;
375: PetscInt in, out;
377: PetscFunctionBegin;
378: PetscCall(TSGetTime(ts, &t));
379: for (i = 0; i < event->nevents; i++) {
380: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
381: if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
382: event->status = TSEVENT_LOCATED_INTERVAL;
383: if (event->monitor) {
384: 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));
385: }
386: continue;
387: }
388: if (PetscAbsScalar(event->fvalue_prev[i]) < event->vtol[i]) continue; /* avoid duplicative detection if the previous endpoint is an event location */
389: fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
390: fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
391: if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign)) {
392: if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
393: event->status = TSEVENT_LOCATED_INTERVAL;
394: 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));
395: }
396: }
397: in = (PetscInt)event->status;
398: PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ts)));
399: event->status = (TSEventStatus)out;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode TSEventLocation(TS ts, PetscReal *dt)
404: {
405: TSEvent event = ts->event;
406: PetscReal diff = PetscAbsReal((event->ptime_right - event->ptime_prev) / 2);
407: PetscInt i;
408: PetscReal t;
409: PetscInt fvalue_sign, fvalueprev_sign;
410: PetscInt rollback = 0, in[2], out[2];
412: PetscFunctionBegin;
413: PetscCall(TSGetTime(ts, &t));
414: event->nevents_zero = 0;
415: for (i = 0; i < event->nevents; i++) {
416: if (event->zerocrossing[i]) {
417: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i] || *dt < event->timestep_min || PetscAbsReal(*dt) < diff * event->vtol[i]) { /* stopping criteria */
418: event->status = TSEVENT_ZERO;
419: event->fvalue_right[i] = event->fvalue[i];
420: continue;
421: }
422: /* Compute new time step */
423: *dt = TSEventComputeStepSize(event->ptime_prev, t, event->ptime_right, event->fvalue_prev[i], event->fvalue[i], event->fvalue_right[i], event->side[i], *dt);
424: fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
425: fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
426: switch (event->direction[i]) {
427: case -1:
428: if (fvalue_sign < 0) {
429: rollback = 1;
430: event->fvalue_right[i] = event->fvalue[i];
431: event->side[i] = 1;
432: }
433: break;
434: case 1:
435: if (fvalue_sign > 0) {
436: rollback = 1;
437: event->fvalue_right[i] = event->fvalue[i];
438: event->side[i] = 1;
439: }
440: break;
441: case 0:
442: if (fvalue_sign != fvalueprev_sign) { /* trigger rollback only when there is a sign change */
443: rollback = 1;
444: event->fvalue_right[i] = event->fvalue[i];
445: event->side[i] = 1;
446: }
447: break;
448: }
449: if (event->status == TSEVENT_PROCESSING) event->side[i] = -1;
450: }
451: }
452: in[0] = (PetscInt)event->status;
453: in[1] = rollback;
454: PetscCall(MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ts)));
455: event->status = (TSEventStatus)out[0];
456: rollback = out[1];
457: /* 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 */
458: if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
459: if (event->status == TSEVENT_ZERO) {
460: for (i = 0; i < event->nevents; i++) {
461: if (event->zerocrossing[i]) {
462: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i] || *dt < event->timestep_min || PetscAbsReal(*dt) < diff * event->vtol[i]) { /* stopping criteria */
463: event->events_zero[event->nevents_zero++] = i;
464: 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));
465: event->zerocrossing[i] = PETSC_FALSE;
466: }
467: }
468: event->side[i] = 0;
469: }
470: }
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: PetscErrorCode TSEventHandler(TS ts)
475: {
476: TSEvent event;
477: PetscReal t;
478: Vec U;
479: PetscInt i;
480: PetscReal dt, dt_min, dt_reset = 0.0;
482: PetscFunctionBegin;
484: if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
485: event = ts->event;
487: PetscCall(TSGetTime(ts, &t));
488: PetscCall(TSGetTimeStep(ts, &dt));
489: PetscCall(TSGetSolution(ts, &U));
491: if (event->status == TSEVENT_NONE) {
492: event->timestep_prev = dt;
493: event->ptime_end = t;
494: }
495: if (event->status == TSEVENT_RESET_NEXTSTEP) {
496: /* user has specified a PostEventInterval dt */
497: dt = event->timestep_posteventinterval;
498: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
499: PetscReal maxdt = ts->max_time - t;
500: dt = dt > maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
501: }
502: PetscCall(TSSetTimeStep(ts, dt));
503: event->status = TSEVENT_NONE;
504: }
506: PetscCall(VecLockReadPush(U));
507: PetscCall((*event->eventhandler)(ts, t, U, event->fvalue, event->ctx));
508: PetscCall(VecLockReadPop(U));
510: /* Detect the events */
511: PetscCall(TSEventDetection(ts));
513: /* Locate the events */
514: if (event->status == TSEVENT_LOCATED_INTERVAL || event->status == TSEVENT_PROCESSING) {
515: /* Approach the zero crosing by setting a new step size */
516: PetscCall(TSEventLocation(ts, &dt));
517: /* Roll back when new events are detected */
518: if (event->status == TSEVENT_LOCATED_INTERVAL) {
519: PetscCall(TSRollBack(ts));
520: PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
521: event->iterctr++;
522: }
523: PetscCall(MPIU_Allreduce(&dt, &dt_min, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
524: if (dt_reset > 0.0 && dt_reset < dt_min) dt_min = dt_reset;
525: PetscCall(TSSetTimeStep(ts, dt_min));
526: /* Found the zero crossing */
527: if (event->status == TSEVENT_ZERO) {
528: PetscCall(TSPostEvent(ts, t, U));
530: dt = event->ptime_end - t;
531: if (PetscAbsReal(dt) < PETSC_SMALL) { /* we hit the event, continue with the candidate time step */
532: dt = event->timestep_prev;
533: event->status = TSEVENT_NONE;
534: }
535: if (event->timestep_postevent) { /* user has specified a PostEvent dt*/
536: dt = event->timestep_postevent;
537: }
538: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
539: PetscReal maxdt = ts->max_time - t;
540: dt = dt > maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
541: }
542: PetscCall(TSSetTimeStep(ts, dt));
543: event->iterctr = 0;
544: }
545: /* Have not found the zero crosing yet */
546: if (event->status == TSEVENT_PROCESSING) {
547: 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));
548: event->iterctr++;
549: }
550: }
551: if (event->status == TSEVENT_LOCATED_INTERVAL) { /* The step has been rolled back */
552: event->status = TSEVENT_PROCESSING;
553: event->ptime_right = t;
554: } else {
555: for (i = 0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];
556: event->ptime_prev = t;
557: }
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: PetscErrorCode TSAdjointEventHandler(TS ts)
562: {
563: TSEvent event;
564: PetscReal t;
565: Vec U;
566: PetscInt ctr;
567: PetscBool forwardsolve = PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
569: PetscFunctionBegin;
571: if (!ts->event) PetscFunctionReturn(PETSC_SUCCESS);
572: event = ts->event;
574: PetscCall(TSGetTime(ts, &t));
575: PetscCall(TSGetSolution(ts, &U));
577: ctr = event->recorder.ctr - 1;
578: if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
579: /* Call the user postevent function */
580: if (event->postevent) {
581: PetscCall((*event->postevent)(ts, event->recorder.nevents[ctr], event->recorder.eventidx[ctr], t, U, forwardsolve, event->ctx));
582: event->recorder.ctr--;
583: }
584: }
586: PetscCall(PetscBarrier((PetscObject)ts));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /*@
591: TSGetNumEvents - Get the numbers of events currently set to be detected
593: Logically Collective
595: Input Parameter:
596: . ts - the `TS` context
598: Output Parameter:
599: . nevents - the number of events
601: Level: intermediate
603: .seealso: [](ch_ts), `TSEvent`, `TSSetEventHandler()`
604: @*/
605: PetscErrorCode TSGetNumEvents(TS ts, PetscInt *nevents)
606: {
607: PetscFunctionBegin;
608: *nevents = ts->event->nevents;
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }