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