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