Actual source code: ex2.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"
9: "Using two event functions = piecewise-polynomials, zeros = 1 (rank-0), 9 (last rank)\n"
10: "Options:\n"
11: "-dir d : zero-crossing direction for events\n"
12: "-flg : additional output in Postevent\n"
13: "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n"
14: "-restart : flag for TSRestartStep() in PostEvent\n"
15: "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n"
16: " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n"
17: " if x == 0, nothing happens\n";
19: #define MAX_NFUNC 100 // max event functions per rank
20: #define MAX_NEV 5000 // max zero crossings for each rank
22: typedef struct {
23: PetscMPIInt rank, size;
24: PetscReal pi;
25: PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals
26: PetscReal evres[MAX_NEV]; // times of found zero-crossings
27: PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking
28: PetscInt cnt; // counter
29: PetscInt cntref; // actual length of 'ref' on the given rank
30: PetscBool flg; // flag for additional print in PostEvent
31: PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default)
32: PetscBool restart; // flag for TSRestartStep() in PostEvent
33: PetscReal dtpost; // post-event step
34: PetscInt postcnt; // counter for PostEvent calls
35: } AppCtx;
37: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx);
38: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx);
40: int main(int argc, char **argv)
41: {
42: TS ts;
43: Mat A;
44: Vec sol;
45: PetscInt n, dir0, m = 0;
46: PetscInt dir[MAX_NFUNC], inds[2];
47: PetscBool term[MAX_NFUNC];
48: PetscScalar *x, vals[4];
49: AppCtx ctx;
51: PetscFunctionBeginUser;
52: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
53: setbuf(stdout, NULL);
54: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
55: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
56: ctx.pi = PetscAcosReal(-1.0);
57: ctx.cnt = 0;
58: ctx.cntref = 0;
59: ctx.flg = PETSC_FALSE;
60: ctx.errtol = 1e-5;
61: ctx.restart = PETSC_FALSE;
62: ctx.dtpost = 0;
63: ctx.postcnt = 0;
65: // The linear problem has a 2*2 matrix. The matrix is constant
66: if (ctx.rank == 0) m = 2;
67: inds[0] = 0;
68: inds[1] = 1;
69: vals[0] = 0;
70: vals[1] = 0.2;
71: vals[2] = -0.2;
72: vals[3] = 0;
73: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A));
74: PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES));
75: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
76: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
77: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
79: PetscCall(MatCreateVecs(A, &sol, NULL));
80: PetscCall(VecGetArray(sol, &x));
81: if (ctx.rank == 0) { // initial conditions
82: x[0] = 0; // sin(0)
83: x[1] = 1; // cos(0)
84: }
85: PetscCall(VecRestoreArray(sol, &x));
87: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
88: PetscCall(TSSetProblemType(ts, TS_LINEAR));
90: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
91: PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL));
93: PetscCall(TSSetTimeStep(ts, 0.1));
94: PetscCall(TSSetType(ts, TSBEULER));
95: PetscCall(TSSetMaxSteps(ts, 10000));
96: PetscCall(TSSetMaxTime(ts, 10.0));
97: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
98: PetscCall(TSSetFromOptions(ts));
100: // Set the event handling
101: dir0 = 0;
102: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction
103: PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output
104: PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events
105: PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep()
106: PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step
108: n = 0; // event counter
109: if (ctx.rank == 0) { // first event -- on rank-0
110: dir[n] = dir0;
111: term[n++] = PETSC_FALSE;
112: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0;
113: }
114: if (ctx.rank == ctx.size - 1) { // second event -- on last rank
115: dir[n] = dir0;
116: term[n++] = PETSC_FALSE;
117: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 9.0;
118: }
119: if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
120: PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
122: // Solution
123: PetscCall(TSSolve(ts, sol));
125: // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
126: for (PetscInt j = 0; j < ctx.cnt; j++) {
127: PetscReal err = 10.0;
128: if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
129: 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"));
130: }
131: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
133: PetscCall(MatDestroy(&A));
134: PetscCall(TSDestroy(&ts));
135: PetscCall(VecDestroy(&sol));
137: PetscCall(PetscFinalize());
138: return 0;
139: }
141: /*
142: User callback for defining the event-functions
143: */
144: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx)
145: {
146: PetscInt n = 0;
147: AppCtx *Ctx = (AppCtx *)ctx;
149: PetscFunctionBeginUser;
150: // for the test purposes, event-functions are defined based on t
151: // first event -- on rank-0
152: if (Ctx->rank == 0) {
153: if (t < 2.0) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.0, 12));
154: else gval[n++] = 0.5;
155: }
157: // second event -- on last rank
158: if (Ctx->rank == Ctx->size - 1) {
159: if (t > 8.0) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.0, 12));
160: else gval[n++] = 0.25;
161: }
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /*
166: User callback for the post-event stuff
167: */
168: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx)
169: {
170: AppCtx *Ctx = (AppCtx *)ctx;
172: PetscFunctionBeginUser;
173: if (Ctx->flg) {
174: PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
175: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
176: for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
177: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
178: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
179: }
181: if (Ctx->cnt + nev_zero < MAX_NEV)
182: for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
184: #ifdef NEW_VERSION
185: Ctx->postcnt++; // sync
186: if (Ctx->dtpost > 0) {
187: if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
188: else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
189: }
190: #endif
192: if (Ctx->restart) PetscCall(TSRestartStep(ts));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
195: /*---------------------------------------------------------------------------------------------*/
196: /*
197: Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
198: which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
199: explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
200: simply remove this option altogether. This will result in using the defaults
201: (which is PETSC_DECIDE).
202: */
203: /*TEST
204: test:
205: suffix: 0s1
206: requires: !single
207: output_file: output/ex2_0s1.out
208: args: -dir 0
209: args: -restart 0
210: args: -dtpost 0.25
211: args: -ts_event_post_event_step {{-1 0.31}}
212: args: -ts_type {{beuler rk}}
213: args: -ts_adapt_type {{none basic}}
214: nsize: 1
216: test:
217: suffix: 0s4
218: requires: !single
219: output_file: output/ex2_0s4.out
220: args: -dir 0
221: args: -restart 1
222: args: -dtpost {{0 0.25}}
223: args: -ts_event_post_event_step -1
224: args: -ts_type {{beuler rk}}
225: args: -ts_adapt_type {{none basic}}
226: nsize: 4
227: filter: sort
228: filter_output: sort
230: test:
231: suffix: pos
232: requires: !single
233: output_file: output/ex2_pos.out
234: args: -dir 1
235: args: -restart {{0 1}}
236: args: -dtpost 0
237: args: -ts_event_post_event_step 0.31005
238: args: -ts_type rk
239: args: -ts_adapt_type {{none basic}}
240: nsize: {{1 4}}
241: filter: sort
242: filter_output: sort
244: test:
245: suffix: ns1
246: requires: !single
247: output_file: output/ex2_ns1.out
248: args: -dir -1
249: args: -restart 0
250: args: -dtpost 0.25
251: args: -ts_event_post_event_step {{-1 0.305}}
252: args: -ts_type {{beuler rk}}
253: args: -ts_adapt_type {{none basic}}
254: nsize: 1
256: test:
257: suffix: ns4
258: requires: !single
259: output_file: output/ex2_ns4.out
260: args: -dir -1
261: args: -restart 1
262: args: -dtpost {{0 0.25}}
263: args: -ts_event_post_event_step -1
264: args: -ts_type {{beuler rk}}
265: args: -ts_adapt_type {{none basic}}
266: nsize: 4
267: filter: sort
268: filter_output: sort
270: test:
271: suffix: 0s1single
272: requires: single
273: output_file: output/ex2_0s1.out
274: args: -dir 0 -ts_event_dt_min 1e-6 -errtol 5e-5
275: args: -restart {{0 1}}
276: args: -dtpost 0
277: args: -ts_event_post_event_step 0.31
278: args: -ts_type {{beuler rk}}
279: args: -ts_adapt_type {{none basic}}
280: nsize: 1
282: test:
283: suffix: 0s4single
284: requires: single
285: output_file: output/ex2_0s4.out
286: args: -dir 0 -ts_event_dt_min 1e-6 -errtol 5e-5
287: args: -restart 0
288: args: -dtpost 0.25
289: args: -ts_event_post_event_step {{-1 0.315}}
290: args: -ts_type {{beuler rk}}
291: args: -ts_adapt_type {{none basic}}
292: nsize: 4
293: filter: sort
294: filter_output: sort
296: test:
297: suffix: possingle
298: requires: single
299: output_file: output/ex2_pos.out
300: args: -dir 1 -ts_event_dt_min 1e-6 -errtol 5e-5
301: args: -restart 1
302: args: -dtpost {{0 0.25}}
303: args: -ts_event_post_event_step -1
304: args: -ts_type {{beuler rk}}
305: args: -ts_adapt_type basic
306: nsize: {{1 4}}
307: filter: sort
308: filter_output: sort
310: test:
311: suffix: ns1single
312: requires: single
313: output_file: output/ex2_ns1.out
314: args: -dir -1 -ts_event_dt_min 1e-6 -errtol 5e-5
315: args: -restart {{0 1}}
316: args: -dtpost 0
317: args: -ts_event_post_event_step 0.30501
318: args: -ts_type {{beuler rk}}
319: args: -ts_adapt_type {{none basic}}
320: nsize: 1
322: test:
323: suffix: ns4single
324: requires: single
325: output_file: output/ex2_ns4.out
326: args: -dir -1 -ts_event_dt_min 1e-6 -errtol 5e-5
327: args: -restart 0
328: args: -dtpost 0.25
329: args: -ts_event_post_event_step {{-1 0.31}}
330: args: -ts_type {{beuler rk}}
331: args: -ts_adapt_type {{none basic}}
332: nsize: 4
333: filter: sort
334: filter_output: sort
335: TEST*/