Actual source code: tsevent.c
petsc-3.6.1 2015-08-06
2: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
6: /*
7: TSEventMonitorInitialize - Initializes TSEvent for TSSolve
8: */
9: PetscErrorCode TSEventMonitorInitialize(TS ts)
10: {
12: PetscReal t;
13: Vec U;
14: TSEvent event=ts->event;
18: TSGetTime(ts,&t);
19: TSGetTimeStep(ts,&event->initial_timestep);
20: TSGetSolution(ts,&U);
21: event->ptime_prev = t;
22: (*event->monitor)(ts,t,U,event->fvalue_prev,event->monitorcontext);
24: /* Initialize the event recorder */
25: event->recorder.ctr = 0;
27: return(0);
28: }
32: /*@
33: TSSetEventTolerances - Set tolerances for event zero crossings when using event handler
35: Logically Collective
37: Input Arguments:
38: + ts - time integration context
39: . tol - scalar tolerance, PETSC_DECIDE to leave current value
40: - vtol - array of tolerances or NULL, used in preference to tol if present
42: - -ts_event_tol <tol> tolerance for event zero crossing
44: Notes:
45: Must call TSSetEventMonitor() before setting the tolerances.
47: The size of vtol is equal to the number of events.
49: Level: beginner
51: .seealso: TS, TSEvent, TSSetEventMonitor()
52: @*/
53: PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal *vtol)
54: {
55: TSEvent event=ts->event;
56: PetscInt i;
59: if(!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventMonitor()");
60: if(vtol) {
61: for(i=0; i < event->nevents; i++) event->vtol[i] = vtol[i];
62: } else {
63: if(tol != PETSC_DECIDE || tol != PETSC_DEFAULT) {
64: for(i=0; i < event->nevents; i++) event->vtol[i] = tol;
65: }
66: }
67: return(0);
68: }
72: /*@C
73: TSSetEventMonitor - Sets a monitoring function used for detecting events
75: Logically Collective on TS
77: Input Parameters:
78: + ts - the TS context obtained from TSCreate()
79: . nevents - number of local events
80: . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction,
81: +1 => Zero crossing in positive direction, 0 => both ways (one for each event)
82: . terminate - flag to indicate whether time stepping should be terminated after
83: event is detected (one for each event)
84: . eventmonitor - event monitoring routine
85: . postevent - [optional] post-event function
86: - mectx - [optional] user-defined context for private data for the
87: event monitor and post event routine (use NULL if no
88: context is desired)
90: Calling sequence of eventmonitor:
91: PetscErrorCode EventMonitor(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void* mectx)
93: Input Parameters:
94: + ts - the TS context
95: . t - current time
96: . U - current iterate
97: - ctx - [optional] context passed with eventmonitor
99: Output parameters:
100: . fvalue - function value of events at time t
101:
102: Calling sequence of postevent:
103: PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero, PetscInt events_zero, PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
105: Input Parameters:
106: + ts - the TS context
107: . nevents_zero - number of local events whose event function is zero
108: . events_zero - indices of local events which have reached zero
109: . t - current time
110: . U - current solution
111: . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0)
112: - ctx - the context passed with eventmonitor
114: Level: intermediate
116: .keywords: TS, event, set, monitor
118: .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason()
119: @*/
120: PetscErrorCode TSSetEventMonitor(TS ts,PetscInt nevents,PetscInt *direction,PetscBool *terminate,PetscErrorCode (*eventmonitor)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,PetscBool,void*),void *mectx)
121: {
123: TSEvent event;
124: PetscInt i;
125: PetscBool flg;
126: PetscReal tol=1e-6;
129: PetscNew(&event);
130: PetscMalloc1(nevents,&event->fvalue);
131: PetscMalloc1(nevents,&event->fvalue_prev);
132: PetscMalloc1(nevents,&event->direction);
133: PetscMalloc1(nevents,&event->terminate);
134: PetscMalloc1(nevents,&event->vtol);
135: for (i=0; i < nevents; i++) {
136: event->direction[i] = direction[i];
137: event->terminate[i] = terminate[i];
138: }
139: PetscMalloc1(nevents,&event->events_zero);
140: event->monitor = eventmonitor;
141: event->postevent = postevent;
142: event->monitorcontext = (void*)mectx;
143: event->nevents = nevents;
145: for(i=0; i < MAXEVENTRECORDERS; i++) {
146: PetscMalloc1(nevents,&event->recorder.eventidx[i]);
147: }
149: PetscOptionsBegin(((PetscObject)ts)->comm,"","TS Event options","");
150: {
151: PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","",tol,&tol,NULL);
152: PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);
153: }
154: PetscOptionsEnd();
156: for(i=0; i < event->nevents; i++) event->vtol[i] = tol;
159: if(flg) {
160: PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->mon);
161: }
162: ts->event = event;
163: return(0);
164: }
168: /*
169: TSPostEvent - Does post event processing by calling the user-defined postevent function
171: Logically Collective on TS
173: Input Parameters:
174: + ts - the TS context
175: . nevents_zero - number of local events whose event function is zero
176: . events_zero - indices of local events which have reached zero
177: . t - current time
178: . U - current solution
179: . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0)
180: - ctx - the context passed with eventmonitor
182: Level: intermediate
184: .keywords: TS, event, set, monitor
186: .seealso: TSSetEventMonitor(),TSEvent
187: */
190: PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void *ctx)
191: {
193: TSEvent event=ts->event;
194: PetscBool terminate=PETSC_FALSE;
195: PetscInt i,ctr,stepnum;
196: PetscBool ts_terminate;
199: if (event->postevent) {
200: (*event->postevent)(ts,nevents_zero,events_zero,t,U,forwardsolve,ctx);
201: }
202: for(i = 0; i < nevents_zero;i++) {
203: terminate = (PetscBool)(terminate || event->terminate[events_zero[i]]);
204: }
205: MPI_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);
206: if (ts_terminate) {
207: TSSetConvergedReason(ts,TS_CONVERGED_EVENT);
208: event->status = TSEVENT_NONE;
209: } else {
210: event->status = TSEVENT_RESET_NEXTSTEP;
211: }
213: /* Record the event in the event recorder */
214: TSGetTimeStepNumber(ts,&stepnum);
215: ctr = event->recorder.ctr;
216: if (ctr == MAXEVENTRECORDERS) {
217: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Exceeded limit (=%d) for number of events recorded",MAXEVENTRECORDERS);
218: }
219: event->recorder.time[ctr] = t;
220: event->recorder.stepnum[ctr] = stepnum;
221: event->recorder.nevents[ctr] = nevents_zero;
222: for(i=0; i < nevents_zero; i++) event->recorder.eventidx[ctr][i] = events_zero[i];
223: event->recorder.ctr++;
225: /* Reset the event residual functions as states might get changed by the postevent callback */
226: (*event->monitor)(ts,t,U,event->fvalue_prev,event->monitorcontext);
227: event->ptime_prev = t;
229: return(0);
230: }
234: PetscErrorCode TSEventMonitorDestroy(TSEvent *event)
235: {
237: PetscInt i;
240: PetscFree((*event)->fvalue);
241: PetscFree((*event)->fvalue_prev);
242: PetscFree((*event)->direction);
243: PetscFree((*event)->terminate);
244: PetscFree((*event)->events_zero);
245: PetscFree((*event)->vtol);
246: for(i=0; i < MAXEVENTRECORDERS; i++) {
247: PetscFree((*event)->recorder.eventidx[i]);
248: }
249: PetscViewerDestroy(&(*event)->mon);
250: PetscFree(*event);
251: *event = NULL;
252: return(0);
253: }
257: PetscErrorCode TSEventMonitor(TS ts)
258: {
260: TSEvent event=ts->event;
261: PetscReal t;
262: Vec U;
263: PetscInt i;
264: PetscReal dt;
265: PetscInt rollback=0,in[2],out[2];
266: PetscBool forwardsolve=PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */
267: PetscInt fvalue_sign,fvalueprev_sign;
271: TSGetTime(ts,&t);
272: TSGetSolution(ts,&U);
274: TSGetTimeStep(ts,&dt);
275: if (event->status == TSEVENT_RESET_NEXTSTEP) {
276: /* Take initial time step */
277: dt = event->initial_timestep;
278: ts->time_step = dt;
279: event->status = TSEVENT_NONE;
280: }
282: if (event->status == TSEVENT_NONE) {
283: event->tstepend = t;
284: }
286: event->nevents_zero = 0;
288: (*event->monitor)(ts,t,U,event->fvalue,event->monitorcontext);
290: for (i = 0; i < event->nevents; i++) {
291: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
292: event->status = TSEVENT_ZERO;
293: continue;
294: }
295: fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
296: fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
297: if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) {
298: switch (event->direction[i]) {
299: case -1:
300: if (fvalue_sign < 0) {
301: rollback = 1;
302: /* Compute linearly interpolated new time step */
303: dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i])));
304: if(event->mon) {
305: PetscViewerASCIIPrintf(event->mon,"TSEvent : Event %D interval located [%g - %g]\n",i,(double)event->ptime_prev,(double)t);
306: }
307: }
308: break;
309: case 1:
310: if (fvalue_sign > 0) {
311: rollback = 1;
312: /* Compute linearly interpolated new time step */
313: dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i])));
314: if(event->mon) {
315: PetscViewerASCIIPrintf(event->mon,"TSEvent : Event %D interval located [%g - %g]\n",i,(double)event->ptime_prev,(double)t);
316: }
317: }
318: break;
319: case 0:
320: rollback = 1;
321: /* Compute linearly interpolated new time step */
322: dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i])));
323: if(event->mon) {
324: PetscViewerASCIIPrintf(event->mon,"TSEvent : Event %D interval located [%g - %g]\n",i,(double)event->ptime_prev,(double)t);
325: }
326: break;
327: }
328: }
329: }
330: if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
331:
332: in[0] = event->status;
333: in[1] = rollback;
334: MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);
335:
336: event->status = (TSEventStatus)out[0];
337: rollback = out[1];
338: if (rollback) {
339: event->status = TSEVENT_LOCATED_INTERVAL;
340: }
342: if (event->status == TSEVENT_ZERO) {
343: for(i=0; i < event->nevents; i++) {
344: if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
345: event->events_zero[event->nevents_zero++] = i;
346: if(event->mon) {
347: PetscViewerASCIIPrintf(event->mon,"TSEvent : Event %D zero crossing at time %g\n",i,(double)t);
348: }
349: }
350: }
352: TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->monitorcontext);
354: dt = event->tstepend-t;
355: if(PetscAbsReal(dt) < PETSC_SMALL) dt += event->initial_timestep;
356: TSSetTimeStep(ts,dt);
357: return(0);
358: }
360: if (event->status == TSEVENT_LOCATED_INTERVAL) {
361: TSRollBack(ts);
362: ts->steps--;
363: ts->total_steps--;
364: ts->reason = TS_CONVERGED_ITERATING;
365: event->status = TSEVENT_PROCESSING;
366: } else {
367: for (i = 0; i < event->nevents; i++) {
368: event->fvalue_prev[i] = event->fvalue[i];
369: }
370: event->ptime_prev = t;
371: if (event->status == TSEVENT_PROCESSING) {
372: dt = event->tstepend - event->ptime_prev;
373: }
374: }
376: MPI_Allreduce(&dt,&(ts->time_step),1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);
377: return(0);
378: }
382: PetscErrorCode TSAdjointEventMonitor(TS ts)
383: {
385: TSEvent event=ts->event;
386: PetscReal t;
387: Vec U;
388: PetscInt ctr;
389: PetscBool forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
393: TSGetTime(ts,&t);
394: TSGetSolution(ts,&U);
396: ctr = event->recorder.ctr-1;
397: if(ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
398: /* Call the user postevent function */
399: if(event->postevent) {
400: (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->monitorcontext);
401: event->recorder.ctr--;
402: }
403: }
405: PetscBarrier((PetscObject)ts);
406: return(0);
407: }