Actual source code: ex3span.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 sin(pi*t), zeros = 1,...,10\n"
13: "TimeSpan = [0.01, 0.21, 1.01, ..., 6.21, 6.99, 7.21,... 9.21] plus the points: {3, 4, 4+D, 5-D, 5, 6-D, 6, 6+D} with user-defined D\n"
15: "Options:\n"
16: "-dir d : zero-crossing direction for events: 0 (default), 1, -1\n"
17: "-flg : additional output in Postevent (default: nothing)\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 (default: no)\n"
20: "-term : flag to terminate at 9.05 event (true by default)\n"
21: "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n"
22: " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n"
23: " if x == 0, nothing happens (default)\n"
24: "-D z : a small real number to define additional TimeSpan points (default = 0.02)\n"
25: "-dt2_at6 t : second time step set after event at t=6 (if nothing is specified, no action is done)\n"
26: "-mult7 m : after event at t=7, the linear system coeffs '0.2' are multiplied by m (default = 1.0)\n";
28: #define MAX_NFUNC 100 // max event functions per rank
29: #define MAX_NEV 5000 // max zero crossings for each rank
31: typedef struct {
32: PetscMPIInt rank, size;
33: PetscReal pi;
34: PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals
35: PetscReal evres[MAX_NEV]; // times of found zero-crossings
36: PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking
37: PetscInt cnt; // counter
38: PetscInt cntref; // actual length of 'ref' on the given rank
39: PetscBool flg; // flag for additional print in PostEvent
40: PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default)
41: PetscBool restart; // flag for TSRestartStep() in PostEvent
42: PetscBool term; // flag to terminate at 9.05 event
43: PetscReal dtpost; // first post-event step
44: PetscReal dt2_at6; // second time step set after event at t=6
45: PetscReal mult7; // multiplier for coeffs at t=7
46: PetscInt postcnt; // counter for PostEvent calls
47: Mat A; // system matrix
48: PetscInt m; // local size of A
49: } AppCtx;
51: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx);
52: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx);
53: PetscErrorCode Fill_mat(PetscReal coeff, PetscInt m, Mat A); // Fills the system matrix (2*2)
55: int main(int argc, char **argv)
56: {
57: TS ts;
58: Vec sol;
59: PetscInt n, dir0;
60: PetscReal tol = 1e-7, D = 0.02;
61: PetscInt dir[MAX_NFUNC];
62: PetscBool term[MAX_NFUNC], match;
63: PetscScalar *x;
64: PetscReal tspan[28], dtlast, tlast, tlast_expected, maxtime;
65: AppCtx ctx;
66: TSConvergedReason reason;
67: TSAdapt adapt;
69: PetscFunctionBeginUser;
70: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
71: setbuf(stdout, NULL);
72: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
73: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
74: ctx.pi = PetscAcosReal(-1.0);
75: ctx.cnt = 0;
76: ctx.cntref = 0;
77: ctx.flg = PETSC_FALSE;
78: ctx.errtol = 1e-5;
79: ctx.restart = PETSC_FALSE;
80: ctx.term = PETSC_TRUE;
81: ctx.dtpost = 0;
82: ctx.dt2_at6 = -2;
83: ctx.mult7 = 1.0;
84: ctx.postcnt = 0;
85: ctx.m = 0;
87: // The linear problem has a 2*2 matrix. The matrix is constant
88: if (ctx.rank == 0) ctx.m = 2;
89: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, ctx.m, ctx.m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &ctx.A));
90: PetscCallBack("Fill_mat", Fill_mat(0.2, ctx.m, ctx.A));
91: PetscCall(MatCreateVecs(ctx.A, &sol, NULL));
92: PetscCall(VecGetArray(sol, &x));
93: if (ctx.rank == 0) { // initial conditions
94: x[0] = 0; // sin(0)
95: x[1] = 1; // cos(0)
96: }
97: PetscCall(VecRestoreArray(sol, &x));
99: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
100: PetscCall(TSSetProblemType(ts, TS_LINEAR));
102: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
103: PetscCall(TSSetRHSJacobian(ts, ctx.A, ctx.A, TSComputeRHSJacobianConstant, NULL));
105: PetscCall(TSSetTimeStep(ts, 0.099));
106: PetscCall(TSSetType(ts, TSBEULER));
107: PetscCall(TSSetMaxSteps(ts, 10000));
108: PetscCall(TSSetMaxTime(ts, 10.0));
109: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
111: // Set the event handling
112: dir0 = 0;
113: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction
114: PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output
115: PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events
116: PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep()
117: PetscCall(PetscOptionsGetBool(NULL, NULL, "-term", &ctx.term, NULL)); // flag to terminate at 9.05 event
118: PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step
119: PetscCall(PetscOptionsGetReal(NULL, NULL, "-dt2_at6", &ctx.dt2_at6, NULL)); // second time step set after event at t=6
120: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mult7", &ctx.mult7, NULL)); // multiplier for coeffs at t=7
121: PetscCall(PetscOptionsGetReal(NULL, NULL, "-D", &D, NULL)); // small number for tspan
123: n = 0; // event counter
124: if (ctx.rank == 0) { // first event -- on rank-0
125: dir[n] = dir0;
126: term[n++] = PETSC_FALSE;
127: if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.05;
128: }
129: if (ctx.rank == ctx.size - 1) { // second event (with optional termination) -- on last rank
130: dir[n] = dir0;
131: term[n++] = ctx.term;
132: if (dir0 <= 0) ctx.ref[ctx.cntref++] = 9.05;
133: }
134: if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size
135: dir[n] = dir0;
136: term[n++] = PETSC_FALSE;
138: for (PetscInt i = 1; i < MAX_NEV - 2; i++) {
139: if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i;
140: if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i;
141: }
142: }
143: if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
144: PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
145: PetscCall(TSSetEventTolerances(ts, tol, NULL));
147: // Set the time span
148: for (PetscInt i = 0; i < 10; i++) {
149: tspan[2 * i] = 0.01 + i + (i == 7 ? -0.02 : 0);
150: tspan[2 * i + 1] = 0.21 + i;
151: }
152: tspan[20] = 3;
153: tspan[21] = 4;
154: tspan[22] = 4 + D;
155: tspan[23] = 5 - D;
156: tspan[24] = 5;
157: tspan[25] = 6 - D;
158: tspan[26] = 6;
159: tspan[27] = 6 + D;
160: PetscCall(PetscSortReal(28, tspan));
161: PetscCall(TSSetTimeSpan(ts, 28, tspan));
162: PetscCall(TSSetFromOptions(ts));
164: // Solution
165: PetscCall(TSSolve(ts, sol));
166: PetscCall(TSGetConvergedReason(ts, &reason));
167: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT));
169: // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
170: for (PetscInt j = 0; j < ctx.cnt; j++) {
171: PetscReal err = 10.0;
172: if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
173: 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"));
174: }
175: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
177: // print the final time and step
178: PetscCall(TSGetTime(ts, &tlast));
179: PetscCall(TSGetTimeStep(ts, &dtlast));
180: PetscCall(TSGetAdapt(ts, &adapt));
181: PetscCall(PetscObjectTypeCompare((PetscObject)adapt, TSADAPTNONE, &match));
183: PetscCall(TSGetMaxTime(ts, &maxtime));
184: tlast_expected = ((dir0 == 1 || !ctx.term) ? maxtime : PetscMin(maxtime, 9.05));
185: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final time = %g, max time = %g, %s\n", (double)tlast, (double)maxtime, PetscAbsReal(tlast - tlast_expected) < ctx.errtol ? "pass" : "fail"));
187: if (match) {
188: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Adapt = none\n"));
189: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Last dt = %g\n", (double)dtlast));
190: }
192: PetscCall(MatDestroy(&ctx.A));
193: PetscCall(TSDestroy(&ts));
194: PetscCall(VecDestroy(&sol));
196: PetscCall(PetscFinalize());
197: return 0;
198: }
200: /*
201: User callback for defining the event-functions
202: */
203: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx)
204: {
205: PetscInt n = 0;
206: AppCtx *Ctx = (AppCtx *)ctx;
208: PetscFunctionBeginUser;
209: // for the test purposes, event-functions are defined based on t
210: // first event -- on rank-0
211: if (Ctx->rank == 0) {
212: if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12));
213: else gval[n++] = 0.5;
214: }
216: // second event -- on last rank
217: if (Ctx->rank == Ctx->size - 1) {
218: if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12));
219: else gval[n++] = 0.25;
220: }
222: // third event -- on rank = 1%ctx.size
223: if (Ctx->rank == 1 % Ctx->size) { gval[n++] = PetscSinReal(Ctx->pi * t); }
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: /*
228: User callback for the post-event stuff
229: */
230: PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx)
231: {
232: AppCtx *Ctx = (AppCtx *)ctx;
233: PetscBool mat_changed = PETSC_FALSE;
235: PetscFunctionBeginUser;
236: if (Ctx->flg) {
237: PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
238: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
239: for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
240: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
241: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
242: }
244: if (Ctx->cnt + nev_zero < MAX_NEV)
245: for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
247: #ifdef NEW_VERSION
248: Ctx->postcnt++; // sync
249: if (Ctx->dtpost > 0) {
250: if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
251: else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
252: }
253: #endif
255: // t==6: set the second post-event step
256: if (PetscAbsReal(t - (PetscReal)6.0) < 0.01 && Ctx->dt2_at6 != -2) PetscCall(TSSetPostEventSecondStep(ts, Ctx->dt2_at6));
258: // t==7: change the system matrix
259: if (PetscAbsReal(t - 7) < 0.01 && Ctx->mult7 != 1) {
260: PetscCallBack("Fill_mat", Fill_mat(0.2 * Ctx->mult7, Ctx->m, Ctx->A));
261: PetscCall(TSSetRHSJacobian(ts, Ctx->A, Ctx->A, TSComputeRHSJacobianConstant, NULL));
262: mat_changed = PETSC_TRUE;
263: }
265: if (Ctx->restart || mat_changed) PetscCall(TSRestartStep(ts));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*
270: Fills the system matrix (2*2)
271: */
272: PetscErrorCode Fill_mat(PetscReal coeff, PetscInt m, Mat A)
273: {
274: PetscInt inds[2];
275: PetscScalar vals[4];
277: PetscFunctionBeginUser;
278: inds[0] = 0;
279: inds[1] = 1;
280: vals[0] = 0;
281: vals[1] = coeff;
282: vals[2] = -coeff;
283: vals[3] = 0;
284: PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES));
285: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
286: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
287: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
290: /*---------------------------------------------------------------------------------------------*/
291: /*
292: Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
293: which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
294: explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
295: simply remove this option altogether. This will result in using the defaults
296: (which is PETSC_DECIDE).
297: */
298: /*TEST
299: test:
300: suffix: 1
301: requires: !single
302: output_file: output/ex3span_1.out
303: args: -ts_monitor -ts_adapt_type none -restart
304: args: -dtpost 0.1127 -D 0.0015 -dir 0 -ts_max_time 9.8 -ts_dt 0.18
305: nsize: 1
307: test:
308: suffix: 1single
309: requires: single
310: output_file: output/ex3span_1single.out
311: args: -ts_monitor -ts_adapt_type none -restart -ts_event_dt_min 1e-6
312: args: -dtpost 0.1127 -D 0.0015 -dir 0 -ts_max_time 9.8 -ts_dt 0.18
313: nsize: 1
315: test:
316: suffix: 2
317: output_file: output/ex3span_2.out
318: args: -ts_event_dt_min 1e-6 -dtpost 1 -term 0 -ts_max_time 9.61
319: nsize: 1
321: test:
322: suffix: 3none
323: output_file: output/ex3span_3none.out
324: args: -ts_event_dt_min 1e-6 -ts_adapt_type none -dir 0
325: args: -ts_event_post_event_step {{-1 0.11}}
326: args: -ts_event_post_event_second_step 0.12
327: args: -dt2_at6 {{-2 0.08 0.15}}
328: nsize: 3
330: test:
331: suffix: 3basic
332: output_file: output/ex3span_3basic.out
333: args: -ts_event_dt_min 1e-6 -ts_adapt_type basic -dir 0
334: args: -ts_event_post_event_step {{-1 0.11}}
335: args: -ts_event_post_event_second_step 0.12
336: args: -dt2_at6 {{-2 0.08 0.15}}
337: args: -mult7 {{1 2}}
338: nsize: 2
340: test:
341: suffix: fin
342: output_file: output/ex3span_fin.out
343: args: -ts_max_time {{8.21 8.99 9 9.04 9.05 9.06 9.21 9.99 12}}
344: args: -ts_event_dt_min 1e-6
345: args: -ts_adapt_type {{none basic}}
346: args: -dtpost 0.1125
347: args: -D 0.0025
348: args: -dir {{0 -1 1}}
349: args: -ts_dt 0.3025
350: args: -ts_type {{rk bdf}}
351: filter: grep "Final time ="
352: filter_output: grep "Final time ="
353: nsize: 2
355: test:
356: suffix: adaptmonitor
357: requires: !single
358: output_file: output/ex3span_adaptmonitor.out
359: args: -ts_adapt_monitor -dir 1
360: nsize: 1
361: TEST*/