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