Actual source code: tsevent.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
5: /*
6: TSEventInitialize - Initializes TSEvent for TSSolve
7: */
8: PetscErrorCode TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U)
9: {
13: if (!event) return(0);
17: event->ptime_prev = t;
18: event->iterctr = 0;
19: (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);
20: /* Initialize the event recorder */
21: event->recorder.ctr = 0;
22: return(0);
23: }
27: PetscErrorCode TSEventDestroy(TSEvent *event)
28: {
30: PetscInt i;
34: if (!*event) return(0);
35: PetscFree((*event)->fvalue);
36: PetscFree((*event)->fvalue_prev);
37: PetscFree((*event)->fvalue_right);
38: PetscFree((*event)->zerocrossing);
39: PetscFree((*event)->side);
40: PetscFree((*event)->direction);
41: PetscFree((*event)->terminate);
42: PetscFree((*event)->events_zero);
43: PetscFree((*event)->vtol);
45: for (i=0; i < (*event)->recsize; i++) {
46: PetscFree((*event)->recorder.eventidx[i]);
47: }
48: PetscFree((*event)->recorder.eventidx);
49: PetscFree((*event)->recorder.nevents);
50: PetscFree((*event)->recorder.stepnum);
51: PetscFree((*event)->recorder.time);
53: PetscViewerDestroy(&(*event)->monitor);
54: PetscFree(*event);
55: return(0);
56: }
60: /*@
61: TSSetEventTolerances - Set tolerances for event zero crossings when using event handler
63: Logically Collective
65: Input Arguments:
66: + ts - time integration context
67: . tol - scalar tolerance, PETSC_DECIDE to leave current value
68: - vtol - array of tolerances or NULL, used in preference to tol if present
70: Options Database Keys:
71: . -ts_event_tol <tol> tolerance for event zero crossing
73: Notes:
74: Must call TSSetEventHandler() before setting the tolerances.
76: The size of vtol is equal to the number of events.
78: Level: beginner
80: .seealso: TS, TSEvent, TSSetEventHandler()
81: @*/
82: PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[])
83: {
84: TSEvent event;
85: PetscInt i;
90: if (!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()");
92: event = ts->event;
93: if (vtol) {
94: for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i];
95: } else {
96: if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) {
97: for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
98: }
99: }
100: return(0);
101: }
105: /*@C
106: TSSetEventHandler - Sets a monitoring function used for detecting events
108: Logically Collective on TS
110: Input Parameters:
111: + ts - the TS context obtained from TSCreate()
112: . nevents - number of local events
113: . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction,
114: +1 => Zero crossing in positive direction, 0 => both ways (one for each event)
115: . terminate - flag to indicate whether time stepping should be terminated after
116: event is detected (one for each event)
117: . eventhandler - event monitoring routine
118: . postevent - [optional] post-event function
119: - ctx - [optional] user-defined context for private data for the
120: event monitor and post event routine (use NULL if no
121: context is desired)
123: Calling sequence of eventhandler:
124: PetscErrorCode EventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx)
126: Input Parameters:
127: + ts - the TS context
128: . t - current time
129: . U - current iterate
130: - ctx - [optional] context passed with eventhandler
132: Output parameters:
133: . fvalue - function value of events at time t
135: Calling sequence of postevent:
136: PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
138: Input Parameters:
139: + ts - the TS context
140: . nevents_zero - number of local events whose event function is zero
141: . events_zero - indices of local events which have reached zero
142: . t - current time
143: . U - current solution
144: . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0)
145: - ctx - the context passed with eventhandler
147: Level: intermediate
149: .keywords: TS, event, set
151: .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason()
152: @*/
153: 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)
154: {
156: TSEvent event;
157: PetscInt i;
158: PetscBool flg;
159: #if defined PETSC_USE_REAL_SINGLE
160: PetscReal tol=1e-4;
161: #else
162: PetscReal tol=1e-6;
163: #endif
170: PetscNewLog(ts,&event);
171: PetscMalloc1(nevents,&event->fvalue);
172: PetscMalloc1(nevents,&event->fvalue_prev);
173: PetscMalloc1(nevents,&event->fvalue_right);
174: PetscMalloc1(nevents,&event->zerocrossing);
175: PetscMalloc1(nevents,&event->side);
176: PetscMalloc1(nevents,&event->direction);
177: PetscMalloc1(nevents,&event->terminate);
178: PetscMalloc1(nevents,&event->vtol);
179: for (i=0; i < nevents; i++) {
180: event->direction[i] = direction[i];
181: event->terminate[i] = terminate[i];
182: event->zerocrossing[i] = PETSC_FALSE;
183: event->side[i] = 0;
184: }
185: PetscMalloc1(nevents,&event->events_zero);
186: event->nevents = nevents;
187: event->eventhandler = eventhandler;
188: event->postevent = postevent;
189: event->ctx = ctx;
191: event->recsize = 8; /* Initial size of the recorder */
192: PetscOptionsBegin(((PetscObject)ts)->comm,"","TS Event options","");
193: {
194: PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","TSSetEventTolerances",tol,&tol,NULL);
195: PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);
196: PetscOptionsInt("-ts_event_recorder_initial_size","Initial size of event recorder","",event->recsize,&event->recsize,NULL);
197: }
198: PetscOptionsEnd();
200: PetscMalloc1(event->recsize,&event->recorder.time);
201: PetscMalloc1(event->recsize,&event->recorder.stepnum);
202: PetscMalloc1(event->recsize,&event->recorder.nevents);
203: PetscMalloc1(event->recsize,&event->recorder.eventidx);
204: for (i=0; i < event->recsize; i++) {
205: PetscMalloc1(event->nevents,&event->recorder.eventidx[i]);
206: }
208: for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
209: if (flg) {PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);}
211: TSEventDestroy(&ts->event);
212: ts->event = event;
213: return(0);
214: }
218: /*
219: TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
220: is reached.
221: */
222: static PetscErrorCode TSEventRecorderResize(TSEvent event)
223: {
225: PetscReal *time;
226: PetscInt *stepnum;
227: PetscInt *nevents;
228: PetscInt **eventidx;
229: PetscInt i,fact=2;
233: /* Create large arrays */
234: PetscMalloc1(fact*event->recsize,&time);
235: PetscMalloc1(fact*event->recsize,&stepnum);
236: PetscMalloc1(fact*event->recsize,&nevents);
237: PetscMalloc1(fact*event->recsize,&eventidx);
238: for (i=0; i < fact*event->recsize; i++) {
239: PetscMalloc1(event->nevents,&eventidx[i]);
240: }
242: /* Copy over data */
243: PetscMemcpy(time,event->recorder.time,event->recsize*sizeof(PetscReal));
244: PetscMemcpy(stepnum,event->recorder.stepnum,event->recsize*sizeof(PetscInt));
245: PetscMemcpy(nevents,event->recorder.nevents,event->recsize*sizeof(PetscInt));
246: for (i=0; i < event->recsize; i++) {
247: PetscMemcpy(eventidx[i],event->recorder.eventidx[i],event->recorder.nevents[i]*sizeof(PetscInt));
248: }
250: /* Destroy old arrays */
251: for (i=0; i < event->recsize; i++) {
252: PetscFree(event->recorder.eventidx[i]);
253: }
254: PetscFree(event->recorder.eventidx);
255: PetscFree(event->recorder.nevents);
256: PetscFree(event->recorder.stepnum);
257: PetscFree(event->recorder.time);
259: /* Set pointers */
260: event->recorder.time = time;
261: event->recorder.stepnum = stepnum;
262: event->recorder.nevents = nevents;
263: event->recorder.eventidx = eventidx;
265: /* Double size */
266: event->recsize *= fact;
268: return(0);
269: }
271: /*
272: Helper rutine to handle user postenvents and recording
273: */
276: static PetscErrorCode TSPostEvent(TS ts,PetscReal t,Vec U)
277: {
279: TSEvent event = ts->event;
280: PetscBool terminate = PETSC_FALSE;
281: PetscBool restart = PETSC_FALSE;
282: PetscInt i,ctr,stepnum;
283: PetscBool ts_terminate,ts_restart;
284: PetscBool forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */
287: if (event->postevent) {
288: PetscObjectState state_prev,state_post;
289: PetscObjectStateGet((PetscObject)U,&state_prev);
290: (*event->postevent)(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->ctx);
291: PetscObjectStateGet((PetscObject)U,&state_post);
292: if (state_prev != state_post) restart = PETSC_TRUE;
293: }
295: /* Handle termination events and step restart */
296: for (i=0; i<event->nevents_zero; i++) if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
297: MPIU_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);
298: if (ts_terminate) {TSSetConvergedReason(ts,TS_CONVERGED_EVENT);}
299: event->status = ts_terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;
300: MPIU_Allreduce(&restart,&ts_restart,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);
301: if (ts_restart) ts->steprestart = PETSC_TRUE;
303: event->ptime_prev = t;
304: /* Reset event residual functions as states might get changed by the postevent callback */
305: if (event->postevent) {(*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);}
306: /* Cache current event residual functions */
307: for (i=0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];
309: /* Record the event in the event recorder */
310: TSGetTimeStepNumber(ts,&stepnum);
311: ctr = event->recorder.ctr;
312: if (ctr == event->recsize) {
313: TSEventRecorderResize(event);
314: }
315: event->recorder.time[ctr] = t;
316: event->recorder.stepnum[ctr] = stepnum;
317: event->recorder.nevents[ctr] = event->nevents_zero;
318: for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
319: event->recorder.ctr++;
320: return(0);
321: }
323: /* Uses Anderson-Bjorck variant of regula falsi method */
324: PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt)
325: {
326: PetscReal new_dt, scal = 1.0;
327: if (PetscRealPart(fleft)*PetscRealPart(f) < 0) {
328: if (side == 1) {
329: scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright);
330: if (scal < PETSC_SMALL) scal = 0.5;
331: }
332: new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
333: } else {
334: if (side == -1) {
335: scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft);
336: if (scal < PETSC_SMALL) scal = 0.5;
337: }
338: new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t;
339: }
340: return PetscMin(dt,new_dt);
341: }
346: PetscErrorCode TSEventHandler(TS ts)
347: {
349: TSEvent event;
350: PetscReal t;
351: Vec U;
352: PetscInt i;
353: PetscReal dt,dt_min;
354: PetscInt rollback=0,in[2],out[2];
355: PetscInt fvalue_sign,fvalueprev_sign;
359: if (!ts->event) return(0);
360: event = ts->event;
362: TSGetTime(ts,&t);
363: TSGetTimeStep(ts,&dt);
364: TSGetSolution(ts,&U);
366: if (event->status == TSEVENT_NONE) {
367: if (ts->steps == 1) /* After first successful step */
368: event->timestep_orig = ts->ptime - ts->ptime_prev;
369: event->timestep_prev = dt;
370: }
372: if (event->status == TSEVENT_RESET_NEXTSTEP) {
373: /* Restore time step */
374: dt = PetscMin(event->timestep_orig,event->timestep_prev);
375: TSSetTimeStep(ts,dt);
376: event->status = TSEVENT_NONE;
377: }
379: if (event->status == TSEVENT_NONE) {
380: event->ptime_end = t;
381: }
383: (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);
385: for (i=0; i < event->nevents; i++) {
386: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
387: event->status = TSEVENT_ZERO;
388: event->fvalue_right[i] = event->fvalue[i];
389: continue;
390: }
391: fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
392: fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
393: if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) {
394: switch (event->direction[i]) {
395: case -1:
396: if (fvalue_sign < 0) {
397: rollback = 1;
399: /* Compute new time step */
400: dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
402: if (event->monitor) {
403: PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);
404: }
405: event->fvalue_right[i] = event->fvalue[i];
406: event->side[i] = 1;
408: if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
409: event->status = TSEVENT_LOCATED_INTERVAL;
410: }
411: break;
412: case 1:
413: if (fvalue_sign > 0) {
414: rollback = 1;
416: /* Compute new time step */
417: dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
419: if (event->monitor) {
420: PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);
421: }
422: event->fvalue_right[i] = event->fvalue[i];
423: event->side[i] = 1;
425: if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
426: event->status = TSEVENT_LOCATED_INTERVAL;
427: }
428: break;
429: case 0:
430: rollback = 1;
432: /* Compute new time step */
433: dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
435: if (event->monitor) {
436: PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);
437: }
438: event->fvalue_right[i] = event->fvalue[i];
439: event->side[i] = 1;
441: if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
442: event->status = TSEVENT_LOCATED_INTERVAL;
444: break;
445: }
446: }
447: }
449: in[0] = event->status; in[1] = rollback;
450: MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));
451: event->status = (TSEventStatus)out[0]; rollback = out[1];
452: if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
454: event->nevents_zero = 0;
455: if (event->status == TSEVENT_ZERO) {
456: for (i=0; i < event->nevents; i++) {
457: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
458: event->events_zero[event->nevents_zero++] = i;
459: if (event->monitor) {
460: PetscViewerASCIIPrintf(event->monitor,"TSEvent: Event %D zero crossing at time %g located in %D iterations\n",i,(double)t,event->iterctr);
461: }
462: event->zerocrossing[i] = PETSC_FALSE;
463: }
464: event->side[i] = 0;
465: }
466: TSPostEvent(ts,t,U);
468: dt = event->ptime_end - t;
469: if (PetscAbsReal(dt) < PETSC_SMALL) dt += PetscMin(event->timestep_orig,event->timestep_prev); /* XXX Should be done better */
470: TSSetTimeStep(ts,dt);
471: event->iterctr = 0;
472: return(0);
473: }
475: if (event->status == TSEVENT_LOCATED_INTERVAL) {
476: TSRollBack(ts);
477: TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);
478: event->status = TSEVENT_PROCESSING;
479: event->ptime_right = t;
480: } else {
481: for (i=0; i < event->nevents; i++) {
482: if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) {
483: /* Compute new time step */
484: dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
485: event->side[i] = -1;
486: }
487: event->fvalue_prev[i] = event->fvalue[i];
488: }
489: if (event->monitor && event->status == TSEVENT_PROCESSING) {
490: PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Stepping forward as no event detected in interval [%g - %g]\n",event->iterctr,(double)event->ptime_prev,(double)t);
491: }
492: event->ptime_prev = t;
493: }
495: if (event->status == TSEVENT_PROCESSING) event->iterctr++;
497: MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));
498: TSSetTimeStep(ts,dt_min);
499: return(0);
500: }
504: PetscErrorCode TSAdjointEventHandler(TS ts)
505: {
507: TSEvent event;
508: PetscReal t;
509: Vec U;
510: PetscInt ctr;
511: PetscBool forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
515: if (!ts->event) return(0);
516: event = ts->event;
518: TSGetTime(ts,&t);
519: TSGetSolution(ts,&U);
521: ctr = event->recorder.ctr-1;
522: if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
523: /* Call the user postevent function */
524: if (event->postevent) {
525: (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);
526: event->recorder.ctr--;
527: }
528: }
530: PetscBarrier((PetscObject)ts);
531: return(0);
532: }