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