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