xref: /petsc/src/ts/event/tsevent.c (revision 458122a4cabef3ef10f725d657ce3f74347b02da)
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 /*@
55*458122a4SShri Abhyankar   TSSetPostEventIntervalStep - Set the time-step immediately following the event interval
56*458122a4SShri Abhyankar 
57*458122a4SShri Abhyankar   Logically Collective
58*458122a4SShri Abhyankar 
59*458122a4SShri Abhyankar   Input Arguments:
60*458122a4SShri Abhyankar + ts - time integration context
61*458122a4SShri Abhyankar - dt - post event interval step
62*458122a4SShri Abhyankar 
63*458122a4SShri Abhyankar   Options Database Keys:
64*458122a4SShri Abhyankar . -ts_event_post_eventinterval_step <dt> time-step after event interval
65*458122a4SShri Abhyankar 
66*458122a4SShri Abhyankar   Notes:
67*458122a4SShri Abhyankar   TSSetPostEventIntervalStep allows to set a time-step immediately following an event interval.
68*458122a4SShri Abhyankar 
69*458122a4SShri Abhyankar   This function should be called from the postevent function set with TSSetEventHandler().
70*458122a4SShri Abhyankar 
71*458122a4SShri Abhyankar   The post event interval time-step should be selected based on the dynamics following the event.
72*458122a4SShri Abhyankar   If the dynamics are stiff, a conservative (small) step should be used.
73*458122a4SShri Abhyankar   If not, then a larger time-step can be used.
74*458122a4SShri Abhyankar 
75*458122a4SShri Abhyankar   Level: Advanced
76*458122a4SShri Abhyankar   .seealso: TS, TSEvent, TSSetEventHandler()
77*458122a4SShri Abhyankar @*/
78*458122a4SShri Abhyankar PetscErrorCode TSSetPostEventIntervalStep(TS ts,PetscReal dt)
79*458122a4SShri Abhyankar {
80*458122a4SShri Abhyankar   PetscFunctionBegin;
81*458122a4SShri Abhyankar   ts->event->timestep_posteventinterval = dt;
82*458122a4SShri Abhyankar   PetscFunctionReturn(0);
83*458122a4SShri Abhyankar }
84*458122a4SShri Abhyankar 
85*458122a4SShri Abhyankar /*@
86e3005195SShri Abhyankar    TSSetEventTolerances - Set tolerances for event zero crossings when using event handler
87e3005195SShri Abhyankar 
88e3005195SShri Abhyankar    Logically Collective
89e3005195SShri Abhyankar 
90e3005195SShri Abhyankar    Input Arguments:
91e3005195SShri Abhyankar +  ts - time integration context
92e3005195SShri Abhyankar .  tol - scalar tolerance, PETSC_DECIDE to leave current value
93e3005195SShri Abhyankar -  vtol - array of tolerances or NULL, used in preference to tol if present
94e3005195SShri Abhyankar 
952589fa24SLisandro Dalcin    Options Database Keys:
962589fa24SLisandro Dalcin .  -ts_event_tol <tol> tolerance for event zero crossing
97e3005195SShri Abhyankar 
98e3005195SShri Abhyankar    Notes:
99f25fe674SLisandro Dalcin    Must call TSSetEventHandler() before setting the tolerances.
100e3005195SShri Abhyankar 
101e3005195SShri Abhyankar    The size of vtol is equal to the number of events.
102e3005195SShri Abhyankar 
103e3005195SShri Abhyankar    Level: beginner
104e3005195SShri Abhyankar 
105f25fe674SLisandro Dalcin .seealso: TS, TSEvent, TSSetEventHandler()
106e3005195SShri Abhyankar @*/
1076427ac75SLisandro Dalcin PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[])
108e3005195SShri Abhyankar {
1096427ac75SLisandro Dalcin   TSEvent        event;
110e3005195SShri Abhyankar   PetscInt       i;
111e3005195SShri Abhyankar 
112e3005195SShri Abhyankar   PetscFunctionBegin;
1136427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1146427ac75SLisandro Dalcin   if (vtol) PetscValidRealPointer(vtol,3);
115f25fe674SLisandro Dalcin   if (!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()");
1166427ac75SLisandro Dalcin 
1176427ac75SLisandro Dalcin   event = ts->event;
118e3005195SShri Abhyankar   if (vtol) {
119e3005195SShri Abhyankar     for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i];
120e3005195SShri Abhyankar   } else {
121e3005195SShri Abhyankar     if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) {
122e3005195SShri Abhyankar       for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
123e3005195SShri Abhyankar     }
124e3005195SShri Abhyankar   }
125e3005195SShri Abhyankar   PetscFunctionReturn(0);
126e3005195SShri Abhyankar }
127e3005195SShri Abhyankar 
1282dc7a7e3SShri Abhyankar /*@C
129f25fe674SLisandro Dalcin    TSSetEventHandler - Sets a monitoring function used for detecting events
1302dc7a7e3SShri Abhyankar 
1312dc7a7e3SShri Abhyankar    Logically Collective on TS
1322dc7a7e3SShri Abhyankar 
1332dc7a7e3SShri Abhyankar    Input Parameters:
1342dc7a7e3SShri Abhyankar +  ts - the TS context obtained from TSCreate()
1352dc7a7e3SShri Abhyankar .  nevents - number of local events
136d94325d3SShri Abhyankar .  direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction,
137d94325d3SShri Abhyankar                +1 => Zero crossing in positive direction, 0 => both ways (one for each event)
138d94325d3SShri Abhyankar .  terminate - flag to indicate whether time stepping should be terminated after
139d94325d3SShri Abhyankar                event is detected (one for each event)
1406427ac75SLisandro Dalcin .  eventhandler - event monitoring routine
1412dc7a7e3SShri Abhyankar .  postevent - [optional] post-event function
1422589fa24SLisandro Dalcin -  ctx       - [optional] user-defined context for private data for the
1432dc7a7e3SShri Abhyankar                event monitor and post event routine (use NULL if no
1442dc7a7e3SShri Abhyankar                context is desired)
1452dc7a7e3SShri Abhyankar 
1466427ac75SLisandro Dalcin    Calling sequence of eventhandler:
1473a88037aSBarry Smith    PetscErrorCode PetscEventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx)
1482dc7a7e3SShri Abhyankar 
1492dc7a7e3SShri Abhyankar    Input Parameters:
1502dc7a7e3SShri Abhyankar +  ts  - the TS context
1512dc7a7e3SShri Abhyankar .  t   - current time
1522dc7a7e3SShri Abhyankar .  U   - current iterate
1536427ac75SLisandro Dalcin -  ctx - [optional] context passed with eventhandler
1542dc7a7e3SShri Abhyankar 
1552dc7a7e3SShri Abhyankar    Output parameters:
156d94325d3SShri Abhyankar .  fvalue    - function value of events at time t
1572dc7a7e3SShri Abhyankar 
1582dc7a7e3SShri Abhyankar    Calling sequence of postevent:
1592589fa24SLisandro Dalcin    PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
1602dc7a7e3SShri Abhyankar 
1612dc7a7e3SShri Abhyankar    Input Parameters:
1622dc7a7e3SShri Abhyankar +  ts - the TS context
1632dc7a7e3SShri Abhyankar .  nevents_zero - number of local events whose event function is zero
1642dc7a7e3SShri Abhyankar .  events_zero  - indices of local events which have reached zero
1652dc7a7e3SShri Abhyankar .  t            - current time
1662dc7a7e3SShri Abhyankar .  U            - current solution
167031fbad4SShri Abhyankar .  forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0)
1686427ac75SLisandro Dalcin -  ctx          - the context passed with eventhandler
1692dc7a7e3SShri Abhyankar 
1702dc7a7e3SShri Abhyankar    Level: intermediate
1712dc7a7e3SShri Abhyankar 
17202749585SLisandro Dalcin .keywords: TS, event, set
1732dc7a7e3SShri Abhyankar 
1742dc7a7e3SShri Abhyankar .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason()
1752dc7a7e3SShri Abhyankar @*/
1762589fa24SLisandro 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)
1772dc7a7e3SShri Abhyankar {
1782dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
1792dc7a7e3SShri Abhyankar   TSEvent        event;
180d94325d3SShri Abhyankar   PetscInt       i;
181006e6a18SShri Abhyankar   PetscBool      flg;
182a6c783d2SShri Abhyankar #if defined PETSC_USE_REAL_SINGLE
183a6c783d2SShri Abhyankar   PetscReal      tol=1e-4;
184a6c783d2SShri Abhyankar #else
185d569cc17SSatish Balay   PetscReal      tol=1e-6;
186a6c783d2SShri Abhyankar #endif
1872dc7a7e3SShri Abhyankar 
1882dc7a7e3SShri Abhyankar   PetscFunctionBegin;
1896427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1900a82b154SShri   if (nevents) {
1916427ac75SLisandro Dalcin     PetscValidIntPointer(direction,2);
1926427ac75SLisandro Dalcin     PetscValidIntPointer(terminate,3);
1930a82b154SShri   }
1946427ac75SLisandro Dalcin 
1956427ac75SLisandro Dalcin   ierr = PetscNewLog(ts,&event);CHKERRQ(ierr);
196854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->fvalue);CHKERRQ(ierr);
197854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->fvalue_prev);CHKERRQ(ierr);
198e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->fvalue_right);CHKERRQ(ierr);
199e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->zerocrossing);CHKERRQ(ierr);
200e2cdd850SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->side);CHKERRQ(ierr);
201854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->direction);CHKERRQ(ierr);
202854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->terminate);CHKERRQ(ierr);
203e3005195SShri Abhyankar   ierr = PetscMalloc1(nevents,&event->vtol);CHKERRQ(ierr);
204d94325d3SShri Abhyankar   for (i=0; i < nevents; i++) {
205d94325d3SShri Abhyankar     event->direction[i] = direction[i];
206d94325d3SShri Abhyankar     event->terminate[i] = terminate[i];
207e2cdd850SShri Abhyankar     event->zerocrossing[i] = PETSC_FALSE;
208e2cdd850SShri Abhyankar     event->side[i] = 0;
209d94325d3SShri Abhyankar   }
210854ce69bSBarry Smith   ierr = PetscMalloc1(nevents,&event->events_zero);CHKERRQ(ierr);
2112589fa24SLisandro Dalcin   event->nevents = nevents;
2126427ac75SLisandro Dalcin   event->eventhandler = eventhandler;
2132dc7a7e3SShri Abhyankar   event->postevent = postevent;
2146427ac75SLisandro Dalcin   event->ctx = ctx;
215*458122a4SShri Abhyankar   event->timestep_posteventinterval = ts->time_step;
2162dc7a7e3SShri Abhyankar 
21702749585SLisandro Dalcin   event->recsize = 8;  /* Initial size of the recorder */
2185fbebc36SLisandro Dalcin   ierr = PetscOptionsBegin(((PetscObject)ts)->comm,((PetscObject)ts)->prefix,"TS Event options","TS");CHKERRQ(ierr);
2192dc7a7e3SShri Abhyankar   {
22002749585SLisandro Dalcin     ierr = PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","TSSetEventTolerances",tol,&tol,NULL);CHKERRQ(ierr);
221006e6a18SShri Abhyankar     ierr = PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);CHKERRQ(ierr);
222a4ffd976SShri Abhyankar     ierr = PetscOptionsInt("-ts_event_recorder_initial_size","Initial size of event recorder","",event->recsize,&event->recsize,NULL);CHKERRQ(ierr);
223*458122a4SShri Abhyankar     ierr = PetscOptionsReal("-ts_event_post_eventinterval_step","Time step after event interval","",event->timestep_posteventinterval,&event->timestep_posteventinterval,NULL);CHKERRQ(ierr);
2242dc7a7e3SShri Abhyankar   }
225b6518c14SStefano Zampini   ierr = PetscOptionsEnd();CHKERRQ(ierr);
226a4ffd976SShri Abhyankar 
227a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.time);CHKERRQ(ierr);
228a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.stepnum);CHKERRQ(ierr);
229a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.nevents);CHKERRQ(ierr);
230a4ffd976SShri Abhyankar   ierr = PetscMalloc1(event->recsize,&event->recorder.eventidx);CHKERRQ(ierr);
231a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
232a4ffd976SShri Abhyankar     ierr = PetscMalloc1(event->nevents,&event->recorder.eventidx[i]);CHKERRQ(ierr);
233a4ffd976SShri Abhyankar   }
234e7069c78SShri   /* Initialize the event recorder */
235e7069c78SShri   event->recorder.ctr = 0;
236a4ffd976SShri Abhyankar 
237e3005195SShri Abhyankar   for (i=0; i < event->nevents; i++) event->vtol[i] = tol;
2386427ac75SLisandro Dalcin   if (flg) {ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);CHKERRQ(ierr);}
2399e12be75SShri Abhyankar 
2406427ac75SLisandro Dalcin   ierr = TSEventDestroy(&ts->event);CHKERRQ(ierr);
241d94325d3SShri Abhyankar   ts->event = event;
242e7069c78SShri   ts->event->refct = 1;
2432dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
2442dc7a7e3SShri Abhyankar }
2452dc7a7e3SShri Abhyankar 
246a4ffd976SShri Abhyankar /*
247a4ffd976SShri Abhyankar   TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize)
248a4ffd976SShri Abhyankar                           is reached.
249a4ffd976SShri Abhyankar */
25002749585SLisandro Dalcin static PetscErrorCode TSEventRecorderResize(TSEvent event)
251a4ffd976SShri Abhyankar {
252a4ffd976SShri Abhyankar   PetscErrorCode ierr;
253a4ffd976SShri Abhyankar   PetscReal      *time;
254a4ffd976SShri Abhyankar   PetscInt       *stepnum;
255a4ffd976SShri Abhyankar   PetscInt       *nevents;
256a4ffd976SShri Abhyankar   PetscInt       **eventidx;
257a4ffd976SShri Abhyankar   PetscInt       i,fact=2;
258a4ffd976SShri Abhyankar 
259a4ffd976SShri Abhyankar   PetscFunctionBegin;
260a4ffd976SShri Abhyankar 
261a4ffd976SShri Abhyankar   /* Create large arrays */
262a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&time);CHKERRQ(ierr);
263a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&stepnum);CHKERRQ(ierr);
264a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&nevents);CHKERRQ(ierr);
265a4ffd976SShri Abhyankar   ierr = PetscMalloc1(fact*event->recsize,&eventidx);CHKERRQ(ierr);
266a4ffd976SShri Abhyankar   for (i=0; i < fact*event->recsize; i++) {
267a4ffd976SShri Abhyankar     ierr = PetscMalloc1(event->nevents,&eventidx[i]);CHKERRQ(ierr);
268a4ffd976SShri Abhyankar   }
269a4ffd976SShri Abhyankar 
270a4ffd976SShri Abhyankar   /* Copy over data */
271a4ffd976SShri Abhyankar   ierr = PetscMemcpy(time,event->recorder.time,event->recsize*sizeof(PetscReal));CHKERRQ(ierr);
272a4ffd976SShri Abhyankar   ierr = PetscMemcpy(stepnum,event->recorder.stepnum,event->recsize*sizeof(PetscInt));CHKERRQ(ierr);
273a4ffd976SShri Abhyankar   ierr = PetscMemcpy(nevents,event->recorder.nevents,event->recsize*sizeof(PetscInt));CHKERRQ(ierr);
274a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
275a4ffd976SShri Abhyankar     ierr = PetscMemcpy(eventidx[i],event->recorder.eventidx[i],event->recorder.nevents[i]*sizeof(PetscInt));CHKERRQ(ierr);
276a4ffd976SShri Abhyankar   }
277a4ffd976SShri Abhyankar 
278a4ffd976SShri Abhyankar   /* Destroy old arrays */
279a4ffd976SShri Abhyankar   for (i=0; i < event->recsize; i++) {
280a4ffd976SShri Abhyankar     ierr = PetscFree(event->recorder.eventidx[i]);CHKERRQ(ierr);
281a4ffd976SShri Abhyankar   }
282a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.eventidx);CHKERRQ(ierr);
283a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.nevents);CHKERRQ(ierr);
284a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.stepnum);CHKERRQ(ierr);
285a4ffd976SShri Abhyankar   ierr = PetscFree(event->recorder.time);CHKERRQ(ierr);
286a4ffd976SShri Abhyankar 
287a4ffd976SShri Abhyankar   /* Set pointers */
288a4ffd976SShri Abhyankar   event->recorder.time = time;
289a4ffd976SShri Abhyankar   event->recorder.stepnum = stepnum;
290a4ffd976SShri Abhyankar   event->recorder.nevents = nevents;
291a4ffd976SShri Abhyankar   event->recorder.eventidx = eventidx;
292a4ffd976SShri Abhyankar 
293a4ffd976SShri Abhyankar   /* Double size */
294a4ffd976SShri Abhyankar   event->recsize *= fact;
295a4ffd976SShri Abhyankar 
296a4ffd976SShri Abhyankar   PetscFunctionReturn(0);
297a4ffd976SShri Abhyankar }
298a4ffd976SShri Abhyankar 
299031fbad4SShri Abhyankar /*
3004597913aSLisandro Dalcin    Helper rutine to handle user postenvents and recording
301031fbad4SShri Abhyankar */
3024597913aSLisandro Dalcin static PetscErrorCode TSPostEvent(TS ts,PetscReal t,Vec U)
3032dc7a7e3SShri Abhyankar {
3042dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
3052dc7a7e3SShri Abhyankar   TSEvent        event = ts->event;
3062dc7a7e3SShri Abhyankar   PetscBool      terminate = PETSC_FALSE;
30728d5b5d6SLisandro Dalcin   PetscBool      restart = PETSC_FALSE;
308d0578d90SShri Abhyankar   PetscInt       i,ctr,stepnum;
3097324a0ffSLisandro Dalcin   PetscBool      inflag[2],outflag[2];
3104597913aSLisandro Dalcin   PetscBool      forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */
3112dc7a7e3SShri Abhyankar 
3122dc7a7e3SShri Abhyankar   PetscFunctionBegin;
3132dc7a7e3SShri Abhyankar   if (event->postevent) {
31428d5b5d6SLisandro Dalcin     PetscObjectState state_prev,state_post;
31528d5b5d6SLisandro Dalcin     ierr = PetscObjectStateGet((PetscObject)U,&state_prev);CHKERRQ(ierr);
3164597913aSLisandro Dalcin     ierr = (*event->postevent)(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
31728d5b5d6SLisandro Dalcin     ierr = PetscObjectStateGet((PetscObject)U,&state_post);CHKERRQ(ierr);
31828d5b5d6SLisandro Dalcin     if (state_prev != state_post) restart = PETSC_TRUE;
3192dc7a7e3SShri Abhyankar   }
3204597913aSLisandro Dalcin 
32128d5b5d6SLisandro Dalcin   /* Handle termination events and step restart */
32228d5b5d6SLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE;
3237324a0ffSLisandro Dalcin   inflag[0] = restart; inflag[1] = terminate;
3247324a0ffSLisandro Dalcin   ierr = MPIU_Allreduce(inflag,outflag,2,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr);
3257324a0ffSLisandro Dalcin   restart = outflag[0]; terminate = outflag[1];
3267324a0ffSLisandro Dalcin   if (restart) {ierr = TSRestartStep(ts);CHKERRQ(ierr);}
3277324a0ffSLisandro Dalcin   if (terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);}
3287324a0ffSLisandro Dalcin   event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP;
329f7aea88cSShri Abhyankar 
3304597913aSLisandro Dalcin   /* Reset event residual functions as states might get changed by the postevent callback */
331f443add6SLisandro Dalcin   if (event->postevent) {
332f443add6SLisandro Dalcin     ierr = VecLockPush(U);CHKERRQ(ierr);
333f443add6SLisandro Dalcin     ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);
334f443add6SLisandro Dalcin     ierr = VecLockPop(U);CHKERRQ(ierr);
335f443add6SLisandro Dalcin   }
336f443add6SLisandro Dalcin 
337f443add6SLisandro Dalcin   /* Cache current time and event residual functions */
338f443add6SLisandro Dalcin   event->ptime_prev = t;
339f443add6SLisandro Dalcin   for (i=0; i<event->nevents; i++)
340f443add6SLisandro Dalcin     event->fvalue_prev[i] = event->fvalue[i];
3414597913aSLisandro Dalcin 
342d0578d90SShri Abhyankar   /* Record the event in the event recorder */
34380275a0aSLisandro Dalcin   ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr);
344f7aea88cSShri Abhyankar   ctr = event->recorder.ctr;
345a4ffd976SShri Abhyankar   if (ctr == event->recsize) {
346a4ffd976SShri Abhyankar     ierr = TSEventRecorderResize(event);CHKERRQ(ierr);
347a4ffd976SShri Abhyankar   }
348f7aea88cSShri Abhyankar   event->recorder.time[ctr] = t;
349d0578d90SShri Abhyankar   event->recorder.stepnum[ctr] = stepnum;
3504597913aSLisandro Dalcin   event->recorder.nevents[ctr] = event->nevents_zero;
3514597913aSLisandro Dalcin   for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i];
352f7aea88cSShri Abhyankar   event->recorder.ctr++;
3532dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
3542dc7a7e3SShri Abhyankar }
3552dc7a7e3SShri Abhyankar 
35602749585SLisandro Dalcin /* Uses Anderson-Bjorck variant of regula falsi method */
357a3a8645aSShri Abhyankar PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt)
35838bf2713SShri Abhyankar {
35902749585SLisandro Dalcin   PetscReal new_dt, scal = 1.0;
360e2cdd850SShri Abhyankar   if (PetscRealPart(fleft)*PetscRealPart(f) < 0) {
361e2cdd850SShri Abhyankar     if (side == 1) {
362a6c783d2SShri Abhyankar       scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright);
363a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
364e2cdd850SShri Abhyankar     }
36502749585SLisandro Dalcin     new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft;
366e2cdd850SShri Abhyankar   } else {
367e2cdd850SShri Abhyankar     if (side == -1) {
368a6c783d2SShri Abhyankar       scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft);
369a6c783d2SShri Abhyankar       if (scal < PETSC_SMALL) scal = 0.5;
370e2cdd850SShri Abhyankar     }
37102749585SLisandro Dalcin     new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t;
372e2cdd850SShri Abhyankar   }
37302749585SLisandro Dalcin   return PetscMin(dt,new_dt);
37438bf2713SShri Abhyankar }
375e2cdd850SShri Abhyankar 
37638bf2713SShri Abhyankar 
3776427ac75SLisandro Dalcin PetscErrorCode TSEventHandler(TS ts)
3782dc7a7e3SShri Abhyankar {
3792dc7a7e3SShri Abhyankar   PetscErrorCode ierr;
3806427ac75SLisandro Dalcin   TSEvent        event;
3812dc7a7e3SShri Abhyankar   PetscReal      t;
3822dc7a7e3SShri Abhyankar   Vec            U;
3832dc7a7e3SShri Abhyankar   PetscInt       i;
3847dbe0728SLisandro Dalcin   PetscReal      dt,dt_min;
3852dc7a7e3SShri Abhyankar   PetscInt       rollback=0,in[2],out[2];
3869e12be75SShri Abhyankar   PetscInt       fvalue_sign,fvalueprev_sign;
3872dc7a7e3SShri Abhyankar 
3882dc7a7e3SShri Abhyankar   PetscFunctionBegin;
3896427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
3906427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
3916427ac75SLisandro Dalcin   event = ts->event;
3922dc7a7e3SShri Abhyankar 
3932dc7a7e3SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
3947dbe0728SLisandro Dalcin   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
3952dc7a7e3SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
3962dc7a7e3SShri Abhyankar 
3977dbe0728SLisandro Dalcin   if (event->status == TSEVENT_NONE) {
3987dbe0728SLisandro Dalcin     event->timestep_prev = dt;
3997dbe0728SLisandro Dalcin   }
4007dbe0728SLisandro Dalcin 
4012dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_RESET_NEXTSTEP) {
402*458122a4SShri Abhyankar     dt = event->timestep_posteventinterval;
4037dbe0728SLisandro Dalcin     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
4042dc7a7e3SShri Abhyankar     event->status = TSEVENT_NONE;
4052dc7a7e3SShri Abhyankar   }
4062dc7a7e3SShri Abhyankar 
4072dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_NONE) {
4087dbe0728SLisandro Dalcin     event->ptime_end = t;
4092dc7a7e3SShri Abhyankar   }
4102dc7a7e3SShri Abhyankar 
411f443add6SLisandro Dalcin   ierr = VecLockPush(U);CHKERRQ(ierr);
4126427ac75SLisandro Dalcin   ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);
413f443add6SLisandro Dalcin   ierr = VecLockPop(U);CHKERRQ(ierr);
4149e12be75SShri Abhyankar 
4152dc7a7e3SShri Abhyankar   for (i=0; i < event->nevents; i++) {
416e3005195SShri Abhyankar     if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
4172dc7a7e3SShri Abhyankar       event->status = TSEVENT_ZERO;
418e2cdd850SShri Abhyankar       event->fvalue_right[i] = event->fvalue[i];
4199e12be75SShri Abhyankar       continue;
420006e6a18SShri Abhyankar     }
421d94325d3SShri Abhyankar     fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i]));
422d94325d3SShri Abhyankar     fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i]));
423e3005195SShri Abhyankar     if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) {
424d94325d3SShri Abhyankar       switch (event->direction[i]) {
425d94325d3SShri Abhyankar       case -1:
426d94325d3SShri Abhyankar         if (fvalue_sign < 0) {
4272dc7a7e3SShri Abhyankar           rollback = 1;
428e2cdd850SShri Abhyankar 
429e2cdd850SShri Abhyankar           /* Compute new time step */
430e2cdd850SShri 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);
431e2cdd850SShri Abhyankar 
4326427ac75SLisandro Dalcin           if (event->monitor) {
43338bf2713SShri 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);
434006e6a18SShri Abhyankar           }
435e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
436e2cdd850SShri Abhyankar           event->side[i] = 1;
437e2cdd850SShri Abhyankar 
438e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
439e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4402dc7a7e3SShri Abhyankar         }
441d94325d3SShri Abhyankar         break;
442d94325d3SShri Abhyankar       case 1:
443d94325d3SShri Abhyankar         if (fvalue_sign > 0) {
444d94325d3SShri Abhyankar           rollback = 1;
445e2cdd850SShri Abhyankar 
446e2cdd850SShri Abhyankar           /* Compute new time step */
447e2cdd850SShri 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);
448e2cdd850SShri Abhyankar 
4496427ac75SLisandro Dalcin           if (event->monitor) {
45038bf2713SShri 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);
451006e6a18SShri Abhyankar           }
452e2cdd850SShri Abhyankar           event->fvalue_right[i] = event->fvalue[i];
453e2cdd850SShri Abhyankar           event->side[i] = 1;
454e2cdd850SShri Abhyankar 
455e2cdd850SShri Abhyankar           if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
456e2cdd850SShri Abhyankar           event->status = TSEVENT_LOCATED_INTERVAL;
4572dc7a7e3SShri Abhyankar         }
458d94325d3SShri Abhyankar         break;
459d94325d3SShri Abhyankar       case 0:
460d94325d3SShri Abhyankar         rollback = 1;
461e2cdd850SShri Abhyankar 
462e2cdd850SShri Abhyankar         /* Compute new time step */
463e2cdd850SShri 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);
464e2cdd850SShri Abhyankar 
4656427ac75SLisandro Dalcin         if (event->monitor) {
46638bf2713SShri 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);
467006e6a18SShri Abhyankar         }
468e2cdd850SShri Abhyankar         event->fvalue_right[i] = event->fvalue[i];
469e2cdd850SShri Abhyankar         event->side[i] = 1;
470e2cdd850SShri Abhyankar 
471e2cdd850SShri Abhyankar         if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE;
472e2cdd850SShri Abhyankar         event->status = TSEVENT_LOCATED_INTERVAL;
473e2cdd850SShri Abhyankar 
474d94325d3SShri Abhyankar         break;
475d94325d3SShri Abhyankar       }
476d94325d3SShri Abhyankar     }
477d94325d3SShri Abhyankar   }
4787dbe0728SLisandro Dalcin 
4792589fa24SLisandro Dalcin   in[0] = event->status; in[1] = rollback;
4807dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
4812589fa24SLisandro Dalcin   event->status = (TSEventStatus)out[0]; rollback = out[1];
4822589fa24SLisandro Dalcin   if (rollback) event->status = TSEVENT_LOCATED_INTERVAL;
4832dc7a7e3SShri Abhyankar 
4847dbe0728SLisandro Dalcin   event->nevents_zero = 0;
4859e12be75SShri Abhyankar   if (event->status == TSEVENT_ZERO) {
4869e12be75SShri Abhyankar     for (i=0; i < event->nevents; i++) {
4879e12be75SShri Abhyankar       if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) {
4889e12be75SShri Abhyankar         event->events_zero[event->nevents_zero++] = i;
4896427ac75SLisandro Dalcin         if (event->monitor) {
49038bf2713SShri 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);
4919e12be75SShri Abhyankar         }
492e2cdd850SShri Abhyankar         event->zerocrossing[i] = PETSC_FALSE;
4939e12be75SShri Abhyankar       }
494e2cdd850SShri Abhyankar       event->side[i] = 0;
4959e12be75SShri Abhyankar     }
4964597913aSLisandro Dalcin     ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr);
4979e12be75SShri Abhyankar 
498e2cdd850SShri Abhyankar     dt = event->ptime_end - t;
499ad508104SStefano Zampini     if (PetscAbsReal(dt) < PETSC_SMALL) { /* we hit the event, continue with the candidate time step */
500ad508104SStefano Zampini       dt = event->timestep_prev;
501ad508104SStefano Zampini       event->status = TSEVENT_NONE;
502ad508104SStefano Zampini     }
5039e12be75SShri Abhyankar     ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
50438bf2713SShri Abhyankar     event->iterctr = 0;
5059e12be75SShri Abhyankar     PetscFunctionReturn(0);
5069e12be75SShri Abhyankar   }
5079e12be75SShri Abhyankar 
5082dc7a7e3SShri Abhyankar   if (event->status == TSEVENT_LOCATED_INTERVAL) {
5092dc7a7e3SShri Abhyankar     ierr = TSRollBack(ts);CHKERRQ(ierr);
5102589fa24SLisandro Dalcin     ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr);
5112dc7a7e3SShri Abhyankar     event->status = TSEVENT_PROCESSING;
512e2cdd850SShri Abhyankar     event->ptime_right = t;
5132dc7a7e3SShri Abhyankar   } else {
5142dc7a7e3SShri Abhyankar     for (i=0; i < event->nevents; i++) {
515e2cdd850SShri Abhyankar       if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) {
516e2cdd850SShri Abhyankar         /* Compute new time step */
517e2cdd850SShri 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);
518e2cdd850SShri Abhyankar         event->side[i] = -1;
519e2cdd850SShri Abhyankar       }
5202dc7a7e3SShri Abhyankar       event->fvalue_prev[i] = event->fvalue[i];
5212dc7a7e3SShri Abhyankar     }
522e2cdd850SShri Abhyankar     if (event->monitor && event->status == TSEVENT_PROCESSING) {
52338bf2713SShri 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);
5242dc7a7e3SShri Abhyankar     }
52538bf2713SShri Abhyankar     event->ptime_prev = t;
52638bf2713SShri Abhyankar   }
52738bf2713SShri Abhyankar 
52838bf2713SShri Abhyankar   if (event->status == TSEVENT_PROCESSING) event->iterctr++;
5299e12be75SShri Abhyankar 
5307dbe0728SLisandro Dalcin   ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
5317dbe0728SLisandro Dalcin   ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr);
5322dc7a7e3SShri Abhyankar   PetscFunctionReturn(0);
5332dc7a7e3SShri Abhyankar }
5342dc7a7e3SShri Abhyankar 
5356427ac75SLisandro Dalcin PetscErrorCode TSAdjointEventHandler(TS ts)
536d0578d90SShri Abhyankar {
537d0578d90SShri Abhyankar   PetscErrorCode ierr;
5386427ac75SLisandro Dalcin   TSEvent        event;
539d0578d90SShri Abhyankar   PetscReal      t;
540d0578d90SShri Abhyankar   Vec            U;
541d0578d90SShri Abhyankar   PetscInt       ctr;
542d0578d90SShri Abhyankar   PetscBool      forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */
543d0578d90SShri Abhyankar 
544d0578d90SShri Abhyankar   PetscFunctionBegin;
5456427ac75SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
5466427ac75SLisandro Dalcin   if (!ts->event) PetscFunctionReturn(0);
5476427ac75SLisandro Dalcin   event = ts->event;
548d0578d90SShri Abhyankar 
549d0578d90SShri Abhyankar   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
550d0578d90SShri Abhyankar   ierr = TSGetSolution(ts,&U);CHKERRQ(ierr);
551d0578d90SShri Abhyankar 
552d0578d90SShri Abhyankar   ctr = event->recorder.ctr-1;
553bcbf8bb3SShri Abhyankar   if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) {
554d0578d90SShri Abhyankar     /* Call the user postevent function */
555d0578d90SShri Abhyankar     if (event->postevent) {
5566427ac75SLisandro Dalcin       ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr);
557d0578d90SShri Abhyankar       event->recorder.ctr--;
558d0578d90SShri Abhyankar     }
559d0578d90SShri Abhyankar   }
560d0578d90SShri Abhyankar 
561d0578d90SShri Abhyankar   PetscBarrier((PetscObject)ts);
562d0578d90SShri Abhyankar   PetscFunctionReturn(0);
563d0578d90SShri Abhyankar }
564