1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "TSEventInitialize" 5 /* 6 TSEventInitialize - Initializes TSEvent for TSSolve 7 */ 8 PetscErrorCode TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U) 9 { 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 if (!event) PetscFunctionReturn(0); 14 PetscValidPointer(event,1); 15 PetscValidHeaderSpecific(ts,TS_CLASSID,2); 16 PetscValidHeaderSpecific(U,VEC_CLASSID,4); 17 18 ierr = TSGetTimeStep(ts,&event->initial_timestep);CHKERRQ(ierr); 19 event->ptime_prev = t; 20 ierr = (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);CHKERRQ(ierr); 21 /* Initialize the event recorder */ 22 event->recorder.ctr = 0; 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "TSSetEventTolerances" 28 /*@ 29 TSSetEventTolerances - Set tolerances for event zero crossings when using event handler 30 31 Logically Collective 32 33 Input Arguments: 34 + ts - time integration context 35 . tol - scalar tolerance, PETSC_DECIDE to leave current value 36 - vtol - array of tolerances or NULL, used in preference to tol if present 37 38 - -ts_event_tol <tol> tolerance for event zero crossing 39 40 Notes: 41 Must call TSSetEventHandler() before setting the tolerances. 42 43 The size of vtol is equal to the number of events. 44 45 Level: beginner 46 47 .seealso: TS, TSEvent, TSSetEventHandler() 48 @*/ 49 PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[]) 50 { 51 TSEvent event; 52 PetscInt i; 53 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 56 if (vtol) PetscValidRealPointer(vtol,3); 57 if(!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()"); 58 59 event = ts->event; 60 if (vtol) { 61 for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i]; 62 } else { 63 if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) { 64 for (i=0; i < event->nevents; i++) event->vtol[i] = tol; 65 } 66 } 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "TSSetEventHandler" 72 /*@C 73 TSSetEventHandler - Sets a monitoring function used for detecting events 74 75 Logically Collective on TS 76 77 Input Parameters: 78 + ts - the TS context obtained from TSCreate() 79 . nevents - number of local events 80 . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction, 81 +1 => Zero crossing in positive direction, 0 => both ways (one for each event) 82 . terminate - flag to indicate whether time stepping should be terminated after 83 event is detected (one for each event) 84 . eventhandler - event monitoring routine 85 . postevent - [optional] post-event function 86 - mectx - [optional] user-defined context for private data for the 87 event monitor and post event routine (use NULL if no 88 context is desired) 89 90 Calling sequence of eventhandler: 91 PetscErrorCode EventHandler(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void* mectx) 92 93 Input Parameters: 94 + ts - the TS context 95 . t - current time 96 . U - current iterate 97 - ctx - [optional] context passed with eventhandler 98 99 Output parameters: 100 . fvalue - function value of events at time t 101 102 Calling sequence of postevent: 103 PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero, PetscInt events_zero, PetscReal t,Vec U,PetscBool forwardsolve,void* ctx) 104 105 Input Parameters: 106 + ts - the TS context 107 . nevents_zero - number of local events whose event function is zero 108 . events_zero - indices of local events which have reached zero 109 . t - current time 110 . U - current solution 111 . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0) 112 - ctx - the context passed with eventhandler 113 114 Level: intermediate 115 116 .keywords: TS, event, set, monitor 117 118 .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason() 119 @*/ 120 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) 121 { 122 PetscErrorCode ierr; 123 TSEvent event; 124 PetscInt i; 125 PetscBool flg; 126 PetscReal tol=1e-6; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 130 PetscValidIntPointer(direction,2); 131 PetscValidIntPointer(terminate,3); 132 133 ierr = PetscNewLog(ts,&event);CHKERRQ(ierr); 134 ierr = PetscMalloc1(nevents,&event->fvalue);CHKERRQ(ierr); 135 ierr = PetscMalloc1(nevents,&event->fvalue_prev);CHKERRQ(ierr); 136 ierr = PetscMalloc1(nevents,&event->direction);CHKERRQ(ierr); 137 ierr = PetscMalloc1(nevents,&event->terminate);CHKERRQ(ierr); 138 ierr = PetscMalloc1(nevents,&event->vtol);CHKERRQ(ierr); 139 for (i=0; i < nevents; i++) { 140 event->direction[i] = direction[i]; 141 event->terminate[i] = terminate[i]; 142 } 143 ierr = PetscMalloc1(nevents,&event->events_zero);CHKERRQ(ierr); 144 event->eventhandler = eventhandler; 145 event->postevent = postevent; 146 event->ctx = ctx; 147 event->nevents = nevents; 148 149 for (i=0; i < MAXEVENTRECORDERS; i++) { 150 ierr = PetscMalloc1(nevents,&event->recorder.eventidx[i]);CHKERRQ(ierr); 151 } 152 153 ierr = PetscOptionsBegin(((PetscObject)ts)->comm,"","TS Event options","");CHKERRQ(ierr); 154 { 155 ierr = PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","",tol,&tol,NULL);CHKERRQ(ierr); 156 ierr = PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);CHKERRQ(ierr); 157 } 158 PetscOptionsEnd(); 159 for (i=0; i < event->nevents; i++) event->vtol[i] = tol; 160 if (flg) {ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);CHKERRQ(ierr);} 161 162 ierr = TSEventDestroy(&ts->event);CHKERRQ(ierr); 163 ts->event = event; 164 PetscFunctionReturn(0); 165 } 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "TSPostEvent" 169 /* 170 TSPostEvent - Does post event processing by calling the user-defined postevent function 171 172 Logically Collective on TS 173 174 Input Parameters: 175 + ts - the TS context 176 . nevents_zero - number of local events whose event function is zero 177 . events_zero - indices of local events which have reached zero 178 . t - current time 179 . U - current solution 180 - forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0) 181 182 Level: intermediate 183 184 .keywords: TS, event, set, monitor 185 186 .seealso: TSSetEventHandler(), TSEvent 187 */ 188 #undef __FUNCT__ 189 #define __FUNCT__ "TSPostEvent" 190 PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve) 191 { 192 PetscErrorCode ierr; 193 TSEvent event=ts->event; 194 PetscBool terminate=PETSC_FALSE; 195 PetscInt i,ctr,stepnum; 196 PetscBool ts_terminate; 197 198 PetscFunctionBegin; 199 if (event->postevent) { 200 ierr = (*event->postevent)(ts,nevents_zero,events_zero,t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 201 } 202 for(i = 0; i < nevents_zero;i++) { 203 terminate = (PetscBool)(terminate || event->terminate[events_zero[i]]); 204 } 205 ierr = MPIU_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr); 206 if (ts_terminate) { 207 ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr); 208 event->status = TSEVENT_NONE; 209 } else { 210 event->status = TSEVENT_RESET_NEXTSTEP; 211 } 212 213 /* Record the event in the event recorder */ 214 ierr = TSGetTimeStepNumber(ts,&stepnum);CHKERRQ(ierr); 215 ctr = event->recorder.ctr; 216 if (ctr == MAXEVENTRECORDERS) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Exceeded limit (=%d) for number of events recorded",MAXEVENTRECORDERS); 217 event->recorder.time[ctr] = t; 218 event->recorder.stepnum[ctr] = stepnum; 219 event->recorder.nevents[ctr] = nevents_zero; 220 for(i=0; i < nevents_zero; i++) event->recorder.eventidx[ctr][i] = events_zero[i]; 221 event->recorder.ctr++; 222 223 /* Reset the event residual functions as states might get changed by the postevent callback */ 224 ierr = (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);CHKERRQ(ierr); 225 event->ptime_prev = t; 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "TSEventDestroy" 231 PetscErrorCode TSEventDestroy(TSEvent *event) 232 { 233 PetscErrorCode ierr; 234 PetscInt i; 235 236 PetscFunctionBegin; 237 PetscValidPointer(event,1); 238 if (!*event) PetscFunctionReturn(0); 239 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 240 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 241 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 242 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 243 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 244 ierr = PetscFree((*event)->vtol);CHKERRQ(ierr); 245 for(i=0; i < MAXEVENTRECORDERS; i++) { 246 ierr = PetscFree((*event)->recorder.eventidx[i]);CHKERRQ(ierr); 247 } 248 ierr = PetscViewerDestroy(&(*event)->monitor);CHKERRQ(ierr); 249 ierr = PetscFree(*event);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "TSEventHandler" 255 PetscErrorCode TSEventHandler(TS ts) 256 { 257 PetscErrorCode ierr; 258 TSEvent event; 259 PetscReal t; 260 Vec U; 261 PetscInt i; 262 PetscReal dt; 263 PetscInt rollback=0,in[2],out[2]; 264 PetscBool forwardsolve=PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */ 265 PetscInt fvalue_sign,fvalueprev_sign; 266 267 PetscFunctionBegin; 268 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 269 if (!ts->event) PetscFunctionReturn(0); 270 event = ts->event; 271 272 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 273 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 274 275 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 276 if (event->status == TSEVENT_RESET_NEXTSTEP) { 277 /* Take initial time step */ 278 dt = event->initial_timestep; 279 ts->time_step = dt; 280 event->status = TSEVENT_NONE; 281 } 282 283 if (event->status == TSEVENT_NONE) { 284 event->tstepend = t; 285 } 286 287 event->nevents_zero = 0; 288 289 ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr); 290 291 for (i = 0; i < event->nevents; i++) { 292 if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 293 event->status = TSEVENT_ZERO; 294 continue; 295 } 296 fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 297 fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 298 if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign) && (PetscAbsScalar(event->fvalue_prev[i]) > event->vtol[i])) { 299 switch (event->direction[i]) { 300 case -1: 301 if (fvalue_sign < 0) { 302 rollback = 1; 303 /* Compute linearly interpolated new time step */ 304 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 305 if(event->monitor) { 306 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent : Event %D interval located [%g - %g]\n",i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 307 } 308 } 309 break; 310 case 1: 311 if (fvalue_sign > 0) { 312 rollback = 1; 313 /* Compute linearly interpolated new time step */ 314 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 315 if(event->monitor) { 316 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent : Event %D interval located [%g - %g]\n",i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 317 } 318 } 319 break; 320 case 0: 321 rollback = 1; 322 /* Compute linearly interpolated new time step */ 323 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 324 if(event->monitor) { 325 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent : Event %D interval located [%g - %g]\n",i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 326 } 327 break; 328 } 329 } 330 } 331 if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 332 333 in[0] = event->status; 334 in[1] = rollback; 335 ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 336 337 event->status = (TSEventStatus)out[0]; 338 rollback = out[1]; 339 if (rollback) { 340 event->status = TSEVENT_LOCATED_INTERVAL; 341 } 342 343 if (event->status == TSEVENT_ZERO) { 344 for(i=0; i < event->nevents; i++) { 345 if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 346 event->events_zero[event->nevents_zero++] = i; 347 if(event->monitor) { 348 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent : Event %D zero crossing at time %g\n",i,(double)t);CHKERRQ(ierr); 349 } 350 } 351 } 352 353 ierr = TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve);CHKERRQ(ierr); 354 355 dt = event->tstepend-t; 356 if(PetscAbsReal(dt) < PETSC_SMALL) dt += event->initial_timestep; 357 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 358 PetscFunctionReturn(0); 359 } 360 361 if (event->status == TSEVENT_LOCATED_INTERVAL) { 362 ierr = TSRollBack(ts);CHKERRQ(ierr); 363 ts->steps--; 364 ts->total_steps--; 365 ts->reason = TS_CONVERGED_ITERATING; 366 event->status = TSEVENT_PROCESSING; 367 } else { 368 for (i = 0; i < event->nevents; i++) { 369 event->fvalue_prev[i] = event->fvalue[i]; 370 } 371 event->ptime_prev = t; 372 if (event->status == TSEVENT_PROCESSING) { 373 dt = event->tstepend - event->ptime_prev; 374 } 375 } 376 377 ierr = MPIU_Allreduce(&dt,&(ts->time_step),1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 378 PetscFunctionReturn(0); 379 } 380 381 #undef __FUNCT__ 382 #define __FUNCT__ "TSAdjointEventHandler" 383 PetscErrorCode TSAdjointEventHandler(TS ts) 384 { 385 PetscErrorCode ierr; 386 TSEvent event; 387 PetscReal t; 388 Vec U; 389 PetscInt ctr; 390 PetscBool forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */ 391 392 PetscFunctionBegin; 393 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 394 if (!ts->event) PetscFunctionReturn(0); 395 event = ts->event; 396 397 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 398 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 399 400 ctr = event->recorder.ctr-1; 401 if(ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) { 402 /* Call the user postevent function */ 403 if(event->postevent) { 404 ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 405 event->recorder.ctr--; 406 } 407 } 408 409 PetscBarrier((PetscObject)ts); 410 PetscFunctionReturn(0); 411 } 412