xref: /petsc/src/ts/event/tsevent.c (revision 7324a0ffd8f54dc9690180a83c9ccd78651131f8)
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;
276*7324a0ffSLisandro 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;
290*7324a0ffSLisandro Dalcin   inflag[0] = restart; inflag[1] = terminate;
291*7324a0ffSLisandro Dalcin   ierr = MPIU_Allreduce(inflag,outflag,2,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr);
292*7324a0ffSLisandro Dalcin   restart = outflag[0]; terminate = outflag[1];
293*7324a0ffSLisandro Dalcin   if (restart) {ierr = TSRestartStep(ts);CHKERRQ(ierr);}
294*7324a0ffSLisandro Dalcin   if (terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);}
295*7324a0ffSLisandro Dalcin   event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;
296f7aea88cSShri Abhyankar 
2974597913aSLisandro Dalcin   event->ptime_prev = t;
2984597913aSLisandro Dalcin   /* Reset event residual functions as states might get changed by the postevent callback */
2994597913aSLisandro Dalcin   if (event->postevent) {ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);}
3004597913aSLisandro Dalcin   /* Cache current event residual functions */
3014597913aSLisandro Dalcin   for (i=0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];
3024597913aSLisandro Dalcin 
303d0578d90SShri Abhyankar   /* Record the event in the event recorder */
30480275a0aSLisandro Dalcin   ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr);
305f7aea88cSShri Abhyankar   ctr = event->recorder.ctr;
306a4ffd976SShri Abhyankar   if (ctr == event->recsize) {
307a4ffd976SShri Abhyankar     ierr = TSEventRecorderResize(event);CHKERRQ(ierr);
308a4ffd976SShri Abhyankar   }
309f7aea88cSShri Abhyankar   event->recorder.time[ctr] = t;
310d0578d90SShri Abhyankar   event->recorder.stepnum[ctr] = stepnum;
3114597913aSLisandro Dalcin   event->recorder.nevents[ctr] = event->nevents_zero;
3124597913aSLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
313f7aea88cSShri Abhyankar   event->recorder.ctr++;
3142dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
3152dc7a7e3SShri Abhyankar }
3162dc7a7e3SShri Abhyankar 
31702749585SLisandro Dalcin /* Uses Anderson-Bjorck variant of regula falsi method */
318a3a8645aSShri Abhyankar PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt)
31938bf2713SShri Abhyankar {
32002749585SLisandro Dalcin   PetscReal new_dt, scal = 1.0;
321e2cdd850SShri Abhyankar   if (PetscRealPart(fleft)*PetscRealPart(f) < 0) {
322e2cdd850SShri Abhyankar     if (side == 1) {
323a6c783d2SShri Abhyankar       scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright);
324a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
325e2cdd850SShri Abhyankar     }
32602749585SLisandro Dalcin     new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
327e2cdd850SShri Abhyankar   } else {
328e2cdd850SShri Abhyankar     if (side == -1) {
329a6c783d2SShri Abhyankar       scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft);
330a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
331e2cdd850SShri Abhyankar     }
33202749585SLisandro Dalcin     new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t;
333e2cdd850SShri Abhyankar   }
33402749585SLisandro Dalcin   return PetscMin(dt,new_dt);
33538bf2713SShri Abhyankar }
336e2cdd850SShri Abhyankar 
33738bf2713SShri Abhyankar 
3386427ac75SLisandro Dalcin PetscErrorCode TSEventHandler(TS ts)
3392dc7a7e3SShri Abhyankar {
3402dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
3416427ac75SLisandro Dalcin   TSEvent        event;
3422dc7a7e3SShri Abhyankar   PetscReal      t;
3432dc7a7e3SShri Abhyankar   Vec            U;
3442dc7a7e3SShri Abhyankar   PetscInt       i;
3457dbe0728SLisandro Dalcin   PetscReal      dt,dt_min;
3462dc7a7e3SShri Abhyankar   PetscInt       rollback=0,in[2],out[2];
3479e12be75SShri Abhyankar   PetscInt       fvalue_sign,fvalueprev_sign;
3482dc7a7e3SShri Abhyankar 
3492dc7a7e3SShri Abhyankar   PetscFunctionBegin;
3506427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
3516427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
3526427ac75SLisandro Dalcin   event = ts->event;
3532dc7a7e3SShri Abhyankar 
3542dc7a7e3SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
3557dbe0728SLisandro Dalcin   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
3562dc7a7e3SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
3572dc7a7e3SShri Abhyankar 
3587dbe0728SLisandro Dalcin   if (event->status == TSEVENT_NONE) {
3592d5bd56dSLisandro Dalcin     if (ts->steps == 1) /* After first successful step */
3602d5bd56dSLisandro Dalcin       event->timestep_orig = ts->ptime - ts->ptime_prev;
3617dbe0728SLisandro Dalcin     event->timestep_prev = dt;
3627dbe0728SLisandro Dalcin   }
3637dbe0728SLisandro Dalcin 
3642dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_RESET_NEXTSTEP) {
3652d5bd56dSLisandro Dalcin     /* Restore time step */
3662d5bd56dSLisandro Dalcin     dt = PetscMin(event->timestep_orig,event->timestep_prev);
3677dbe0728SLisandro Dalcin     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
3682dc7a7e3SShri Abhyankar     event->status = TSEVENT_NONE;
3692dc7a7e3SShri Abhyankar   }
3702dc7a7e3SShri Abhyankar 
3712dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_NONE) {
3727dbe0728SLisandro Dalcin     event->ptime_end = t;
3732dc7a7e3SShri Abhyankar   }
3742dc7a7e3SShri Abhyankar 
3756427ac75SLisandro Dalcin   ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);
3769e12be75SShri Abhyankar 
3772dc7a7e3SShri Abhyankar   for (i=0; i < event->nevents; i++) {
378e3005195SShri Abhyankar     if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
3792dc7a7e3SShri Abhyankar       event->status = TSEVENT_ZERO;
380e2cdd850SShri Abhyankar       event->fvalue_right[i] = event->fvalue[i];
3819e12be75SShri Abhyankar       continue;
382006e6a18SShri Abhyankar     }
383d94325d3SShri Abhyankar     fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
384d94325d3SShri Abhyankar     fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
385e3005195SShri Abhyankar     if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) {
386d94325d3SShri Abhyankar       switch (event->direction[i]) {
387d94325d3SShri Abhyankar       case -1:
388d94325d3SShri Abhyankar         if (fvalue_sign < 0) {
3892dc7a7e3SShri Abhyankar           rollback = 1;
390e2cdd850SShri Abhyankar 
391e2cdd850SShri Abhyankar           /* Compute new time step */
392e2cdd850SShri 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);
393e2cdd850SShri Abhyankar 
3946427ac75SLisandro Dalcin           if (event->monitor) {
39538bf2713SShri 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);
396006e6a18SShri Abhyankar           }
397e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
398e2cdd850SShri Abhyankar           event->side[i] = 1;
399e2cdd850SShri Abhyankar 
400e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
401e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4022dc7a7e3SShri Abhyankar         }
403d94325d3SShri Abhyankar         break;
404d94325d3SShri Abhyankar       case 1:
405d94325d3SShri Abhyankar         if (fvalue_sign > 0) {
406d94325d3SShri Abhyankar           rollback = 1;
407e2cdd850SShri Abhyankar 
408e2cdd850SShri Abhyankar           /* Compute new time step */
409e2cdd850SShri 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);
410e2cdd850SShri Abhyankar 
4116427ac75SLisandro Dalcin           if (event->monitor) {
41238bf2713SShri 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);
413006e6a18SShri Abhyankar           }
414e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
415e2cdd850SShri Abhyankar           event->side[i] = 1;
416e2cdd850SShri Abhyankar 
417e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
418e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4192dc7a7e3SShri Abhyankar         }
420d94325d3SShri Abhyankar         break;
421d94325d3SShri Abhyankar       case 0:
422d94325d3SShri Abhyankar         rollback = 1;
423e2cdd850SShri Abhyankar 
424e2cdd850SShri Abhyankar         /* Compute new time step */
425e2cdd850SShri 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);
426e2cdd850SShri Abhyankar 
4276427ac75SLisandro Dalcin         if (event->monitor) {
42838bf2713SShri 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);
429006e6a18SShri Abhyankar         }
430e2cdd850SShri Abhyankar         event->fvalue_right[i] = event->fvalue[i];
431e2cdd850SShri Abhyankar         event->side[i] = 1;
432e2cdd850SShri Abhyankar 
433e2cdd850SShri Abhyankar         if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
434e2cdd850SShri Abhyankar         event->status = TSEVENT_LOCATED_INTERVAL;
435e2cdd850SShri Abhyankar 
436d94325d3SShri Abhyankar         break;
437d94325d3SShri Abhyankar       }
438d94325d3SShri Abhyankar     }
439d94325d3SShri Abhyankar   }
4407dbe0728SLisandro Dalcin 
4412589fa24SLisandro Dalcin   in[0] = event->status; in[1] = rollback;
4427dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
4432589fa24SLisandro Dalcin   event->status = (TSEventStatus)out[0]; rollback = out[1];
4442589fa24SLisandro Dalcin   if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
4452dc7a7e3SShri Abhyankar 
4467dbe0728SLisandro Dalcin   event->nevents_zero = 0;
4479e12be75SShri Abhyankar   if (event->status == TSEVENT_ZERO) {
4489e12be75SShri Abhyankar     for (i=0; i < event->nevents; i++) {
4499e12be75SShri Abhyankar       if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
4509e12be75SShri Abhyankar         event->events_zero[event->nevents_zero++] = i;
4516427ac75SLisandro Dalcin         if (event->monitor) {
45238bf2713SShri 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);
4539e12be75SShri Abhyankar         }
454e2cdd850SShri Abhyankar         event->zerocrossing[i] = PETSC_FALSE;
4559e12be75SShri Abhyankar       }
456e2cdd850SShri Abhyankar       event->side[i] = 0;
4579e12be75SShri Abhyankar     }
4584597913aSLisandro Dalcin     ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr);
4599e12be75SShri Abhyankar 
460e2cdd850SShri Abhyankar     dt = event->ptime_end - t;
4612d5bd56dSLisandro Dalcin     if (PetscAbsReal(dt) < PETSC_SMALL) dt += PetscMin(event->timestep_orig,event->timestep_prev); /* XXX Should be done better */
4629e12be75SShri Abhyankar     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
46338bf2713SShri Abhyankar     event->iterctr = 0;
4649e12be75SShri Abhyankar     PetscFunctionReturn(0);
4659e12be75SShri Abhyankar   }
4669e12be75SShri Abhyankar 
4672dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_LOCATED_INTERVAL) {
4682dc7a7e3SShri Abhyankar     ierr = TSRollBack(ts);CHKERRQ(ierr);
4692589fa24SLisandro Dalcin     ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr);
4702dc7a7e3SShri Abhyankar     event->status = TSEVENT_PROCESSING;
471e2cdd850SShri Abhyankar     event->ptime_right = t;
4722dc7a7e3SShri Abhyankar   } else {
4732dc7a7e3SShri Abhyankar     for (i=0; i < event->nevents; i++) {
474e2cdd850SShri Abhyankar       if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) {
475e2cdd850SShri Abhyankar         /* Compute new time step */
476e2cdd850SShri 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);
477e2cdd850SShri Abhyankar         event->side[i] = -1;
478e2cdd850SShri Abhyankar       }
4792dc7a7e3SShri Abhyankar       event->fvalue_prev[i] = event->fvalue[i];
4802dc7a7e3SShri Abhyankar     }
481e2cdd850SShri Abhyankar     if (event->monitor && event->status == TSEVENT_PROCESSING) {
48238bf2713SShri 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);
4832dc7a7e3SShri Abhyankar     }
48438bf2713SShri Abhyankar     event->ptime_prev = t;
48538bf2713SShri Abhyankar   }
48638bf2713SShri Abhyankar 
48738bf2713SShri Abhyankar   if (event->status == TSEVENT_PROCESSING) event->iterctr++;
4889e12be75SShri Abhyankar 
4897dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
4907dbe0728SLisandro Dalcin   ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr);
4912dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
4922dc7a7e3SShri Abhyankar }
4932dc7a7e3SShri Abhyankar 
4946427ac75SLisandro Dalcin PetscErrorCode TSAdjointEventHandler(TS ts)
495d0578d90SShri Abhyankar {
496d0578d90SShri Abhyankar   PetscErrorCode ierr;
4976427ac75SLisandro Dalcin   TSEvent        event;
498d0578d90SShri Abhyankar   PetscReal      t;
499d0578d90SShri Abhyankar   Vec            U;
500d0578d90SShri Abhyankar   PetscInt       ctr;
501d0578d90SShri Abhyankar   PetscBool      forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
502d0578d90SShri Abhyankar 
503d0578d90SShri Abhyankar   PetscFunctionBegin;
5046427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
5056427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
5066427ac75SLisandro Dalcin   event = ts->event;
507d0578d90SShri Abhyankar 
508d0578d90SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
509d0578d90SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
510d0578d90SShri Abhyankar 
511d0578d90SShri Abhyankar   ctr = event->recorder.ctr-1;
512bcbf8bb3SShri Abhyankar   if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
513d0578d90SShri Abhyankar     /* Call the user postevent function */
514d0578d90SShri Abhyankar     if (event->postevent) {
5156427ac75SLisandro Dalcin       ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
516d0578d90SShri Abhyankar       event->recorder.ctr--;
517d0578d90SShri Abhyankar     }
518d0578d90SShri Abhyankar   }
519d0578d90SShri Abhyankar 
520d0578d90SShri Abhyankar   PetscBarrier((PetscObject)ts);
521d0578d90SShri Abhyankar   PetscFunctionReturn(0);
522d0578d90SShri Abhyankar }
523