1af0996ceSBarry Smith #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 22dc7a7e3SShri Abhyankar 36fea3669SShri Abhyankar /* 46427ac75SLisandro Dalcin TSEventInitialize - Initializes TSEvent for TSSolve 56fea3669SShri Abhyankar */ 66427ac75SLisandro Dalcin PetscErrorCode TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U) 76fea3669SShri Abhyankar { 86fea3669SShri Abhyankar PetscErrorCode ierr; 96fea3669SShri Abhyankar 106fea3669SShri Abhyankar PetscFunctionBegin; 116427ac75SLisandro Dalcin if (!event) PetscFunctionReturn(0); 126427ac75SLisandro Dalcin PetscValidPointer(event,1); 136427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,2); 146427ac75SLisandro Dalcin PetscValidHeaderSpecific(U,VEC_CLASSID,4); 156fea3669SShri Abhyankar event->ptime_prev = t; 1638bf2713SShri Abhyankar event->iterctr = 0; 176427ac75SLisandro Dalcin ierr = (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);CHKERRQ(ierr); 186fea3669SShri Abhyankar PetscFunctionReturn(0); 196fea3669SShri Abhyankar } 206fea3669SShri Abhyankar 217dbe0728SLisandro Dalcin PetscErrorCode TSEventDestroy(TSEvent *event) 227dbe0728SLisandro Dalcin { 237dbe0728SLisandro Dalcin PetscErrorCode ierr; 247dbe0728SLisandro Dalcin PetscInt i; 257dbe0728SLisandro Dalcin 267dbe0728SLisandro Dalcin PetscFunctionBegin; 277dbe0728SLisandro Dalcin PetscValidPointer(event,1); 287dbe0728SLisandro Dalcin if (!*event) PetscFunctionReturn(0); 29e7069c78SShri if (--(*event)->refct > 0) {*event = 0; PetscFunctionReturn(0);} 30e7069c78SShri 317dbe0728SLisandro Dalcin ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 327dbe0728SLisandro Dalcin ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 33e2cdd850SShri Abhyankar ierr = PetscFree((*event)->fvalue_right);CHKERRQ(ierr); 34e2cdd850SShri Abhyankar ierr = PetscFree((*event)->zerocrossing);CHKERRQ(ierr); 35e2cdd850SShri Abhyankar ierr = PetscFree((*event)->side);CHKERRQ(ierr); 367dbe0728SLisandro Dalcin ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 377dbe0728SLisandro Dalcin ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 387dbe0728SLisandro Dalcin ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 397dbe0728SLisandro Dalcin ierr = PetscFree((*event)->vtol);CHKERRQ(ierr); 40a4ffd976SShri Abhyankar 41a4ffd976SShri Abhyankar for (i=0; i < (*event)->recsize; i++) { 427dbe0728SLisandro Dalcin ierr = PetscFree((*event)->recorder.eventidx[i]);CHKERRQ(ierr); 437dbe0728SLisandro Dalcin } 44a4ffd976SShri Abhyankar ierr = PetscFree((*event)->recorder.eventidx);CHKERRQ(ierr); 45a4ffd976SShri Abhyankar ierr = PetscFree((*event)->recorder.nevents);CHKERRQ(ierr); 46a4ffd976SShri Abhyankar ierr = PetscFree((*event)->recorder.stepnum);CHKERRQ(ierr); 47a4ffd976SShri Abhyankar ierr = PetscFree((*event)->recorder.time);CHKERRQ(ierr); 48a4ffd976SShri Abhyankar 497dbe0728SLisandro Dalcin ierr = PetscViewerDestroy(&(*event)->monitor);CHKERRQ(ierr); 507dbe0728SLisandro Dalcin ierr = PetscFree(*event);CHKERRQ(ierr); 517dbe0728SLisandro Dalcin PetscFunctionReturn(0); 527dbe0728SLisandro Dalcin } 537dbe0728SLisandro Dalcin 54e3005195SShri Abhyankar /*@ 55e3005195SShri Abhyankar TSSetEventTolerances - Set tolerances for event zero crossings when using event handler 56e3005195SShri Abhyankar 57e3005195SShri Abhyankar Logically Collective 58e3005195SShri Abhyankar 59e3005195SShri Abhyankar Input Arguments: 60e3005195SShri Abhyankar + ts - time integration context 61e3005195SShri Abhyankar . tol - scalar tolerance, PETSC_DECIDE to leave current value 62e3005195SShri Abhyankar - vtol - array of tolerances or NULL, used in preference to tol if present 63e3005195SShri Abhyankar 642589fa24SLisandro Dalcin Options Database Keys: 652589fa24SLisandro Dalcin . -ts_event_tol <tol> tolerance for event zero crossing 66e3005195SShri Abhyankar 67e3005195SShri Abhyankar Notes: 68f25fe674SLisandro Dalcin Must call TSSetEventHandler() before setting the tolerances. 69e3005195SShri Abhyankar 70e3005195SShri Abhyankar The size of vtol is equal to the number of events. 71e3005195SShri Abhyankar 72e3005195SShri Abhyankar Level: beginner 73e3005195SShri Abhyankar 74f25fe674SLisandro Dalcin .seealso: TS, TSEvent, TSSetEventHandler() 75e3005195SShri Abhyankar @*/ 766427ac75SLisandro Dalcin PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[]) 77e3005195SShri Abhyankar { 786427ac75SLisandro Dalcin TSEvent event; 79e3005195SShri Abhyankar PetscInt i; 80e3005195SShri Abhyankar 81e3005195SShri Abhyankar PetscFunctionBegin; 826427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 836427ac75SLisandro Dalcin if (vtol) PetscValidRealPointer(vtol,3); 84f25fe674SLisandro Dalcin if (!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()"); 856427ac75SLisandro Dalcin 866427ac75SLisandro Dalcin event = ts->event; 87e3005195SShri Abhyankar if (vtol) { 88e3005195SShri Abhyankar for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i]; 89e3005195SShri Abhyankar } else { 90e3005195SShri Abhyankar if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) { 91e3005195SShri Abhyankar for (i=0; i < event->nevents; i++) event->vtol[i] = tol; 92e3005195SShri Abhyankar } 93e3005195SShri Abhyankar } 94e3005195SShri Abhyankar PetscFunctionReturn(0); 95e3005195SShri Abhyankar } 96e3005195SShri Abhyankar 972dc7a7e3SShri Abhyankar /*@C 98f25fe674SLisandro Dalcin TSSetEventHandler - Sets a monitoring function used for detecting events 992dc7a7e3SShri Abhyankar 1002dc7a7e3SShri Abhyankar Logically Collective on TS 1012dc7a7e3SShri Abhyankar 1022dc7a7e3SShri Abhyankar Input Parameters: 1032dc7a7e3SShri Abhyankar + ts - the TS context obtained from TSCreate() 1042dc7a7e3SShri Abhyankar . nevents - number of local events 105d94325d3SShri Abhyankar . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction, 106d94325d3SShri Abhyankar +1 => Zero crossing in positive direction, 0 => both ways (one for each event) 107d94325d3SShri Abhyankar . terminate - flag to indicate whether time stepping should be terminated after 108d94325d3SShri Abhyankar event is detected (one for each event) 1096427ac75SLisandro Dalcin . eventhandler - event monitoring routine 1102dc7a7e3SShri Abhyankar . postevent - [optional] post-event function 1112589fa24SLisandro Dalcin - ctx - [optional] user-defined context for private data for the 1122dc7a7e3SShri Abhyankar event monitor and post event routine (use NULL if no 1132dc7a7e3SShri Abhyankar context is desired) 1142dc7a7e3SShri Abhyankar 1156427ac75SLisandro Dalcin Calling sequence of eventhandler: 1163a88037aSBarry Smith PetscErrorCode PetscEventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx) 1172dc7a7e3SShri Abhyankar 1182dc7a7e3SShri Abhyankar Input Parameters: 1192dc7a7e3SShri Abhyankar + ts - the TS context 1202dc7a7e3SShri Abhyankar . t - current time 1212dc7a7e3SShri Abhyankar . U - current iterate 1226427ac75SLisandro Dalcin - ctx - [optional] context passed with eventhandler 1232dc7a7e3SShri Abhyankar 1242dc7a7e3SShri Abhyankar Output parameters: 125d94325d3SShri Abhyankar . fvalue - function value of events at time t 1262dc7a7e3SShri Abhyankar 1272dc7a7e3SShri Abhyankar Calling sequence of postevent: 1282589fa24SLisandro Dalcin PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx) 1292dc7a7e3SShri Abhyankar 1302dc7a7e3SShri Abhyankar Input Parameters: 1312dc7a7e3SShri Abhyankar + ts - the TS context 1322dc7a7e3SShri Abhyankar . nevents_zero - number of local events whose event function is zero 1332dc7a7e3SShri Abhyankar . events_zero - indices of local events which have reached zero 1342dc7a7e3SShri Abhyankar . t - current time 1352dc7a7e3SShri Abhyankar . U - current solution 136031fbad4SShri Abhyankar . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0) 1376427ac75SLisandro Dalcin - ctx - the context passed with eventhandler 1382dc7a7e3SShri Abhyankar 1392dc7a7e3SShri Abhyankar Level: intermediate 1402dc7a7e3SShri Abhyankar 14102749585SLisandro Dalcin .keywords: TS, event, set 1422dc7a7e3SShri Abhyankar 1432dc7a7e3SShri Abhyankar .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason() 1442dc7a7e3SShri Abhyankar @*/ 1452589fa24SLisandro Dalcin PetscErrorCode TSSetEventHandler(TS ts,PetscInt nevents,PetscInt direction[],PetscBool terminate[],PetscErrorCode (*eventhandler)(TS,PetscReal,Vec,PetscScalar[],void*),PetscErrorCode (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,PetscBool,void*),void *ctx) 1462dc7a7e3SShri Abhyankar { 1472dc7a7e3SShri Abhyankar PetscErrorCode ierr; 1482dc7a7e3SShri Abhyankar TSEvent event; 149d94325d3SShri Abhyankar PetscInt i; 150006e6a18SShri Abhyankar PetscBool flg; 151a6c783d2SShri Abhyankar #if defined PETSC_USE_REAL_SINGLE 152a6c783d2SShri Abhyankar PetscReal tol=1e-4; 153a6c783d2SShri Abhyankar #else 154d569cc17SSatish Balay PetscReal tol=1e-6; 155a6c783d2SShri Abhyankar #endif 1562dc7a7e3SShri Abhyankar 1572dc7a7e3SShri Abhyankar PetscFunctionBegin; 1586427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1590a82b154SShri if(nevents) { 1606427ac75SLisandro Dalcin PetscValidIntPointer(direction,2); 1616427ac75SLisandro Dalcin PetscValidIntPointer(terminate,3); 1620a82b154SShri } 1636427ac75SLisandro Dalcin 1646427ac75SLisandro Dalcin ierr = PetscNewLog(ts,&event);CHKERRQ(ierr); 165854ce69bSBarry Smith ierr = PetscMalloc1(nevents,&event->fvalue);CHKERRQ(ierr); 166854ce69bSBarry Smith ierr = PetscMalloc1(nevents,&event->fvalue_prev);CHKERRQ(ierr); 167e2cdd850SShri Abhyankar ierr = PetscMalloc1(nevents,&event->fvalue_right);CHKERRQ(ierr); 168e2cdd850SShri Abhyankar ierr = PetscMalloc1(nevents,&event->zerocrossing);CHKERRQ(ierr); 169e2cdd850SShri Abhyankar ierr = PetscMalloc1(nevents,&event->side);CHKERRQ(ierr); 170854ce69bSBarry Smith ierr = PetscMalloc1(nevents,&event->direction);CHKERRQ(ierr); 171854ce69bSBarry Smith ierr = PetscMalloc1(nevents,&event->terminate);CHKERRQ(ierr); 172e3005195SShri Abhyankar ierr = PetscMalloc1(nevents,&event->vtol);CHKERRQ(ierr); 173d94325d3SShri Abhyankar for (i=0; i < nevents; i++) { 174d94325d3SShri Abhyankar event->direction[i] = direction[i]; 175d94325d3SShri Abhyankar event->terminate[i] = terminate[i]; 176e2cdd850SShri Abhyankar event->zerocrossing[i] = PETSC_FALSE; 177e2cdd850SShri Abhyankar event->side[i] = 0; 178d94325d3SShri Abhyankar } 179854ce69bSBarry Smith ierr = PetscMalloc1(nevents,&event->events_zero);CHKERRQ(ierr); 1802589fa24SLisandro Dalcin event->nevents = nevents; 1816427ac75SLisandro Dalcin event->eventhandler = eventhandler; 1822dc7a7e3SShri Abhyankar event->postevent = postevent; 1836427ac75SLisandro Dalcin event->ctx = ctx; 1842dc7a7e3SShri Abhyankar 18502749585SLisandro Dalcin event->recsize = 8; /* Initial size of the recorder */ 186a9514d71SShri Abhyankar ierr = PetscOptionsBegin(((PetscObject)ts)->comm,"","TS Event options","");CHKERRQ(ierr); 1872dc7a7e3SShri Abhyankar { 18802749585SLisandro Dalcin ierr = PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","TSSetEventTolerances",tol,&tol,NULL);CHKERRQ(ierr); 189006e6a18SShri Abhyankar ierr = PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);CHKERRQ(ierr); 190a4ffd976SShri Abhyankar ierr = PetscOptionsInt("-ts_event_recorder_initial_size","Initial size of event recorder","",event->recsize,&event->recsize,NULL);CHKERRQ(ierr); 1912dc7a7e3SShri Abhyankar } 1929e12be75SShri Abhyankar PetscOptionsEnd(); 193a4ffd976SShri Abhyankar 194a4ffd976SShri Abhyankar ierr = PetscMalloc1(event->recsize,&event->recorder.time);CHKERRQ(ierr); 195a4ffd976SShri Abhyankar ierr = PetscMalloc1(event->recsize,&event->recorder.stepnum);CHKERRQ(ierr); 196a4ffd976SShri Abhyankar ierr = PetscMalloc1(event->recsize,&event->recorder.nevents);CHKERRQ(ierr); 197a4ffd976SShri Abhyankar ierr = PetscMalloc1(event->recsize,&event->recorder.eventidx);CHKERRQ(ierr); 198a4ffd976SShri Abhyankar for (i=0; i < event->recsize; i++) { 199a4ffd976SShri Abhyankar ierr = PetscMalloc1(event->nevents,&event->recorder.eventidx[i]);CHKERRQ(ierr); 200a4ffd976SShri Abhyankar } 201e7069c78SShri /* Initialize the event recorder */ 202e7069c78SShri event->recorder.ctr = 0; 203a4ffd976SShri Abhyankar 204e3005195SShri Abhyankar for (i=0; i < event->nevents; i++) event->vtol[i] = tol; 2056427ac75SLisandro Dalcin if (flg) {ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);CHKERRQ(ierr);} 2069e12be75SShri Abhyankar 2076427ac75SLisandro Dalcin ierr = TSEventDestroy(&ts->event);CHKERRQ(ierr); 208d94325d3SShri Abhyankar ts->event = event; 209e7069c78SShri ts->event->refct = 1; 2102dc7a7e3SShri Abhyankar PetscFunctionReturn(0); 2112dc7a7e3SShri Abhyankar } 2122dc7a7e3SShri Abhyankar 213a4ffd976SShri Abhyankar /* 214a4ffd976SShri Abhyankar TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize) 215a4ffd976SShri Abhyankar is reached. 216a4ffd976SShri Abhyankar */ 21702749585SLisandro Dalcin static PetscErrorCode TSEventRecorderResize(TSEvent event) 218a4ffd976SShri Abhyankar { 219a4ffd976SShri Abhyankar PetscErrorCode ierr; 220a4ffd976SShri Abhyankar PetscReal *time; 221a4ffd976SShri Abhyankar PetscInt *stepnum; 222a4ffd976SShri Abhyankar PetscInt *nevents; 223a4ffd976SShri Abhyankar PetscInt **eventidx; 224a4ffd976SShri Abhyankar PetscInt i,fact=2; 225a4ffd976SShri Abhyankar 226a4ffd976SShri Abhyankar PetscFunctionBegin; 227a4ffd976SShri Abhyankar 228a4ffd976SShri Abhyankar /* Create large arrays */ 229a4ffd976SShri Abhyankar ierr = PetscMalloc1(fact*event->recsize,&time);CHKERRQ(ierr); 230a4ffd976SShri Abhyankar ierr = PetscMalloc1(fact*event->recsize,&stepnum);CHKERRQ(ierr); 231a4ffd976SShri Abhyankar ierr = PetscMalloc1(fact*event->recsize,&nevents);CHKERRQ(ierr); 232a4ffd976SShri Abhyankar ierr = PetscMalloc1(fact*event->recsize,&eventidx);CHKERRQ(ierr); 233a4ffd976SShri Abhyankar for (i=0; i < fact*event->recsize; i++) { 234a4ffd976SShri Abhyankar ierr = PetscMalloc1(event->nevents,&eventidx[i]);CHKERRQ(ierr); 235a4ffd976SShri Abhyankar } 236a4ffd976SShri Abhyankar 237a4ffd976SShri Abhyankar /* Copy over data */ 238a4ffd976SShri Abhyankar ierr = PetscMemcpy(time,event->recorder.time,event->recsize*sizeof(PetscReal));CHKERRQ(ierr); 239a4ffd976SShri Abhyankar ierr = PetscMemcpy(stepnum,event->recorder.stepnum,event->recsize*sizeof(PetscInt));CHKERRQ(ierr); 240a4ffd976SShri Abhyankar ierr = PetscMemcpy(nevents,event->recorder.nevents,event->recsize*sizeof(PetscInt));CHKERRQ(ierr); 241a4ffd976SShri Abhyankar for (i=0; i < event->recsize; i++) { 242a4ffd976SShri Abhyankar ierr = PetscMemcpy(eventidx[i],event->recorder.eventidx[i],event->recorder.nevents[i]*sizeof(PetscInt));CHKERRQ(ierr); 243a4ffd976SShri Abhyankar } 244a4ffd976SShri Abhyankar 245a4ffd976SShri Abhyankar /* Destroy old arrays */ 246a4ffd976SShri Abhyankar for (i=0; i < event->recsize; i++) { 247a4ffd976SShri Abhyankar ierr = PetscFree(event->recorder.eventidx[i]);CHKERRQ(ierr); 248a4ffd976SShri Abhyankar } 249a4ffd976SShri Abhyankar ierr = PetscFree(event->recorder.eventidx);CHKERRQ(ierr); 250a4ffd976SShri Abhyankar ierr = PetscFree(event->recorder.nevents);CHKERRQ(ierr); 251a4ffd976SShri Abhyankar ierr = PetscFree(event->recorder.stepnum);CHKERRQ(ierr); 252a4ffd976SShri Abhyankar ierr = PetscFree(event->recorder.time);CHKERRQ(ierr); 253a4ffd976SShri Abhyankar 254a4ffd976SShri Abhyankar /* Set pointers */ 255a4ffd976SShri Abhyankar event->recorder.time = time; 256a4ffd976SShri Abhyankar event->recorder.stepnum = stepnum; 257a4ffd976SShri Abhyankar event->recorder.nevents = nevents; 258a4ffd976SShri Abhyankar event->recorder.eventidx = eventidx; 259a4ffd976SShri Abhyankar 260a4ffd976SShri Abhyankar /* Double size */ 261a4ffd976SShri Abhyankar event->recsize *= fact; 262a4ffd976SShri Abhyankar 263a4ffd976SShri Abhyankar PetscFunctionReturn(0); 264a4ffd976SShri Abhyankar } 265a4ffd976SShri Abhyankar 266031fbad4SShri Abhyankar /* 2674597913aSLisandro Dalcin Helper rutine to handle user postenvents and recording 268031fbad4SShri Abhyankar */ 2694597913aSLisandro Dalcin static PetscErrorCode TSPostEvent(TS ts,PetscReal t,Vec U) 2702dc7a7e3SShri Abhyankar { 2712dc7a7e3SShri Abhyankar PetscErrorCode ierr; 2722dc7a7e3SShri Abhyankar TSEvent event = ts->event; 2732dc7a7e3SShri Abhyankar PetscBool terminate = PETSC_FALSE; 27428d5b5d6SLisandro Dalcin PetscBool restart = PETSC_FALSE; 275d0578d90SShri Abhyankar PetscInt i,ctr,stepnum; 276*7324a0ffSLisandro Dalcin PetscBool inflag[2],outflag[2]; 2774597913aSLisandro Dalcin PetscBool forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */ 2782dc7a7e3SShri Abhyankar 2792dc7a7e3SShri Abhyankar PetscFunctionBegin; 2802dc7a7e3SShri Abhyankar if (event->postevent) { 28128d5b5d6SLisandro Dalcin PetscObjectState state_prev,state_post; 28228d5b5d6SLisandro Dalcin ierr = PetscObjectStateGet((PetscObject)U,&state_prev);CHKERRQ(ierr); 2834597913aSLisandro Dalcin ierr = (*event->postevent)(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 28428d5b5d6SLisandro Dalcin ierr = PetscObjectStateGet((PetscObject)U,&state_post);CHKERRQ(ierr); 28528d5b5d6SLisandro Dalcin if (state_prev != state_post) restart = PETSC_TRUE; 2862dc7a7e3SShri Abhyankar } 2874597913aSLisandro Dalcin 28828d5b5d6SLisandro Dalcin /* Handle termination events and step restart */ 28928d5b5d6SLisandro Dalcin for (i=0; i<event->nevents_zero; i++) if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE; 290*7324a0ffSLisandro Dalcin inflag[0] = restart; inflag[1] = terminate; 291*7324a0ffSLisandro Dalcin ierr = MPIU_Allreduce(inflag,outflag,2,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr); 292*7324a0ffSLisandro Dalcin restart = outflag[0]; terminate = outflag[1]; 293*7324a0ffSLisandro Dalcin if (restart) {ierr = TSRestartStep(ts);CHKERRQ(ierr);} 294*7324a0ffSLisandro Dalcin if (terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);} 295*7324a0ffSLisandro Dalcin event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP; 296f7aea88cSShri Abhyankar 2974597913aSLisandro Dalcin event->ptime_prev = t; 2984597913aSLisandro Dalcin /* Reset event residual functions as states might get changed by the postevent callback */ 2994597913aSLisandro Dalcin if (event->postevent) {ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr);} 3004597913aSLisandro Dalcin /* Cache current event residual functions */ 3014597913aSLisandro Dalcin for (i=0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i]; 3024597913aSLisandro Dalcin 303d0578d90SShri Abhyankar /* Record the event in the event recorder */ 30480275a0aSLisandro Dalcin ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr); 305f7aea88cSShri Abhyankar ctr = event->recorder.ctr; 306a4ffd976SShri Abhyankar if (ctr == event->recsize) { 307a4ffd976SShri Abhyankar ierr = TSEventRecorderResize(event);CHKERRQ(ierr); 308a4ffd976SShri Abhyankar } 309f7aea88cSShri Abhyankar event->recorder.time[ctr] = t; 310d0578d90SShri Abhyankar event->recorder.stepnum[ctr] = stepnum; 3114597913aSLisandro Dalcin event->recorder.nevents[ctr] = event->nevents_zero; 3124597913aSLisandro Dalcin for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i]; 313f7aea88cSShri Abhyankar event->recorder.ctr++; 3142dc7a7e3SShri Abhyankar PetscFunctionReturn(0); 3152dc7a7e3SShri Abhyankar } 3162dc7a7e3SShri Abhyankar 31702749585SLisandro Dalcin /* Uses Anderson-Bjorck variant of regula falsi method */ 318a3a8645aSShri Abhyankar PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt) 31938bf2713SShri Abhyankar { 32002749585SLisandro Dalcin PetscReal new_dt, scal = 1.0; 321e2cdd850SShri Abhyankar if (PetscRealPart(fleft)*PetscRealPart(f) < 0) { 322e2cdd850SShri Abhyankar if (side == 1) { 323a6c783d2SShri Abhyankar scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright); 324a6c783d2SShri Abhyankar if (scal < PETSC_SMALL) scal = 0.5; 325e2cdd850SShri Abhyankar } 32602749585SLisandro Dalcin new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft; 327e2cdd850SShri Abhyankar } else { 328e2cdd850SShri Abhyankar if (side == -1) { 329a6c783d2SShri Abhyankar scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft); 330a6c783d2SShri Abhyankar if (scal < PETSC_SMALL) scal = 0.5; 331e2cdd850SShri Abhyankar } 33202749585SLisandro Dalcin new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t; 333e2cdd850SShri Abhyankar } 33402749585SLisandro Dalcin return PetscMin(dt,new_dt); 33538bf2713SShri Abhyankar } 336e2cdd850SShri Abhyankar 33738bf2713SShri Abhyankar 3386427ac75SLisandro Dalcin PetscErrorCode TSEventHandler(TS ts) 3392dc7a7e3SShri Abhyankar { 3402dc7a7e3SShri Abhyankar PetscErrorCode ierr; 3416427ac75SLisandro Dalcin TSEvent event; 3422dc7a7e3SShri Abhyankar PetscReal t; 3432dc7a7e3SShri Abhyankar Vec U; 3442dc7a7e3SShri Abhyankar PetscInt i; 3457dbe0728SLisandro Dalcin PetscReal dt,dt_min; 3462dc7a7e3SShri Abhyankar PetscInt rollback=0,in[2],out[2]; 3479e12be75SShri Abhyankar PetscInt fvalue_sign,fvalueprev_sign; 3482dc7a7e3SShri Abhyankar 3492dc7a7e3SShri Abhyankar PetscFunctionBegin; 3506427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3516427ac75SLisandro Dalcin if (!ts->event) PetscFunctionReturn(0); 3526427ac75SLisandro Dalcin event = ts->event; 3532dc7a7e3SShri Abhyankar 3542dc7a7e3SShri Abhyankar ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 3557dbe0728SLisandro Dalcin ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 3562dc7a7e3SShri Abhyankar ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 3572dc7a7e3SShri Abhyankar 3587dbe0728SLisandro Dalcin if (event->status == TSEVENT_NONE) { 3592d5bd56dSLisandro Dalcin if (ts->steps == 1) /* After first successful step */ 3602d5bd56dSLisandro Dalcin event->timestep_orig = ts->ptime - ts->ptime_prev; 3617dbe0728SLisandro Dalcin event->timestep_prev = dt; 3627dbe0728SLisandro Dalcin } 3637dbe0728SLisandro Dalcin 3642dc7a7e3SShri Abhyankar if (event->status == TSEVENT_RESET_NEXTSTEP) { 3652d5bd56dSLisandro Dalcin /* Restore time step */ 3662d5bd56dSLisandro Dalcin dt = PetscMin(event->timestep_orig,event->timestep_prev); 3677dbe0728SLisandro Dalcin ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 3682dc7a7e3SShri Abhyankar event->status = TSEVENT_NONE; 3692dc7a7e3SShri Abhyankar } 3702dc7a7e3SShri Abhyankar 3712dc7a7e3SShri Abhyankar if (event->status == TSEVENT_NONE) { 3727dbe0728SLisandro Dalcin event->ptime_end = t; 3732dc7a7e3SShri Abhyankar } 3742dc7a7e3SShri Abhyankar 3756427ac75SLisandro Dalcin ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr); 3769e12be75SShri Abhyankar 3772dc7a7e3SShri Abhyankar for (i=0; i < event->nevents; i++) { 378e3005195SShri Abhyankar if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 3792dc7a7e3SShri Abhyankar event->status = TSEVENT_ZERO; 380e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 3819e12be75SShri Abhyankar continue; 382006e6a18SShri Abhyankar } 383d94325d3SShri Abhyankar fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 384d94325d3SShri Abhyankar fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 385e3005195SShri Abhyankar if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) { 386d94325d3SShri Abhyankar switch (event->direction[i]) { 387d94325d3SShri Abhyankar case -1: 388d94325d3SShri Abhyankar if (fvalue_sign < 0) { 3892dc7a7e3SShri Abhyankar rollback = 1; 390e2cdd850SShri Abhyankar 391e2cdd850SShri Abhyankar /* Compute new time step */ 392e2cdd850SShri Abhyankar dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt); 393e2cdd850SShri Abhyankar 3946427ac75SLisandro Dalcin if (event->monitor) { 39538bf2713SShri Abhyankar ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 396006e6a18SShri Abhyankar } 397e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 398e2cdd850SShri Abhyankar event->side[i] = 1; 399e2cdd850SShri Abhyankar 400e2cdd850SShri Abhyankar if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 401e2cdd850SShri Abhyankar event->status = TSEVENT_LOCATED_INTERVAL; 4022dc7a7e3SShri Abhyankar } 403d94325d3SShri Abhyankar break; 404d94325d3SShri Abhyankar case 1: 405d94325d3SShri Abhyankar if (fvalue_sign > 0) { 406d94325d3SShri Abhyankar rollback = 1; 407e2cdd850SShri Abhyankar 408e2cdd850SShri Abhyankar /* Compute new time step */ 409e2cdd850SShri Abhyankar dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt); 410e2cdd850SShri Abhyankar 4116427ac75SLisandro Dalcin if (event->monitor) { 41238bf2713SShri Abhyankar ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 413006e6a18SShri Abhyankar } 414e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 415e2cdd850SShri Abhyankar event->side[i] = 1; 416e2cdd850SShri Abhyankar 417e2cdd850SShri Abhyankar if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 418e2cdd850SShri Abhyankar event->status = TSEVENT_LOCATED_INTERVAL; 4192dc7a7e3SShri Abhyankar } 420d94325d3SShri Abhyankar break; 421d94325d3SShri Abhyankar case 0: 422d94325d3SShri Abhyankar rollback = 1; 423e2cdd850SShri Abhyankar 424e2cdd850SShri Abhyankar /* Compute new time step */ 425e2cdd850SShri Abhyankar dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt); 426e2cdd850SShri Abhyankar 4276427ac75SLisandro Dalcin if (event->monitor) { 42838bf2713SShri Abhyankar ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 429006e6a18SShri Abhyankar } 430e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 431e2cdd850SShri Abhyankar event->side[i] = 1; 432e2cdd850SShri Abhyankar 433e2cdd850SShri Abhyankar if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 434e2cdd850SShri Abhyankar event->status = TSEVENT_LOCATED_INTERVAL; 435e2cdd850SShri Abhyankar 436d94325d3SShri Abhyankar break; 437d94325d3SShri Abhyankar } 438d94325d3SShri Abhyankar } 439d94325d3SShri Abhyankar } 4407dbe0728SLisandro Dalcin 4412589fa24SLisandro Dalcin in[0] = event->status; in[1] = rollback; 4427dbe0728SLisandro Dalcin ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4432589fa24SLisandro Dalcin event->status = (TSEventStatus)out[0]; rollback = out[1]; 4442589fa24SLisandro Dalcin if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 4452dc7a7e3SShri Abhyankar 4467dbe0728SLisandro Dalcin event->nevents_zero = 0; 4479e12be75SShri Abhyankar if (event->status == TSEVENT_ZERO) { 4489e12be75SShri Abhyankar for (i=0; i < event->nevents; i++) { 4499e12be75SShri Abhyankar if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 4509e12be75SShri Abhyankar event->events_zero[event->nevents_zero++] = i; 4516427ac75SLisandro Dalcin if (event->monitor) { 45238bf2713SShri Abhyankar ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: Event %D zero crossing at time %g located in %D iterations\n",i,(double)t,event->iterctr);CHKERRQ(ierr); 4539e12be75SShri Abhyankar } 454e2cdd850SShri Abhyankar event->zerocrossing[i] = PETSC_FALSE; 4559e12be75SShri Abhyankar } 456e2cdd850SShri Abhyankar event->side[i] = 0; 4579e12be75SShri Abhyankar } 4584597913aSLisandro Dalcin ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr); 4599e12be75SShri Abhyankar 460e2cdd850SShri Abhyankar dt = event->ptime_end - t; 4612d5bd56dSLisandro Dalcin if (PetscAbsReal(dt) < PETSC_SMALL) dt += PetscMin(event->timestep_orig,event->timestep_prev); /* XXX Should be done better */ 4629e12be75SShri Abhyankar ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 46338bf2713SShri Abhyankar event->iterctr = 0; 4649e12be75SShri Abhyankar PetscFunctionReturn(0); 4659e12be75SShri Abhyankar } 4669e12be75SShri Abhyankar 4672dc7a7e3SShri Abhyankar if (event->status == TSEVENT_LOCATED_INTERVAL) { 4682dc7a7e3SShri Abhyankar ierr = TSRollBack(ts);CHKERRQ(ierr); 4692589fa24SLisandro Dalcin ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr); 4702dc7a7e3SShri Abhyankar event->status = TSEVENT_PROCESSING; 471e2cdd850SShri Abhyankar event->ptime_right = t; 4722dc7a7e3SShri Abhyankar } else { 4732dc7a7e3SShri Abhyankar for (i=0; i < event->nevents; i++) { 474e2cdd850SShri Abhyankar if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) { 475e2cdd850SShri Abhyankar /* Compute new time step */ 476e2cdd850SShri Abhyankar dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],dt); 477e2cdd850SShri Abhyankar event->side[i] = -1; 478e2cdd850SShri Abhyankar } 4792dc7a7e3SShri Abhyankar event->fvalue_prev[i] = event->fvalue[i]; 4802dc7a7e3SShri Abhyankar } 481e2cdd850SShri Abhyankar if (event->monitor && event->status == TSEVENT_PROCESSING) { 48238bf2713SShri Abhyankar ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Stepping forward as no event detected in interval [%g - %g]\n",event->iterctr,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 4832dc7a7e3SShri Abhyankar } 48438bf2713SShri Abhyankar event->ptime_prev = t; 48538bf2713SShri Abhyankar } 48638bf2713SShri Abhyankar 48738bf2713SShri Abhyankar if (event->status == TSEVENT_PROCESSING) event->iterctr++; 4889e12be75SShri Abhyankar 4897dbe0728SLisandro Dalcin ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4907dbe0728SLisandro Dalcin ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr); 4912dc7a7e3SShri Abhyankar PetscFunctionReturn(0); 4922dc7a7e3SShri Abhyankar } 4932dc7a7e3SShri Abhyankar 4946427ac75SLisandro Dalcin PetscErrorCode TSAdjointEventHandler(TS ts) 495d0578d90SShri Abhyankar { 496d0578d90SShri Abhyankar PetscErrorCode ierr; 4976427ac75SLisandro Dalcin TSEvent event; 498d0578d90SShri Abhyankar PetscReal t; 499d0578d90SShri Abhyankar Vec U; 500d0578d90SShri Abhyankar PetscInt ctr; 501d0578d90SShri Abhyankar PetscBool forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */ 502d0578d90SShri Abhyankar 503d0578d90SShri Abhyankar PetscFunctionBegin; 5046427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5056427ac75SLisandro Dalcin if (!ts->event) PetscFunctionReturn(0); 5066427ac75SLisandro Dalcin event = ts->event; 507d0578d90SShri Abhyankar 508d0578d90SShri Abhyankar ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 509d0578d90SShri Abhyankar ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 510d0578d90SShri Abhyankar 511d0578d90SShri Abhyankar ctr = event->recorder.ctr-1; 512bcbf8bb3SShri Abhyankar if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) { 513d0578d90SShri Abhyankar /* Call the user postevent function */ 514d0578d90SShri Abhyankar if (event->postevent) { 5156427ac75SLisandro Dalcin ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 516d0578d90SShri Abhyankar event->recorder.ctr--; 517d0578d90SShri Abhyankar } 518d0578d90SShri Abhyankar } 519d0578d90SShri Abhyankar 520d0578d90SShri Abhyankar PetscBarrier((PetscObject)ts); 521d0578d90SShri Abhyankar PetscFunctionReturn(0); 522d0578d90SShri Abhyankar } 523