1 static char help[] = "Serial bouncing ball example to test TS event feature.\n";
2
3 /*
4 The dynamics of the bouncing ball is described by the ODE
5 u1_t = u2
6 u2_t = -9.8
7
8 There is one event set in this example, which checks for the ball hitting the
9 ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
10 a factor of 0.9. On reaching the limit on the number of ball bounces,
11 the TS run is requested to terminate from the PostEvent() callback.
12 */
13
14 #include <petscts.h>
15
16 typedef struct {
17 PetscInt maxbounces;
18 PetscInt nbounces;
19 } AppCtx;
20
EventFunction(TS ts,PetscReal t,Vec U,PetscReal * fvalue,PetscCtx ctx)21 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, PetscCtx ctx)
22 {
23 const PetscScalar *u;
24
25 PetscFunctionBeginUser;
26 /* Event for ball height */
27 PetscCall(VecGetArrayRead(U, &u));
28 fvalue[0] = PetscRealPart(u[0]);
29 PetscCall(VecRestoreArrayRead(U, &u));
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,PetscCtx ctx)33 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx)
34 {
35 AppCtx *app = (AppCtx *)ctx;
36 PetscScalar *u;
37
38 PetscFunctionBeginUser;
39 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds\n", (double)t));
40 /* Set new initial conditions with .9 attenuation */
41 PetscCall(VecGetArray(U, &u));
42 u[0] = 0.0;
43 u[1] = -0.9 * u[1];
44 PetscCall(VecRestoreArray(U, &u));
45 app->nbounces++;
46 if (app->nbounces >= app->maxbounces) {
47 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball bounced %" PetscInt_FMT " times\n", app->nbounces));
48 PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER)); // request TS to terminate; since the program is serial, no need to sync this call
49 }
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52
53 /*
54 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
55 */
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,PetscCtx ctx)56 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx)
57 {
58 PetscScalar *f;
59 const PetscScalar *u;
60
61 PetscFunctionBeginUser;
62 /* The following lines allow us to access the entries of the vectors directly */
63 PetscCall(VecGetArrayRead(U, &u));
64 PetscCall(VecGetArray(F, &f));
65
66 f[0] = u[1];
67 f[1] = -9.8;
68
69 PetscCall(VecRestoreArrayRead(U, &u));
70 PetscCall(VecRestoreArray(F, &f));
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
74 /*
75 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
76 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,PetscCtx ctx)77 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, PetscCtx ctx)
78 {
79 PetscInt rowcol[] = {0, 1};
80 PetscScalar J[2][2];
81 const PetscScalar *u;
82
83 PetscFunctionBeginUser;
84 PetscCall(VecGetArrayRead(U, &u));
85
86 J[0][0] = 0.0;
87 J[0][1] = 1.0;
88 J[1][0] = 0.0;
89 J[1][1] = 0.0;
90 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
91
92 PetscCall(VecRestoreArrayRead(U, &u));
93
94 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
95 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
96 if (A != B) {
97 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
98 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
99 }
100 PetscFunctionReturn(PETSC_SUCCESS);
101 }
102
103 /*
104 Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
105 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)106 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
107 {
108 PetscScalar *f;
109 const PetscScalar *u, *udot;
110
111 PetscFunctionBeginUser;
112 /* The next three lines allow us to access the entries of the vectors directly */
113 PetscCall(VecGetArrayRead(U, &u));
114 PetscCall(VecGetArrayRead(Udot, &udot));
115 PetscCall(VecGetArray(F, &f));
116
117 f[0] = udot[0] - u[1];
118 f[1] = udot[1] + 9.8;
119
120 PetscCall(VecRestoreArrayRead(U, &u));
121 PetscCall(VecRestoreArrayRead(Udot, &udot));
122 PetscCall(VecRestoreArray(F, &f));
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
126 /*
127 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
128 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,PetscCtx ctx)129 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
130 {
131 PetscInt rowcol[] = {0, 1};
132 PetscScalar J[2][2];
133 const PetscScalar *u, *udot;
134
135 PetscFunctionBeginUser;
136 PetscCall(VecGetArrayRead(U, &u));
137 PetscCall(VecGetArrayRead(Udot, &udot));
138
139 J[0][0] = a;
140 J[0][1] = -1.0;
141 J[1][0] = 0.0;
142 J[1][1] = a;
143 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
144
145 PetscCall(VecRestoreArrayRead(U, &u));
146 PetscCall(VecRestoreArrayRead(Udot, &udot));
147
148 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
149 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
150 if (A != B) {
151 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
152 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
153 }
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
main(int argc,char ** argv)157 int main(int argc, char **argv)
158 {
159 TS ts; /* ODE integrator */
160 Vec U; /* solution will be stored here */
161 PetscMPIInt size;
162 PetscInt n = 2;
163 PetscScalar *u;
164 AppCtx app;
165 PetscInt direction[1];
166 PetscBool terminate[1];
167 PetscBool rhs_form = PETSC_FALSE, hist = PETSC_TRUE;
168 TSAdapt adapt;
169
170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 Initialize program
172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173 PetscFunctionBeginUser;
174 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
175 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
176 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
177
178 app.nbounces = 0;
179 app.maxbounces = 10;
180 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex40 options", "");
181 PetscCall(PetscOptionsInt("-maxbounces", "", "", app.maxbounces, &app.maxbounces, NULL));
182 PetscCall(PetscOptionsBool("-test_adapthistory", "", "", hist, &hist, NULL));
183 PetscOptionsEnd();
184
185 Mat A; /* Jacobian matrix */
186 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
187 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
188 PetscCall(MatSetType(A, MATDENSE));
189 PetscCall(MatSetFromOptions(A));
190 PetscCall(MatSetUp(A));
191 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192 Create timestepping solver context
193 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
195 PetscCall(TSSetType(ts, TSROSW));
196
197 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198 Set ODE routines
199 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
201 /* Users are advised against the following branching and code duplication.
202 For problems without a mass matrix like the one at hand, the RHSFunction
203 (and companion RHSJacobian) interface is enough to support both explicit
204 and implicit timesteppers. This tutorial example also deals with the
205 IFunction/IJacobian interface for demonstration and testing purposes. */
206 PetscCall(PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL));
207 if (rhs_form) {
208 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
209 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
210 } else {
211 PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL));
212 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, NULL));
213 }
214
215 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216 Set initial conditions
217 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218 PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
219 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
220 PetscCall(VecSetUp(U));
221 PetscCall(VecGetArray(U, &u));
222 u[0] = 0.0;
223 u[1] = 20.0;
224 PetscCall(VecRestoreArray(U, &u));
225 PetscCall(TSSetSolution(ts, U));
226
227 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228 Set solver options
229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230 if (hist) PetscCall(TSSetSaveTrajectory(ts));
231 PetscCall(TSSetMaxTime(ts, 30.0));
232 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
233 PetscCall(TSSetTimeStep(ts, 0.1));
234 /* The adaptive time step controller could take very large timesteps
235 jumping over the next event zero-crossing point. A maximum step size
236 limit is enforced here to avoid this issue. */
237 PetscCall(TSGetAdapt(ts, &adapt));
238 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));
239
240 /* Set direction and terminate flag for the event */
241 direction[0] = -1;
242 terminate[0] = PETSC_FALSE;
243 PetscCall(TSSetEventHandler(ts, 1, direction, terminate, EventFunction, PostEventFunction, (void *)&app));
244
245 PetscCall(TSSetFromOptions(ts));
246
247 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248 Run timestepping solver
249 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
250 PetscCall(TSSolve(ts, U));
251
252 if (hist) { /* replay following history */
253 TSTrajectory tj;
254 PetscReal tf, t0, dt;
255
256 app.nbounces = 0;
257 PetscCall(TSGetTime(ts, &tf));
258 PetscCall(TSSetMaxTime(ts, tf));
259 PetscCall(TSSetStepNumber(ts, 0));
260 PetscCall(TSRestartStep(ts));
261 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
262 PetscCall(TSSetFromOptions(ts));
263 PetscCall(TSGetAdapt(ts, &adapt));
264 PetscCall(TSAdaptSetType(adapt, TSADAPTHISTORY));
265 PetscCall(TSGetTrajectory(ts, &tj));
266 PetscCall(TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE));
267 PetscCall(TSAdaptHistoryGetStep(adapt, 0, &t0, &dt));
268 /* this example fails with single (or smaller) precision */
269 #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
270 /*
271 In the first TSSolve() the final time 'tf' is the event location found after a few event handler iterations.
272 If 'tf' is set as the max time for the second run, the TS solver may approach this point by
273 slightly different steps, resulting in a slightly different solution and fvalue[] at 'tf',
274 so that the event may not be triggered at 'tf' anymore. Fix: apply safety factor 1.05
275 */
276 PetscCall(TSSetMaxTime(ts, tf * 1.05));
277 PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC));
278 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));
279 PetscCall(TSSetFromOptions(ts));
280 #endif
281 PetscCall(TSSetTime(ts, t0));
282 PetscCall(TSSetTimeStep(ts, dt));
283 PetscCall(TSResetTrajectory(ts));
284 PetscCall(VecGetArray(U, &u));
285 u[0] = 0.0;
286 u[1] = 20.0;
287 PetscCall(VecRestoreArray(U, &u));
288 PetscCall(TSSolve(ts, U));
289 }
290 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291 Free work space. All PETSc objects should be destroyed when they are no longer needed.
292 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
293 PetscCall(MatDestroy(&A));
294 PetscCall(VecDestroy(&U));
295 PetscCall(TSDestroy(&ts));
296
297 PetscCall(PetscFinalize());
298 return 0;
299 }
300
301 /*TEST
302
303 test:
304 suffix: a
305 args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
306 output_file: output/ex40.out
307
308 test:
309 suffix: b
310 args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
311 output_file: output/ex40.out
312
313 test:
314 suffix: c
315 args: -snes_mf_operator -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
316 output_file: output/ex40.out
317
318 test:
319 suffix: cr
320 args: -rhs-form -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_cr_dir
321 output_file: output/ex40.out
322
323 test:
324 suffix: crmf
325 args: -rhs-form -snes_mf_operator -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_crmf_dir
326 output_file: output/ex40.out
327
328 test:
329 suffix: d
330 args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
331 output_file: output/ex40.out
332
333 test:
334 suffix: e
335 args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
336 output_file: output/ex40.out
337
338 test:
339 suffix: f
340 args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
341 output_file: output/ex40.out
342
343 test:
344 suffix: g
345 args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
346 output_file: output/ex40.out
347
348 test:
349 suffix: h
350 args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
351 output_file: output/ex40.out
352
353 test:
354 suffix: i
355 args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
356 output_file: output/ex40.out
357
358 test:
359 suffix: j
360 args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
361 output_file: output/ex40.out
362
363 test:
364 suffix: k
365 args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir
366 output_file: output/ex40.out
367
368 test:
369 suffix: l
370 args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir
371 args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01
372 output_file: output/ex40.out
373
374 test:
375 suffix: m
376 args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false
377 args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}}
378
379 test:
380 requires: !single
381 suffix: n
382 args: -test_adapthistory false
383 args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view
384 args: -ts_time_step 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
385 args: -ts_max_steps 1 -ts_max_step_rejections {{0 1 2}separate_output} -ts_error_if_step_fails false
386
387 test:
388 requires: !single
389 suffix: o
390 args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir
391 output_file: output/ex40.out
392 TEST*/
393