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