Actual source code: tsevent.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  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: }