Actual source code: ex40.c
1: static char help[] = "Serial bouncing ball example to test TS event feature.\n";
3: /*
4: The dynamics of the bouncing ball is described by the ODE
5: u1_t = u2
6: u2_t = -9.8
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: */
14: #include <petscts.h>
16: typedef struct {
17: PetscInt maxbounces;
18: PetscInt nbounces;
19: } AppCtx;
21: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx)
22: {
23: const PetscScalar *u;
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: }
33: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
34: {
35: AppCtx *app = (AppCtx *)ctx;
36: PetscScalar *u;
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: }
53: /*
54: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
55: */
56: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
57: {
58: PetscScalar *f;
59: const PetscScalar *u;
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));
66: f[0] = u[1];
67: f[1] = -9.8;
69: PetscCall(VecRestoreArrayRead(U, &u));
70: PetscCall(VecRestoreArray(F, &f));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*
75: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
76: */
77: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
78: {
79: PetscInt rowcol[] = {0, 1};
80: PetscScalar J[2][2];
81: const PetscScalar *u;
83: PetscFunctionBeginUser;
84: PetscCall(VecGetArrayRead(U, &u));
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));
92: PetscCall(VecRestoreArrayRead(U, &u));
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: }
103: /*
104: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
105: */
106: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
107: {
108: PetscScalar *f;
109: const PetscScalar *u, *udot;
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));
117: f[0] = udot[0] - u[1];
118: f[1] = udot[1] + 9.8;
120: PetscCall(VecRestoreArrayRead(U, &u));
121: PetscCall(VecRestoreArrayRead(Udot, &udot));
122: PetscCall(VecRestoreArray(F, &f));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*
127: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
128: */
129: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
130: {
131: PetscInt rowcol[] = {0, 1};
132: PetscScalar J[2][2];
133: const PetscScalar *u, *udot;
135: PetscFunctionBeginUser;
136: PetscCall(VecGetArrayRead(U, &u));
137: PetscCall(VecGetArrayRead(Udot, &udot));
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));
145: PetscCall(VecRestoreArrayRead(U, &u));
146: PetscCall(VecRestoreArrayRead(Udot, &udot));
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: }
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;
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Initialize program
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: PetscFunctionBeginUser;
174: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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");
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();
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));
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: }
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));
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));
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));
245: PetscCall(TSSetFromOptions(ts));
247: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248: Run timestepping solver
249: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
250: PetscCall(TSSolve(ts, U));
252: if (hist) { /* replay following history */
253: TSTrajectory tj;
254: PetscReal tf, t0, dt;
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));
297: PetscCall(PetscFinalize());
298: return 0;
299: }
301: /*TEST
303: test:
304: suffix: a
305: args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
306: output_file: output/ex40.out
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
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
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
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
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
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
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
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
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
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
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
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
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
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}}
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_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
385: args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false
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*/