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