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
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 V*sin(pi*t), zeros = 1,...,10\n"
13 "After each event location the tolerance for the sin() event is multiplied by 4\n"
14
15 "Options:\n"
16 "-dir d : zero-crossing direction for events: 0, 1, -1\n"
17 "-flg : additional output in Postevent\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\n"
20 "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n"
21 " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n"
22 " if x == 0, nothing happens\n"
23 "-v {float}: scaling of the sin() event function; for small v this event is triggered by the function values,\n"
24 " for large v the event is triggered by the small step size\n"
25 "-change5 : flag to change the state vector at t=5 PostEvent\n";
26
27 #define MAX_NFUNC 100 // max event functions per rank
28 #define MAX_NEV 5000 // max zero crossings for each rank
29
30 typedef struct {
31 PetscMPIInt rank, size;
32 PetscReal pi;
33 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals
34 PetscReal evres[MAX_NEV]; // times of found zero-crossings
35 PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking
36 PetscInt cnt; // counter
37 PetscInt cntref; // actual length of 'ref' on the given rank
38 PetscBool flg; // flag for additional print in PostEvent
39 PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default)
40 PetscBool restart; // flag for TSRestartStep() in PostEvent
41 PetscReal dtpost; // post-event step
42 PetscInt postcnt; // counter for PostEvent calls
43 PetscReal V; // vertical scaling for sin()
44 PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage
45 PetscBool change5; // flag to change the state vector at t=5 PostEvent
46 } AppCtx;
47
48 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx);
49 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx);
50
main(int argc,char ** argv)51 int main(int argc, char **argv)
52 {
53 TS ts;
54 Mat A;
55 Vec sol;
56 PetscInt n, dir0, m = 0;
57 PetscReal tol = 1e-7;
58 PetscInt dir[MAX_NFUNC], inds[2];
59 PetscBool term[MAX_NFUNC];
60 PetscScalar *x, vals[4];
61 AppCtx ctx;
62 TSConvergedReason reason;
63
64 PetscFunctionBeginUser;
65 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
66 setbuf(stdout, NULL);
67 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
68 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
69 ctx.pi = PetscAcosReal(-1.0);
70 ctx.cnt = 0;
71 ctx.cntref = 0;
72 ctx.flg = PETSC_FALSE;
73 ctx.errtol = 1e-5;
74 ctx.restart = PETSC_FALSE;
75 ctx.dtpost = 0;
76 ctx.postcnt = 0;
77 ctx.V = 1.0;
78 ctx.change5 = PETSC_FALSE;
79
80 // The linear problem has a 2*2 matrix. The matrix is constant
81 if (ctx.rank == 0) m = 2;
82 inds[0] = 0;
83 inds[1] = 1;
84 vals[0] = 0;
85 vals[1] = 0.2;
86 vals[2] = -0.2;
87 vals[3] = 0;
88 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A));
89 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES));
90 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
91 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
92 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
93
94 PetscCall(MatCreateVecs(A, &sol, NULL));
95 PetscCall(VecGetArray(sol, &x));
96 if (ctx.rank == 0) { // initial conditions
97 x[0] = 0; // sin(0)
98 x[1] = 1; // cos(0)
99 }
100 PetscCall(VecRestoreArray(sol, &x));
101
102 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
103 PetscCall(TSSetProblemType(ts, TS_LINEAR));
104
105 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
106 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL));
107
108 PetscCall(TSSetTimeStep(ts, 0.1));
109 PetscCall(TSSetType(ts, TSBEULER));
110 PetscCall(TSSetMaxSteps(ts, 10000));
111 PetscCall(TSSetMaxTime(ts, 10.0));
112 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
113 PetscCall(TSSetFromOptions(ts));
114
115 // Set the event handling
116 dir0 = 0;
117 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction
118 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output
119 PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events
120 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep()
121 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step
122 PetscCall(PetscOptionsGetReal(NULL, NULL, "-v", &ctx.V, NULL));
123 PetscCall(PetscOptionsGetBool(NULL, NULL, "-change5", &ctx.change5, NULL)); // flag to change the state vector at t=5 PostEvent
124
125 n = 0; // event counter
126 if (ctx.rank == 0) { // first event -- on rank-0
127 ctx.vtol[n] = tol * 10;
128 dir[n] = dir0;
129 term[n++] = PETSC_FALSE;
130 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.05;
131 }
132 if (ctx.rank == ctx.size - 1) { // second event (with termination) -- on last rank
133 ctx.vtol[n] = tol * 10;
134 dir[n] = dir0;
135 term[n++] = PETSC_TRUE;
136 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 9.05;
137 }
138 if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size
139 ctx.vtol[n] = tol;
140 dir[n] = dir0;
141 term[n++] = PETSC_FALSE;
142
143 for (PetscInt i = 1; i < MAX_NEV - 2; i++) {
144 if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i;
145 if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i;
146 }
147 }
148 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
149 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
150 PetscCall(TSSetEventTolerances(ts, tol, ctx.vtol));
151
152 // Solution
153 PetscCall(TSSolve(ts, sol));
154 PetscCall(TSGetConvergedReason(ts, &reason));
155 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT));
156
157 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
158 for (PetscInt j = 0; j < ctx.cnt; j++) {
159 PetscReal err = 10.0;
160 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
161 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"));
162 }
163 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
164
165 PetscCall(MatDestroy(&A));
166 PetscCall(TSDestroy(&ts));
167 PetscCall(VecDestroy(&sol));
168
169 PetscCall(PetscFinalize());
170 return 0;
171 }
172
173 /*
174 User callback for defining the event-functions
175 */
EventFunction(TS ts,PetscReal t,Vec U,PetscReal gval[],PetscCtx ctx)176 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx)
177 {
178 PetscInt n = 0;
179 AppCtx *Ctx = (AppCtx *)ctx;
180
181 PetscFunctionBeginUser;
182 // for the test purposes, event-functions are defined based on t
183 // first event -- on rank-0
184 if (Ctx->rank == 0) {
185 if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12));
186 else gval[n++] = 0.5;
187 }
188
189 // second event -- on last rank
190 if (Ctx->rank == Ctx->size - 1) {
191 if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12));
192 else gval[n++] = 0.25;
193 }
194
195 // third event -- on rank = 1%ctx.size
196 if (Ctx->rank == 1 % Ctx->size) gval[n++] = Ctx->V * PetscSinReal(Ctx->pi * t);
197 PetscFunctionReturn(PETSC_SUCCESS);
198 }
199
200 /*
201 User callback for the post-event stuff
202 */
Postevent(TS ts,PetscInt nev_zero,PetscInt evs_zero[],PetscReal t,Vec U,PetscBool fwd,PetscCtx ctx)203 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx)
204 {
205 PetscInt n = 0;
206 PetscScalar *x;
207 AppCtx *Ctx = (AppCtx *)ctx;
208
209 PetscFunctionBeginUser;
210 if (Ctx->flg) {
211 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
212 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
213 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
214 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
215 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
216 }
217
218 // change the state vector near t=5.0
219 if (PetscAbsReal(t - (PetscReal)5.0) < 0.01 && Ctx->change5) {
220 PetscCall(VecGetArray(U, &x));
221 if (Ctx->rank == 0) x[1] = -x[1];
222 PetscCall(VecRestoreArray(U, &x));
223 }
224
225 // update vtol's
226 if (Ctx->rank == 0) n++; // first event -- on rank-0
227 if (Ctx->rank == Ctx->size - 1) n++; // second event -- on last rank
228 if (Ctx->rank == 1 % Ctx->size) { // third event -- on rank = 1%ctx.size
229 if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "vtol for sin: %g -> ", (double)Ctx->vtol[n]));
230 Ctx->vtol[n] *= 4;
231 if (PetscAbsReal(t - (PetscReal)5.0) < 0.01) Ctx->vtol[n] /= 100; // one-off decrease
232 if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g\n", (double)Ctx->vtol[n]));
233 n++;
234 }
235 PetscCall(TSSetEventTolerances(ts, 0, Ctx->vtol));
236
237 if (Ctx->cnt + nev_zero < MAX_NEV)
238 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
239
240 #ifdef NEW_VERSION
241 Ctx->postcnt++; // sync
242 if (Ctx->dtpost > 0) {
243 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
244 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
245 }
246 #endif
247
248 if (Ctx->restart) PetscCall(TSRestartStep(ts));
249 PetscFunctionReturn(PETSC_SUCCESS);
250 }
251 /*---------------------------------------------------------------------------------------------*/
252 /*
253 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
254 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
255 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
256 simply remove this option altogether. This will result in using the defaults
257 (which is PETSC_DECIDE).
258 */
259 /*TEST
260 test:
261 suffix: V
262 output_file: output/ex3_V.out
263 args: -ts_type beuler
264 args: -ts_adapt_type basic
265 args: -v {{1e2 1e5 1e8}}
266 args: -ts_adapt_dt_min 1e-6
267 args: -change5 {{0 1}}
268 nsize: 1
269
270 test:
271 suffix: neu1
272 output_file: output/ex3_neu1.out
273 args: -dir 0
274 args: -v 1e5
275 args: -ts_adapt_dt_min 1e-6
276 args: -restart 1
277 args: -dtpost 0.24
278 args: -ts_event_post_event_step 0.31
279 args: -ts_type {{beuler rk}}
280 args: -ts_adapt_type {{none basic}}
281 nsize: 1
282
283 test:
284 suffix: neu2
285 output_file: output/ex3_neu2.out
286 args: -dir 0
287 args: -v 1e5
288 args: -ts_adapt_dt_min 1e-6
289 args: -restart 1
290 args: -dtpost 0
291 args: -ts_event_post_event_step {{-1 0.31}}
292 args: -ts_type rk
293 args: -ts_adapt_type {{none basic}}
294 nsize: 2
295 filter: sort
296 filter_output: sort
297
298 test:
299 suffix: neu4
300 output_file: output/ex3_neu4.out
301 args: -dir 0
302 args: -v 1e5
303 args: -ts_adapt_dt_min 1e-6
304 args: -restart {{0 1}}
305 args: -dtpost 0.24
306 args: -ts_event_post_event_step 0.21
307 args: -ts_type beuler
308 args: -ts_adapt_type {{none basic}}
309 nsize: 4
310 filter: sort
311 filter_output: sort
312
313 test:
314 suffix: pos1
315 output_file: output/ex3_pos1.out
316 args: -dir 1
317 args: -v 1e5
318 args: -ts_adapt_dt_min 1e-6
319 args: -restart 0
320 args: -dtpost 0.24
321 args: -ts_type {{beuler rk}}
322 args: -ts_adapt_type {{none basic}}
323 nsize: 1
324
325 test:
326 suffix: pos2
327 output_file: output/ex3_pos2.out
328 args: -dir 1
329 args: -v 1e5
330 args: -ts_adapt_dt_min 1e-6
331 args: -restart 1
332 args: -dtpost {{0 0.24}}
333 args: -ts_type rk
334 args: -ts_adapt_type {{none basic}}
335 nsize: 2
336 filter: sort
337 filter_output: sort
338
339 test:
340 suffix: pos4
341 output_file: output/ex3_pos4.out
342 args: -dir 1
343 args: -v 1e9
344 args: -ts_adapt_dt_min 1e-6
345 args: -restart 0
346 args: -dtpost 0
347 args: -ts_event_post_event_step {{-1 0.32}}
348 args: -ts_type beuler
349 args: -ts_adapt_type {{none basic}}
350 args: -change5 1
351 nsize: 4
352 filter: sort
353 filter_output: sort
354
355 test:
356 suffix: neg1
357 output_file: output/ex3_neg1.out
358 args: -dir -1
359 args: -v 1e5
360 args: -ts_adapt_dt_min 1e-6
361 args: -restart 1
362 args: -dtpost {{0 0.24}}
363 args: -ts_type {{beuler rk}}
364 args: -ts_adapt_type basic
365 nsize: 1
366
367 test:
368 suffix: neg2
369 output_file: output/ex3_neg2.out
370 args: -dir -1
371 args: -v 1e5
372 args: -ts_adapt_dt_min 1e-6
373 args: -restart 0
374 args: -dtpost {{0 0.24}}
375 args: -ts_type rk
376 args: -ts_adapt_type {{none basic}}
377 nsize: 2
378 filter: sort
379 filter_output: sort
380
381 test:
382 suffix: neg4
383 output_file: output/ex3_neg4.out
384 args: -dir -1
385 args: -v 1e5
386 args: -ts_adapt_dt_min 1e-6
387 args: -restart 0
388 args: -dtpost {{0 0.24}}
389 args: -ts_event_post_event_step 0.3
390 args: -ts_type beuler
391 args: -ts_adapt_type {{none basic}}
392 nsize: 4
393 filter: sort
394 filter_output: sort
395 TEST*/
396