Actual source code: ex44.c
1: static char help[] = "Parallel bouncing ball example formulated as a second-order system to test TS event feature.\n";
3: /*
4: The dynamics of the bouncing ball with drag coefficient Cd is described by the ODE
6: u_tt = -9.8 - 1/2 Cd (u_t)^2 sign(u_t)
8: There are two events set in this example. The first one checks for the ball hitting the
9: ground (u = 0). Every time the ball hits the ground, its velocity u_t is attenuated by
10: a restitution coefficient Cr. The second event sets a limit on the number of ball bounces.
11: */
13: #include <petscts.h>
15: typedef struct {
16: PetscReal Cd; /* drag coefficient */
17: PetscReal Cr; /* restitution coefficient */
18: PetscInt bounces;
19: PetscInt maxbounces;
20: } AppCtx;
22: static PetscErrorCode Event(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx)
23: {
24: AppCtx *app = (AppCtx *)ctx;
25: Vec V;
26: const PetscScalar *u, *v;
29: /* Event for ball height */
30: TS2GetSolution(ts, &U, &V);
31: VecGetArrayRead(U, &u);
32: VecGetArrayRead(V, &v);
33: fvalue[0] = u[0];
34: /* Event for number of bounces */
35: fvalue[1] = app->maxbounces - app->bounces;
36: VecRestoreArrayRead(U, &u);
37: VecRestoreArrayRead(V, &v);
38: return 0;
39: }
41: static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
42: {
43: AppCtx *app = (AppCtx *)ctx;
44: Vec V;
45: PetscScalar *u, *v;
46: PetscMPIInt rank;
49: if (!nevents) return 0;
50: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
51: if (event_list[0] == 0) {
52: PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball hit the ground at t = %5.2f seconds\n", rank, (double)t);
53: /* Set new initial conditions with .9 attenuation */
54: TS2GetSolution(ts, &U, &V);
55: VecGetArray(U, &u);
56: VecGetArray(V, &v);
57: u[0] = 0.0;
58: v[0] = -app->Cr * v[0];
59: VecRestoreArray(U, &u);
60: VecRestoreArray(V, &v);
61: app->bounces++;
62: } else if (event_list[0] == 1) {
63: PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces);
64: }
65: return 0;
66: }
68: static PetscErrorCode I2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F, void *ctx)
69: {
70: AppCtx *app = (AppCtx *)ctx;
71: const PetscScalar *u, *v, *a;
72: PetscScalar Res, *f;
75: VecGetArrayRead(U, &u);
76: VecGetArrayRead(V, &v);
77: VecGetArrayRead(A, &a);
78: Res = a[0] + 9.8 + 0.5 * app->Cd * v[0] * v[0] * PetscSignReal(PetscRealPart(v[0]));
79: VecRestoreArrayRead(U, &u);
80: VecRestoreArrayRead(V, &v);
81: VecRestoreArrayRead(A, &a);
83: VecGetArray(F, &f);
84: f[0] = Res;
85: VecRestoreArray(F, &f);
86: return 0;
87: }
89: static PetscErrorCode I2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx)
90: {
91: AppCtx *app = (AppCtx *)ctx;
92: const PetscScalar *u, *v, *a;
93: PetscInt i;
94: PetscScalar Jac;
97: VecGetArrayRead(U, &u);
98: VecGetArrayRead(V, &v);
99: VecGetArrayRead(A, &a);
100: Jac = shiftA + shiftV * app->Cd * v[0];
101: VecRestoreArrayRead(U, &u);
102: VecRestoreArrayRead(V, &v);
103: VecRestoreArrayRead(A, &a);
105: MatGetOwnershipRange(P, &i, NULL);
106: MatSetValue(P, i, i, Jac, INSERT_VALUES);
107: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
108: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
109: if (J != P) {
110: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
112: }
113: return 0;
114: }
116: int main(int argc, char **argv)
117: {
118: TS ts; /* ODE integrator */
119: Vec U, V; /* solution will be stored here */
120: Vec F; /* residual vector */
121: Mat J; /* Jacobian matrix */
122: PetscMPIInt rank;
123: PetscScalar *u, *v;
124: AppCtx app;
125: PetscInt direction[2];
126: PetscBool terminate[2];
127: TSAdapt adapt;
130: PetscInitialize(&argc, &argv, NULL, help);
131: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
133: app.Cd = 0.0;
134: app.Cr = 0.9;
135: app.bounces = 0;
136: app.maxbounces = 10;
137: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex44 options", "");
138: PetscOptionsReal("-Cd", "Drag coefficient", "", app.Cd, &app.Cd, NULL);
139: PetscOptionsReal("-Cr", "Restitution coefficient", "", app.Cr, &app.Cr, NULL);
140: PetscOptionsInt("-maxbounces", "Maximum number of bounces", "", app.maxbounces, &app.maxbounces, NULL);
141: PetscOptionsEnd();
143: TSCreate(PETSC_COMM_WORLD, &ts);
144: /*TSSetSaveTrajectory(ts);*/
145: TSSetProblemType(ts, TS_NONLINEAR);
146: TSSetType(ts, TSALPHA2);
148: TSSetMaxTime(ts, PETSC_INFINITY);
149: TSSetTimeStep(ts, 0.1);
150: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
151: TSGetAdapt(ts, &adapt);
152: TSAdaptSetStepLimits(adapt, 0.0, 0.5);
154: direction[0] = -1;
155: terminate[0] = PETSC_FALSE;
156: direction[1] = -1;
157: terminate[1] = PETSC_TRUE;
158: TSSetEventHandler(ts, 2, direction, terminate, Event, PostEvent, &app);
160: MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &J);
161: MatSetFromOptions(J);
162: MatSetUp(J);
163: MatCreateVecs(J, NULL, &F);
164: TSSetI2Function(ts, F, I2Function, &app);
165: TSSetI2Jacobian(ts, J, J, I2Jacobian, &app);
166: VecDestroy(&F);
167: MatDestroy(&J);
169: TSGetI2Jacobian(ts, &J, NULL, NULL, NULL);
170: MatCreateVecs(J, &U, NULL);
171: MatCreateVecs(J, &V, NULL);
172: VecGetArray(U, &u);
173: VecGetArray(V, &v);
174: u[0] = 5.0 * rank;
175: v[0] = 20.0;
176: VecRestoreArray(U, &u);
177: VecRestoreArray(V, &v);
179: TS2SetSolution(ts, U, V);
180: TSSetFromOptions(ts);
181: TSSolve(ts, NULL);
183: VecDestroy(&U);
184: VecDestroy(&V);
185: TSDestroy(&ts);
187: PetscFinalize();
188: return 0;
189: }
191: /*TEST
193: test:
194: suffix: a
195: args: -ts_alpha_radius {{1.0 0.5}}
196: output_file: output/ex44.out
198: test:
199: suffix: b
200: args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic
201: output_file: output/ex44.out
203: test:
204: suffix: 2
205: nsize: 2
206: args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic
207: output_file: output/ex44_2.out
208: filter: sort -b
209: filter_output: sort -b
211: test:
212: requires: !single
213: args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
214: args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false
216: TEST*/