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 . eventmonitor - event monitoring routine 15 . postevent - [optional] post-event function 16 - mectx - [optional] user-defined context for private data for the 17 event monitor and post event routine (use NULL if no 18 context is desired) 19 20 Calling sequence of eventmonitor: 21 PetscErrorCode EventMonitor(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,PetscInt *direction,PetscBool *terminate,void* mectx) 22 23 Input Parameters: 24 + ts - the TS context 25 . t - current time 26 . U - current iterate 27 - ctx - [optional] context passed with eventmonitor 28 29 Output parameters: 30 + fvalue - function value of events at time t 31 . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction, 32 +1 => Zero crossing in positive direction, 0 => both ways 33 - terminate - terminate time stepping after event is detected. 34 35 Calling sequence of postevent: 36 PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero, PetscInt events_zero, PetscReal t,Vec U,void* ctx) 37 38 Input Parameters: 39 + ts - the TS context 40 . nevents_zero - number of local events whose event function is zero 41 . events_zero - indices of local events which have reached zero 42 . t - current time 43 . U - current solution 44 - ctx - the context passed with eventmonitor 45 46 Level: intermediate 47 48 .keywords: TS, event, set, monitor 49 50 .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason() 51 @*/ 52 PetscErrorCode TSSetEventMonitor(TS ts,PetscInt nevents,PetscErrorCode (*eventmonitor)(TS,PetscReal,Vec,PetscScalar*,PetscInt*,PetscBool*,void*),PetscErrorCode (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,void*),void *mectx) 53 { 54 PetscErrorCode ierr; 55 PetscReal t; 56 Vec U; 57 TSEvent event; 58 59 PetscFunctionBegin; 60 ierr = PetscNew(struct _p_TSEvent,&ts->event);CHKERRQ(ierr); 61 event = ts->event; 62 ierr = PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue);CHKERRQ(ierr); 63 ierr = PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue_prev);CHKERRQ(ierr); 64 ierr = PetscMalloc(nevents*sizeof(PetscBool),&event->terminate);CHKERRQ(ierr); 65 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->direction);CHKERRQ(ierr); 66 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->events_zero);CHKERRQ(ierr); 67 event->monitor = eventmonitor; 68 event->postevent = postevent; 69 event->monitorcontext = (void*)mectx; 70 event->nevents = nevents; 71 72 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 73 ierr = TSGetTimeStep(ts,&event->initial_timestep);CHKERRQ(ierr); 74 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 75 event->ptime_prev = t; 76 ierr = (*event->monitor)(ts,t,U,event->fvalue_prev,event->direction,event->terminate,mectx);CHKERRQ(ierr); 77 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TS Event options","");CHKERRQ(ierr); 78 { 79 event->tol = 1.0e-6; 80 ierr = PetscOptionsReal("-ts_event_tol","","",event->tol,&event->tol,NULL);CHKERRQ(ierr); 81 } 82 ierr = PetscOptionsEnd();CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "TSPostEvent" 88 PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,void *ctx) 89 { 90 PetscErrorCode ierr; 91 TSEvent event=ts->event; 92 PetscBool terminate=PETSC_FALSE; 93 PetscInt i; 94 PetscBool ts_terminate; 95 96 PetscFunctionBegin; 97 if (event->postevent) { 98 ierr = (*event->postevent)(ts,nevents_zero,events_zero,t,U,ctx);CHKERRQ(ierr); 99 } 100 for(i = 0; i < nevents_zero;i++) { 101 terminate = terminate || event->terminate[events_zero[i]]; 102 } 103 ierr = MPI_Allreduce(&terminate,&ts_terminate,1,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 104 if (terminate) { 105 ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr); 106 event->status = TSEVENT_NONE; 107 } else { 108 event->status = TSEVENT_RESET_NEXTSTEP; 109 } 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "TSEventMonitorDestroy" 115 PetscErrorCode TSEventMonitorDestroy(TSEvent *event) 116 { 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 121 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 122 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 123 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 124 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 125 ierr = PetscFree(*event);CHKERRQ(ierr); 126 *event = NULL; 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "TSEventMonitor" 132 PetscErrorCode TSEventMonitor(TS ts) 133 { 134 PetscErrorCode ierr; 135 TSEvent event=ts->event; 136 PetscReal t; 137 Vec U; 138 PetscInt i; 139 PetscReal dt; 140 PetscInt status=event->status; 141 PetscInt rollback=0,in[2],out[2]; 142 143 PetscFunctionBegin; 144 145 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 146 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 147 148 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 149 if (event->status == TSEVENT_RESET_NEXTSTEP) { 150 /* Take initial time step */ 151 dt = event->initial_timestep; 152 ts->time_step = dt; 153 event->status = TSEVENT_NONE; 154 } 155 156 if (event->status == TSEVENT_NONE) { 157 event->tstepend = t; 158 } 159 160 event->nevents_zero = 0; 161 162 ierr = (*event->monitor)(ts,t,U,event->fvalue,event->direction,event->terminate,event->monitorcontext);CHKERRQ(ierr); 163 if (event->status != TSEVENT_NONE) { 164 for (i=0; i < event->nevents; i++) { 165 if (PetscAbs(event->fvalue[i]) < event->tol) { 166 event->status = TSEVENT_ZERO; 167 event->events_zero[event->nevents_zero++] = i; 168 } 169 } 170 } 171 172 status = event->status; 173 ierr = MPI_Allreduce(&status,&event->status,1,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 174 175 if (event->status == TSEVENT_ZERO) { 176 dt = event->tstepend-t; 177 ts->time_step = dt; 178 ierr = TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,event->monitorcontext);CHKERRQ(ierr); 179 for (i = 0; i < event->nevents; i++) { 180 event->fvalue_prev[i] = event->fvalue[i]; 181 } 182 event->ptime_prev = t; 183 PetscFunctionReturn(0); 184 } 185 186 for (i = 0; i < event->nevents; i++) { 187 if ((event->direction[i] < 0 && PetscSign(event->fvalue[i]) <= 0 && PetscSign(event->fvalue_prev[i]) >= 0) || \ 188 (event->direction[i] > 0 && PetscSign(event->fvalue[i]) >= 0 && PetscSign(event->fvalue_prev[i]) <= 0) || \ 189 (event->direction[i] == 0 && PetscSign(event->fvalue[i])*PetscSign(event->fvalue_prev[i]) <= 0)) { 190 191 event->status = TSEVENT_LOCATED_INTERVAL; 192 rollback = 1; 193 /* Compute linearly interpolated new time step */ 194 dt = PetscMin(dt,-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i])); 195 } 196 } 197 in[0] = event->status; 198 in[1] = rollback; 199 ierr = MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 200 201 rollback = out[1]; 202 if (rollback) { 203 event->status = TSEVENT_LOCATED_INTERVAL; 204 } 205 206 if (event->status == TSEVENT_LOCATED_INTERVAL) { 207 ierr = TSRollBack(ts);CHKERRQ(ierr); 208 event->status = TSEVENT_PROCESSING; 209 } else { 210 for (i = 0; i < event->nevents; i++) { 211 event->fvalue_prev[i] = event->fvalue[i]; 212 } 213 event->ptime_prev = t; 214 if (event->status == TSEVENT_PROCESSING) { 215 dt = event->tstepend - event->ptime_prev; 216 } 217 } 218 PetscReal time_step; 219 ierr = MPI_Allreduce(&dt,&time_step,1,MPIU_REAL,MPI_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 220 ts->time_step = time_step; 221 PetscFunctionReturn(0); 222 } 223 224