Actual source code: ex4.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 one or several event functions (on rank-0)\n"
10: "This program is mostly intended to test the Anderson-Bjorck iteration with challenging event-functions\n"
11: "Options:\n"
12: "-dir d : zero-crossing direction for events\n"
13: "-flg : additional output in Postevent\n"
14: "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n"
15: "-restart : flag for TSRestartStep() in PostEvent\n"
16: "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n"
17: " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n"
18: " if x == 0, nothing happens\n"
19: "-func F : selects the event function [0, ..., 11], if F == -1 (default) is set, all event functions are taken\n";
21: #define MAX_NFUNC 100 // max event functions per rank
22: #define MAX_NEV 5000 // max zero crossings for each rank
24: typedef struct {
25: PetscMPIInt rank, size;
26: PetscReal pi;
27: PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals
28: PetscReal evres[MAX_NEV]; // times of found zero-crossings
29: PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking
30: PetscInt cnt; // counter
31: PetscInt cntref; // actual length of 'ref' on the given rank
32: PetscBool flg; // flag for additional print in PostEvent
33: PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default)
34: PetscBool restart; // flag for TSRestartStep() in PostEvent
35: PetscReal dtpost; // post-event step
36: PetscInt postcnt; // counter for PostEvent calls
37: PetscInt F; // event-function index
38: PetscInt Fnum; // total available event functions
39: } AppCtx;
41: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx);
42: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx);
44: int main(int argc, char **argv)
45: {
46: TS ts;
47: Mat A;
48: Vec sol;
49: PetscInt n, dir0, m = 0;
50: PetscInt dir[MAX_NFUNC], inds[2];
51: PetscBool term[MAX_NFUNC];
52: PetscBool pass = PETSC_TRUE;
53: PetscScalar *x, vals[4];
54: AppCtx ctx;
56: PetscFunctionBeginUser;
57: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
58: setbuf(stdout, NULL);
59: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
60: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
61: ctx.pi = PetscAcosReal(-1.0);
62: ctx.cnt = 0;
63: ctx.cntref = 0;
64: ctx.flg = PETSC_FALSE;
65: ctx.errtol = 1e-5;
66: ctx.restart = PETSC_FALSE;
67: ctx.dtpost = 0;
68: ctx.postcnt = 0;
69: ctx.F = -1;
70: ctx.Fnum = 12;
72: // The linear problem has a 2*2 matrix. The matrix is constant
73: if (ctx.rank == 0) m = 2;
74: inds[0] = 0;
75: inds[1] = 1;
76: vals[0] = 0;
77: vals[1] = 0.2;
78: vals[2] = -0.2;
79: vals[3] = 0;
80: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A));
81: PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES));
82: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
84: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
86: PetscCall(MatCreateVecs(A, &sol, NULL));
87: PetscCall(VecGetArray(sol, &x));
88: if (ctx.rank == 0) { // initial conditions
89: x[0] = 0; // sin(0)
90: x[1] = 1; // cos(0)
91: }
92: PetscCall(VecRestoreArray(sol, &x));
94: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
95: PetscCall(TSSetProblemType(ts, TS_LINEAR));
97: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
98: PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL));
100: PetscCall(TSSetTime(ts, 0.03));
101: PetscCall(TSSetTimeStep(ts, 0.1));
102: PetscCall(TSSetType(ts, TSBEULER));
103: PetscCall(TSSetMaxSteps(ts, 10000));
104: PetscCall(TSSetMaxTime(ts, 4.0));
105: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
106: PetscCall(TSSetFromOptions(ts));
108: // Set the event handling
109: dir0 = 0;
110: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &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
115: PetscCall(PetscOptionsGetInt(NULL, NULL, "-F", &ctx.F, NULL)); // event-function index
116: PetscCheck(ctx.F >= -1 && ctx.F < ctx.Fnum, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Value of 'F' is out of range");
118: n = 0; // event counter
119: if (ctx.rank == 0) { // all events -- on rank-0
120: if (ctx.F == -1)
121: for (n = 0; n < ctx.Fnum; n++) { // all event-functions
122: dir[n] = dir0;
123: term[n] = PETSC_FALSE;
124: }
125: else { // single event-function
126: dir[n] = dir0;
127: term[n++] = PETSC_FALSE;
128: }
130: // set the reference values
131: if (ctx.F == 0 || ctx.F == -1) {
132: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 2.0;
133: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 4.0;
134: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 1.0;
135: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 3.0;
136: }
137: if (ctx.F == 1 || ctx.F == -1)
138: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0;
139: if (ctx.F == 2 || ctx.F == -1)
140: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0;
141: if (ctx.F == 3 || ctx.F == -1)
142: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.696812386809752;
143: if (ctx.F == 4 || ctx.F == -1)
144: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.1;
145: if (ctx.F == 5 || ctx.F == -1)
146: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.02;
147: if (ctx.F == 6 || ctx.F == -1)
148: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 0.90478821787302;
149: if (ctx.F == 7 || ctx.F == -1)
150: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.05;
151: if (ctx.F == 8 || ctx.F == -1)
152: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.552704666678489;
153: if (ctx.F == 9 || ctx.F == -1)
154: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.399422291710969;
155: if (ctx.F == 10 || ctx.F == -1) {
156: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.874511220376091;
157: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 2.388359335869107;
158: }
159: if (ctx.F == 11 || ctx.F == -1)
160: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 2.5;
161: }
162: if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
163: PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
165: // Solution
166: PetscCall(TSSolve(ts, sol));
168: // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
169: for (PetscInt j = 0; j < ctx.cnt; j++) {
170: PetscReal err = 10.0;
171: if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
172: 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"));
173: pass = (pass && err < ctx.errtol ? PETSC_TRUE : PETSC_FALSE);
174: }
175: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
176: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "This test: %s\n", pass ? "PASSED" : "FAILED"));
178: PetscCall(MatDestroy(&A));
179: PetscCall(TSDestroy(&ts));
180: PetscCall(VecDestroy(&sol));
182: PetscCall(PetscFinalize());
183: return 0;
184: }
186: /*
187: User callback for defining the event-functions
188: */
189: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx)
190: {
191: PetscInt n = 0;
192: AppCtx *Ctx = (AppCtx *)ctx;
194: PetscFunctionBeginUser;
195: // for the test purposes, event-functions are defined based on t
196: // all events -- on rank-0
197: if (Ctx->rank == 0) {
198: if (Ctx->F == 0 || Ctx->F == -1) gval[n++] = PetscSinReal(Ctx->pi * t) / Ctx->pi; // FUNC-0, roots 1, 2, 3, 4
199: if (Ctx->F == 1 || Ctx->F == -1) gval[n++] = PetscLogReal(t); // FUNC-2, root 1
200: if (Ctx->F == 2 || Ctx->F == -1) { // FUNC-3, root 1
201: if (t < 2) gval[n++] = (1 - PetscPowReal(t - 2, 12)) / 12.0;
202: else gval[n++] = 1 / 12.0;
203: }
204: if (Ctx->F == 3 || Ctx->F == -1) gval[n++] = t - PetscExpReal(PetscSinReal(t)) + 1; // FUNC-5, root 1.69681
205: if (Ctx->F == 4 || Ctx->F == -1) gval[n++] = (1e10 * PetscPowReal(t, 1 / t) - 1) / 100; // FUNC-6, root 0.1
206: if (Ctx->F == 5 || Ctx->F == -1) gval[n++] = PetscLogReal(t - 0.02) * PetscLogReal(t - 0.02) * PetscSignReal(t - 1.02) * 1e7; // FUNC-7, root 1.02
207: if (Ctx->F == 6 || Ctx->F == -1) gval[n++] = 4 * PetscCosReal(t) - PetscExpReal(t); // FUNC-14, root 0.904788
208: if (Ctx->F == 7 || Ctx->F == -1) gval[n++] = (20.0 * t - 1) / (19.0 * t) / 10; // FUNC-15, root 0.05
209: if (Ctx->F == 8 || Ctx->F == -1) gval[n++] = ((t - 1) * PetscExpReal(-20 * t) + PetscPowReal(t, 20)) * 1e4; // FUNC-16, root 0.5527
210: if (Ctx->F == 9 || Ctx->F == -1) gval[n++] = (t * t * (t * t / 3.0 + PetscSqrtReal(2.0) * PetscSinReal(t)) - PetscSqrtReal(3.0) / 18) * 10; // FUNC-17, root 0.399
211: if (Ctx->F == 10 || Ctx->F == -1) gval[n++] = ((t * t + 1) * PetscSinReal(t) - PetscExpReal(PetscSqrtReal(t)) * (t - 1) * (t * t - 5)) / 10; // FUNC-18, roots 0.87, 2.388
212: if (Ctx->F == 11 || Ctx->F == -1) gval[n++] = 2 * t - 5; // FUNC-21, root 2.5
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /*
218: User callback for the post-event stuff
219: */
220: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx)
221: {
222: AppCtx *Ctx = (AppCtx *)ctx;
224: PetscFunctionBeginUser;
225: if (Ctx->flg) {
226: PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
227: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
228: for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
229: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
230: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
231: }
233: if (Ctx->cnt + nev_zero < MAX_NEV)
234: for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
236: #ifdef NEW_VERSION
237: Ctx->postcnt++; // sync
238: if (Ctx->dtpost > 0) {
239: if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
240: else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
241: }
242: #endif
244: if (Ctx->restart) PetscCall(TSRestartStep(ts));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
247: /*---------------------------------------------------------------------------------------------*/
248: /*
249: Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
250: which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
251: explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
252: simply remove this option altogether. This will result in using the defaults
253: (which is PETSC_DECIDE).
254: */
255: /*TEST
256: test:
257: suffix: 0
258: requires: !single
259: output_file: output/ex4_0.out
260: args: -dir 0
261: args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-10
262: args: -ts_dt 0.25
263: args: -restart 0
264: args: -ts_event_tol {{1e-8 1e-15}}
265: args: -errtol 1e-7
266: args: -ts_adapt_type {{none basic}}
267: args: -dtpost 0
268: args: -ts_event_post_event_step -1
269: args: -ts_type rk
270: nsize: 2
271: filter: sort
272: filter_output: sort
274: test:
275: suffix: 0single
276: requires: single
277: output_file: output/ex4_0single.out
278: args: -dir 0
279: args: -ts_adapt_dt_min 1e-6 -ts_event_dt_min 1e-6
280: args: -ts_dt 0.3
281: args: -ts_event_tol {{1e-7 1e-10}}
282: args: -ts_adapt_type {{none basic}}
283: args: -dtpost 0.23
284: args: -ts_event_post_event_step -1
285: args: -ts_type beuler
286: nsize: 3
287: filter: sort
288: filter_output: sort
290: test:
291: suffix: F7
292: output_file: output/ex4_F7.out
293: args: -dir 0
294: args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6
295: args: -ts_dt 0.4
296: args: -F 7
297: args: -ts_event_tol {{1e-8 1e-15}}
298: args: -ts_adapt_type {{none basic}}
299: args: -ts_type rk
300: nsize: 1
302: test:
303: suffix: F7revisit
304: output_file: output/ex4_F7revisit.out
305: args: -ts_event_monitor -F 7 -ts_dt 0.04 -ts_event_dt_min 0.016 -errtol 0.005
306: nsize: 1
308: test:
309: suffix: 2all
310: output_file: output/ex4_2.out
311: args: -dir 0
312: args: -F {{-1 0 1 2 3 4 5 6 7 8 9 10 11}}
313: args: -ts_event_dt_min 1e-6 -ts_dt 0.4 -ts_event_tol 1e-8
314: args: -ts_adapt_type {{none basic}}
315: args: -dtpost 0.35
316: args: -ts_type rk
317: filter: grep "This test"
318: nsize: 1
320: test:
321: suffix: 2pos
322: output_file: output/ex4_2.out
323: args: -dir 1
324: args: -F {{-1 0 1 2 3 4 5 7 8 9 10 11}}
325: args: -ts_event_dt_min 1e-6 -ts_dt 0.4 -ts_event_tol 1e-8
326: args: -ts_adapt_type none
327: args: -dtpost 0.34
328: args: -ts_type beuler
329: filter: grep "This test"
330: nsize: 1
332: test:
333: suffix: 2neg
334: output_file: output/ex4_2.out
335: args: -dir -1
336: args: -F {{-1 0 6 10}}
337: args: -ts_event_dt_min 1e-6 -ts_dt 0.4 -ts_event_tol 1e-8
338: args: -ts_adapt_type {{none basic}}
339: args: -dtpost 0.33
340: args: -ts_type rk
341: filter: grep "This test"
342: nsize: 1
344: test:
345: suffix: pos
346: output_file: output/ex4_pos.out
347: args: -dir 1
348: args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6
349: args: -ts_dt 0.4
350: args: -restart 0
351: args: -ts_event_tol {{1e-8 1e-15}}
352: args: -ts_adapt_type {{none basic}}
353: args: -dtpost 0.25
354: args: -ts_event_post_event_step -1
355: args: -ts_type {{beuler rk}}
356: nsize: 1
358: test:
359: suffix: neg
360: output_file: output/ex4_neg.out
361: args: -dir -1
362: args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6
363: args: -ts_dt 0.4
364: args: -restart 1
365: args: -ts_event_tol {{1e-8 1e-15}}
366: args: -ts_adapt_type {{none basic}}
367: args: -dtpost 0.25
368: args: -ts_event_post_event_step 0.35
369: args: -ts_type rk
370: nsize: 2
371: filter: sort
372: filter_output: sort
373: TEST*/