Actual source code: ex41.c
1: static char help[] = "Parallel 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: Each processor is assigned one ball.
10: The event function routine checks for the ball hitting the
11: ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
12: a factor of 0.9 and its height set to 1.0*rank.
13: */
15: #include <petscts.h>
17: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx)
18: {
19: const PetscScalar *u;
22: /* Event for ball height */
23: VecGetArrayRead(U, &u);
24: fvalue[0] = u[0];
25: VecRestoreArrayRead(U, &u);
26: return 0;
27: }
29: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
30: {
31: PetscScalar *u;
32: PetscMPIInt rank;
35: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
36: if (nevents) {
37: PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds -> Processor[%d]\n", (double)t, rank);
38: /* Set new initial conditions with .9 attenuation */
39: VecGetArray(U, &u);
40: u[0] = 1.0 * rank;
41: u[1] = -0.9 * u[1];
42: VecRestoreArray(U, &u);
43: }
44: return 0;
45: }
47: /*
48: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
49: */
50: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
51: {
52: PetscScalar *f;
53: const PetscScalar *u;
56: /* The next three lines allow us to access the entries of the vectors directly */
57: VecGetArrayRead(U, &u);
58: VecGetArray(F, &f);
60: f[0] = u[1];
61: f[1] = -9.8;
63: VecRestoreArrayRead(U, &u);
64: VecRestoreArray(F, &f);
65: return 0;
66: }
68: /*
69: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian.
70: */
71: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
72: {
73: PetscInt rowcol[2], rstart;
74: PetscScalar J[2][2];
75: const PetscScalar *u;
78: VecGetArrayRead(U, &u);
80: MatGetOwnershipRange(B, &rstart, NULL);
81: rowcol[0] = rstart;
82: rowcol[1] = rstart + 1;
84: J[0][0] = 0.0;
85: J[0][1] = 1.0;
86: J[1][0] = 0.0;
87: J[1][1] = 0.0;
88: MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
90: VecRestoreArrayRead(U, &u);
91: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
93: if (A != B) {
94: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
96: }
97: return 0;
98: }
100: /*
101: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
102: */
103: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
104: {
105: PetscScalar *f;
106: const PetscScalar *u, *udot;
109: /* The next three lines allow us to access the entries of the vectors directly */
110: VecGetArrayRead(U, &u);
111: VecGetArrayRead(Udot, &udot);
112: VecGetArray(F, &f);
114: f[0] = udot[0] - u[1];
115: f[1] = udot[1] + 9.8;
117: VecRestoreArrayRead(U, &u);
118: VecRestoreArrayRead(Udot, &udot);
119: VecRestoreArray(F, &f);
120: return 0;
121: }
123: /*
124: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
125: */
126: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
127: {
128: PetscInt rowcol[2], rstart;
129: PetscScalar J[2][2];
130: const PetscScalar *u, *udot;
133: VecGetArrayRead(U, &u);
134: VecGetArrayRead(Udot, &udot);
136: MatGetOwnershipRange(B, &rstart, NULL);
137: rowcol[0] = rstart;
138: rowcol[1] = rstart + 1;
140: J[0][0] = a;
141: J[0][1] = -1.0;
142: J[1][0] = 0.0;
143: J[1][1] = a;
144: MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
146: VecRestoreArrayRead(U, &u);
147: VecRestoreArrayRead(Udot, &udot);
149: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
150: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
151: if (A != B) {
152: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
154: }
155: return 0;
156: }
158: int main(int argc, char **argv)
159: {
160: TS ts; /* ODE integrator */
161: Vec U; /* solution will be stored here */
162: PetscMPIInt rank;
163: PetscInt n = 2;
164: PetscScalar *u;
165: PetscInt direction = -1;
166: PetscBool terminate = PETSC_FALSE;
167: PetscBool rhs_form = PETSC_FALSE, hist = PETSC_TRUE;
168: TSAdapt adapt;
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Initialize program
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: PetscInitialize(&argc, &argv, (char *)0, help);
175: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Create timestepping solver context
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: TSCreate(PETSC_COMM_WORLD, &ts);
181: TSSetType(ts, TSROSW);
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Set ODE routines
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: TSSetProblemType(ts, TS_NONLINEAR);
187: /* Users are advised against the following branching and code duplication.
188: For problems without a mass matrix like the one at hand, the RHSFunction
189: (and companion RHSJacobian) interface is enough to support both explicit
190: and implicit timesteppers. This tutorial example also deals with the
191: IFunction/IJacobian interface for demonstration and testing purposes. */
192: PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL);
193: if (rhs_form) {
194: TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
195: TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL);
196: } else {
197: TSSetIFunction(ts, NULL, IFunction, NULL);
198: TSSetIJacobian(ts, NULL, NULL, IJacobian, NULL);
199: }
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Set initial conditions
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: VecCreate(PETSC_COMM_WORLD, &U);
205: VecSetSizes(U, n, PETSC_DETERMINE);
206: VecSetUp(U);
207: VecGetArray(U, &u);
208: u[0] = 1.0 * rank;
209: u[1] = 20.0;
210: VecRestoreArray(U, &u);
211: TSSetSolution(ts, U);
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Set solver options
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: TSSetSaveTrajectory(ts);
217: TSSetMaxTime(ts, 30.0);
218: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
219: TSSetTimeStep(ts, 0.1);
220: /* The adaptive time step controller could take very large timesteps resulting in
221: the same event occurring multiple times in the same interval. A maximum step size
222: limit is enforced here to avoid this issue. */
223: TSGetAdapt(ts, &adapt);
224: TSAdaptSetType(adapt, TSADAPTBASIC);
225: TSAdaptSetStepLimits(adapt, 0.0, 0.5);
227: /* Set direction and terminate flag for the event */
228: TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL);
230: TSSetFromOptions(ts);
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Run timestepping solver
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: TSSolve(ts, U);
237: if (hist) { /* replay following history */
238: TSTrajectory tj;
239: PetscReal tf, t0, dt;
241: TSGetTime(ts, &tf);
242: TSSetMaxTime(ts, tf);
243: TSSetStepNumber(ts, 0);
244: TSRestartStep(ts);
245: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
246: TSSetFromOptions(ts);
247: TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL);
248: TSGetAdapt(ts, &adapt);
249: TSAdaptSetType(adapt, TSADAPTHISTORY);
250: TSGetTrajectory(ts, &tj);
251: TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE);
252: TSAdaptHistoryGetStep(adapt, 0, &t0, &dt);
253: /* this example fails with single (or smaller) precision */
254: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
255: TSAdaptSetType(adapt, TSADAPTBASIC);
256: TSAdaptSetStepLimits(adapt, 0.0, 0.5);
257: TSSetFromOptions(ts);
258: #endif
259: TSSetTime(ts, t0);
260: TSSetTimeStep(ts, dt);
261: TSResetTrajectory(ts);
262: VecGetArray(U, &u);
263: u[0] = 1.0 * rank;
264: u[1] = 20.0;
265: VecRestoreArray(U, &u);
266: TSSolve(ts, U);
267: }
268: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: Free work space. All PETSc objects should be destroyed when they are no longer needed.
270: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: VecDestroy(&U);
272: TSDestroy(&ts);
274: PetscFinalize();
275: return 0;
276: }
278: /*TEST
280: test:
281: suffix: a
282: nsize: 2
283: args: -ts_trajectory_type memory -snes_stol 1e-4
284: filter: sort -b
286: test:
287: suffix: b
288: nsize: 2
289: args: -ts_trajectory_type memory -ts_type arkimex -snes_stol 1e-4
290: filter: sort -b
292: test:
293: suffix: c
294: nsize: 2
295: args: -ts_trajectory_type memory -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
296: filter: sort -b
298: test:
299: suffix: d
300: nsize: 2
301: args: -ts_trajectory_type memory -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
302: filter: sort -b
304: test:
305: suffix: e
306: nsize: 2
307: args: -ts_trajectory_type memory -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
308: filter: sort -b
310: test:
311: suffix: f
312: nsize: 2
313: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 3bs
314: filter: sort -b
316: test:
317: suffix: g
318: nsize: 2
319: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 5bs
320: filter: sort -b
322: test:
323: suffix: h
324: nsize: 2
325: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 6vr
326: filter: sort -b
327: output_file: output/ex41_g.out
329: test:
330: suffix: i
331: nsize: 2
332: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 7vr
333: filter: sort -b
334: output_file: output/ex41_g.out
336: test:
337: suffix: j
338: nsize: 2
339: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 8vr
340: filter: sort -b
341: output_file: output/ex41_g.out
343: TEST*/