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,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 - ctx - the context passed with eventmonitor 46 47 Level: intermediate 48 49 .keywords: TS, event, set, monitor 50 51 .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason() 52 @*/ 53 PetscErrorCode TSSetEventMonitor(TS ts,PetscInt nevents,PetscInt *direction,PetscBool *terminate,PetscErrorCode (*eventmonitor)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,void*),void *mectx) 54 { 55 PetscErrorCode ierr; 56 PetscReal t; 57 Vec U; 58 TSEvent event; 59 PetscInt i; 60 61 PetscFunctionBegin; 62 ierr = PetscNew(&event);CHKERRQ(ierr); 63 ierr = PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue);CHKERRQ(ierr); 64 ierr = PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue_prev);CHKERRQ(ierr); 65 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->direction);CHKERRQ(ierr); 66 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->terminate);CHKERRQ(ierr); 67 for (i=0; i < nevents; i++) { 68 event->direction[i] = direction[i]; 69 event->terminate[i] = terminate[i]; 70 } 71 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->events_zero);CHKERRQ(ierr); 72 event->monitor = eventmonitor; 73 event->postevent = postevent; 74 event->monitorcontext = (void*)mectx; 75 event->nevents = nevents; 76 77 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 78 ierr = TSGetTimeStep(ts,&event->initial_timestep);CHKERRQ(ierr); 79 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 80 event->ptime_prev = t; 81 ierr = (*event->monitor)(ts,t,U,event->fvalue_prev,mectx);CHKERRQ(ierr); 82 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TS Event options","");CHKERRQ(ierr); 83 { 84 event->tol = 1.0e-6; 85 ierr = PetscOptionsReal("-ts_event_tol","","",event->tol,&event->tol,NULL);CHKERRQ(ierr); 86 } 87 ierr = PetscOptionsEnd();CHKERRQ(ierr); 88 ts->event = event; 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "TSPostEvent" 94 PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,void *ctx) 95 { 96 PetscErrorCode ierr; 97 TSEvent event=ts->event; 98 PetscBool terminate=PETSC_FALSE; 99 PetscInt i; 100 PetscBool ts_terminate; 101 102 PetscFunctionBegin; 103 if (event->postevent) { 104 ierr = (*event->postevent)(ts,nevents_zero,events_zero,t,U,ctx);CHKERRQ(ierr); 105 } 106 for(i = 0; i < nevents_zero;i++) { 107 terminate = (PetscBool)(terminate || event->terminate[events_zero[i]]); 108 } 109 ierr = MPI_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);CHKERRQ(ierr); 110 if (terminate) { 111 ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr); 112 event->status = TSEVENT_NONE; 113 } else { 114 event->status = TSEVENT_RESET_NEXTSTEP; 115 } 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "TSEventMonitorDestroy" 121 PetscErrorCode TSEventMonitorDestroy(TSEvent *event) 122 { 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 127 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 128 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 129 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 130 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 131 ierr = PetscFree(*event);CHKERRQ(ierr); 132 *event = NULL; 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "TSEventMonitor" 138 PetscErrorCode TSEventMonitor(TS ts) 139 { 140 PetscErrorCode ierr; 141 TSEvent event=ts->event; 142 PetscReal t; 143 Vec U; 144 PetscInt i; 145 PetscReal dt; 146 TSEventStatus status = event->status; 147 PetscInt rollback=0,in[2],out[2]; 148 149 PetscFunctionBegin; 150 151 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 152 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 153 154 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 155 if (event->status == TSEVENT_RESET_NEXTSTEP) { 156 /* Take initial time step */ 157 dt = event->initial_timestep; 158 ts->time_step = dt; 159 event->status = TSEVENT_NONE; 160 } 161 162 if (event->status == TSEVENT_NONE) { 163 event->tstepend = t; 164 } 165 166 event->nevents_zero = 0; 167 168 ierr = (*event->monitor)(ts,t,U,event->fvalue,event->monitorcontext);CHKERRQ(ierr); 169 if (event->status != TSEVENT_NONE) { 170 for (i=0; i < event->nevents; i++) { 171 if (PetscAbsScalar(event->fvalue[i]) < event->tol) { 172 event->status = TSEVENT_ZERO; 173 event->events_zero[event->nevents_zero++] = i; 174 } 175 } 176 } 177 178 status = event->status; 179 ierr = MPI_Allreduce((PetscEnum*)&status,(PetscEnum*)&event->status,1,MPIU_ENUM,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 180 181 if (event->status == TSEVENT_ZERO) { 182 dt = event->tstepend-t; 183 ts->time_step = dt; 184 ierr = TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,event->monitorcontext);CHKERRQ(ierr); 185 for (i = 0; i < event->nevents; i++) { 186 event->fvalue_prev[i] = event->fvalue[i]; 187 } 188 event->ptime_prev = t; 189 PetscFunctionReturn(0); 190 } 191 192 for (i = 0; i < event->nevents; i++) { 193 PetscInt fvalue_sign,fvalueprev_sign; 194 fvalue_sign = PetscSign(PetscRealPart(event->fvalue[i])); 195 fvalueprev_sign = PetscSign(PetscRealPart(event->fvalue_prev[i])); 196 if (fvalueprev_sign != 0 && (fvalue_sign != fvalueprev_sign)) { 197 switch (event->direction[i]) { 198 case -1: 199 if (fvalue_sign < 0) { 200 rollback = 1; 201 /* Compute linearly interpolated new time step */ 202 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 203 } 204 break; 205 case 1: 206 if (fvalue_sign > 0) { 207 rollback = 1; 208 /* Compute linearly interpolated new time step */ 209 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 210 } 211 break; 212 case 0: 213 rollback = 1; 214 /* Compute linearly interpolated new time step */ 215 dt = PetscMin(dt,PetscRealPart(-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i]))); 216 break; 217 } 218 } 219 } 220 if (rollback) event->status = TSEVENT_LOCATED_INTERVAL; 221 222 in[0] = event->status; 223 in[1] = rollback; 224 ierr = MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 225 226 rollback = out[1]; 227 if (rollback) { 228 event->status = TSEVENT_LOCATED_INTERVAL; 229 } 230 231 if (event->status == TSEVENT_LOCATED_INTERVAL) { 232 ierr = TSRollBack(ts);CHKERRQ(ierr); 233 event->status = TSEVENT_PROCESSING; 234 } else { 235 for (i = 0; i < event->nevents; i++) { 236 event->fvalue_prev[i] = event->fvalue[i]; 237 } 238 event->ptime_prev = t; 239 if (event->status == TSEVENT_PROCESSING) { 240 dt = event->tstepend - event->ptime_prev; 241 } 242 } 243 ierr = MPI_Allreduce(&dt,&(ts->time_step),1,MPIU_REAL,MPI_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247