1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 /* 4 TSEventInitialize - Initializes TSEvent for TSSolve 5 */ 6 PetscErrorCode TSEventInitialize(TSEvent event,TS ts,PetscReal t,Vec U) 7 { 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 if (!event) PetscFunctionReturn(0); 12 PetscValidPointer(event,1); 13 PetscValidHeaderSpecific(ts,TS_CLASSID,2); 14 PetscValidHeaderSpecific(U,VEC_CLASSID,4); 15 event->ptime_prev = t; 16 event->iterctr = 0; 17 ierr = (*event->eventhandler)(ts,t,U,event->fvalue_prev,event->ctx);CHKERRQ(ierr); 18 PetscFunctionReturn(0); 19 } 20 21 PetscErrorCode TSEventDestroy(TSEvent *event) 22 { 23 PetscErrorCode ierr; 24 PetscInt i; 25 26 PetscFunctionBegin; 27 PetscValidPointer(event,1); 28 if (!*event) PetscFunctionReturn(0); 29 if (--(*event)->refct > 0) {*event = NULL; PetscFunctionReturn(0);} 30 31 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 32 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 33 ierr = PetscFree((*event)->fvalue_right);CHKERRQ(ierr); 34 ierr = PetscFree((*event)->zerocrossing);CHKERRQ(ierr); 35 ierr = PetscFree((*event)->side);CHKERRQ(ierr); 36 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 37 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 38 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 39 ierr = PetscFree((*event)->vtol);CHKERRQ(ierr); 40 41 for (i=0; i < (*event)->recsize; i++) { 42 ierr = PetscFree((*event)->recorder.eventidx[i]);CHKERRQ(ierr); 43 } 44 ierr = PetscFree((*event)->recorder.eventidx);CHKERRQ(ierr); 45 ierr = PetscFree((*event)->recorder.nevents);CHKERRQ(ierr); 46 ierr = PetscFree((*event)->recorder.stepnum);CHKERRQ(ierr); 47 ierr = PetscFree((*event)->recorder.time);CHKERRQ(ierr); 48 49 ierr = PetscViewerDestroy(&(*event)->monitor);CHKERRQ(ierr); 50 ierr = PetscFree(*event);CHKERRQ(ierr); 51 PetscFunctionReturn(0); 52 } 53 54 /*@ 55 TSSetPostEventIntervalStep - Set the time-step used immediately following the event interval 56 57 Logically Collective 58 59 Input Parameters: 60 + ts - time integration context 61 - dt - post event interval step 62 63 Options Database Keys: 64 . -ts_event_post_eventinterval_step <dt> time-step after event interval 65 66 Notes: 67 TSSetPostEventIntervalStep allows one to set a time-step that is used immediately following an event interval. 68 69 This function should be called from the postevent function set with TSSetEventHandler(). 70 71 The post event interval time-step should be selected based on the dynamics following the event. 72 If the dynamics are stiff, a conservative (small) step should be used. 73 If not, then a larger time-step can be used. 74 75 Level: Advanced 76 .seealso: TS, TSEvent, TSSetEventHandler() 77 @*/ 78 PetscErrorCode TSSetPostEventIntervalStep(TS ts,PetscReal dt) 79 { 80 PetscFunctionBegin; 81 ts->event->timestep_posteventinterval = dt; 82 PetscFunctionReturn(0); 83 } 84 85 /*@ 86 TSSetEventTolerances - Set tolerances for event zero crossings when using event handler 87 88 Logically Collective 89 90 Input Parameters: 91 + ts - time integration context 92 . tol - scalar tolerance, PETSC_DECIDE to leave current value 93 - vtol - array of tolerances or NULL, used in preference to tol if present 94 95 Options Database Keys: 96 . -ts_event_tol <tol> tolerance for event zero crossing 97 98 Notes: 99 Must call TSSetEventHandler() before setting the tolerances. 100 101 The size of vtol is equal to the number of events. 102 103 Level: beginner 104 105 .seealso: TS, TSEvent, TSSetEventHandler() 106 @*/ 107 PetscErrorCode TSSetEventTolerances(TS ts,PetscReal tol,PetscReal vtol[]) 108 { 109 TSEvent event; 110 PetscInt i; 111 112 PetscFunctionBegin; 113 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 114 if (vtol) PetscValidRealPointer(vtol,3); 115 if (!ts->event) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set the events first by calling TSSetEventHandler()"); 116 117 event = ts->event; 118 if (vtol) { 119 for (i=0; i < event->nevents; i++) event->vtol[i] = vtol[i]; 120 } else { 121 if (tol != PETSC_DECIDE || tol != PETSC_DEFAULT) { 122 for (i=0; i < event->nevents; i++) event->vtol[i] = tol; 123 } 124 } 125 PetscFunctionReturn(0); 126 } 127 128 /*@C 129 TSSetEventHandler - Sets a function used for detecting events 130 131 Logically Collective on TS 132 133 Input Parameters: 134 + ts - the TS context obtained from TSCreate() 135 . nevents - number of local events 136 . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction, 137 +1 => Zero crossing in positive direction, 0 => both ways (one for each event) 138 . terminate - flag to indicate whether time stepping should be terminated after 139 event is detected (one for each event) 140 . eventhandler - event monitoring routine 141 . postevent - [optional] post-event function 142 - ctx - [optional] user-defined context for private data for the 143 event monitor and post event routine (use NULL if no 144 context is desired) 145 146 Calling sequence of eventhandler: 147 PetscErrorCode PetscEventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx) 148 149 Input Parameters: 150 + ts - the TS context 151 . t - current time 152 . U - current iterate 153 - ctx - [optional] context passed with eventhandler 154 155 Output parameters: 156 . fvalue - function value of events at time t 157 158 Calling sequence of postevent: 159 PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx) 160 161 Input Parameters: 162 + ts - the TS context 163 . nevents_zero - number of local events whose event function is zero 164 . events_zero - indices of local events which have reached zero 165 . t - current time 166 . U - current solution 167 . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0) 168 - ctx - the context passed with eventhandler 169 170 Level: intermediate 171 172 .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason() 173 @*/ 174 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) 175 { 176 PetscErrorCode ierr; 177 TSAdapt adapt; 178 PetscReal hmin; 179 TSEvent event; 180 PetscInt i; 181 PetscBool flg; 182 #if defined PETSC_USE_REAL_SINGLE 183 PetscReal tol=1e-4; 184 #else 185 PetscReal tol=1e-6; 186 #endif 187 188 PetscFunctionBegin; 189 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 190 if (nevents) { 191 PetscValidIntPointer(direction,3); 192 PetscValidBoolPointer(terminate,4); 193 } 194 195 ierr = PetscNewLog(ts,&event);CHKERRQ(ierr); 196 ierr = PetscMalloc1(nevents,&event->fvalue);CHKERRQ(ierr); 197 ierr = PetscMalloc1(nevents,&event->fvalue_prev);CHKERRQ(ierr); 198 ierr = PetscMalloc1(nevents,&event->fvalue_right);CHKERRQ(ierr); 199 ierr = PetscMalloc1(nevents,&event->zerocrossing);CHKERRQ(ierr); 200 ierr = PetscMalloc1(nevents,&event->side);CHKERRQ(ierr); 201 ierr = PetscMalloc1(nevents,&event->direction);CHKERRQ(ierr); 202 ierr = PetscMalloc1(nevents,&event->terminate);CHKERRQ(ierr); 203 ierr = PetscMalloc1(nevents,&event->vtol);CHKERRQ(ierr); 204 for (i=0; i < nevents; i++) { 205 event->direction[i] = direction[i]; 206 event->terminate[i] = terminate[i]; 207 event->zerocrossing[i] = PETSC_FALSE; 208 event->side[i] = 0; 209 } 210 ierr = PetscMalloc1(nevents,&event->events_zero);CHKERRQ(ierr); 211 event->nevents = nevents; 212 event->eventhandler = eventhandler; 213 event->postevent = postevent; 214 event->ctx = ctx; 215 event->timestep_posteventinterval = ts->time_step; 216 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 217 ierr = TSAdaptGetStepLimits(adapt,&hmin,NULL);CHKERRQ(ierr); 218 event->timestep_min = hmin; 219 220 event->recsize = 8; /* Initial size of the recorder */ 221 ierr = PetscOptionsBegin(((PetscObject)ts)->comm,((PetscObject)ts)->prefix,"TS Event options","TS");CHKERRQ(ierr); 222 { 223 ierr = PetscOptionsReal("-ts_event_tol","Scalar event tolerance for zero crossing check","TSSetEventTolerances",tol,&tol,NULL);CHKERRQ(ierr); 224 ierr = PetscOptionsName("-ts_event_monitor","Print choices made by event handler","",&flg);CHKERRQ(ierr); 225 ierr = PetscOptionsInt("-ts_event_recorder_initial_size","Initial size of event recorder","",event->recsize,&event->recsize,NULL);CHKERRQ(ierr); 226 ierr = PetscOptionsReal("-ts_event_post_eventinterval_step","Time step after event interval","",event->timestep_posteventinterval,&event->timestep_posteventinterval,NULL);CHKERRQ(ierr); 227 ierr = PetscOptionsReal("-ts_event_post_event_step","Time step after event","",event->timestep_postevent,&event->timestep_postevent,NULL);CHKERRQ(ierr); 228 ierr = PetscOptionsReal("-ts_event_dt_min","Minimum time step considered for TSEvent","",event->timestep_min,&event->timestep_min,NULL);CHKERRQ(ierr); 229 } 230 ierr = PetscOptionsEnd();CHKERRQ(ierr); 231 232 ierr = PetscMalloc1(event->recsize,&event->recorder.time);CHKERRQ(ierr); 233 ierr = PetscMalloc1(event->recsize,&event->recorder.stepnum);CHKERRQ(ierr); 234 ierr = PetscMalloc1(event->recsize,&event->recorder.nevents);CHKERRQ(ierr); 235 ierr = PetscMalloc1(event->recsize,&event->recorder.eventidx);CHKERRQ(ierr); 236 for (i=0; i < event->recsize; i++) { 237 ierr = PetscMalloc1(event->nevents,&event->recorder.eventidx[i]);CHKERRQ(ierr); 238 } 239 /* Initialize the event recorder */ 240 event->recorder.ctr = 0; 241 242 for (i=0; i < event->nevents; i++) event->vtol[i] = tol; 243 if (flg) {ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,"stdout",&event->monitor);CHKERRQ(ierr);} 244 245 ierr = TSEventDestroy(&ts->event);CHKERRQ(ierr); 246 ts->event = event; 247 ts->event->refct = 1; 248 PetscFunctionReturn(0); 249 } 250 251 /* 252 TSEventRecorderResize - Resizes (2X) the event recorder arrays whenever the recording limit (event->recsize) 253 is reached. 254 */ 255 static PetscErrorCode TSEventRecorderResize(TSEvent event) 256 { 257 PetscErrorCode ierr; 258 PetscReal *time; 259 PetscInt *stepnum; 260 PetscInt *nevents; 261 PetscInt **eventidx; 262 PetscInt i,fact=2; 263 264 PetscFunctionBegin; 265 266 /* Create large arrays */ 267 ierr = PetscMalloc1(fact*event->recsize,&time);CHKERRQ(ierr); 268 ierr = PetscMalloc1(fact*event->recsize,&stepnum);CHKERRQ(ierr); 269 ierr = PetscMalloc1(fact*event->recsize,&nevents);CHKERRQ(ierr); 270 ierr = PetscMalloc1(fact*event->recsize,&eventidx);CHKERRQ(ierr); 271 for (i=0; i < fact*event->recsize; i++) { 272 ierr = PetscMalloc1(event->nevents,&eventidx[i]);CHKERRQ(ierr); 273 } 274 275 /* Copy over data */ 276 ierr = PetscArraycpy(time,event->recorder.time,event->recsize);CHKERRQ(ierr); 277 ierr = PetscArraycpy(stepnum,event->recorder.stepnum,event->recsize);CHKERRQ(ierr); 278 ierr = PetscArraycpy(nevents,event->recorder.nevents,event->recsize);CHKERRQ(ierr); 279 for (i=0; i < event->recsize; i++) { 280 ierr = PetscArraycpy(eventidx[i],event->recorder.eventidx[i],event->recorder.nevents[i]);CHKERRQ(ierr); 281 } 282 283 /* Destroy old arrays */ 284 for (i=0; i < event->recsize; i++) { 285 ierr = PetscFree(event->recorder.eventidx[i]);CHKERRQ(ierr); 286 } 287 ierr = PetscFree(event->recorder.eventidx);CHKERRQ(ierr); 288 ierr = PetscFree(event->recorder.nevents);CHKERRQ(ierr); 289 ierr = PetscFree(event->recorder.stepnum);CHKERRQ(ierr); 290 ierr = PetscFree(event->recorder.time);CHKERRQ(ierr); 291 292 /* Set pointers */ 293 event->recorder.time = time; 294 event->recorder.stepnum = stepnum; 295 event->recorder.nevents = nevents; 296 event->recorder.eventidx = eventidx; 297 298 /* Double size */ 299 event->recsize *= fact; 300 301 PetscFunctionReturn(0); 302 } 303 304 /* 305 Helper routine to handle user postevents and recording 306 */ 307 static PetscErrorCode TSPostEvent(TS ts,PetscReal t,Vec U) 308 { 309 PetscErrorCode ierr; 310 TSEvent event = ts->event; 311 PetscBool terminate = PETSC_FALSE; 312 PetscBool restart = PETSC_FALSE; 313 PetscInt i,ctr,stepnum; 314 PetscBool inflag[2],outflag[2]; 315 PetscBool forwardsolve = PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */ 316 317 PetscFunctionBegin; 318 if (event->postevent) { 319 PetscObjectState state_prev,state_post; 320 ierr = PetscObjectStateGet((PetscObject)U,&state_prev);CHKERRQ(ierr); 321 ierr = (*event->postevent)(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 322 ierr = PetscObjectStateGet((PetscObject)U,&state_post);CHKERRQ(ierr); 323 if (state_prev != state_post) restart = PETSC_TRUE; 324 } 325 326 /* Handle termination events and step restart */ 327 for (i=0; i<event->nevents_zero; i++) if (event->terminate[event->events_zero[i]]) terminate = PETSC_TRUE; 328 inflag[0] = restart; inflag[1] = terminate; 329 ierr = MPIU_Allreduce(inflag,outflag,2,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRMPI(ierr); 330 restart = outflag[0]; terminate = outflag[1]; 331 if (restart) {ierr = TSRestartStep(ts);CHKERRQ(ierr);} 332 if (terminate) {ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr);} 333 event->status = terminate ? TSEVENT_NONE : TSEVENT_RESET_NEXTSTEP; 334 335 /* Reset event residual functions as states might get changed by the postevent callback */ 336 if (event->postevent) { 337 ierr = VecLockReadPush(U);CHKERRQ(ierr); 338 ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr); 339 ierr = VecLockReadPop(U);CHKERRQ(ierr); 340 } 341 342 /* Cache current time and event residual functions */ 343 event->ptime_prev = t; 344 for (i=0; i<event->nevents; i++) 345 event->fvalue_prev[i] = event->fvalue[i]; 346 347 /* Record the event in the event recorder */ 348 ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr); 349 ctr = event->recorder.ctr; 350 if (ctr == event->recsize) { 351 ierr = TSEventRecorderResize(event);CHKERRQ(ierr); 352 } 353 event->recorder.time[ctr] = t; 354 event->recorder.stepnum[ctr] = stepnum; 355 event->recorder.nevents[ctr] = event->nevents_zero; 356 for (i=0; i<event->nevents_zero; i++) event->recorder.eventidx[ctr][i] = event->events_zero[i]; 357 event->recorder.ctr++; 358 PetscFunctionReturn(0); 359 } 360 361 /* Uses Anderson-Bjorck variant of regula falsi method */ 362 PETSC_STATIC_INLINE PetscReal TSEventComputeStepSize(PetscReal tleft,PetscReal t,PetscReal tright,PetscScalar fleft,PetscScalar f,PetscScalar fright,PetscInt side,PetscReal dt) 363 { 364 PetscReal new_dt, scal = 1.0; 365 if (PetscRealPart(fleft)*PetscRealPart(f) < 0) { 366 if (side == 1) { 367 scal = (PetscRealPart(fright) - PetscRealPart(f))/PetscRealPart(fright); 368 if (scal < PETSC_SMALL) scal = 0.5; 369 } 370 new_dt = (scal*PetscRealPart(fleft)*t - PetscRealPart(f)*tleft)/(scal*PetscRealPart(fleft) - PetscRealPart(f)) - tleft; 371 } else { 372 if (side == -1) { 373 scal = (PetscRealPart(fleft) - PetscRealPart(f))/PetscRealPart(fleft); 374 if (scal < PETSC_SMALL) scal = 0.5; 375 } 376 new_dt = (PetscRealPart(f)*tright - scal*PetscRealPart(fright)*t)/(PetscRealPart(f) - scal*PetscRealPart(fright)) - t; 377 } 378 return PetscMin(dt,new_dt); 379 } 380 381 static PetscErrorCode TSEventDetection(TS ts) 382 { 383 PetscErrorCode ierr; 384 TSEvent event = ts->event; 385 PetscReal t; 386 PetscInt i; 387 PetscInt fvalue_sign,fvalueprev_sign; 388 PetscInt in,out; 389 390 PetscFunctionBegin; 391 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 392 for (i=0; i < event->nevents; i++) { 393 if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i]) { 394 if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 395 event->status = TSEVENT_LOCATED_INTERVAL; 396 if (event->monitor) { 397 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected due to zero value (tol=%g) [%g - %g]\n",event->iterctr,i,(double)event->vtol[i],(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 398 } 399 continue; 400 } 401 if (PetscAbsScalar(event->fvalue_prev[i]) < event->vtol[i]) continue; /* avoid duplicative detection if the previous endpoint is an event location */ 402 fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 403 fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 404 if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign)) { 405 if (!event->iterctr) event->zerocrossing[i] = PETSC_TRUE; 406 event->status = TSEVENT_LOCATED_INTERVAL; 407 if (event->monitor) { 408 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D interval detected due to sign change [%g - %g]\n",event->iterctr,i,(double)event->ptime_prev,(double)t);CHKERRQ(ierr); 409 } 410 } 411 } 412 in = event->status; 413 ierr = MPIU_Allreduce(&in,&out,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRMPI(ierr); 414 event->status = (TSEventStatus)out; 415 PetscFunctionReturn(0); 416 } 417 418 static PetscErrorCode TSEventLocation(TS ts,PetscReal *dt) 419 { 420 PetscErrorCode ierr; 421 TSEvent event = ts->event; 422 PetscInt i; 423 PetscReal t; 424 PetscInt fvalue_sign,fvalueprev_sign; 425 PetscInt rollback=0,in[2],out[2]; 426 427 PetscFunctionBegin; 428 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 429 event->nevents_zero = 0; 430 for (i=0; i < event->nevents; i++) { 431 if (event->zerocrossing[i]) { 432 if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i] || *dt < event->timestep_min || PetscAbsReal((*dt)/((event->ptime_right-event->ptime_prev)/2)) < event->vtol[i]) { /* stopping criteria */ 433 event->status = TSEVENT_ZERO; 434 event->fvalue_right[i] = event->fvalue[i]; 435 continue; 436 } 437 /* Compute new time step */ 438 *dt = TSEventComputeStepSize(event->ptime_prev,t,event->ptime_right,event->fvalue_prev[i],event->fvalue[i],event->fvalue_right[i],event->side[i],*dt); 439 fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 440 fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 441 switch (event->direction[i]) { 442 case -1: 443 if (fvalue_sign < 0) { 444 rollback = 1; 445 event->fvalue_right[i] = event->fvalue[i]; 446 event->side[i] = 1; 447 } 448 break; 449 case 1: 450 if (fvalue_sign > 0) { 451 rollback = 1; 452 event->fvalue_right[i] = event->fvalue[i]; 453 event->side[i] = 1; 454 } 455 break; 456 case 0: 457 if (fvalue_sign != fvalueprev_sign) { /* trigger rollback only when there is a sign change */ 458 rollback = 1; 459 event->fvalue_right[i] = event->fvalue[i]; 460 event->side[i] = 1; 461 } 462 break; 463 } 464 if (event->status == TSEVENT_PROCESSING) event->side[i] = -1; 465 } 466 } 467 in[0] = event->status; in[1] = rollback; 468 ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRMPI(ierr); 469 event->status = (TSEventStatus)out[0]; rollback = out[1]; 470 /* If rollback is true, the status will be overwritten so that an event at the endtime of current time step will be postponed to guarantee corret order */ 471 if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 472 if (event->status == TSEVENT_ZERO) { 473 for (i=0; i < event->nevents; i++) { 474 if (event->zerocrossing[i]) { 475 if (PetscAbsScalar(event->fvalue[i]) < event->vtol[i] || *dt < event->timestep_min || PetscAbsReal((*dt)/((event->ptime_right-event->ptime_prev)/2)) < event->vtol[i]) { /* stopping criteria */ 476 event->events_zero[event->nevents_zero++] = i; 477 if (event->monitor) { 478 ierr = PetscViewerASCIIPrintf(event->monitor,"TSEvent: iter %D - Event %D zero crossing located at time %g\n",event->iterctr,i,(double)t);CHKERRQ(ierr); 479 } 480 event->zerocrossing[i] = PETSC_FALSE; 481 } 482 } 483 event->side[i] = 0; 484 } 485 } 486 PetscFunctionReturn(0); 487 } 488 489 PetscErrorCode TSEventHandler(TS ts) 490 { 491 PetscErrorCode ierr; 492 TSEvent event; 493 PetscReal t; 494 Vec U; 495 PetscInt i; 496 PetscReal dt,dt_min,dt_reset = 0.0; 497 498 PetscFunctionBegin; 499 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 500 if (!ts->event) PetscFunctionReturn(0); 501 event = ts->event; 502 503 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 504 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 505 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 506 507 if (event->status == TSEVENT_NONE) { 508 event->timestep_prev = dt; 509 event->ptime_end = t; 510 } 511 if (event->status == TSEVENT_RESET_NEXTSTEP) { 512 /* user has specified a PostEventInterval dt */ 513 dt = event->timestep_posteventinterval; 514 if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { 515 PetscReal maxdt = ts->max_time-t; 516 dt = dt > maxdt ? maxdt : (PetscIsCloseAtTol(dt,maxdt,10*PETSC_MACHINE_EPSILON,0) ? maxdt : dt); 517 } 518 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 519 event->status = TSEVENT_NONE; 520 } 521 522 ierr = VecLockReadPush(U);CHKERRQ(ierr); 523 ierr = (*event->eventhandler)(ts,t,U,event->fvalue,event->ctx);CHKERRQ(ierr); 524 ierr = VecLockReadPop(U);CHKERRQ(ierr); 525 526 /* Detect the events */ 527 ierr = TSEventDetection(ts);CHKERRQ(ierr); 528 529 /* Locate the events */ 530 if (event->status == TSEVENT_LOCATED_INTERVAL || event->status == TSEVENT_PROCESSING) { 531 /* Approach the zero crosing by setting a new step size */ 532 ierr = TSEventLocation(ts,&dt);CHKERRQ(ierr); 533 /* Roll back when new events are detected */ 534 if (event->status == TSEVENT_LOCATED_INTERVAL) { 535 ierr = TSRollBack(ts);CHKERRQ(ierr); 536 ierr = TSSetConvergedReason(ts,TS_CONVERGED_ITERATING);CHKERRQ(ierr); 537 event->iterctr++; 538 } 539 ierr = MPIU_Allreduce(&dt,&dt_min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRMPI(ierr); 540 if (dt_reset > 0.0 && dt_reset < dt_min) dt_min = dt_reset; 541 ierr = TSSetTimeStep(ts,dt_min);CHKERRQ(ierr); 542 /* Found the zero crossing */ 543 if (event->status == TSEVENT_ZERO) { 544 ierr = TSPostEvent(ts,t,U);CHKERRQ(ierr); 545 546 dt = event->ptime_end - t; 547 if (PetscAbsReal(dt) < PETSC_SMALL) { /* we hit the event, continue with the candidate time step */ 548 dt = event->timestep_prev; 549 event->status = TSEVENT_NONE; 550 } 551 if (event->timestep_postevent) { /* user has specified a PostEvent dt*/ 552 dt = event->timestep_postevent; 553 } 554 if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { 555 PetscReal maxdt = ts->max_time-t; 556 dt = dt > maxdt ? maxdt : (PetscIsCloseAtTol(dt,maxdt,10*PETSC_MACHINE_EPSILON,0) ? maxdt : dt); 557 } 558 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 559 event->iterctr = 0; 560 } 561 /* Have not found the zero crosing yet */ 562 if (event->status == TSEVENT_PROCESSING) { 563 if (event->monitor) { 564 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); 565 } 566 event->iterctr++; 567 } 568 } 569 if (event->status == TSEVENT_LOCATED_INTERVAL) { /* The step has been rolled back */ 570 event->status = TSEVENT_PROCESSING; 571 event->ptime_right = t; 572 } else { 573 for (i=0; i < event->nevents; i++) event->fvalue_prev[i] = event->fvalue[i]; 574 event->ptime_prev = t; 575 } 576 PetscFunctionReturn(0); 577 } 578 579 PetscErrorCode TSAdjointEventHandler(TS ts) 580 { 581 PetscErrorCode ierr; 582 TSEvent event; 583 PetscReal t; 584 Vec U; 585 PetscInt ctr; 586 PetscBool forwardsolve=PETSC_FALSE; /* Flag indicating that TS is doing an adjoint solve */ 587 588 PetscFunctionBegin; 589 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 590 if (!ts->event) PetscFunctionReturn(0); 591 event = ts->event; 592 593 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 594 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 595 596 ctr = event->recorder.ctr-1; 597 if (ctr >= 0 && PetscAbsReal(t - event->recorder.time[ctr]) < PETSC_SMALL) { 598 /* Call the user postevent function */ 599 if (event->postevent) { 600 ierr = (*event->postevent)(ts,event->recorder.nevents[ctr],event->recorder.eventidx[ctr],t,U,forwardsolve,event->ctx);CHKERRQ(ierr); 601 event->recorder.ctr--; 602 } 603 } 604 605 PetscBarrier((PetscObject)ts); 606 PetscFunctionReturn(0); 607 } 608 609 /*@ 610 TSGetNumEvents - Get the numbers of events set 611 612 Logically Collective 613 614 Input Parameter: 615 . ts - the TS context 616 617 Output Parameter: 618 . nevents - number of events 619 620 Level: intermediate 621 622 .seealso: TSSetEventHandler() 623 624 @*/ 625 PetscErrorCode TSGetNumEvents(TS ts,PetscInt * nevents) 626 { 627 PetscFunctionBegin; 628 *nevents = ts->event->nevents; 629 PetscFunctionReturn(0); 630 } 631