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 16 event functions:\n"
10 "7 polynomials (dir=+1) with zeros: 1+2^i, i=-3,...3, on ranks=(i+3)%size\n"
11 "7 polynomials (dir=-1) with zeros: 1+(8-2^i), i=-3,...3, on ranks=(i+3)%size\n"
12 "(t-5)^2 * sin(pi*t), with zeros = 1,2,...10, on rank-0\n"
13 " 0.5 * cos(pi*t), with zeros = 0.5,1.5,...9.5, on last-rank\n"
14 "Options:\n"
15 "-dir d : zero-crossing direction for events\n"
16 "-flg : additional output in Postevent\n"
17 "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n"
18 "-restart : flag for TSRestartStep() in PostEvent\n"
19 "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n"
20 " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n"
21 " if x == 0, nothing happens\n";
22
23 #define MAX_NFUNC 100 // max event functions per rank
24 #define MAX_NEV 5000 // max zero crossings for each rank
25
26 typedef struct {
27 PetscMPIInt rank, size;
28 PetscReal pi;
29 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals
30 PetscReal evres[MAX_NEV]; // times of found zero-crossings
31 PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking
32 PetscInt cnt; // counter
33 PetscInt cntref; // actual length of 'ref' on the given rank
34 PetscBool flg; // flag for additional print in PostEvent
35 PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default)
36 PetscBool restart; // flag for TSRestartStep() in PostEvent
37 PetscReal dtpost; // post-event step
38 PetscInt postcnt; // counter for PostEvent calls
39 PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage
40 PetscInt dir0; // desired zero-crossing direction
41 } AppCtx;
42
43 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx);
44 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx);
45 static inline void SetVtols(PetscMPIInt rank, PetscMPIInt size, PetscReal tol0, PetscReal tolsin, PetscReal *vtol); // helper function to fill vtol[]
46
main(int argc,char ** argv)47 int main(int argc, char **argv)
48 {
49 TS ts;
50 Mat A;
51 Vec sol;
52 PetscInt n, m = 0;
53 PetscInt dir[MAX_NFUNC], inds[2];
54 PetscBool term[MAX_NFUNC];
55 PetscScalar *x, vals[4];
56 PetscReal aux;
57 AppCtx ctx;
58
59 PetscFunctionBeginUser;
60 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
61 setbuf(stdout, NULL);
62 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
63 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
64 ctx.pi = PetscAcosReal(-1.0);
65 ctx.cnt = 0;
66 ctx.cntref = 0;
67 ctx.flg = PETSC_FALSE;
68 ctx.errtol = 1e-5;
69 ctx.restart = PETSC_FALSE;
70 ctx.dtpost = 0;
71 ctx.postcnt = 0;
72
73 // The linear problem has a 2*2 matrix. The matrix is constant
74 if (ctx.rank == 0) m = 2;
75 inds[0] = 0;
76 inds[1] = 1;
77 vals[0] = 0;
78 vals[1] = 0.2;
79 vals[2] = -0.2;
80 vals[3] = 0;
81 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A));
82 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES));
83 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
84 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
85 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
86
87 PetscCall(MatCreateVecs(A, &sol, NULL));
88 PetscCall(VecGetArray(sol, &x));
89 if (ctx.rank == 0) { // initial conditions
90 x[0] = 0; // sin(0)
91 x[1] = 1; // cos(0)
92 }
93 PetscCall(VecRestoreArray(sol, &x));
94
95 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
96 PetscCall(TSSetProblemType(ts, TS_LINEAR));
97
98 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL));
99 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL));
100
101 PetscCall(TSSetTimeStep(ts, 0.1));
102 PetscCall(TSSetType(ts, TSBEULER));
103 PetscCall(TSSetMaxSteps(ts, 10000));
104 PetscCall(TSSetMaxTime(ts, 10.0));
105 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
106 PetscCall(TSSetFromOptions(ts));
107
108 // Set the event handling
109 ctx.dir0 = 0;
110 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &ctx.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
116 n = 0; // event counter
117 aux = 1.0 / 8.0;
118 for (PetscInt i = -3; i <= 3; i++) { // pos-polynomials
119 if (ctx.rank == (i + 3) % ctx.size) {
120 dir[n] = ctx.dir0;
121 term[n++] = PETSC_FALSE;
122 if (ctx.dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0 + aux;
123 }
124 aux *= 2;
125 }
126 aux = 1.0 / 8.0;
127 for (PetscInt i = -3; i <= 3; i++) { // neg-polynomials
128 if (ctx.rank == (i + 3) % ctx.size) {
129 dir[n] = ctx.dir0;
130 term[n++] = PETSC_FALSE;
131 if (ctx.dir0 <= 0) ctx.ref[ctx.cntref++] = 9.0 - aux;
132 }
133 aux *= 2;
134 }
135 if (ctx.rank == 0) { // sin-event -- on rank-0
136 dir[n] = ctx.dir0;
137 term[n++] = PETSC_FALSE;
138 for (PetscInt i = 1; i < MAX_NEV / 2 - 10; i++) {
139 if (i % 2 == 1 && ctx.dir0 <= 0) ctx.ref[ctx.cntref++] = i;
140 if (i % 2 == 0 && ctx.dir0 >= 0) ctx.ref[ctx.cntref++] = i;
141 }
142 }
143 if (ctx.rank == ctx.size - 1) { // cos-event -- on last rank
144 dir[n] = ctx.dir0;
145 term[n++] = PETSC_FALSE;
146 for (PetscInt i = 1; i < MAX_NEV / 2 - 10; i++) {
147 if (i % 2 == 1 && ctx.dir0 <= 0) ctx.ref[ctx.cntref++] = i - 0.5;
148 if (i % 2 == 0 && ctx.dir0 >= 0) ctx.ref[ctx.cntref++] = i - 0.5;
149 }
150 }
151 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
152 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
153 SetVtols(ctx.rank, ctx.size, 1e-8, 1e-8, ctx.vtol);
154 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, ctx.vtol));
155
156 // Solution
157 PetscCall(TSSolve(ts, sol));
158
159 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
160 for (PetscInt j = 0; j < ctx.cnt; j++) {
161 PetscReal err = 10.0;
162 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
163 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"));
164 }
165 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
166
167 PetscCall(MatDestroy(&A));
168 PetscCall(TSDestroy(&ts));
169 PetscCall(VecDestroy(&sol));
170
171 PetscCall(PetscFinalize());
172 return 0;
173 }
174
175 /*
176 User callback for defining the event-functions
177 */
EventFunction(TS ts,PetscReal t,Vec U,PetscReal gval[],PetscCtx ctx)178 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx)
179 {
180 PetscInt n = 0;
181 AppCtx *Ctx = (AppCtx *)ctx;
182 PetscReal P;
183
184 PetscFunctionBeginUser;
185 // for the test purposes, event-functions are defined based on t
186 for (PetscInt i = -3; i <= 3; i++) { // pos-polynomials
187 if (Ctx->rank == (i + 3) % Ctx->size) {
188 P = PetscPowReal(2.0, i);
189 if (t < 2 + P) gval[n++] = 1 - PetscPowReal(2 + P - t, i + 5);
190 else gval[n++] = 1;
191 }
192 }
193 for (PetscInt i = -3; i <= 3; i++) { // neg-polynomials
194 if (Ctx->rank == (i + 3) % Ctx->size) {
195 P = PetscPowReal(2.0, i);
196 if (t > 8 - P) gval[n++] = 1 - PetscPowReal(t - 8 + P, i + 5);
197 else gval[n++] = 1;
198 }
199 }
200 if (Ctx->rank == 0) gval[n++] = (t - 5) * (t - 5) * PetscSinReal(Ctx->pi * t); // sin-event -- on rank-0
201 if (Ctx->rank == Ctx->size - 1) gval[n++] = 0.5 * PetscCosReal(Ctx->pi * t); // cos-event -- on last rank
202 PetscFunctionReturn(PETSC_SUCCESS);
203 }
204
205 /*
206 User callback for the post-event stuff
207 */
Postevent(TS ts,PetscInt nev_zero,PetscInt evs_zero[],PetscReal t,Vec U,PetscBool fwd,PetscCtx ctx)208 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx)
209 {
210 AppCtx *Ctx = (AppCtx *)ctx;
211
212 PetscFunctionBeginUser;
213 if (Ctx->flg) {
214 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
215 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
216 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
217 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
218 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
219 }
220
221 if (Ctx->cnt + nev_zero < MAX_NEV)
222 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
223
224 #ifdef NEW_VERSION
225 Ctx->postcnt++; // sync
226 if (Ctx->dtpost > 0) {
227 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
228 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
229 }
230 #endif
231
232 if ((Ctx->dir0 == 0 && PetscAbsReal(t - (PetscReal)4.0) < 0.01) || (Ctx->dir0 == -1 && PetscAbsReal(t - (PetscReal)3.0) < 0.01)) {
233 SetVtols(Ctx->rank, Ctx->size, 1e-8, 1e-26, Ctx->vtol); // for better resolution of sin-event at t=5.0
234 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, Ctx->vtol));
235 }
236 if (PetscAbsReal(t - (PetscReal)5.0) < 0.01) {
237 SetVtols(Ctx->rank, Ctx->size, 1e-8, 1e-8, Ctx->vtol); // back to normal
238 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, Ctx->vtol));
239 }
240
241 if (Ctx->restart) PetscCall(TSRestartStep(ts));
242 PetscFunctionReturn(PETSC_SUCCESS);
243 }
244
245 // helper function to fill vtol[]
SetVtols(PetscMPIInt rank,PetscMPIInt size,PetscReal tol0,PetscReal tolsin,PetscReal * vtol)246 static inline void SetVtols(PetscMPIInt rank, PetscMPIInt size, PetscReal tol0, PetscReal tolsin, PetscReal *vtol)
247 {
248 PetscInt n = 0;
249 for (PetscInt i = -3; i <= 3; i++)
250 if (rank == (i + 3) % size) vtol[n++] = tol0; // pos-polynomials
251 for (PetscInt i = -3; i <= 3; i++)
252 if (rank == (i + 3) % size) vtol[n++] = tol0; // neg-polynomials
253 if (rank == 0) vtol[n++] = tolsin; // sin-event -- on rank-0
254 if (rank == size - 1) vtol[n++] = tol0; // cos-event -- on last rank
255 }
256 /*---------------------------------------------------------------------------------------------*/
257 /*
258 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
259 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
260 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
261 simply remove this option altogether. This will result in using the defaults
262 (which is PETSC_DECIDE).
263 */
264 /*TEST
265 test:
266 suffix: pos1
267 output_file: output/ex5_pos1.out
268 args: -dir 1 -ts_event_dt_min 1e-6
269 args: -restart 1
270 args: -dtpost {{0 0.25}}
271 args: -ts_event_post_event_step 0.31
272 args: -ts_type rk
273 args: -ts_adapt_type {{none basic}}
274 nsize: 1
275
276 test:
277 suffix: pos4
278 output_file: output/ex5_pos4.out
279 args: -dir 1 -ts_event_dt_min 1e-6 -ts_time_step 0.25
280 args: -restart 0
281 args: -dtpost 0
282 args: -ts_event_post_event_step -1
283 args: -ts_type {{beuler rk}}
284 args: -ts_adapt_type {{none basic}}
285 nsize: 4
286 filter: sort
287 filter_output: sort
288
289 test:
290 suffix: neu1
291 output_file: output/ex5_neu1.out
292 args: -dir 0 -ts_event_dt_min 1e-6
293 args: -restart 0
294 args: -dtpost {{0 0.25}}
295 args: -ts_event_post_event_step -1
296 args: -ts_type rk
297 args: -ts_adapt_type {{none basic}}
298 nsize: 1
299
300 test:
301 suffix: neu4
302 output_file: output/ex5_neu4.out
303 args: -dir 0 -ts_event_dt_min 1e-6 -ts_time_step 0.25
304 args: -dtpost 0
305 args: -ts_event_post_event_step {{-1 0.29}}
306 args: -ts_event_post_event_second_step {{-1 0.31}}
307 args: -ts_type rk
308 args: -ts_adapt_type {{none basic}}
309 nsize: 4
310 filter: sort
311 filter_output: sort
312
313 test:
314 suffix: neg2
315 output_file: output/ex5_neg2.out
316 args: -dir -1 -ts_event_dt_min 1e-6
317 args: -restart 1
318 args: -dtpost {{0 0.25}}
319 args: -ts_event_post_event_step 0.31
320 args: -ts_type beuler
321 args: -ts_adapt_type {{none basic}}
322 nsize: 2
323 filter: sort
324 filter_output: sort
325
326 test:
327 suffix: neg4
328 output_file: output/ex5_neg4.out
329 args: -dir -1 -ts_event_dt_min 1e-6 -ts_time_step 0.25
330 args: -restart 0
331 args: -dtpost 0
332 args: -ts_event_post_event_step -1
333 args: -ts_type {{beuler rk}}
334 args: -ts_adapt_type {{none basic}}
335 nsize: 4
336 filter: sort
337 filter_output: sort
338
339 TEST*/
340