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 event function = sin(pi*t), zeros = 1,...,10\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
113 for (PetscInt i = 1; i < MAX_NEV; i++) {
114 if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i;
115 if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i;
116 }
117 }
118 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref));
119 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx));
120
121 // Solution
122 PetscCall(TSSolve(ts, sol));
123
124 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"]
125 for (PetscInt j = 0; j < ctx.cnt; j++) {
126 PetscReal err = 10.0;
127 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]);
128 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"));
129 }
130 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
131
132 PetscCall(MatDestroy(&A));
133 PetscCall(TSDestroy(&ts));
134 PetscCall(VecDestroy(&sol));
135
136 PetscCall(PetscFinalize());
137 return 0;
138 }
139
140 /*
141 User callback for defining the event-functions
142 */
EventFunction(TS ts,PetscReal t,Vec U,PetscReal gval[],PetscCtx ctx)143 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], PetscCtx ctx)
144 {
145 PetscInt n = 0;
146 AppCtx *Ctx = (AppCtx *)ctx;
147
148 PetscFunctionBeginUser;
149 // for the test purposes, event-functions are defined based on t
150 // first event -- on rank-0
151 if (Ctx->rank == 0) gval[n++] = PetscSinReal(Ctx->pi * t);
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
155 /*
156 User callback for the post-event stuff
157 */
Postevent(TS ts,PetscInt nev_zero,PetscInt evs_zero[],PetscReal t,Vec U,PetscBool fwd,PetscCtx ctx)158 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, PetscCtx ctx)
159 {
160 AppCtx *Ctx = (AppCtx *)ctx;
161
162 PetscFunctionBeginUser;
163 if (Ctx->flg) {
164 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx));
165 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero));
166 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]]));
167 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
168 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
169 }
170
171 if (Ctx->cnt + nev_zero < MAX_NEV)
172 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing
173
174 #ifdef NEW_VERSION
175 Ctx->postcnt++; // sync
176 if (Ctx->dtpost > 0) {
177 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost));
178 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE));
179 }
180 #endif
181
182 if (Ctx->restart) PetscCall(TSRestartStep(ts));
183 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185 /*---------------------------------------------------------------------------------------------*/
186 /*
187 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1,
188 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to
189 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed,
190 simply remove this option altogether. This will result in using the defaults
191 (which is PETSC_DECIDE).
192 */
193 /*TEST
194 test:
195 suffix: 0
196 requires: !single
197 output_file: output/ex1sin_0.out
198 args: -dir 0
199 args: -restart 0
200 args: -dtpost {{0 0.25}}
201 args: -ts_event_post_event_step -1
202 args: -ts_type {{beuler rk}}
203 args: -ts_adapt_type {{none basic}}
204 nsize: 1
205
206 test:
207 suffix: p
208 requires: !single
209 output_file: output/ex1sin_p.out
210 args: -dir 1
211 args: -restart 1
212 args: -dtpost {{0 0.31}}
213 args: -ts_event_post_event_step {{-1 0.25}}
214 args: -ts_type rk
215 args: -ts_adapt_type {{none basic}}
216 nsize: 2
217 filter: sort
218 filter_output: sort
219
220 test:
221 suffix: n
222 requires: !single
223 output_file: output/ex1sin_n.out
224 args: -dir -1
225 args: -restart {{0 1}}
226 args: -dtpost {{0 0.25}}
227 args: -ts_event_post_event_step 0.31
228 args: -ts_type {{beuler rk}}
229 args: -ts_adapt_type basic
230 nsize: 4
231 filter: sort
232 filter_output: sort
233
234 test:
235 suffix: 0single
236 requires: single
237 output_file: output/ex1sin_0.out
238 args: -dir 0 -ts_event_dt_min 1e-6 -errtol 5e-5
239 args: -restart 1
240 args: -dtpost {{0 0.25}}
241 args: -ts_event_post_event_step {{-1 0.31}}
242 args: -ts_type beuler
243 args: -ts_adapt_type {{none basic}}
244 nsize: 4
245 filter: sort
246 filter_output: sort
247
248 test:
249 suffix: psingle
250 requires: single
251 output_file: output/ex1sin_p.out
252 args: -dir 1 -ts_event_dt_min 1e-6 -errtol 5e-5
253 args: -restart 0
254 args: -dtpost {{0 0.31}}
255 args: -ts_event_post_event_step 0.25
256 args: -ts_type {{beuler rk}}
257 args: -ts_adapt_type {{none basic}}
258 nsize: 1
259
260 test:
261 suffix: nsingle
262 requires: single
263 output_file: output/ex1sin_n.out
264 args: -dir -1 -ts_event_dt_min 1e-6 -errtol 5e-5
265 args: -restart 1
266 args: -dtpost {{0 0.25}}
267 args: -ts_event_post_event_step -1
268 args: -ts_type {{beuler rk}}
269 args: -ts_adapt_type {{none basic}}
270 nsize: 2
271 filter: sort
272 filter_output: sort
273 TEST*/
274