Actual source code: tsevent.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
6: /*@C
7: TSSetEventMonitor - Sets a monitoring function used for detecting events
9: Logically Collective on TS
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)
24: Calling sequence of eventmonitor:
25: PetscErrorCode EventMonitor(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void* mectx)
27: Input Parameters:
28: + ts - the TS context
29: . t - current time
30: . U - current iterate
31: - ctx - [optional] context passed with eventmonitor
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)
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
47: Level: intermediate
49: .keywords: TS, event, set, monitor
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: {
56: PetscReal t;
57: Vec U;
58: TSEvent event;
59: PetscInt i;
62: PetscNew(&event);
63: PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue);
64: PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue_prev);
65: PetscMalloc(nevents*sizeof(PetscInt),&event->direction);
66: PetscMalloc(nevents*sizeof(PetscInt),&event->terminate);
67: for (i=0; i < nevents; i++) {
68: event->direction[i] = direction[i];
69: event->terminate[i] = terminate[i];
70: }
71: PetscMalloc(nevents*sizeof(PetscInt),&event->events_zero);
72: event->monitor = eventmonitor;
73: event->postevent = postevent;
74: event->monitorcontext = (void*)mectx;
75: event->nevents = nevents;
77: TSGetTime(ts,&t);
78: TSGetTimeStep(ts,&event->initial_timestep);
79: TSGetSolution(ts,&U);
80: event->ptime_prev = t;
81: (*event->monitor)(ts,t,U,event->fvalue_prev,mectx);
82: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TS Event options","");
83: {
84: event->tol = 1.0e-6;
85: PetscOptionsReal("-ts_event_tol","","",event->tol,&event->tol,NULL);
86: }
87: PetscOptionsEnd();
88: ts->event = event;
89: return(0);
90: }
94: PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,void *ctx)
95: {
97: TSEvent event=ts->event;
98: PetscBool terminate=PETSC_FALSE;
99: PetscInt i;
100: PetscBool ts_terminate;
103: if (event->postevent) {
104: (*event->postevent)(ts,nevents_zero,events_zero,t,U,ctx);
105: }
106: for(i = 0; i < nevents_zero;i++) {
107: terminate = (PetscBool)(terminate || event->terminate[events_zero[i]]);
108: }
109: MPI_Allreduce(&terminate,&ts_terminate,1,MPIU_BOOL,MPI_LOR,((PetscObject)ts)->comm);
110: if (terminate) {
111: TSSetConvergedReason(ts,TS_CONVERGED_EVENT);
112: event->status = TSEVENT_NONE;
113: } else {
114: event->status = TSEVENT_RESET_NEXTSTEP;
115: }
116: return(0);
117: }
121: PetscErrorCode TSEventMonitorDestroy(TSEvent *event)
122: {
126: PetscFree((*event)->fvalue);
127: PetscFree((*event)->fvalue_prev);
128: PetscFree((*event)->direction);
129: PetscFree((*event)->terminate);
130: PetscFree((*event)->events_zero);
131: PetscFree(*event);
132: *event = NULL;
133: return(0);
134: }
138: PetscErrorCode TSEventMonitor(TS ts)
139: {
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];
151: TSGetTime(ts,&t);
152: TSGetSolution(ts,&U);
154: TSGetTimeStep(ts,&dt);
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: }
162: if (event->status == TSEVENT_NONE) {
163: event->tstepend = t;
164: }
166: event->nevents_zero = 0;
168: (*event->monitor)(ts,t,U,event->fvalue,event->monitorcontext);
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: }
178: status = event->status;
179: MPI_Allreduce((PetscEnum*)&status,(PetscEnum*)&event->status,1,MPIU_ENUM,MPI_MAX,((PetscObject)ts)->comm);
181: if (event->status == TSEVENT_ZERO) {
182: dt = event->tstepend-t;
183: ts->time_step = dt;
184: TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,event->monitorcontext);
185: for (i = 0; i < event->nevents; i++) {
186: event->fvalue_prev[i] = event->fvalue[i];
187: }
188: event->ptime_prev = t;
189: return(0);
190: }
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: MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);
225:
226: rollback = out[1];
227: if (rollback) {
228: event->status = TSEVENT_LOCATED_INTERVAL;
229: }
231: if (event->status == TSEVENT_LOCATED_INTERVAL) {
232: TSRollBack(ts);
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: MPI_Allreduce(&dt,&(ts->time_step),1,MPIU_REAL,MPI_MIN,((PetscObject)ts)->comm);
244: return(0);
245: }