xref: /petsc/src/ts/event/tsevent.c (revision e7069c78e74748c563594f592e4924a0e908de35)
1af0996ceSBarry Smith #include <petsc/private/tsimpl.h> /*I  "petscts.h" I*/
22dc7a7e3SShri Abhyankar 
32dc7a7e3SShri Abhyankar #undef __FUNCT__
46427ac75SLisandro Dalcin #define __FUNCT__ "TSEventInitialize"
56fea3669SShri Abhyankar /*
66427ac75SLisandro Dalcin   TSEventInitialize - Initializes TSEvent for TSSolve
76fea3669SShri Abhyankar */
86427ac75SLisandro Dalcin PetscErrorCode TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U)
96fea3669SShri Abhyankar {
106fea3669SShri Abhyankar   PetscErrorCode ierr;
116fea3669SShri Abhyankar 
126fea3669SShri Abhyankar   PetscFunctionBegin;
136427ac75SLisandro Dalcin   if (!event) PetscFunctionReturn(0);
146427ac75SLisandro Dalcin   PetscValidPointer(event,1);
156427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,2);
166427ac75SLisandro Dalcin   PetscValidHeaderSpecific(U,VEC_CLASSID,4);
176fea3669SShri Abhyankar   event->ptime_prev = t;
1838bf2713SShri Abhyankar   event->iterctr = 0;
196427ac75SLisandro Dalcin   ierr = (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);CHKERRQ(ierr);
206fea3669SShri Abhyankar   PetscFunctionReturn(0);
216fea3669SShri Abhyankar }
226fea3669SShri Abhyankar 
236fea3669SShri Abhyankar #undef __FUNCT__
247dbe0728SLisandro Dalcin #define __FUNCT__ "TSEventDestroy"
257dbe0728SLisandro Dalcin PetscErrorCode TSEventDestroy(TSEvent *event)
267dbe0728SLisandro Dalcin {
277dbe0728SLisandro Dalcin   PetscErrorCode ierr;
287dbe0728SLisandro Dalcin   PetscInt       i;
297dbe0728SLisandro Dalcin 
307dbe0728SLisandro Dalcin   PetscFunctionBegin;
317dbe0728SLisandro Dalcin   PetscValidPointer(event,1);
327dbe0728SLisandro Dalcin   if (!*event) PetscFunctionReturn(0);
33*e7069c78SShri   if (--(*event)->refct > 0) {*event = 0; PetscFunctionReturn(0);}
34*e7069c78SShri 
357dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr);
367dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr);
37e2cdd850SShri Abhyankar   ierr = PetscFree((*event)->fvalue_right);CHKERRQ(ierr);
38e2cdd850SShri Abhyankar   ierr = PetscFree((*event)->zerocrossing);CHKERRQ(ierr);
39e2cdd850SShri Abhyankar   ierr = PetscFree((*event)->side);CHKERRQ(ierr);
407dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->direction);CHKERRQ(ierr);
417dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->terminate);CHKERRQ(ierr);
427dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr);
437dbe0728SLisandro Dalcin   ierr = PetscFree((*event)->vtol);CHKERRQ(ierr);
44a4ffd976SShri Abhyankar 
45a4ffd976SShri Abhyankar   for (i=0; i < (*event)->recsize; i++) {
467dbe0728SLisandro Dalcin     ierr = PetscFree((*event)->recorder.eventidx[i]);CHKERRQ(ierr);
477dbe0728SLisandro Dalcin   }
48a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.eventidx);CHKERRQ(ierr);
49a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.nevents);CHKERRQ(ierr);
50a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.stepnum);CHKERRQ(ierr);
51a4ffd976SShri Abhyankar   ierr = PetscFree((*event)->recorder.time);CHKERRQ(ierr);
52a4ffd976SShri Abhyankar 
537dbe0728SLisandro Dalcin   ierr = PetscViewerDestroy(&(*event)->monitor);CHKERRQ(ierr);
547dbe0728SLisandro Dalcin   ierr = PetscFree(*event);CHKERRQ(ierr);
557dbe0728SLisandro Dalcin   PetscFunctionReturn(0);
567dbe0728SLisandro Dalcin }
577dbe0728SLisandro Dalcin 
587dbe0728SLisandro Dalcin #undef __FUNCT__
59e3005195SShri Abhyankar #define __FUNCT__ "TSSetEventTolerances"
60e3005195SShri Abhyankar /*@
61e3005195SShri Abhyankar    TSSetEventTolerances - Set tolerances for event zero crossings when using event handler
62e3005195SShri Abhyankar 
63e3005195SShri Abhyankar    Logically Collective
64e3005195SShri Abhyankar 
65e3005195SShri Abhyankar    Input Arguments:
66e3005195SShri Abhyankar +  ts - time integration context
67e3005195SShri Abhyankar .  tol - scalar tolerance, PETSC_DECIDE to leave current value
68e3005195SShri Abhyankar -  vtol - array of tolerances or NULL, used in preference to tol if present
69e3005195SShri Abhyankar 
702589fa24SLisandro Dalcin    Options Database Keys:
712589fa24SLisandro Dalcin .  -ts_event_tol <tol> tolerance for event zero crossing
72e3005195SShri Abhyankar 
73e3005195SShri Abhyankar    Notes:
74f25fe674SLisandro Dalcin    Must call TSSetEventHandler() before setting the tolerances.
75e3005195SShri Abhyankar 
76e3005195SShri Abhyankar    The size of vtol is equal to the number of events.
77e3005195SShri Abhyankar 
78e3005195SShri Abhyankar    Level: beginner
79e3005195SShri Abhyankar 
80f25fe674SLisandro Dalcin .seealso: TS, TSEvent, TSSetEventHandler()
81e3005195SShri Abhyankar @*/
826427ac75SLisandro Dalcin PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[])
83e3005195SShri Abhyankar {
846427ac75SLisandro Dalcin   TSEvent        event;
85e3005195SShri Abhyankar   PetscInt       i;
86e3005195SShri Abhyankar 
87e3005195SShri Abhyankar   PetscFunctionBegin;
886427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
896427ac75SLisandro Dalcin   if (vtol) PetscValidRealPointer(vtol,3);
90f25fe674SLisandro Dalcin   if (!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()");
916427ac75SLisandro Dalcin 
926427ac75SLisandro Dalcin   event = ts->event;
93e3005195SShri Abhyankar   if (vtol) {
94e3005195SShri Abhyankar     for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i];
95e3005195SShri Abhyankar   } else {
96e3005195SShri Abhyankar     if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) {
97e3005195SShri Abhyankar       for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
98e3005195SShri Abhyankar     }
99e3005195SShri Abhyankar   }
100e3005195SShri Abhyankar   PetscFunctionReturn(0);
101e3005195SShri Abhyankar }
102e3005195SShri Abhyankar 
103e3005195SShri Abhyankar #undef __FUNCT__
104f25fe674SLisandro Dalcin #define __FUNCT__ "TSSetEventHandler"
1052dc7a7e3SShri Abhyankar /*@C
106f25fe674SLisandro Dalcin    TSSetEventHandler - Sets a monitoring function used for detecting events
1072dc7a7e3SShri Abhyankar 
1082dc7a7e3SShri Abhyankar    Logically Collective on TS
1092dc7a7e3SShri Abhyankar 
1102dc7a7e3SShri Abhyankar    Input Parameters:
1112dc7a7e3SShri Abhyankar +  ts - the TS context obtained from TSCreate()
1122dc7a7e3SShri Abhyankar .  nevents - number of local events
113d94325d3SShri Abhyankar .  direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction,
114d94325d3SShri Abhyankar                +1 => Zero crossing in positive direction, 0 => both ways (one for each event)
115d94325d3SShri Abhyankar .  terminate - flag to indicate whether time stepping should be terminated after
116d94325d3SShri Abhyankar                event is detected (one for each event)
1176427ac75SLisandro Dalcin .  eventhandler - event monitoring routine
1182dc7a7e3SShri Abhyankar .  postevent - [optional] post-event function
1192589fa24SLisandro Dalcin -  ctx       - [optional] user-defined context for private data for the
1202dc7a7e3SShri Abhyankar                event monitor and post event routine (use NULL if no
1212dc7a7e3SShri Abhyankar                context is desired)
1222dc7a7e3SShri Abhyankar 
1236427ac75SLisandro Dalcin    Calling sequence of eventhandler:
1243a88037aSBarry Smith    PetscErrorCode PetscEventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx)
1252dc7a7e3SShri Abhyankar 
1262dc7a7e3SShri Abhyankar    Input Parameters:
1272dc7a7e3SShri Abhyankar +  ts  - the TS context
1282dc7a7e3SShri Abhyankar .  t   - current time
1292dc7a7e3SShri Abhyankar .  U   - current iterate
1306427ac75SLisandro Dalcin -  ctx - [optional] context passed with eventhandler
1312dc7a7e3SShri Abhyankar 
1322dc7a7e3SShri Abhyankar    Output parameters:
133d94325d3SShri Abhyankar .  fvalue    - function value of events at time t
1342dc7a7e3SShri Abhyankar 
1352dc7a7e3SShri Abhyankar    Calling sequence of postevent:
1362589fa24SLisandro Dalcin    PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
1372dc7a7e3SShri Abhyankar 
1382dc7a7e3SShri Abhyankar    Input Parameters:
1392dc7a7e3SShri Abhyankar +  ts - the TS context
1402dc7a7e3SShri Abhyankar .  nevents_zero - number of local events whose event function is zero
1412dc7a7e3SShri Abhyankar .  events_zero  - indices of local events which have reached zero
1422dc7a7e3SShri Abhyankar .  t            - current time
1432dc7a7e3SShri Abhyankar .  U            - current solution
144031fbad4SShri Abhyankar .  forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0)
1456427ac75SLisandro Dalcin -  ctx          - the context passed with eventhandler
1462dc7a7e3SShri Abhyankar 
1472dc7a7e3SShri Abhyankar    Level: intermediate
1482dc7a7e3SShri Abhyankar 
14902749585SLisandro Dalcin .keywords: TS, event, set
1502dc7a7e3SShri Abhyankar 
1512dc7a7e3SShri Abhyankar .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason()
1522dc7a7e3SShri Abhyankar @*/
1532589fa24SLisandro 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)
1542dc7a7e3SShri Abhyankar {
1552dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
1562dc7a7e3SShri Abhyankar   TSEvent        event;
157d94325d3SShri Abhyankar   PetscInt       i;
158006e6a18SShri Abhyankar   PetscBool      flg;
159a6c783d2SShri Abhyankar #if defined PETSC_USE_REAL_SINGLE
160a6c783d2SShri Abhyankar   PetscReal      tol=1e-4;
161a6c783d2SShri Abhyankar #else
162d569cc17SSatish Balay   PetscReal      tol=1e-6;
163a6c783d2SShri Abhyankar #endif
1642dc7a7e3SShri Abhyankar 
1652dc7a7e3SShri Abhyankar   PetscFunctionBegin;
1666427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1676427ac75SLisandro Dalcin   PetscValidIntPointer(direction,2);
1686427ac75SLisandro Dalcin   PetscValidIntPointer(terminate,3);
1696427ac75SLisandro Dalcin 
1706427ac75SLisandro Dalcin   ierr = PetscNewLog(ts,&event);CHKERRQ(ierr);
171854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->fvalue);CHKERRQ(ierr);
172854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->fvalue_prev);CHKERRQ(ierr);
173e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->fvalue_right);CHKERRQ(ierr);
174e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->zerocrossing);CHKERRQ(ierr);
175e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->side);CHKERRQ(ierr);
176854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->direction);CHKERRQ(ierr);
177854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->terminate);CHKERRQ(ierr);
178e3005195SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->vtol);CHKERRQ(ierr);
179d94325d3SShri Abhyankar   for (i=0; i < nevents; i++) {
180d94325d3SShri Abhyankar     event->direction[i] = direction[i];
181d94325d3SShri Abhyankar     event->terminate[i] = terminate[i];
182e2cdd850SShri Abhyankar     event->zerocrossing[i] = PETSC_FALSE;
183e2cdd850SShri Abhyankar     event->side[i] = 0;
184d94325d3SShri Abhyankar   }
185854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->events_zero);CHKERRQ(ierr);
1862589fa24SLisandro Dalcin   event->nevents = nevents;
1876427ac75SLisandro Dalcin   event->eventhandler = eventhandler;
1882dc7a7e3SShri Abhyankar   event->postevent = postevent;
1896427ac75SLisandro Dalcin   event->ctx = ctx;
1902dc7a7e3SShri Abhyankar 
19102749585SLisandro Dalcin   event->recsize = 8;  /* Initial size of the recorder */
192a9514d71SShri Abhyankar   ierr = PetscOptionsBegin(((PetscObject)ts)->comm,"","TS Event options","");CHKERRQ(ierr);
1932dc7a7e3SShri Abhyankar   {
19402749585SLisandro Dalcin     ierr = PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","TSSetEventTolerances",tol,&tol,NULL);CHKERRQ(ierr);
195006e6a18SShri Abhyankar     ierr = PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);CHKERRQ(ierr);
196a4ffd976SShri Abhyankar     ierr = PetscOptionsInt("-ts_event_recorder_initial_size","Initial size of event recorder","",event->recsize,&event->recsize,NULL);CHKERRQ(ierr);
1972dc7a7e3SShri Abhyankar   }
1989e12be75SShri Abhyankar   PetscOptionsEnd();
199a4ffd976SShri Abhyankar 
200a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.time);CHKERRQ(ierr);
201a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.stepnum);CHKERRQ(ierr);
202a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.nevents);CHKERRQ(ierr);
203a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.eventidx);CHKERRQ(ierr);
204a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
205a4ffd976SShri Abhyankar     ierr = PetscMalloc1(event->nevents,&event->recorder.eventidx[i]);CHKERRQ(ierr);
206a4ffd976SShri Abhyankar   }
207*e7069c78SShri   /* Initialize the event recorder */
208*e7069c78SShri   event->recorder.ctr = 0;
209a4ffd976SShri Abhyankar 
210e3005195SShri Abhyankar   for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
2116427ac75SLisandro Dalcin   if (flg) {ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);CHKERRQ(ierr);}
2129e12be75SShri Abhyankar 
2136427ac75SLisandro Dalcin   ierr = TSEventDestroy(&ts->event);CHKERRQ(ierr);
214d94325d3SShri Abhyankar   ts->event = event;
215*e7069c78SShri   ts->event->refct = 1;
2162dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
2172dc7a7e3SShri Abhyankar }
2182dc7a7e3SShri Abhyankar 
219a4ffd976SShri Abhyankar #undef __FUNCT__
220a4ffd976SShri Abhyankar #define __FUNCT__ "TSEventRecorderResize"
221a4ffd976SShri Abhyankar /*
222a4ffd976SShri Abhyankar   TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
223a4ffd976SShri Abhyankar                           is reached.
224a4ffd976SShri Abhyankar */
22502749585SLisandro Dalcin static PetscErrorCode TSEventRecorderResize(TSEvent event)
226a4ffd976SShri Abhyankar {
227a4ffd976SShri Abhyankar   PetscErrorCode ierr;
228a4ffd976SShri Abhyankar   PetscReal      *time;
229a4ffd976SShri Abhyankar   PetscInt       *stepnum;
230a4ffd976SShri Abhyankar   PetscInt       *nevents;
231a4ffd976SShri Abhyankar   PetscInt       **eventidx;
232a4ffd976SShri Abhyankar   PetscInt       i,fact=2;
233a4ffd976SShri Abhyankar 
234a4ffd976SShri Abhyankar   PetscFunctionBegin;
235a4ffd976SShri Abhyankar 
236a4ffd976SShri Abhyankar   /* Create large arrays */
237a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&time);CHKERRQ(ierr);
238a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&stepnum);CHKERRQ(ierr);
239a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&nevents);CHKERRQ(ierr);
240a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&eventidx);CHKERRQ(ierr);
241a4ffd976SShri Abhyankar   for (i=0; i < fact*event->recsize; i++) {
242a4ffd976SShri Abhyankar     ierr = PetscMalloc1(event->nevents,&eventidx[i]);CHKERRQ(ierr);
243a4ffd976SShri Abhyankar   }
244a4ffd976SShri Abhyankar 
245a4ffd976SShri Abhyankar   /* Copy over data */
246a4ffd976SShri Abhyankar   ierr = PetscMemcpy(time,event->recorder.time,event->recsize*sizeof(PetscReal));CHKERRQ(ierr);
247a4ffd976SShri Abhyankar   ierr = PetscMemcpy(stepnum,event->recorder.stepnum,event->recsize*sizeof(PetscInt));CHKERRQ(ierr);
248a4ffd976SShri Abhyankar   ierr = PetscMemcpy(nevents,event->recorder.nevents,event->recsize*sizeof(PetscInt));CHKERRQ(ierr);
249a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
250a4ffd976SShri Abhyankar     ierr = PetscMemcpy(eventidx[i],event->recorder.eventidx[i],event->recorder.nevents[i]*sizeof(PetscInt));CHKERRQ(ierr);
251a4ffd976SShri Abhyankar   }
252a4ffd976SShri Abhyankar 
253a4ffd976SShri Abhyankar   /* Destroy old arrays */
254a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
255a4ffd976SShri Abhyankar     ierr = PetscFree(event->recorder.eventidx[i]);CHKERRQ(ierr);
256a4ffd976SShri Abhyankar   }
257a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.eventidx);CHKERRQ(ierr);
258a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.nevents);CHKERRQ(ierr);
259a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.stepnum);CHKERRQ(ierr);
260a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.time);CHKERRQ(ierr);
261a4ffd976SShri Abhyankar 
262a4ffd976SShri Abhyankar   /* Set pointers */
263a4ffd976SShri Abhyankar   event->recorder.time = time;
264a4ffd976SShri Abhyankar   event->recorder.stepnum = stepnum;
265a4ffd976SShri Abhyankar   event->recorder.nevents = nevents;
266a4ffd976SShri Abhyankar   event->recorder.eventidx = eventidx;
267a4ffd976SShri Abhyankar 
268a4ffd976SShri Abhyankar   /* Double size */
269a4ffd976SShri Abhyankar   event->recsize *= fact;
270a4ffd976SShri Abhyankar 
271a4ffd976SShri Abhyankar   PetscFunctionReturn(0);
272a4ffd976SShri Abhyankar }
273a4ffd976SShri Abhyankar 
274031fbad4SShri Abhyankar /*
2754597913aSLisandro Dalcin    Helper rutine to handle user postenvents and recording
276031fbad4SShri Abhyankar */
277031fbad4SShri Abhyankar #undef __FUNCT__
278031fbad4SShri Abhyankar #define __FUNCT__ "TSPostEvent"
2794597913aSLisandro Dalcin static PetscErrorCode TSPostEvent(TS ts,PetscReal t,Vec U)
2802dc7a7e3SShri Abhyankar {
2812dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
2822dc7a7e3SShri Abhyankar   TSEvent        event = ts->event;
2832dc7a7e3SShri Abhyankar   PetscBool      terminate = PETSC_FALSE;
28428d5b5d6SLisandro Dalcin   PetscBool      restart = PETSC_FALSE;
285d0578d90SShri Abhyankar   PetscInt       i,ctr,stepnum;
28628d5b5d6SLisandro Dalcin   PetscBool      ts_terminate,ts_restart;
2874597913aSLisandro Dalcin   PetscBool      forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */
2882dc7a7e3SShri Abhyankar 
2892dc7a7e3SShri Abhyankar   PetscFunctionBegin;
2902dc7a7e3SShri Abhyankar   if (event->postevent) {
29128d5b5d6SLisandro Dalcin     PetscObjectState state_prev,state_post;
29228d5b5d6SLisandro Dalcin     ierr = PetscObjectStateGet((PetscObject)U,&state_prev);CHKERRQ(ierr);
2934597913aSLisandro Dalcin     ierr = (*event->postevent)(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
29428d5b5d6SLisandro Dalcin     ierr = PetscObjectStateGet((PetscObject)U,&state_post);CHKERRQ(ierr);
29528d5b5d6SLisandro Dalcin     if (state_prev != state_post) restart = PETSC_TRUE;
2962dc7a7e3SShri Abhyankar   }
2974597913aSLisandro Dalcin 
29828d5b5d6SLisandro Dalcin   /* Handle termination events and step restart */
29928d5b5d6SLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
300b2566f29SBarry Smith   ierr = MPIU_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr);
30128d5b5d6SLisandro Dalcin   if (ts_terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);}
30228d5b5d6SLisandro Dalcin   event->status = ts_terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;
30328d5b5d6SLisandro Dalcin   ierr = MPIU_Allreduce(&restart,&ts_restart,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr);
30428d5b5d6SLisandro Dalcin   if (ts_restart) ts->steprestart = PETSC_TRUE;
305f7aea88cSShri Abhyankar 
3064597913aSLisandro Dalcin   event->ptime_prev = t;
3074597913aSLisandro Dalcin   /* Reset event residual functions as states might get changed by the postevent callback */
3084597913aSLisandro Dalcin   if (event->postevent) {ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);}
3094597913aSLisandro Dalcin   /* Cache current event residual functions */
3104597913aSLisandro Dalcin   for (i=0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i];
3114597913aSLisandro Dalcin 
312d0578d90SShri Abhyankar   /* Record the event in the event recorder */
313d0578d90SShri Abhyankar   ierr = TSGetTimeStepNumber(ts,&stepnum);CHKERRQ(ierr);
314f7aea88cSShri Abhyankar   ctr = event->recorder.ctr;
315a4ffd976SShri Abhyankar   if (ctr == event->recsize) {
316a4ffd976SShri Abhyankar     ierr = TSEventRecorderResize(event);CHKERRQ(ierr);
317a4ffd976SShri Abhyankar   }
318f7aea88cSShri Abhyankar   event->recorder.time[ctr] = t;
319d0578d90SShri Abhyankar   event->recorder.stepnum[ctr] = stepnum;
3204597913aSLisandro Dalcin   event->recorder.nevents[ctr] = event->nevents_zero;
3214597913aSLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
322f7aea88cSShri Abhyankar   event->recorder.ctr++;
3232dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
3242dc7a7e3SShri Abhyankar }
3252dc7a7e3SShri Abhyankar 
32602749585SLisandro Dalcin /* Uses Anderson-Bjorck variant of regula falsi method */
327a3a8645aSShri Abhyankar PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt)
32838bf2713SShri Abhyankar {
32902749585SLisandro Dalcin   PetscReal new_dt, scal = 1.0;
330e2cdd850SShri Abhyankar   if (PetscRealPart(fleft)*PetscRealPart(f) < 0) {
331e2cdd850SShri Abhyankar     if (side == 1) {
332a6c783d2SShri Abhyankar       scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright);
333a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
334e2cdd850SShri Abhyankar     }
33502749585SLisandro Dalcin     new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
336e2cdd850SShri Abhyankar   } else {
337e2cdd850SShri Abhyankar     if (side == -1) {
338a6c783d2SShri Abhyankar       scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft);
339a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
340e2cdd850SShri Abhyankar     }
34102749585SLisandro Dalcin     new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t;
342e2cdd850SShri Abhyankar   }
34302749585SLisandro Dalcin   return PetscMin(dt,new_dt);
34438bf2713SShri Abhyankar }
345e2cdd850SShri Abhyankar 
34638bf2713SShri Abhyankar 
3472dc7a7e3SShri Abhyankar #undef __FUNCT__
3486427ac75SLisandro Dalcin #define __FUNCT__ "TSEventHandler"
3496427ac75SLisandro Dalcin PetscErrorCode TSEventHandler(TS ts)
3502dc7a7e3SShri Abhyankar {
3512dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
3526427ac75SLisandro Dalcin   TSEvent        event;
3532dc7a7e3SShri Abhyankar   PetscReal      t;
3542dc7a7e3SShri Abhyankar   Vec            U;
3552dc7a7e3SShri Abhyankar   PetscInt       i;
3567dbe0728SLisandro Dalcin   PetscReal      dt,dt_min;
3572dc7a7e3SShri Abhyankar   PetscInt       rollback=0,in[2],out[2];
3589e12be75SShri Abhyankar   PetscInt       fvalue_sign,fvalueprev_sign;
3592dc7a7e3SShri Abhyankar 
3602dc7a7e3SShri Abhyankar   PetscFunctionBegin;
3616427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
3626427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
3636427ac75SLisandro Dalcin   event = ts->event;
3642dc7a7e3SShri Abhyankar 
3652dc7a7e3SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
3667dbe0728SLisandro Dalcin   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
3672dc7a7e3SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
3682dc7a7e3SShri Abhyankar 
3697dbe0728SLisandro Dalcin   if (event->status == TSEVENT_NONE) {
3702d5bd56dSLisandro Dalcin     if (ts->steps == 1) /* After first successful step */
3712d5bd56dSLisandro Dalcin       event->timestep_orig = ts->ptime - ts->ptime_prev;
3727dbe0728SLisandro Dalcin     event->timestep_prev = dt;
3737dbe0728SLisandro Dalcin   }
3747dbe0728SLisandro Dalcin 
3752dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_RESET_NEXTSTEP) {
3762d5bd56dSLisandro Dalcin     /* Restore time step */
3772d5bd56dSLisandro Dalcin     dt = PetscMin(event->timestep_orig,event->timestep_prev);
3787dbe0728SLisandro Dalcin     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
3792dc7a7e3SShri Abhyankar     event->status = TSEVENT_NONE;
3802dc7a7e3SShri Abhyankar   }
3812dc7a7e3SShri Abhyankar 
3822dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_NONE) {
3837dbe0728SLisandro Dalcin     event->ptime_end = t;
3842dc7a7e3SShri Abhyankar   }
3852dc7a7e3SShri Abhyankar 
3866427ac75SLisandro Dalcin   ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);
3879e12be75SShri Abhyankar 
3882dc7a7e3SShri Abhyankar   for (i=0; i < event->nevents; i++) {
389e3005195SShri Abhyankar     if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
3902dc7a7e3SShri Abhyankar       event->status = TSEVENT_ZERO;
391e2cdd850SShri Abhyankar       event->fvalue_right[i] = event->fvalue[i];
3929e12be75SShri Abhyankar       continue;
393006e6a18SShri Abhyankar     }
394d94325d3SShri Abhyankar     fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
395d94325d3SShri Abhyankar     fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
396e3005195SShri Abhyankar     if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) {
397d94325d3SShri Abhyankar       switch (event->direction[i]) {
398d94325d3SShri Abhyankar       case -1:
399d94325d3SShri Abhyankar         if (fvalue_sign < 0) {
4002dc7a7e3SShri Abhyankar           rollback = 1;
401e2cdd850SShri Abhyankar 
402e2cdd850SShri Abhyankar           /* Compute new time step */
403e2cdd850SShri 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);
404e2cdd850SShri Abhyankar 
4056427ac75SLisandro Dalcin           if (event->monitor) {
40638bf2713SShri 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);
407006e6a18SShri Abhyankar           }
408e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
409e2cdd850SShri Abhyankar           event->side[i] = 1;
410e2cdd850SShri Abhyankar 
411e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
412e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4132dc7a7e3SShri Abhyankar         }
414d94325d3SShri Abhyankar         break;
415d94325d3SShri Abhyankar       case 1:
416d94325d3SShri Abhyankar         if (fvalue_sign > 0) {
417d94325d3SShri Abhyankar           rollback = 1;
418e2cdd850SShri Abhyankar 
419e2cdd850SShri Abhyankar           /* Compute new time step */
420e2cdd850SShri 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);
421e2cdd850SShri Abhyankar 
4226427ac75SLisandro Dalcin           if (event->monitor) {
42338bf2713SShri 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);
424006e6a18SShri Abhyankar           }
425e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
426e2cdd850SShri Abhyankar           event->side[i] = 1;
427e2cdd850SShri Abhyankar 
428e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
429e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4302dc7a7e3SShri Abhyankar         }
431d94325d3SShri Abhyankar         break;
432d94325d3SShri Abhyankar       case 0:
433d94325d3SShri Abhyankar         rollback = 1;
434e2cdd850SShri Abhyankar 
435e2cdd850SShri Abhyankar         /* Compute new time step */
436e2cdd850SShri 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);
437e2cdd850SShri Abhyankar 
4386427ac75SLisandro Dalcin         if (event->monitor) {
43938bf2713SShri 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);
440006e6a18SShri Abhyankar         }
441e2cdd850SShri Abhyankar         event->fvalue_right[i] = event->fvalue[i];
442e2cdd850SShri Abhyankar         event->side[i] = 1;
443e2cdd850SShri Abhyankar 
444e2cdd850SShri Abhyankar         if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
445e2cdd850SShri Abhyankar         event->status = TSEVENT_LOCATED_INTERVAL;
446e2cdd850SShri Abhyankar 
447d94325d3SShri Abhyankar         break;
448d94325d3SShri Abhyankar       }
449d94325d3SShri Abhyankar     }
450d94325d3SShri Abhyankar   }
4517dbe0728SLisandro Dalcin 
4522589fa24SLisandro Dalcin   in[0] = event->status; in[1] = rollback;
4537dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
4542589fa24SLisandro Dalcin   event->status = (TSEventStatus)out[0]; rollback = out[1];
4552589fa24SLisandro Dalcin   if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
4562dc7a7e3SShri Abhyankar 
4577dbe0728SLisandro Dalcin   event->nevents_zero = 0;
4589e12be75SShri Abhyankar   if (event->status == TSEVENT_ZERO) {
4599e12be75SShri Abhyankar     for (i=0; i < event->nevents; i++) {
4609e12be75SShri Abhyankar       if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
4619e12be75SShri Abhyankar         event->events_zero[event->nevents_zero++] = i;
4626427ac75SLisandro Dalcin         if (event->monitor) {
46338bf2713SShri 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);
4649e12be75SShri Abhyankar         }
465e2cdd850SShri Abhyankar         event->zerocrossing[i] = PETSC_FALSE;
4669e12be75SShri Abhyankar       }
467e2cdd850SShri Abhyankar       event->side[i] = 0;
4689e12be75SShri Abhyankar     }
4694597913aSLisandro Dalcin     ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr);
4709e12be75SShri Abhyankar 
471e2cdd850SShri Abhyankar     dt = event->ptime_end - t;
4722d5bd56dSLisandro Dalcin     if (PetscAbsReal(dt) < PETSC_SMALL) dt += PetscMin(event->timestep_orig,event->timestep_prev); /* XXX Should be done better */
4739e12be75SShri Abhyankar     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
47438bf2713SShri Abhyankar     event->iterctr = 0;
4759e12be75SShri Abhyankar     PetscFunctionReturn(0);
4769e12be75SShri Abhyankar   }
4779e12be75SShri Abhyankar 
4782dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_LOCATED_INTERVAL) {
4792dc7a7e3SShri Abhyankar     ierr = TSRollBack(ts);CHKERRQ(ierr);
4802589fa24SLisandro Dalcin     ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr);
4812dc7a7e3SShri Abhyankar     event->status = TSEVENT_PROCESSING;
482e2cdd850SShri Abhyankar     event->ptime_right = t;
4832dc7a7e3SShri Abhyankar   } else {
4842dc7a7e3SShri Abhyankar     for (i=0; i < event->nevents; i++) {
485e2cdd850SShri Abhyankar       if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) {
486e2cdd850SShri Abhyankar         /* Compute new time step */
487e2cdd850SShri 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);
488e2cdd850SShri Abhyankar         event->side[i] = -1;
489e2cdd850SShri Abhyankar       }
4902dc7a7e3SShri Abhyankar       event->fvalue_prev[i] = event->fvalue[i];
4912dc7a7e3SShri Abhyankar     }
492e2cdd850SShri Abhyankar     if (event->monitor && event->status == TSEVENT_PROCESSING) {
49338bf2713SShri 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);
4942dc7a7e3SShri Abhyankar     }
49538bf2713SShri Abhyankar     event->ptime_prev = t;
49638bf2713SShri Abhyankar   }
49738bf2713SShri Abhyankar 
49838bf2713SShri Abhyankar   if (event->status == TSEVENT_PROCESSING) event->iterctr++;
4999e12be75SShri Abhyankar 
5007dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
5017dbe0728SLisandro Dalcin   ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr);
5022dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
5032dc7a7e3SShri Abhyankar }
5042dc7a7e3SShri Abhyankar 
505d0578d90SShri Abhyankar #undef __FUNCT__
5066427ac75SLisandro Dalcin #define __FUNCT__ "TSAdjointEventHandler"
5076427ac75SLisandro Dalcin PetscErrorCode TSAdjointEventHandler(TS ts)
508d0578d90SShri Abhyankar {
509d0578d90SShri Abhyankar   PetscErrorCode ierr;
5106427ac75SLisandro Dalcin   TSEvent        event;
511d0578d90SShri Abhyankar   PetscReal      t;
512d0578d90SShri Abhyankar   Vec            U;
513d0578d90SShri Abhyankar   PetscInt       ctr;
514d0578d90SShri Abhyankar   PetscBool      forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
515d0578d90SShri Abhyankar 
516d0578d90SShri Abhyankar   PetscFunctionBegin;
5176427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
5186427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
5196427ac75SLisandro Dalcin   event = ts->event;
520d0578d90SShri Abhyankar 
521d0578d90SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
522d0578d90SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
523d0578d90SShri Abhyankar 
524d0578d90SShri Abhyankar   ctr = event->recorder.ctr-1;
525bcbf8bb3SShri Abhyankar   if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
526d0578d90SShri Abhyankar     /* Call the user postevent function */
527d0578d90SShri Abhyankar     if (event->postevent) {
5286427ac75SLisandro Dalcin       ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
529d0578d90SShri Abhyankar       event->recorder.ctr--;
530d0578d90SShri Abhyankar     }
531d0578d90SShri Abhyankar   }
532d0578d90SShri Abhyankar 
533d0578d90SShri Abhyankar   PetscBarrier((PetscObject)ts);
534d0578d90SShri Abhyankar   PetscFunctionReturn(0);
535d0578d90SShri Abhyankar }
536