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