Actual source code: ex3.c
1: #include <petscts.h>
2: #include <stdio.h>
4: #define NEW_VERSION // Applicable for the new features; avoid this for the older PETSc versions (without TSSetPostEventStep())
6: static char help[] = "Simple linear problem with events\n"
7: "x_dot = 0.2*y\n"
8: "y_dot = -0.2*x\n"
10: "The following event functions are involved:\n"
11: "- two polynomial event functions on rank-0 and last-rank (with zeros: 1.05, 9.05[terminating])\n"
12: "- one event function on rank = '1%size', equal to V*sin(pi*t), zeros = 1,...,10\n"
13: "After each event location the tolerance for the sin() event is multiplied by 4\n"
15: "Options:\n"
16: "-dir d : zero-crossing direction for events: 0, 1, -1\n"
17: "-flg : additional output in Postevent\n"
18: "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n"
19: "-restart : flag for TSRestartStep() in PostEvent\n"
20: "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n"
21: " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n"
22: " if x == 0, nothing happens\n"
23: "-v {float}: scaling of the sin() event function; for small v this event is triggered by the function values,\n"
24: " for large v the event is triggered by the small step size\n"
25: "-change5 : flag to change the state vector at t=5 PostEvent\n";
27: #define MAX_NFUNC 100 // max event functions per rank
28: #define MAX_NEV 5000 // max zero crossings for each rank
30: typedef struct {
31: PetscMPIInt rank, size;
32: PetscReal pi;
33: PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals
34: PetscReal evres[MAX_NEV]; // times of found zero-crossings
35: PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking
36: PetscInt cnt; // counter
37: PetscInt cntref; // actual length of 'ref' on the given rank
38: PetscBool flg; // flag for additional print in PostEvent
39: PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default)
40: PetscBool restart; // flag for TSRestartStep() in PostEvent
41: PetscReal dtpost; // post-event step
42: PetscInt postcnt; // counter for PostEvent calls
43: PetscReal V; // vertical scaling for sin()
44: PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage
45: PetscBool change5; // flag to change the state vector at t=5 PostEvent
46: } AppCtx;
48: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx);
49: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx);
51: int main(int argc, char **argv)
52: {
53: TS ts;
54: Mat A;
55: Vec sol;
56: PetscInt n, dir0, m = 0;
57: PetscReal tol = 1e-7;
58: PetscInt dir[MAX_NFUNC], inds[2];
59: PetscBool term[MAX_NFUNC];
60: PetscScalar *x, vals[4];
61: AppCtx ctx;
62: TSConvergedReason reason;
64: PetscFunctionBeginUser;
65: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
66: setbuf(stdout, NULL);
67: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
68: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
69: ctx.pi = PetscAcosReal(-1.0);
70: ctx.cnt = 0;
71: ctx.cntref = 0;
72: ctx.flg = PETSC_FALSE;
73: ctx.errtol = 1e-5;
74: ctx.restart = PETSC_FALSE;
75: ctx.dtpost = 0;
76: ctx.postcnt = 0;
77: ctx.V = 1.0;
78: ctx.change5 = PETSC_FALSE;
80: // The linear problem has a 2*2 matrix. The matrix is constant
81: if (ctx.rank == 0) m = 2;
82: inds[0] = 0;
83: inds[1] = 1;
84: vals[0] = 0;
85: vals[1] = 0.2;
86: vals[2] = -0.2;
87: vals[3] = 0;
88: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A));
89: PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES));
90: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
91: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
92: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
94: PetscCall(MatCreateVecs(A, &sol, NULL));
95: PetscCall(VecGetArray(sol, &x));
96: if (ctx.rank == 0) { // initial conditions
97: x[0] = 0; // sin(0)
98: x[1] = 1; // cos(0)
99: }
100: PetscCall(VecRestoreArray(sol, &x));
102: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
103: PetscCall(TSSetProblemType(ts, TS_LINEAR));
105: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
106: PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL));
108: PetscCall(TSSetTimeStep(ts, 0.1));
109: PetscCall(TSSetType(ts, TSBEULER));
110: PetscCall(TSSetMaxSteps(ts, 10000));
111: PetscCall(TSSetMaxTime(ts, 10.0));
112: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
113: PetscCall(TSSetFromOptions(ts));
115: // Set the event handling
116: dir0 = 0;
117: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction
118: PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output
119: PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events
120: PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep()
121: PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step
122: PetscCall(PetscOptionsGetReal(NULL, NULL, "-v", &ctx.V, NULL));
123: PetscCall(PetscOptionsGetBool(NULL, NULL, "-change5", &ctx.change5, NULL)); // flag to change the state vector at t=5 PostEvent
125: n = 0; // event counter
126: if (ctx.rank == 0) { // first event -- on rank-0
127: ctx.vtol[n] = tol * 10;
128: dir[n] = dir0;
129: term[n++] = PETSC_FALSE;
130: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.05;
131: }
132: if (ctx.rank == ctx.size - 1) { // second event (with termination) -- on last rank
133: ctx.vtol[n] = tol * 10;
134: dir[n] = dir0;
135: term[n++] = PETSC_TRUE;
136: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 9.05;
137: }
138: if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size
139: ctx.vtol[n] = tol;
140: dir[n] = dir0;
141: term[n++] = PETSC_FALSE;
143: for (PetscInt i = 1; i < MAX_NEV - 2; i++) {
144: if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i;
145: if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i;
146: }
147: }
148: if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
149: PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
150: PetscCall(TSSetEventTolerances(ts, tol, ctx.vtol));
152: // Solution
153: PetscCall(TSSolve(ts, sol));
154: PetscCall(TSGetConvergedReason(ts, &reason));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT));
157: // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
158: for (PetscInt j = 0; j < ctx.cnt; j++) {
159: PetscReal err = 10.0;
160: if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
161: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d\t%g\t%g\t%s\n", ctx.rank, (double)ctx.evres[j], (double)err, err < ctx.errtol ? "pass" : "fail"));
162: }
163: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
165: PetscCall(MatDestroy(&A));
166: PetscCall(TSDestroy(&ts));
167: PetscCall(VecDestroy(&sol));
169: PetscCall(PetscFinalize());
170: return 0;
171: }
173: /*
174: User callback for defining the event-functions
175: */
176: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx)
177: {
178: PetscInt n = 0;
179: AppCtx *Ctx = (AppCtx *)ctx;
181: PetscFunctionBeginUser;
182: // for the test purposes, event-functions are defined based on t
183: // first event -- on rank-0
184: if (Ctx->rank == 0) {
185: if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12));
186: else gval[n++] = 0.5;
187: }
189: // second event -- on last rank
190: if (Ctx->rank == Ctx->size - 1) {
191: if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12));
192: else gval[n++] = 0.25;
193: }
195: // third event -- on rank = 1%ctx.size
196: if (Ctx->rank == 1 % Ctx->size) { gval[n++] = Ctx->V * PetscSinReal(Ctx->pi * t); }
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /*
201: User callback for the post-event stuff
202: */
203: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx)
204: {
205: PetscInt n = 0;
206: PetscScalar *x;
207: AppCtx *Ctx = (AppCtx *)ctx;
209: PetscFunctionBeginUser;
210: if (Ctx->flg) {
211: PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
212: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
213: for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
214: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
215: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
216: }
218: // change the state vector near t=5.0
219: if (PetscAbsReal(t - (PetscReal)5.0) < 0.01 && Ctx->change5) {
220: PetscCall(VecGetArray(U, &x));
221: if (Ctx->rank == 0) x[1] = -x[1];
222: PetscCall(VecRestoreArray(U, &x));
223: }
225: // update vtol's
226: if (Ctx->rank == 0) n++; // first event -- on rank-0
227: if (Ctx->rank == Ctx->size - 1) n++; // second event -- on last rank
228: if (Ctx->rank == 1 % Ctx->size) { // third event -- on rank = 1%ctx.size
229: if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "vtol for sin: %g -> ", (double)Ctx->vtol[n]));
230: Ctx->vtol[n] *= 4;
231: if (PetscAbsReal(t - (PetscReal)5.0) < 0.01) Ctx->vtol[n] /= 100; // one-off decrease
232: if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g\n", (double)Ctx->vtol[n]));
233: n++;
234: }
235: PetscCall(TSSetEventTolerances(ts, 0, Ctx->vtol));
237: if (Ctx->cnt + nev_zero < MAX_NEV)
238: for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
240: #ifdef NEW_VERSION
241: Ctx->postcnt++; // sync
242: if (Ctx->dtpost > 0) {
243: if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
244: else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
245: }
246: #endif
248: if (Ctx->restart) PetscCall(TSRestartStep(ts));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
251: /*---------------------------------------------------------------------------------------------*/
252: /*
253: Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
254: which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
255: explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
256: simply remove this option altogether. This will result in using the defaults
257: (which is PETSC_DECIDE).
258: */
259: /*TEST
260: test:
261: suffix: V
262: output_file: output/ex3_V.out
263: args: -ts_type beuler
264: args: -ts_adapt_type basic
265: args: -v {{1e2 1e5 1e8}}
266: args: -ts_adapt_dt_min 1e-6
267: args: -change5 {{0 1}}
268: nsize: 1
270: test:
271: suffix: neu1
272: output_file: output/ex3_neu1.out
273: args: -dir 0
274: args: -v 1e5
275: args: -ts_adapt_dt_min 1e-6
276: args: -restart 1
277: args: -dtpost 0.24
278: args: -ts_event_post_event_step 0.31
279: args: -ts_type {{beuler rk}}
280: args: -ts_adapt_type {{none basic}}
281: nsize: 1
283: test:
284: suffix: neu2
285: output_file: output/ex3_neu2.out
286: args: -dir 0
287: args: -v 1e5
288: args: -ts_adapt_dt_min 1e-6
289: args: -restart 1
290: args: -dtpost 0
291: args: -ts_event_post_event_step {{-1 0.31}}
292: args: -ts_type rk
293: args: -ts_adapt_type {{none basic}}
294: nsize: 2
295: filter: sort
296: filter_output: sort
298: test:
299: suffix: neu4
300: output_file: output/ex3_neu4.out
301: args: -dir 0
302: args: -v 1e5
303: args: -ts_adapt_dt_min 1e-6
304: args: -restart {{0 1}}
305: args: -dtpost 0.24
306: args: -ts_event_post_event_step 0.21
307: args: -ts_type beuler
308: args: -ts_adapt_type {{none basic}}
309: nsize: 4
310: filter: sort
311: filter_output: sort
313: test:
314: suffix: pos1
315: output_file: output/ex3_pos1.out
316: args: -dir 1
317: args: -v 1e5
318: args: -ts_adapt_dt_min 1e-6
319: args: -restart 0
320: args: -dtpost 0.24
321: args: -ts_type {{beuler rk}}
322: args: -ts_adapt_type {{none basic}}
323: nsize: 1
325: test:
326: suffix: pos2
327: output_file: output/ex3_pos2.out
328: args: -dir 1
329: args: -v 1e5
330: args: -ts_adapt_dt_min 1e-6
331: args: -restart 1
332: args: -dtpost {{0 0.24}}
333: args: -ts_type rk
334: args: -ts_adapt_type {{none basic}}
335: nsize: 2
336: filter: sort
337: filter_output: sort
339: test:
340: suffix: pos4
341: output_file: output/ex3_pos4.out
342: args: -dir 1
343: args: -v 1e9
344: args: -ts_adapt_dt_min 1e-6
345: args: -restart 0
346: args: -dtpost 0
347: args: -ts_event_post_event_step {{-1 0.32}}
348: args: -ts_type beuler
349: args: -ts_adapt_type {{none basic}}
350: args: -change5 1
351: nsize: 4
352: filter: sort
353: filter_output: sort
355: test:
356: suffix: neg1
357: output_file: output/ex3_neg1.out
358: args: -dir -1
359: args: -v 1e5
360: args: -ts_adapt_dt_min 1e-6
361: args: -restart 1
362: args: -dtpost {{0 0.24}}
363: args: -ts_type {{beuler rk}}
364: args: -ts_adapt_type basic
365: nsize: 1
367: test:
368: suffix: neg2
369: output_file: output/ex3_neg2.out
370: args: -dir -1
371: args: -v 1e5
372: args: -ts_adapt_dt_min 1e-6
373: args: -restart 0
374: args: -dtpost {{0 0.24}}
375: args: -ts_type rk
376: args: -ts_adapt_type {{none basic}}
377: nsize: 2
378: filter: sort
379: filter_output: sort
381: test:
382: suffix: neg4
383: output_file: output/ex3_neg4.out
384: args: -dir -1
385: args: -v 1e5
386: args: -ts_adapt_dt_min 1e-6
387: args: -restart 0
388: args: -dtpost {{0 0.24}}
389: args: -ts_event_post_event_step 0.3
390: args: -ts_type beuler
391: args: -ts_adapt_type {{none basic}}
392: nsize: 4
393: filter: sort
394: filter_output: sort
395: TEST*/