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