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