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