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 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 79 ierr = TSGetTimeStep(ts,&event->initial_timestep);CHKERRQ(ierr); 80 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 81 event->ptime_prev = t; 82 ierr = (*event->monitor)(ts,t,U,event->fvalue_prev,mectx);CHKERRQ(ierr); 83 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TS Event options","");CHKERRQ(ierr); 84 { 85 event->tol = 1.0e-6; 86 ierr = PetscOptionsReal("-ts_event_tol","","",event->tol,&event->tol,NULL);CHKERRQ(ierr); 87 } 88 ierr = PetscOptionsEnd();CHKERRQ(ierr); 89 ts->event = event; 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "TSPostEvent" 95 /* 96 TSPostEvent - Does post event processing by calling the user-defined postevent function 97 98 Logically Collective on TS 99 100 Input Parameters: 101 + ts - the TS context 102 . nevents_zero - number of local events whose event function is zero 103 . events_zero - indices of local events which have reached zero 104 . t - current time 105 . U - current solution 106 . forwardsolve - Flag to indicate whether TS is doing a forward solve (1) or adjoint solve (0) 107 - ctx - the context passed with eventmonitor 108 109 Level: intermediate 110 111 .keywords: TS, event, set, monitor 112 113 .seealso: TSSetEventMonitor(),TSEvent 114 */ 115 #undef __FUNCT__ 116 #define __FUNCT__ "TSPostEvent" 117 PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void *ctx) 118 { 119 PetscErrorCode ierr; 120 TSEvent event=ts->event; 121 PetscBool terminate=PETSC_FALSE; 122 PetscInt i; 123 PetscBool ts_terminate; 124 125 PetscFunctionBegin; 126 if (event->postevent) { 127 ierr = (*event->postevent)(ts,nevents_zero,events_zero,t,U,forwardsolve,ctx);CHKERRQ(ierr); 128 } 129 for(i = 0; i < nevents_zero;i++) { 130 terminate = (PetscBool)(terminate || event->terminate[events_zero[i]]); 131 } 132 ierr = MPI_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr); 133 if (terminate) { 134 ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr); 135 event->status = TSEVENT_NONE; 136 } else { 137 event->status = TSEVENT_RESET_NEXTSTEP; 138 } 139 /* Reset the event residual functions as states might get changed by the postevent callback */ 140 ierr = (*event->monitor)(ts,t,U,event->fvalue_prev,event->monitorcontext);CHKERRQ(ierr); 141 event->ptime_prev = t; 142 143 PetscFunctionReturn(0); 144 } 145 146 #undef __FUNCT__ 147 #define __FUNCT__ "TSEventMonitorDestroy" 148 PetscErrorCode TSEventMonitorDestroy(TSEvent *event) 149 { 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 154 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 155 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 156 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 157 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 158 ierr = PetscFree(*event);CHKERRQ(ierr); 159 *event = NULL; 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "TSEventMonitor" 165 PetscErrorCode TSEventMonitor(TS ts) 166 { 167 PetscErrorCode ierr; 168 TSEvent event=ts->event; 169 PetscReal t; 170 Vec U; 171 PetscInt i; 172 PetscReal dt; 173 TSEventStatus status = event->status; 174 PetscInt rollback=0,in[2],out[2]; 175 PetscBool forwardsolve=PETSC_TRUE; /* Flag indicating that TS is doing a forward solve */ 176 177 PetscFunctionBegin; 178 179 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 180 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 181 182 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 183 if (event->status == TSEVENT_RESET_NEXTSTEP) { 184 /* Take initial time step */ 185 dt = event->initial_timestep; 186 ts->time_step = dt; 187 event->status = TSEVENT_NONE; 188 } 189 190 if (event->status == TSEVENT_NONE) { 191 event->tstepend = t; 192 } 193 194 event->nevents_zero = 0; 195 196 ierr = (*event->monitor)(ts,t,U,event->fvalue,event->monitorcontext);CHKERRQ(ierr); 197 if (event->status != TSEVENT_NONE) { 198 for (i=0; i < event->nevents; i++) { 199 if (PetscAbsScalar(event->fvalue[i]) < event->tol) { 200 event->status = TSEVENT_ZERO; 201 event->events_zero[event->nevents_zero++] = i; 202 } 203 } 204 } 205 206 status = event->status; 207 ierr = MPI_Allreduce((PetscEnum*)&status,(PetscEnum*)&event->status,1,MPIU_ENUM,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 208 209 if (event->status == TSEVENT_ZERO) { 210 dt = event->tstepend-t; 211 ts->time_step = dt; 212 ierr = TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,forwardsolve,event->monitorcontext);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 for (i = 0; i < event->nevents; i++) { 217 PetscInt fvalue_sign,fvalueprev_sign; 218 fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 219 fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 220 if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign)) { 221 switch (event->direction[i]) { 222 case -1: 223 if (fvalue_sign < 0) { 224 rollback = 1; 225 /* Compute linearly interpolated new time step */ 226 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 227 } 228 break; 229 case 1: 230 if (fvalue_sign > 0) { 231 rollback = 1; 232 /* Compute linearly interpolated new time step */ 233 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 234 } 235 break; 236 case 0: 237 rollback = 1; 238 /* Compute linearly interpolated new time step */ 239 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 240 break; 241 } 242 } 243 } 244 if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 245 246 in[0] = event->status; 247 in[1] = rollback; 248 ierr = MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 249 250 rollback = out[1]; 251 if (rollback) { 252 event->status = TSEVENT_LOCATED_INTERVAL; 253 } 254 255 if (event->status == TSEVENT_LOCATED_INTERVAL) { 256 ierr = TSRollBack(ts);CHKERRQ(ierr); 257 event->status = TSEVENT_PROCESSING; 258 } else { 259 for (i = 0; i < event->nevents; i++) { 260 event->fvalue_prev[i] = event->fvalue[i]; 261 } 262 event->ptime_prev = t; 263 if (event->status == TSEVENT_PROCESSING) { 264 dt = event->tstepend - event->ptime_prev; 265 } 266 } 267 ierr = MPI_Allreduce(&dt,&(ts->time_step),1,MPIU_REAL,MPI_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 271