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