xref: /petsc/src/ts/event/tsevent.c (revision f443add64c9a15b09b925d07777e58f4e9aae700)
1af0996ceSBarry Smith #include <petsc/private/tsimpl.h> /*I  "petscts.h" I*/
22dc7a7e3SShri Abhyankar 
36fea3669SShri Abhyankar /*
46427ac75SLisandro Dalcin   TSEventInitialize - Initializes TSEvent for TSSolve
56fea3669SShri Abhyankar */
66427ac75SLisandro Dalcin PetscErrorCode TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U)
76fea3669SShri Abhyankar {
86fea3669SShri Abhyankar   PetscErrorCode ierr;
96fea3669SShri Abhyankar 
106fea3669SShri Abhyankar   PetscFunctionBegin;
116427ac75SLisandro Dalcin   if (!event) PetscFunctionReturn(0);
126427ac75SLisandro Dalcin   PetscValidPointer(event,1);
136427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,2);
146427ac75SLisandro Dalcin   PetscValidHeaderSpecific(U,VEC_CLASSID,4);
156fea3669SShri Abhyankar   event->ptime_prev = t;
1638bf2713SShri Abhyankar   event->iterctr = 0;
176427ac75SLisandro Dalcin   ierr = (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);CHKERRQ(ierr);
186fea3669SShri Abhyankar   PetscFunctionReturn(0);
196fea3669SShri Abhyankar }
206fea3669SShri Abhyankar 
217dbe0728SLisandro Dalcin PetscErrorCode TSEventDestroy(TSEvent *event)
227dbe0728SLisandro Dalcin {
237dbe0728SLisandro Dalcin   PetscErrorCode ierr;
247dbe0728SLisandro Dalcin   PetscInt       i;
257dbe0728SLisandro Dalcin 
267dbe0728SLisandro Dalcin   PetscFunctionBegin;
277dbe0728SLisandro Dalcin   PetscValidPointer(event,1);
287dbe0728SLisandro Dalcin   if (!*event) PetscFunctionReturn(0);
29e7069c78SShri   if (--(*event)->refct > 0) {*event = 0; PetscFunctionReturn(0);}
30e7069c78SShri 
317dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr);
327dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr);
33e2cdd850SShri Abhyankar   ierr = PetscFree((*event)->fvalue_right);CHKERRQ(ierr);
34e2cdd850SShri Abhyankar   ierr = PetscFree((*event)->zerocrossing);CHKERRQ(ierr);
35e2cdd850SShri Abhyankar   ierr = PetscFree((*event)->side);CHKERRQ(ierr);
367dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->direction);CHKERRQ(ierr);
377dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->terminate);CHKERRQ(ierr);
387dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr);
397dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->vtol);CHKERRQ(ierr);
40a4ffd976SShri Abhyankar 
41a4ffd976SShri Abhyankar   for (i=0; i < (*event)->recsize; i++) {
427dbe0728SLisandro Dalcin     ierr = PetscFree((*event)->recorder.eventidx[i]);CHKERRQ(ierr);
437dbe0728SLisandro Dalcin   }
44a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.eventidx);CHKERRQ(ierr);
45a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.nevents);CHKERRQ(ierr);
46a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.stepnum);CHKERRQ(ierr);
47a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.time);CHKERRQ(ierr);
48a4ffd976SShri Abhyankar 
497dbe0728SLisandro Dalcin   ierr = PetscViewerDestroy(&(*event)->monitor);CHKERRQ(ierr);
507dbe0728SLisandro Dalcin   ierr = PetscFree(*event);CHKERRQ(ierr);
517dbe0728SLisandro Dalcin   PetscFunctionReturn(0);
527dbe0728SLisandro Dalcin }
537dbe0728SLisandro Dalcin 
54e3005195SShri Abhyankar /*@
55e3005195SShri Abhyankar    TSSetEventTolerances - Set tolerances for event zero crossings when using event handler
56e3005195SShri Abhyankar 
57e3005195SShri Abhyankar    Logically Collective
58e3005195SShri Abhyankar 
59e3005195SShri Abhyankar    Input Arguments:
60e3005195SShri Abhyankar +  ts - time integration context
61e3005195SShri Abhyankar .  tol - scalar tolerance, PETSC_DECIDE to leave current value
62e3005195SShri Abhyankar -  vtol - array of tolerances or NULL, used in preference to tol if present
63e3005195SShri Abhyankar 
642589fa24SLisandro Dalcin    Options Database Keys:
652589fa24SLisandro Dalcin .  -ts_event_tol <tol> tolerance for event zero crossing
66e3005195SShri Abhyankar 
67e3005195SShri Abhyankar    Notes:
68f25fe674SLisandro Dalcin    Must call TSSetEventHandler() before setting the tolerances.
69e3005195SShri Abhyankar 
70e3005195SShri Abhyankar    The size of vtol is equal to the number of events.
71e3005195SShri Abhyankar 
72e3005195SShri Abhyankar    Level: beginner
73e3005195SShri Abhyankar 
74f25fe674SLisandro Dalcin .seealso: TS, TSEvent, TSSetEventHandler()
75e3005195SShri Abhyankar @*/
766427ac75SLisandro Dalcin PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[])
77e3005195SShri Abhyankar {
786427ac75SLisandro Dalcin   TSEvent        event;
79e3005195SShri Abhyankar   PetscInt       i;
80e3005195SShri Abhyankar 
81e3005195SShri Abhyankar   PetscFunctionBegin;
826427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
836427ac75SLisandro Dalcin   if (vtol) PetscValidRealPointer(vtol,3);
84f25fe674SLisandro Dalcin   if (!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()");
856427ac75SLisandro Dalcin 
866427ac75SLisandro Dalcin   event = ts->event;
87e3005195SShri Abhyankar   if (vtol) {
88e3005195SShri Abhyankar     for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i];
89e3005195SShri Abhyankar   } else {
90e3005195SShri Abhyankar     if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) {
91e3005195SShri Abhyankar       for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
92e3005195SShri Abhyankar     }
93e3005195SShri Abhyankar   }
94e3005195SShri Abhyankar   PetscFunctionReturn(0);
95e3005195SShri Abhyankar }
96e3005195SShri Abhyankar 
972dc7a7e3SShri Abhyankar /*@C
98f25fe674SLisandro Dalcin    TSSetEventHandler - Sets a monitoring function used for detecting events
992dc7a7e3SShri Abhyankar 
1002dc7a7e3SShri Abhyankar    Logically Collective on TS
1012dc7a7e3SShri Abhyankar 
1022dc7a7e3SShri Abhyankar    Input Parameters:
1032dc7a7e3SShri Abhyankar +  ts - the TS context obtained from TSCreate()
1042dc7a7e3SShri Abhyankar .  nevents - number of local events
105d94325d3SShri Abhyankar .  direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction,
106d94325d3SShri Abhyankar                +1 => Zero crossing in positive direction, 0 => both ways (one for each event)
107d94325d3SShri Abhyankar .  terminate - flag to indicate whether time stepping should be terminated after
108d94325d3SShri Abhyankar                event is detected (one for each event)
1096427ac75SLisandro Dalcin .  eventhandler - event monitoring routine
1102dc7a7e3SShri Abhyankar .  postevent - [optional] post-event function
1112589fa24SLisandro Dalcin -  ctx       - [optional] user-defined context for private data for the
1122dc7a7e3SShri Abhyankar                event monitor and post event routine (use NULL if no
1132dc7a7e3SShri Abhyankar                context is desired)
1142dc7a7e3SShri Abhyankar 
1156427ac75SLisandro Dalcin    Calling sequence of eventhandler:
1163a88037aSBarry Smith    PetscErrorCode PetscEventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx)
1172dc7a7e3SShri Abhyankar 
1182dc7a7e3SShri Abhyankar    Input Parameters:
1192dc7a7e3SShri Abhyankar +  ts  - the TS context
1202dc7a7e3SShri Abhyankar .  t   - current time
1212dc7a7e3SShri Abhyankar .  U   - current iterate
1226427ac75SLisandro Dalcin -  ctx - [optional] context passed with eventhandler
1232dc7a7e3SShri Abhyankar 
1242dc7a7e3SShri Abhyankar    Output parameters:
125d94325d3SShri Abhyankar .  fvalue    - function value of events at time t
1262dc7a7e3SShri Abhyankar 
1272dc7a7e3SShri Abhyankar    Calling sequence of postevent:
1282589fa24SLisandro Dalcin    PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
1292dc7a7e3SShri Abhyankar 
1302dc7a7e3SShri Abhyankar    Input Parameters:
1312dc7a7e3SShri Abhyankar +  ts - the TS context
1322dc7a7e3SShri Abhyankar .  nevents_zero - number of local events whose event function is zero
1332dc7a7e3SShri Abhyankar .  events_zero  - indices of local events which have reached zero
1342dc7a7e3SShri Abhyankar .  t            - current time
1352dc7a7e3SShri Abhyankar .  U            - current solution
136031fbad4SShri Abhyankar .  forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0)
1376427ac75SLisandro Dalcin -  ctx          - the context passed with eventhandler
1382dc7a7e3SShri Abhyankar 
1392dc7a7e3SShri Abhyankar    Level: intermediate
1402dc7a7e3SShri Abhyankar 
14102749585SLisandro Dalcin .keywords: TS, event, set
1422dc7a7e3SShri Abhyankar 
1432dc7a7e3SShri Abhyankar .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason()
1442dc7a7e3SShri Abhyankar @*/
1452589fa24SLisandro Dalcin 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)
1462dc7a7e3SShri Abhyankar {
1472dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
1482dc7a7e3SShri Abhyankar   TSEvent        event;
149d94325d3SShri Abhyankar   PetscInt       i;
150006e6a18SShri Abhyankar   PetscBool      flg;
151a6c783d2SShri Abhyankar #if defined PETSC_USE_REAL_SINGLE
152a6c783d2SShri Abhyankar   PetscReal      tol=1e-4;
153a6c783d2SShri Abhyankar #else
154d569cc17SSatish Balay   PetscReal      tol=1e-6;
155a6c783d2SShri Abhyankar #endif
1562dc7a7e3SShri Abhyankar 
1572dc7a7e3SShri Abhyankar   PetscFunctionBegin;
1586427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1590a82b154SShri   if(nevents) {
1606427ac75SLisandro Dalcin     PetscValidIntPointer(direction,2);
1616427ac75SLisandro Dalcin     PetscValidIntPointer(terminate,3);
1620a82b154SShri   }
1636427ac75SLisandro Dalcin 
1646427ac75SLisandro Dalcin   ierr = PetscNewLog(ts,&event);CHKERRQ(ierr);
165854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->fvalue);CHKERRQ(ierr);
166854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->fvalue_prev);CHKERRQ(ierr);
167e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->fvalue_right);CHKERRQ(ierr);
168e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->zerocrossing);CHKERRQ(ierr);
169e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->side);CHKERRQ(ierr);
170854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->direction);CHKERRQ(ierr);
171854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->terminate);CHKERRQ(ierr);
172e3005195SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->vtol);CHKERRQ(ierr);
173d94325d3SShri Abhyankar   for (i=0; i < nevents; i++) {
174d94325d3SShri Abhyankar     event->direction[i] = direction[i];
175d94325d3SShri Abhyankar     event->terminate[i] = terminate[i];
176e2cdd850SShri Abhyankar     event->zerocrossing[i] = PETSC_FALSE;
177e2cdd850SShri Abhyankar     event->side[i] = 0;
178d94325d3SShri Abhyankar   }
179854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->events_zero);CHKERRQ(ierr);
1802589fa24SLisandro Dalcin   event->nevents = nevents;
1816427ac75SLisandro Dalcin   event->eventhandler = eventhandler;
1822dc7a7e3SShri Abhyankar   event->postevent = postevent;
1836427ac75SLisandro Dalcin   event->ctx = ctx;
1842dc7a7e3SShri Abhyankar 
18502749585SLisandro Dalcin   event->recsize = 8;  /* Initial size of the recorder */
186a9514d71SShri Abhyankar   ierr = PetscOptionsBegin(((PetscObject)ts)->comm,"","TS Event options","");CHKERRQ(ierr);
1872dc7a7e3SShri Abhyankar   {
18802749585SLisandro Dalcin     ierr = PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","TSSetEventTolerances",tol,&tol,NULL);CHKERRQ(ierr);
189006e6a18SShri Abhyankar     ierr = PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);CHKERRQ(ierr);
190a4ffd976SShri Abhyankar     ierr = PetscOptionsInt("-ts_event_recorder_initial_size","Initial size of event recorder","",event->recsize,&event->recsize,NULL);CHKERRQ(ierr);
1912dc7a7e3SShri Abhyankar   }
1929e12be75SShri Abhyankar   PetscOptionsEnd();
193a4ffd976SShri Abhyankar 
194a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.time);CHKERRQ(ierr);
195a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.stepnum);CHKERRQ(ierr);
196a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.nevents);CHKERRQ(ierr);
197a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.eventidx);CHKERRQ(ierr);
198a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
199a4ffd976SShri Abhyankar     ierr = PetscMalloc1(event->nevents,&event->recorder.eventidx[i]);CHKERRQ(ierr);
200a4ffd976SShri Abhyankar   }
201e7069c78SShri   /* Initialize the event recorder */
202e7069c78SShri   event->recorder.ctr = 0;
203a4ffd976SShri Abhyankar 
204e3005195SShri Abhyankar   for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
2056427ac75SLisandro Dalcin   if (flg) {ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);CHKERRQ(ierr);}
2069e12be75SShri Abhyankar 
2076427ac75SLisandro Dalcin   ierr = TSEventDestroy(&ts->event);CHKERRQ(ierr);
208d94325d3SShri Abhyankar   ts->event = event;
209e7069c78SShri   ts->event->refct = 1;
2102dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
2112dc7a7e3SShri Abhyankar }
2122dc7a7e3SShri Abhyankar 
213a4ffd976SShri Abhyankar /*
214a4ffd976SShri Abhyankar   TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
215a4ffd976SShri Abhyankar                           is reached.
216a4ffd976SShri Abhyankar */
21702749585SLisandro Dalcin static PetscErrorCode TSEventRecorderResize(TSEvent event)
218a4ffd976SShri Abhyankar {
219a4ffd976SShri Abhyankar   PetscErrorCode ierr;
220a4ffd976SShri Abhyankar   PetscReal      *time;
221a4ffd976SShri Abhyankar   PetscInt       *stepnum;
222a4ffd976SShri Abhyankar   PetscInt       *nevents;
223a4ffd976SShri Abhyankar   PetscInt       **eventidx;
224a4ffd976SShri Abhyankar   PetscInt       i,fact=2;
225a4ffd976SShri Abhyankar 
226a4ffd976SShri Abhyankar   PetscFunctionBegin;
227a4ffd976SShri Abhyankar 
228a4ffd976SShri Abhyankar   /* Create large arrays */
229a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&time);CHKERRQ(ierr);
230a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&stepnum);CHKERRQ(ierr);
231a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&nevents);CHKERRQ(ierr);
232a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&eventidx);CHKERRQ(ierr);
233a4ffd976SShri Abhyankar   for (i=0; i < fact*event->recsize; i++) {
234a4ffd976SShri Abhyankar     ierr = PetscMalloc1(event->nevents,&eventidx[i]);CHKERRQ(ierr);
235a4ffd976SShri Abhyankar   }
236a4ffd976SShri Abhyankar 
237a4ffd976SShri Abhyankar   /* Copy over data */
238a4ffd976SShri Abhyankar   ierr = PetscMemcpy(time,event->recorder.time,event->recsize*sizeof(PetscReal));CHKERRQ(ierr);
239a4ffd976SShri Abhyankar   ierr = PetscMemcpy(stepnum,event->recorder.stepnum,event->recsize*sizeof(PetscInt));CHKERRQ(ierr);
240a4ffd976SShri Abhyankar   ierr = PetscMemcpy(nevents,event->recorder.nevents,event->recsize*sizeof(PetscInt));CHKERRQ(ierr);
241a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
242a4ffd976SShri Abhyankar     ierr = PetscMemcpy(eventidx[i],event->recorder.eventidx[i],event->recorder.nevents[i]*sizeof(PetscInt));CHKERRQ(ierr);
243a4ffd976SShri Abhyankar   }
244a4ffd976SShri Abhyankar 
245a4ffd976SShri Abhyankar   /* Destroy old arrays */
246a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
247a4ffd976SShri Abhyankar     ierr = PetscFree(event->recorder.eventidx[i]);CHKERRQ(ierr);
248a4ffd976SShri Abhyankar   }
249a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.eventidx);CHKERRQ(ierr);
250a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.nevents);CHKERRQ(ierr);
251a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.stepnum);CHKERRQ(ierr);
252a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.time);CHKERRQ(ierr);
253a4ffd976SShri Abhyankar 
254a4ffd976SShri Abhyankar   /* Set pointers */
255a4ffd976SShri Abhyankar   event->recorder.time = time;
256a4ffd976SShri Abhyankar   event->recorder.stepnum = stepnum;
257a4ffd976SShri Abhyankar   event->recorder.nevents = nevents;
258a4ffd976SShri Abhyankar   event->recorder.eventidx = eventidx;
259a4ffd976SShri Abhyankar 
260a4ffd976SShri Abhyankar   /* Double size */
261a4ffd976SShri Abhyankar   event->recsize *= fact;
262a4ffd976SShri Abhyankar 
263a4ffd976SShri Abhyankar   PetscFunctionReturn(0);
264a4ffd976SShri Abhyankar }
265a4ffd976SShri Abhyankar 
266031fbad4SShri Abhyankar /*
2674597913aSLisandro Dalcin    Helper rutine to handle user postenvents and recording
268031fbad4SShri Abhyankar */
2694597913aSLisandro Dalcin static PetscErrorCode TSPostEvent(TS ts,PetscReal t,Vec U)
2702dc7a7e3SShri Abhyankar {
2712dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
2722dc7a7e3SShri Abhyankar   TSEvent        event = ts->event;
2732dc7a7e3SShri Abhyankar   PetscBool      terminate = PETSC_FALSE;
27428d5b5d6SLisandro Dalcin   PetscBool      restart = PETSC_FALSE;
275d0578d90SShri Abhyankar   PetscInt       i,ctr,stepnum;
2767324a0ffSLisandro Dalcin   PetscBool      inflag[2],outflag[2];
2774597913aSLisandro Dalcin   PetscBool      forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */
2782dc7a7e3SShri Abhyankar 
2792dc7a7e3SShri Abhyankar   PetscFunctionBegin;
2802dc7a7e3SShri Abhyankar   if (event->postevent) {
28128d5b5d6SLisandro Dalcin     PetscObjectState state_prev,state_post;
28228d5b5d6SLisandro Dalcin     ierr = PetscObjectStateGet((PetscObject)U,&state_prev);CHKERRQ(ierr);
2834597913aSLisandro Dalcin     ierr = (*event->postevent)(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
28428d5b5d6SLisandro Dalcin     ierr = PetscObjectStateGet((PetscObject)U,&state_post);CHKERRQ(ierr);
28528d5b5d6SLisandro Dalcin     if (state_prev != state_post) restart = PETSC_TRUE;
2862dc7a7e3SShri Abhyankar   }
2874597913aSLisandro Dalcin 
28828d5b5d6SLisandro Dalcin   /* Handle termination events and step restart */
28928d5b5d6SLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
2907324a0ffSLisandro Dalcin   inflag[0] = restart; inflag[1] = terminate;
2917324a0ffSLisandro Dalcin   ierr = MPIU_Allreduce(inflag,outflag,2,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr);
2927324a0ffSLisandro Dalcin   restart = outflag[0]; terminate = outflag[1];
2937324a0ffSLisandro Dalcin   if (restart) {ierr = TSRestartStep(ts);CHKERRQ(ierr);}
2947324a0ffSLisandro Dalcin   if (terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);}
2957324a0ffSLisandro Dalcin   event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;
296f7aea88cSShri Abhyankar 
2974597913aSLisandro Dalcin   /* Reset event residual functions as states might get changed by the postevent callback */
298*f443add6SLisandro Dalcin   if (event->postevent) {
299*f443add6SLisandro Dalcin     ierr = VecLockPush(U);CHKERRQ(ierr);
300*f443add6SLisandro Dalcin     ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);
301*f443add6SLisandro Dalcin     ierr = VecLockPop(U);CHKERRQ(ierr);
302*f443add6SLisandro Dalcin   }
303*f443add6SLisandro Dalcin 
304*f443add6SLisandro Dalcin   /* Cache current time and event residual functions */
305*f443add6SLisandro Dalcin   event->ptime_prev = t;
306*f443add6SLisandro Dalcin   for (i=0; i<event->nevents; i++)
307*f443add6SLisandro Dalcin     event->fvalue_prev[i] = event->fvalue[i];
3084597913aSLisandro Dalcin 
309d0578d90SShri Abhyankar   /* Record the event in the event recorder */
31080275a0aSLisandro Dalcin   ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr);
311f7aea88cSShri Abhyankar   ctr = event->recorder.ctr;
312a4ffd976SShri Abhyankar   if (ctr == event->recsize) {
313a4ffd976SShri Abhyankar     ierr = TSEventRecorderResize(event);CHKERRQ(ierr);
314a4ffd976SShri Abhyankar   }
315f7aea88cSShri Abhyankar   event->recorder.time[ctr] = t;
316d0578d90SShri Abhyankar   event->recorder.stepnum[ctr] = stepnum;
3174597913aSLisandro Dalcin   event->recorder.nevents[ctr] = event->nevents_zero;
3184597913aSLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
319f7aea88cSShri Abhyankar   event->recorder.ctr++;
3202dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
3212dc7a7e3SShri Abhyankar }
3222dc7a7e3SShri Abhyankar 
32302749585SLisandro Dalcin /* Uses Anderson-Bjorck variant of regula falsi method */
324a3a8645aSShri Abhyankar PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt)
32538bf2713SShri Abhyankar {
32602749585SLisandro Dalcin   PetscReal new_dt, scal = 1.0;
327e2cdd850SShri Abhyankar   if (PetscRealPart(fleft)*PetscRealPart(f) < 0) {
328e2cdd850SShri Abhyankar     if (side == 1) {
329a6c783d2SShri Abhyankar       scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright);
330a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
331e2cdd850SShri Abhyankar     }
33202749585SLisandro Dalcin     new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
333e2cdd850SShri Abhyankar   } else {
334e2cdd850SShri Abhyankar     if (side == -1) {
335a6c783d2SShri Abhyankar       scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft);
336a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
337e2cdd850SShri Abhyankar     }
33802749585SLisandro Dalcin     new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t;
339e2cdd850SShri Abhyankar   }
34002749585SLisandro Dalcin   return PetscMin(dt,new_dt);
34138bf2713SShri Abhyankar }
342e2cdd850SShri Abhyankar 
34338bf2713SShri Abhyankar 
3446427ac75SLisandro Dalcin PetscErrorCode TSEventHandler(TS ts)
3452dc7a7e3SShri Abhyankar {
3462dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
3476427ac75SLisandro Dalcin   TSEvent        event;
3482dc7a7e3SShri Abhyankar   PetscReal      t;
3492dc7a7e3SShri Abhyankar   Vec            U;
3502dc7a7e3SShri Abhyankar   PetscInt       i;
3517dbe0728SLisandro Dalcin   PetscReal      dt,dt_min;
3522dc7a7e3SShri Abhyankar   PetscInt       rollback=0,in[2],out[2];
3539e12be75SShri Abhyankar   PetscInt       fvalue_sign,fvalueprev_sign;
3542dc7a7e3SShri Abhyankar 
3552dc7a7e3SShri Abhyankar   PetscFunctionBegin;
3566427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
3576427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
3586427ac75SLisandro Dalcin   event = ts->event;
3592dc7a7e3SShri Abhyankar 
3602dc7a7e3SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
3617dbe0728SLisandro Dalcin   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
3622dc7a7e3SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
3632dc7a7e3SShri Abhyankar 
3647dbe0728SLisandro Dalcin   if (event->status == TSEVENT_NONE) {
3652d5bd56dSLisandro Dalcin     if (ts->steps == 1) /* After first successful step */
3662d5bd56dSLisandro Dalcin       event->timestep_orig = ts->ptime - ts->ptime_prev;
3677dbe0728SLisandro Dalcin     event->timestep_prev = dt;
3687dbe0728SLisandro Dalcin   }
3697dbe0728SLisandro Dalcin 
3702dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_RESET_NEXTSTEP) {
3712d5bd56dSLisandro Dalcin     /* Restore time step */
3722d5bd56dSLisandro Dalcin     dt = PetscMin(event->timestep_orig,event->timestep_prev);
3737dbe0728SLisandro Dalcin     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
3742dc7a7e3SShri Abhyankar     event->status = TSEVENT_NONE;
3752dc7a7e3SShri Abhyankar   }
3762dc7a7e3SShri Abhyankar 
3772dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_NONE) {
3787dbe0728SLisandro Dalcin     event->ptime_end = t;
3792dc7a7e3SShri Abhyankar   }
3802dc7a7e3SShri Abhyankar 
381*f443add6SLisandro Dalcin   ierr = VecLockPush(U);CHKERRQ(ierr);
3826427ac75SLisandro Dalcin   ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);
383*f443add6SLisandro Dalcin   ierr = VecLockPop(U);CHKERRQ(ierr);
3849e12be75SShri Abhyankar 
3852dc7a7e3SShri Abhyankar   for (i=0; i < event->nevents; i++) {
386e3005195SShri Abhyankar     if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
3872dc7a7e3SShri Abhyankar       event->status = TSEVENT_ZERO;
388e2cdd850SShri Abhyankar       event->fvalue_right[i] = event->fvalue[i];
3899e12be75SShri Abhyankar       continue;
390006e6a18SShri Abhyankar     }
391d94325d3SShri Abhyankar     fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
392d94325d3SShri Abhyankar     fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
393e3005195SShri Abhyankar     if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) {
394d94325d3SShri Abhyankar       switch (event->direction[i]) {
395d94325d3SShri Abhyankar       case -1:
396d94325d3SShri Abhyankar         if (fvalue_sign < 0) {
3972dc7a7e3SShri Abhyankar           rollback = 1;
398e2cdd850SShri Abhyankar 
399e2cdd850SShri Abhyankar           /* Compute new time step */
400e2cdd850SShri Abhyankar           dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
401e2cdd850SShri Abhyankar 
4026427ac75SLisandro Dalcin           if (event->monitor) {
40338bf2713SShri Abhyankar             ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr);
404006e6a18SShri Abhyankar           }
405e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
406e2cdd850SShri Abhyankar           event->side[i] = 1;
407e2cdd850SShri Abhyankar 
408e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
409e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4102dc7a7e3SShri Abhyankar         }
411d94325d3SShri Abhyankar         break;
412d94325d3SShri Abhyankar       case 1:
413d94325d3SShri Abhyankar         if (fvalue_sign > 0) {
414d94325d3SShri Abhyankar           rollback = 1;
415e2cdd850SShri Abhyankar 
416e2cdd850SShri Abhyankar           /* Compute new time step */
417e2cdd850SShri Abhyankar           dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
418e2cdd850SShri Abhyankar 
4196427ac75SLisandro Dalcin           if (event->monitor) {
42038bf2713SShri Abhyankar             ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr);
421006e6a18SShri Abhyankar           }
422e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
423e2cdd850SShri Abhyankar           event->side[i] = 1;
424e2cdd850SShri Abhyankar 
425e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
426e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4272dc7a7e3SShri Abhyankar         }
428d94325d3SShri Abhyankar         break;
429d94325d3SShri Abhyankar       case 0:
430d94325d3SShri Abhyankar         rollback = 1;
431e2cdd850SShri Abhyankar 
432e2cdd850SShri Abhyankar         /* Compute new time step */
433e2cdd850SShri Abhyankar         dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
434e2cdd850SShri Abhyankar 
4356427ac75SLisandro Dalcin         if (event->monitor) {
43638bf2713SShri Abhyankar           ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr);
437006e6a18SShri Abhyankar         }
438e2cdd850SShri Abhyankar         event->fvalue_right[i] = event->fvalue[i];
439e2cdd850SShri Abhyankar         event->side[i] = 1;
440e2cdd850SShri Abhyankar 
441e2cdd850SShri Abhyankar         if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
442e2cdd850SShri Abhyankar         event->status = TSEVENT_LOCATED_INTERVAL;
443e2cdd850SShri Abhyankar 
444d94325d3SShri Abhyankar         break;
445d94325d3SShri Abhyankar       }
446d94325d3SShri Abhyankar     }
447d94325d3SShri Abhyankar   }
4487dbe0728SLisandro Dalcin 
4492589fa24SLisandro Dalcin   in[0] = event->status; in[1] = rollback;
4507dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
4512589fa24SLisandro Dalcin   event->status = (TSEventStatus)out[0]; rollback = out[1];
4522589fa24SLisandro Dalcin   if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
4532dc7a7e3SShri Abhyankar 
4547dbe0728SLisandro Dalcin   event->nevents_zero = 0;
4559e12be75SShri Abhyankar   if (event->status == TSEVENT_ZERO) {
4569e12be75SShri Abhyankar     for (i=0; i < event->nevents; i++) {
4579e12be75SShri Abhyankar       if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
4589e12be75SShri Abhyankar         event->events_zero[event->nevents_zero++] = i;
4596427ac75SLisandro Dalcin         if (event->monitor) {
46038bf2713SShri Abhyankar           ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: Event %D zero crossing at time %g located in %D iterations\n",i,(double)t,event->iterctr);CHKERRQ(ierr);
4619e12be75SShri Abhyankar         }
462e2cdd850SShri Abhyankar         event->zerocrossing[i] = PETSC_FALSE;
4639e12be75SShri Abhyankar       }
464e2cdd850SShri Abhyankar       event->side[i] = 0;
4659e12be75SShri Abhyankar     }
4664597913aSLisandro Dalcin     ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr);
4679e12be75SShri Abhyankar 
468e2cdd850SShri Abhyankar     dt = event->ptime_end - t;
4692d5bd56dSLisandro Dalcin     if (PetscAbsReal(dt) < PETSC_SMALL) dt += PetscMin(event->timestep_orig,event->timestep_prev); /* XXX Should be done better */
4709e12be75SShri Abhyankar     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
47138bf2713SShri Abhyankar     event->iterctr = 0;
4729e12be75SShri Abhyankar     PetscFunctionReturn(0);
4739e12be75SShri Abhyankar   }
4749e12be75SShri Abhyankar 
4752dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_LOCATED_INTERVAL) {
4762dc7a7e3SShri Abhyankar     ierr = TSRollBack(ts);CHKERRQ(ierr);
4772589fa24SLisandro Dalcin     ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr);
4782dc7a7e3SShri Abhyankar     event->status = TSEVENT_PROCESSING;
479e2cdd850SShri Abhyankar     event->ptime_right = t;
4802dc7a7e3SShri Abhyankar   } else {
4812dc7a7e3SShri Abhyankar     for (i=0; i < event->nevents; i++) {
482e2cdd850SShri Abhyankar       if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) {
483e2cdd850SShri Abhyankar         /* Compute new time step */
484e2cdd850SShri Abhyankar         dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt);
485e2cdd850SShri Abhyankar         event->side[i] = -1;
486e2cdd850SShri Abhyankar       }
4872dc7a7e3SShri Abhyankar       event->fvalue_prev[i] = event->fvalue[i];
4882dc7a7e3SShri Abhyankar     }
489e2cdd850SShri Abhyankar     if (event->monitor && event->status == TSEVENT_PROCESSING) {
49038bf2713SShri Abhyankar       ierr = 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);CHKERRQ(ierr);
4912dc7a7e3SShri Abhyankar     }
49238bf2713SShri Abhyankar     event->ptime_prev = t;
49338bf2713SShri Abhyankar   }
49438bf2713SShri Abhyankar 
49538bf2713SShri Abhyankar   if (event->status == TSEVENT_PROCESSING) event->iterctr++;
4969e12be75SShri Abhyankar 
4977dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
4987dbe0728SLisandro Dalcin   ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr);
4992dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
5002dc7a7e3SShri Abhyankar }
5012dc7a7e3SShri Abhyankar 
5026427ac75SLisandro Dalcin PetscErrorCode TSAdjointEventHandler(TS ts)
503d0578d90SShri Abhyankar {
504d0578d90SShri Abhyankar   PetscErrorCode ierr;
5056427ac75SLisandro Dalcin   TSEvent        event;
506d0578d90SShri Abhyankar   PetscReal      t;
507d0578d90SShri Abhyankar   Vec            U;
508d0578d90SShri Abhyankar   PetscInt       ctr;
509d0578d90SShri Abhyankar   PetscBool      forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
510d0578d90SShri Abhyankar 
511d0578d90SShri Abhyankar   PetscFunctionBegin;
5126427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
5136427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
5146427ac75SLisandro Dalcin   event = ts->event;
515d0578d90SShri Abhyankar 
516d0578d90SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
517d0578d90SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
518d0578d90SShri Abhyankar 
519d0578d90SShri Abhyankar   ctr = event->recorder.ctr-1;
520bcbf8bb3SShri Abhyankar   if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
521d0578d90SShri Abhyankar     /* Call the user postevent function */
522d0578d90SShri Abhyankar     if (event->postevent) {
5236427ac75SLisandro Dalcin       ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
524d0578d90SShri Abhyankar       event->recorder.ctr--;
525d0578d90SShri Abhyankar     }
526d0578d90SShri Abhyankar   }
527d0578d90SShri Abhyankar 
528d0578d90SShri Abhyankar   PetscBarrier((PetscObject)ts);
529d0578d90SShri Abhyankar   PetscFunctionReturn(0);
530d0578d90SShri Abhyankar }
531