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; 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 events in the event recorder */ 147 ctr = event->recorder.ctr; 148 if (ctr == MAXEVENTRECORDERS) { 149 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Exceeded limit (=%d) for number of events recorded",MAXEVENTRECORDERS); 150 } 151 event->recorder.time[ctr] = t; 152 event->recorder.nevents[ctr] = nevents_zero; 153 for(i=0; i < nevents_zero; i++) event->recorder.eventidx[ctr][i] = events_zero[i]; 154 event->recorder.ctr++; 155 156 /* Reset the event residual functions as states might get changed by the postevent callback */ 157 ierr = (*event->monitor)(ts,t,U,event->fvalue_prev,event->monitorcontext);CHKERRQ(ierr); 158 event->ptime_prev = t; 159 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "TSEventMonitorDestroy" 165 PetscErrorCode TSEventMonitorDestroy(TSEvent *event) 166 { 167 PetscErrorCode ierr; 168 PetscInt i; 169 170 PetscFunctionBegin; 171 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 172 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 173 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 174 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 175 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 176 for(i=0; i < MAXEVENTRECORDERS; i++) { 177 ierr = PetscFree((*event)->recorder.eventidx[i]); 178 } 179 ierr = PetscFree(*event);CHKERRQ(ierr); 180 *event = NULL; 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "TSEventMonitor" 186 PetscErrorCode TSEventMonitor(TS ts) 187 { 188 PetscErrorCode ierr; 189 TSEvent event=ts->event; 190 PetscReal t; 191 Vec U; 192 PetscInt i; 193 PetscReal dt; 194 TSEventStatus status = event->status; 195 PetscInt rollback=0,in[2],out[2]; 196 PetscBool forwardsolve=PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */ 197 198 PetscFunctionBegin; 199 200 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 201 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 202 203 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 204 if (event->status == TSEVENT_RESET_NEXTSTEP) { 205 /* Take initial time step */ 206 dt = event->initial_timestep; 207 ts->time_step = dt; 208 event->status = TSEVENT_NONE; 209 } 210 211 if (event->status == TSEVENT_NONE) { 212 event->tstepend = t; 213 } 214 215 event->nevents_zero = 0; 216 217 ierr = (*event->monitor)(ts,t,U,event->fvalue,event->monitorcontext);CHKERRQ(ierr); 218 if (event->status != TSEVENT_NONE) { 219 for (i=0; i < event->nevents; i++) { 220 if (PetscAbsScalar(event->fvalue[i]) < event->tol) { 221 event->status = TSEVENT_ZERO; 222 event->events_zero[event->nevents_zero++] = i; 223 } 224 } 225 } 226 227 status = event->status; 228 ierr = MPI_Allreduce((PetscEnum*)&status,(PetscEnum*)&event->status,1,MPIU_ENUM,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 229 230 if (event->status == TSEVENT_ZERO) { 231 ierr = TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->monitorcontext);CHKERRQ(ierr); 232 dt = event->tstepend-t; 233 if(PetscAbsScalar(dt) < PETSC_SMALL) dt += event->initial_timestep; 234 ts->time_step = dt; 235 PetscFunctionReturn(0); 236 } 237 238 for (i = 0; i < event->nevents; i++) { 239 PetscInt fvalue_sign,fvalueprev_sign; 240 fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 241 fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 242 if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign)) { 243 switch (event->direction[i]) { 244 case -1: 245 if (fvalue_sign < 0) { 246 rollback = 1; 247 /* Compute linearly interpolated new time step */ 248 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 249 } 250 break; 251 case 1: 252 if (fvalue_sign > 0) { 253 rollback = 1; 254 /* Compute linearly interpolated new time step */ 255 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 256 } 257 break; 258 case 0: 259 rollback = 1; 260 /* Compute linearly interpolated new time step */ 261 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 262 break; 263 } 264 } 265 } 266 if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 267 268 in[0] = event->status; 269 in[1] = rollback; 270 ierr = MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 271 272 rollback = out[1]; 273 if (rollback) { 274 event->status = TSEVENT_LOCATED_INTERVAL; 275 } 276 277 if (event->status == TSEVENT_LOCATED_INTERVAL) { 278 ierr = TSRollBack(ts);CHKERRQ(ierr); 279 ts->steps--; 280 ts->total_steps--; 281 event->status = TSEVENT_PROCESSING; 282 } else { 283 for (i = 0; i < event->nevents; i++) { 284 event->fvalue_prev[i] = event->fvalue[i]; 285 } 286 event->ptime_prev = t; 287 if (event->status == TSEVENT_PROCESSING) { 288 dt = event->tstepend - event->ptime_prev; 289 } 290 } 291 ierr = MPI_Allreduce(&dt,&(ts->time_step),1,MPIU_REAL,MPI_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295