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