1 #include <petscts.h>
2 #include <stdio.h>
3
4 #define NEW_VERSION // Applicable for the new features; avoid this for the older PETSc versions (without TSSetPostEventStep())
5
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";
20
21 #define MAX_NFUNC 100 // max event functions per rank
22 #define MAX_NEV 5000 // max zero crossings for each rank
23
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;
40
41 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx);
42 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx);
43
main(int argc,char ** argv)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;
55
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;
71
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));
85
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));
93
94 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
95 PetscCall(TSSetProblemType(ts, TS_LINEAR));
96
97 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
98 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL));
99
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));
107
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");
117
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 }
129
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));
164
165 // Solution
166 PetscCall(TSSolve(ts, sol));
167
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"));
177
178 PetscCall(MatDestroy(&A));
179 PetscCall(TSDestroy(&ts));
180 PetscCall(VecDestroy(&sol));
181
182 PetscCall(PetscFinalize());
183 return 0;
184 }
185
186 /*
187 User callback for defining the event-functions
188 */
EventFunction(TS ts,PetscReal t,Vec U,PetscReal gval[],PetscCtx ctx)189 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx)
190 {
191 PetscInt n = 0;
192 AppCtx *Ctx = (AppCtx *)ctx;
193
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 }
216
217 /*
218 User callback for the post-event stuff
219 */
Postevent(TS ts,PetscInt nev_zero,PetscInt evs_zero[],PetscReal t,Vec U,PetscBool fwd,PetscCtx ctx)220 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx)
221 {
222 AppCtx *Ctx = (AppCtx *)ctx;
223
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 }
232
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
235
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
243
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_time_step 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
273
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_time_step 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
289
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_time_step 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
301
302 test:
303 suffix: F7revisit
304 output_file: output/ex4_F7revisit.out
305 args: -ts_event_monitor -F 7 -ts_time_step 0.04 -ts_event_dt_min 0.016 -errtol 0.005
306 nsize: 1
307
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_time_step 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
319
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_time_step 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
331
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_time_step 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
343
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_time_step 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
357
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_time_step 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*/
374