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 are two events set in this example. The first one 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. The second event sets a limit on the number of ball bounces.
11: */
13: #include <petscts.h>
15: typedef struct {
16: PetscInt maxbounces;
17: PetscInt nbounces;
18: } AppCtx;
20: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx)
21: {
22: AppCtx *app = (AppCtx *)ctx;
23: const PetscScalar *u;
26: /* Event for ball height */
27: VecGetArrayRead(U, &u);
28: fvalue[0] = u[0];
29: /* Event for number of bounces */
30: fvalue[1] = app->maxbounces - app->nbounces;
31: VecRestoreArrayRead(U, &u);
32: return 0;
33: }
35: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
36: {
37: AppCtx *app = (AppCtx *)ctx;
38: PetscScalar *u;
41: if (event_list[0] == 0) {
42: PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds\n", (double)t);
43: /* Set new initial conditions with .9 attenuation */
44: VecGetArray(U, &u);
45: u[0] = 0.0;
46: u[1] = -0.9 * u[1];
47: VecRestoreArray(U, &u);
48: } else if (event_list[0] == 1) {
49: PetscPrintf(PETSC_COMM_SELF, "Ball bounced %" PetscInt_FMT " times\n", app->nbounces);
50: }
51: app->nbounces++;
52: return 0;
53: }
55: /*
56: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
57: */
58: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
59: {
60: PetscScalar *f;
61: const PetscScalar *u;
64: /* The next three lines allow us to access the entries of the vectors directly */
65: VecGetArrayRead(U, &u);
66: VecGetArray(F, &f);
68: f[0] = u[1];
69: f[1] = -9.8;
71: VecRestoreArrayRead(U, &u);
72: VecRestoreArray(F, &f);
73: return 0;
74: }
76: /*
77: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
78: */
79: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
80: {
81: PetscInt rowcol[] = {0, 1};
82: PetscScalar J[2][2];
83: const PetscScalar *u;
86: VecGetArrayRead(U, &u);
88: J[0][0] = 0.0;
89: J[0][1] = 1.0;
90: J[1][0] = 0.0;
91: J[1][1] = 0.0;
92: MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
94: VecRestoreArrayRead(U, &u);
96: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
98: if (A != B) {
99: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
101: }
102: return 0;
103: }
105: /*
106: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
107: */
108: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
109: {
110: PetscScalar *f;
111: const PetscScalar *u, *udot;
114: /* The next three lines allow us to access the entries of the vectors directly */
115: VecGetArrayRead(U, &u);
116: VecGetArrayRead(Udot, &udot);
117: VecGetArray(F, &f);
119: f[0] = udot[0] - u[1];
120: f[1] = udot[1] + 9.8;
122: VecRestoreArrayRead(U, &u);
123: VecRestoreArrayRead(Udot, &udot);
124: VecRestoreArray(F, &f);
125: return 0;
126: }
128: /*
129: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
130: */
131: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
132: {
133: PetscInt rowcol[] = {0, 1};
134: PetscScalar J[2][2];
135: const PetscScalar *u, *udot;
138: VecGetArrayRead(U, &u);
139: VecGetArrayRead(Udot, &udot);
141: J[0][0] = a;
142: J[0][1] = -1.0;
143: J[1][0] = 0.0;
144: J[1][1] = a;
145: MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
147: VecRestoreArrayRead(U, &u);
148: VecRestoreArrayRead(Udot, &udot);
150: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
151: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
152: if (A != B) {
153: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
154: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
155: }
156: return 0;
157: }
159: int main(int argc, char **argv)
160: {
161: TS ts; /* ODE integrator */
162: Vec U; /* solution will be stored here */
163: PetscMPIInt size;
164: PetscInt n = 2;
165: PetscScalar *u;
166: AppCtx app;
167: PetscInt direction[2];
168: PetscBool terminate[2];
169: PetscBool rhs_form = PETSC_FALSE, hist = PETSC_TRUE;
170: TSAdapt adapt;
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Initialize program
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: PetscInitialize(&argc, &argv, (char *)0, help);
177: MPI_Comm_size(PETSC_COMM_WORLD, &size);
180: app.nbounces = 0;
181: app.maxbounces = 10;
182: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex40 options", "");
183: PetscOptionsInt("-maxbounces", "", "", app.maxbounces, &app.maxbounces, NULL);
184: PetscOptionsBool("-test_adapthistory", "", "", hist, &hist, NULL);
185: PetscOptionsEnd();
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Create timestepping solver context
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190: TSCreate(PETSC_COMM_WORLD, &ts);
191: TSSetType(ts, TSROSW);
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Set ODE routines
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: TSSetProblemType(ts, TS_NONLINEAR);
197: /* Users are advised against the following branching and code duplication.
198: For problems without a mass matrix like the one at hand, the RHSFunction
199: (and companion RHSJacobian) interface is enough to support both explicit
200: and implicit timesteppers. This tutorial example also deals with the
201: IFunction/IJacobian interface for demonstration and testing purposes. */
202: PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL);
203: if (rhs_form) {
204: TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
205: TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL);
206: } else {
207: Mat A; /* Jacobian matrix */
208: MatCreate(PETSC_COMM_WORLD, &A);
209: MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
210: MatSetType(A, MATDENSE);
211: MatSetFromOptions(A);
212: MatSetUp(A);
213: TSSetIFunction(ts, NULL, IFunction, NULL);
214: TSSetIJacobian(ts, A, A, IJacobian, NULL);
215: MatDestroy(&A);
216: }
218: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219: Set initial conditions
220: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221: VecCreate(PETSC_COMM_WORLD, &U);
222: VecSetSizes(U, n, PETSC_DETERMINE);
223: VecSetUp(U);
224: VecGetArray(U, &u);
225: u[0] = 0.0;
226: u[1] = 20.0;
227: VecRestoreArray(U, &u);
228: TSSetSolution(ts, U);
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Set solver options
232: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233: if (hist) TSSetSaveTrajectory(ts);
234: TSSetMaxTime(ts, 30.0);
235: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
236: TSSetTimeStep(ts, 0.1);
237: /* The adaptive time step controller could take very large timesteps resulting in
238: the same event occurring multiple times in the same interval. A maximum step size
239: limit is enforced here to avoid this issue. */
240: TSGetAdapt(ts, &adapt);
241: TSAdaptSetStepLimits(adapt, 0.0, 0.5);
243: /* Set directions and terminate flags for the two events */
244: direction[0] = -1;
245: direction[1] = -1;
246: terminate[0] = PETSC_FALSE;
247: terminate[1] = PETSC_TRUE;
248: TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&app);
250: TSSetFromOptions(ts);
252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: Run timestepping solver
254: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255: TSSolve(ts, U);
257: if (hist) { /* replay following history */
258: TSTrajectory tj;
259: PetscReal tf, t0, dt;
261: app.nbounces = 0;
262: TSGetTime(ts, &tf);
263: TSSetMaxTime(ts, tf);
264: TSSetStepNumber(ts, 0);
265: TSRestartStep(ts);
266: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
267: TSSetFromOptions(ts);
268: TSGetAdapt(ts, &adapt);
269: TSAdaptSetType(adapt, TSADAPTHISTORY);
270: TSGetTrajectory(ts, &tj);
271: TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE);
272: TSAdaptHistoryGetStep(adapt, 0, &t0, &dt);
273: /* this example fails with single (or smaller) precision */
274: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
275: TSAdaptSetType(adapt, TSADAPTBASIC);
276: TSAdaptSetStepLimits(adapt, 0.0, 0.5);
277: TSSetFromOptions(ts);
278: #endif
279: TSSetTime(ts, t0);
280: TSSetTimeStep(ts, dt);
281: TSResetTrajectory(ts);
282: VecGetArray(U, &u);
283: u[0] = 0.0;
284: u[1] = 20.0;
285: VecRestoreArray(U, &u);
286: TSSolve(ts, U);
287: }
288: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: Free work space. All PETSc objects should be destroyed when they are no longer needed.
290: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291: VecDestroy(&U);
292: TSDestroy(&ts);
294: PetscFinalize();
295: return 0;
296: }
298: /*TEST
300: test:
301: suffix: a
302: args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
303: output_file: output/ex40.out
305: test:
306: suffix: b
307: args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
308: output_file: output/ex40.out
310: test:
311: suffix: c
312: args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
313: output_file: output/ex40.out
315: test:
316: suffix: d
317: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
318: output_file: output/ex40.out
320: test:
321: suffix: e
322: args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
323: output_file: output/ex40.out
325: test:
326: suffix: f
327: args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
328: output_file: output/ex40.out
330: test:
331: suffix: g
332: args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
333: output_file: output/ex40.out
335: test:
336: suffix: h
337: args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
338: output_file: output/ex40.out
340: test:
341: suffix: i
342: args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
343: output_file: output/ex40.out
345: test:
346: suffix: j
347: args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
348: output_file: output/ex40.out
350: test:
351: suffix: k
352: args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir
353: output_file: output/ex40.out
355: test:
356: suffix: l
357: args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir
358: args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01
359: output_file: output/ex40.out
361: test:
362: suffix: m
363: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false
364: args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}}
366: test:
367: requires: !single
368: suffix: n
369: args: -test_adapthistory false
370: args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view
371: args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
372: args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false
374: test:
375: requires: !single
376: suffix: o
377: args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir
378: output_file: output/ex40.out
379: TEST*/