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; 2767324a0ffSLisandro 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; 2907324a0ffSLisandro Dalcin inflag[0] = restart; inflag[1] = terminate; 2917324a0ffSLisandro Dalcin ierr = MPIU_Allreduce(inflag,outflag,2,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr); 2927324a0ffSLisandro Dalcin restart = outflag[0]; terminate = outflag[1]; 2937324a0ffSLisandro Dalcin if (restart) {ierr = TSRestartStep(ts);CHKERRQ(ierr);} 2947324a0ffSLisandro Dalcin if (terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);} 2957324a0ffSLisandro Dalcin event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP; 296f7aea88cSShri Abhyankar 2974597913aSLisandro Dalcin /* Reset event residual functions as states might get changed by the postevent callback */ 298*f443add6SLisandro Dalcin if (event->postevent) { 299*f443add6SLisandro Dalcin ierr = VecLockPush(U);CHKERRQ(ierr); 300*f443add6SLisandro Dalcin ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr); 301*f443add6SLisandro Dalcin ierr = VecLockPop(U);CHKERRQ(ierr); 302*f443add6SLisandro Dalcin } 303*f443add6SLisandro Dalcin 304*f443add6SLisandro Dalcin /* Cache current time and event residual functions */ 305*f443add6SLisandro Dalcin event->ptime_prev = t; 306*f443add6SLisandro Dalcin for (i=0; i<event->nevents; i++) 307*f443add6SLisandro Dalcin event->fvalue_prev[i] = event->fvalue[i]; 3084597913aSLisandro Dalcin 309d0578d90SShri Abhyankar /* Record the event in the event recorder */ 31080275a0aSLisandro Dalcin ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr); 311f7aea88cSShri Abhyankar ctr = event->recorder.ctr; 312a4ffd976SShri Abhyankar if (ctr == event->recsize) { 313a4ffd976SShri Abhyankar ierr = TSEventRecorderResize(event);CHKERRQ(ierr); 314a4ffd976SShri Abhyankar } 315f7aea88cSShri Abhyankar event->recorder.time[ctr] = t; 316d0578d90SShri Abhyankar event->recorder.stepnum[ctr] = stepnum; 3174597913aSLisandro Dalcin event->recorder.nevents[ctr] = event->nevents_zero; 3184597913aSLisandro Dalcin for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i]; 319f7aea88cSShri Abhyankar event->recorder.ctr++; 3202dc7a7e3SShri Abhyankar PetscFunctionReturn(0); 3212dc7a7e3SShri Abhyankar } 3222dc7a7e3SShri Abhyankar 32302749585SLisandro Dalcin /* Uses Anderson-Bjorck variant of regula falsi method */ 324a3a8645aSShri Abhyankar PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt) 32538bf2713SShri Abhyankar { 32602749585SLisandro Dalcin PetscReal new_dt, scal = 1.0; 327e2cdd850SShri Abhyankar if (PetscRealPart(fleft)*PetscRealPart(f) < 0) { 328e2cdd850SShri Abhyankar if (side == 1) { 329a6c783d2SShri Abhyankar scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright); 330a6c783d2SShri Abhyankar if (scal < PETSC_SMALL) scal = 0.5; 331e2cdd850SShri Abhyankar } 33202749585SLisandro Dalcin new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft; 333e2cdd850SShri Abhyankar } else { 334e2cdd850SShri Abhyankar if (side == -1) { 335a6c783d2SShri Abhyankar scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft); 336a6c783d2SShri Abhyankar if (scal < PETSC_SMALL) scal = 0.5; 337e2cdd850SShri Abhyankar } 33802749585SLisandro Dalcin new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t; 339e2cdd850SShri Abhyankar } 34002749585SLisandro Dalcin return PetscMin(dt,new_dt); 34138bf2713SShri Abhyankar } 342e2cdd850SShri Abhyankar 34338bf2713SShri Abhyankar 3446427ac75SLisandro Dalcin PetscErrorCode TSEventHandler(TS ts) 3452dc7a7e3SShri Abhyankar { 3462dc7a7e3SShri Abhyankar PetscErrorCode ierr; 3476427ac75SLisandro Dalcin TSEvent event; 3482dc7a7e3SShri Abhyankar PetscReal t; 3492dc7a7e3SShri Abhyankar Vec U; 3502dc7a7e3SShri Abhyankar PetscInt i; 3517dbe0728SLisandro Dalcin PetscReal dt,dt_min; 3522dc7a7e3SShri Abhyankar PetscInt rollback=0,in[2],out[2]; 3539e12be75SShri Abhyankar PetscInt fvalue_sign,fvalueprev_sign; 3542dc7a7e3SShri Abhyankar 3552dc7a7e3SShri Abhyankar PetscFunctionBegin; 3566427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3576427ac75SLisandro Dalcin if (!ts->event) PetscFunctionReturn(0); 3586427ac75SLisandro Dalcin event = ts->event; 3592dc7a7e3SShri Abhyankar 3602dc7a7e3SShri Abhyankar ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 3617dbe0728SLisandro Dalcin ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 3622dc7a7e3SShri Abhyankar ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 3632dc7a7e3SShri Abhyankar 3647dbe0728SLisandro Dalcin if (event->status == TSEVENT_NONE) { 3652d5bd56dSLisandro Dalcin if (ts->steps == 1) /* After first successful step */ 3662d5bd56dSLisandro Dalcin event->timestep_orig = ts->ptime - ts->ptime_prev; 3677dbe0728SLisandro Dalcin event->timestep_prev = dt; 3687dbe0728SLisandro Dalcin } 3697dbe0728SLisandro Dalcin 3702dc7a7e3SShri Abhyankar if (event->status == TSEVENT_RESET_NEXTSTEP) { 3712d5bd56dSLisandro Dalcin /* Restore time step */ 3722d5bd56dSLisandro Dalcin dt = PetscMin(event->timestep_orig,event->timestep_prev); 3737dbe0728SLisandro Dalcin ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 3742dc7a7e3SShri Abhyankar event->status = TSEVENT_NONE; 3752dc7a7e3SShri Abhyankar } 3762dc7a7e3SShri Abhyankar 3772dc7a7e3SShri Abhyankar if (event->status == TSEVENT_NONE) { 3787dbe0728SLisandro Dalcin event->ptime_end = t; 3792dc7a7e3SShri Abhyankar } 3802dc7a7e3SShri Abhyankar 381*f443add6SLisandro Dalcin ierr = VecLockPush(U);CHKERRQ(ierr); 3826427ac75SLisandro Dalcin ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr); 383*f443add6SLisandro Dalcin ierr = VecLockPop(U);CHKERRQ(ierr); 3849e12be75SShri Abhyankar 3852dc7a7e3SShri Abhyankar for (i=0; i < event->nevents; i++) { 386e3005195SShri Abhyankar if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 3872dc7a7e3SShri Abhyankar event->status = TSEVENT_ZERO; 388e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 3899e12be75SShri Abhyankar continue; 390006e6a18SShri Abhyankar } 391d94325d3SShri Abhyankar fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 392d94325d3SShri Abhyankar fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 393e3005195SShri Abhyankar if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) { 394d94325d3SShri Abhyankar switch (event->direction[i]) { 395d94325d3SShri Abhyankar case -1: 396d94325d3SShri Abhyankar if (fvalue_sign < 0) { 3972dc7a7e3SShri Abhyankar rollback = 1; 398e2cdd850SShri Abhyankar 399e2cdd850SShri Abhyankar /* Compute new time step */ 400e2cdd850SShri 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); 401e2cdd850SShri Abhyankar 4026427ac75SLisandro Dalcin if (event->monitor) { 40338bf2713SShri 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); 404006e6a18SShri Abhyankar } 405e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 406e2cdd850SShri Abhyankar event->side[i] = 1; 407e2cdd850SShri Abhyankar 408e2cdd850SShri Abhyankar if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 409e2cdd850SShri Abhyankar event->status = TSEVENT_LOCATED_INTERVAL; 4102dc7a7e3SShri Abhyankar } 411d94325d3SShri Abhyankar break; 412d94325d3SShri Abhyankar case 1: 413d94325d3SShri Abhyankar if (fvalue_sign > 0) { 414d94325d3SShri Abhyankar rollback = 1; 415e2cdd850SShri Abhyankar 416e2cdd850SShri Abhyankar /* Compute new time step */ 417e2cdd850SShri 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); 418e2cdd850SShri Abhyankar 4196427ac75SLisandro Dalcin if (event->monitor) { 42038bf2713SShri 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); 421006e6a18SShri Abhyankar } 422e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 423e2cdd850SShri Abhyankar event->side[i] = 1; 424e2cdd850SShri Abhyankar 425e2cdd850SShri Abhyankar if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 426e2cdd850SShri Abhyankar event->status = TSEVENT_LOCATED_INTERVAL; 4272dc7a7e3SShri Abhyankar } 428d94325d3SShri Abhyankar break; 429d94325d3SShri Abhyankar case 0: 430d94325d3SShri Abhyankar rollback = 1; 431e2cdd850SShri Abhyankar 432e2cdd850SShri Abhyankar /* Compute new time step */ 433e2cdd850SShri 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); 434e2cdd850SShri Abhyankar 4356427ac75SLisandro Dalcin if (event->monitor) { 43638bf2713SShri 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); 437006e6a18SShri Abhyankar } 438e2cdd850SShri Abhyankar event->fvalue_right[i] = event->fvalue[i]; 439e2cdd850SShri Abhyankar event->side[i] = 1; 440e2cdd850SShri Abhyankar 441e2cdd850SShri Abhyankar if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 442e2cdd850SShri Abhyankar event->status = TSEVENT_LOCATED_INTERVAL; 443e2cdd850SShri Abhyankar 444d94325d3SShri Abhyankar break; 445d94325d3SShri Abhyankar } 446d94325d3SShri Abhyankar } 447d94325d3SShri Abhyankar } 4487dbe0728SLisandro Dalcin 4492589fa24SLisandro Dalcin in[0] = event->status; in[1] = rollback; 4507dbe0728SLisandro Dalcin ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4512589fa24SLisandro Dalcin event->status = (TSEventStatus)out[0]; rollback = out[1]; 4522589fa24SLisandro Dalcin if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 4532dc7a7e3SShri Abhyankar 4547dbe0728SLisandro Dalcin event->nevents_zero = 0; 4559e12be75SShri Abhyankar if (event->status == TSEVENT_ZERO) { 4569e12be75SShri Abhyankar for (i=0; i < event->nevents; i++) { 4579e12be75SShri Abhyankar if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 4589e12be75SShri Abhyankar event->events_zero[event->nevents_zero++] = i; 4596427ac75SLisandro Dalcin if (event->monitor) { 46038bf2713SShri 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); 4619e12be75SShri Abhyankar } 462e2cdd850SShri Abhyankar event->zerocrossing[i] = PETSC_FALSE; 4639e12be75SShri Abhyankar } 464e2cdd850SShri Abhyankar event->side[i] = 0; 4659e12be75SShri Abhyankar } 4664597913aSLisandro Dalcin ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr); 4679e12be75SShri Abhyankar 468e2cdd850SShri Abhyankar dt = event->ptime_end - t; 4692d5bd56dSLisandro Dalcin if (PetscAbsReal(dt) < PETSC_SMALL) dt += PetscMin(event->timestep_orig,event->timestep_prev); /* XXX Should be done better */ 4709e12be75SShri Abhyankar ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 47138bf2713SShri Abhyankar event->iterctr = 0; 4729e12be75SShri Abhyankar PetscFunctionReturn(0); 4739e12be75SShri Abhyankar } 4749e12be75SShri Abhyankar 4752dc7a7e3SShri Abhyankar if (event->status == TSEVENT_LOCATED_INTERVAL) { 4762dc7a7e3SShri Abhyankar ierr = TSRollBack(ts);CHKERRQ(ierr); 4772589fa24SLisandro Dalcin ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr); 4782dc7a7e3SShri Abhyankar event->status = TSEVENT_PROCESSING; 479e2cdd850SShri Abhyankar event->ptime_right = t; 4802dc7a7e3SShri Abhyankar } else { 4812dc7a7e3SShri Abhyankar for (i=0; i < event->nevents; i++) { 482e2cdd850SShri Abhyankar if (event->status == TSEVENT_PROCESSING && event->zerocrossing[i]) { 483e2cdd850SShri Abhyankar /* Compute new time step */ 484e2cdd850SShri 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); 485e2cdd850SShri Abhyankar event->side[i] = -1; 486e2cdd850SShri Abhyankar } 4872dc7a7e3SShri Abhyankar event->fvalue_prev[i] = event->fvalue[i]; 4882dc7a7e3SShri Abhyankar } 489e2cdd850SShri Abhyankar if (event->monitor && event->status == TSEVENT_PROCESSING) { 49038bf2713SShri 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); 4912dc7a7e3SShri Abhyankar } 49238bf2713SShri Abhyankar event->ptime_prev = t; 49338bf2713SShri Abhyankar } 49438bf2713SShri Abhyankar 49538bf2713SShri Abhyankar if (event->status == TSEVENT_PROCESSING) event->iterctr++; 4969e12be75SShri Abhyankar 4977dbe0728SLisandro Dalcin ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4987dbe0728SLisandro Dalcin ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr); 4992dc7a7e3SShri Abhyankar PetscFunctionReturn(0); 5002dc7a7e3SShri Abhyankar } 5012dc7a7e3SShri Abhyankar 5026427ac75SLisandro Dalcin PetscErrorCode TSAdjointEventHandler(TS ts) 503d0578d90SShri Abhyankar { 504d0578d90SShri Abhyankar PetscErrorCode ierr; 5056427ac75SLisandro Dalcin TSEvent event; 506d0578d90SShri Abhyankar PetscReal t; 507d0578d90SShri Abhyankar Vec U; 508d0578d90SShri Abhyankar PetscInt ctr; 509d0578d90SShri Abhyankar PetscBool forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */ 510d0578d90SShri Abhyankar 511d0578d90SShri Abhyankar PetscFunctionBegin; 5126427ac75SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5136427ac75SLisandro Dalcin if (!ts->event) PetscFunctionReturn(0); 5146427ac75SLisandro Dalcin event = ts->event; 515d0578d90SShri Abhyankar 516d0578d90SShri Abhyankar ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 517d0578d90SShri Abhyankar ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 518d0578d90SShri Abhyankar 519d0578d90SShri Abhyankar ctr = event->recorder.ctr-1; 520bcbf8bb3SShri Abhyankar if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) { 521d0578d90SShri Abhyankar /* Call the user postevent function */ 522d0578d90SShri Abhyankar if (event->postevent) { 5236427ac75SLisandro Dalcin ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 524d0578d90SShri Abhyankar event->recorder.ctr--; 525d0578d90SShri Abhyankar } 526d0578d90SShri Abhyankar } 527d0578d90SShri Abhyankar 528d0578d90SShri Abhyankar PetscBarrier((PetscObject)ts); 529d0578d90SShri Abhyankar PetscFunctionReturn(0); 530d0578d90SShri Abhyankar } 531