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