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;
28: PetscFunctionBeginUser;
29: /* Event for ball height */
30: PetscCall(TS2GetSolution(ts, &U, &V));
31: PetscCall(VecGetArrayRead(U, &u));
32: PetscCall(VecGetArrayRead(V, &v));
33: fvalue[0] = u[0];
34: /* Event for number of bounces */
35: fvalue[1] = app->maxbounces - app->bounces;
36: PetscCall(VecRestoreArrayRead(U, &u));
37: PetscCall(VecRestoreArrayRead(V, &v));
38: PetscFunctionReturn(PETSC_SUCCESS);
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;
48: PetscFunctionBeginUser;
49: if (!nevents) PetscFunctionReturn(PETSC_SUCCESS);
50: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
51: if (event_list[0] == 0) {
52: PetscCall(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: PetscCall(TS2GetSolution(ts, &U, &V));
55: PetscCall(VecGetArray(U, &u));
56: PetscCall(VecGetArray(V, &v));
57: u[0] = 0.0;
58: v[0] = -app->Cr * v[0];
59: PetscCall(VecRestoreArray(U, &u));
60: PetscCall(VecRestoreArray(V, &v));
61: app->bounces++;
62: } else if (event_list[0] == 1) {
63: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces));
64: }
65: PetscFunctionReturn(PETSC_SUCCESS);
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;
74: PetscFunctionBeginUser;
75: PetscCall(VecGetArrayRead(U, &u));
76: PetscCall(VecGetArrayRead(V, &v));
77: PetscCall(VecGetArrayRead(A, &a));
78: Res = a[0] + 9.8 + 0.5 * app->Cd * v[0] * v[0] * PetscSignReal(PetscRealPart(v[0]));
79: PetscCall(VecRestoreArrayRead(U, &u));
80: PetscCall(VecRestoreArrayRead(V, &v));
81: PetscCall(VecRestoreArrayRead(A, &a));
83: PetscCall(VecGetArray(F, &f));
84: f[0] = Res;
85: PetscCall(VecRestoreArray(F, &f));
86: PetscFunctionReturn(PETSC_SUCCESS);
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;
96: PetscFunctionBeginUser;
97: PetscCall(VecGetArrayRead(U, &u));
98: PetscCall(VecGetArrayRead(V, &v));
99: PetscCall(VecGetArrayRead(A, &a));
100: Jac = shiftA + shiftV * app->Cd * v[0];
101: PetscCall(VecRestoreArrayRead(U, &u));
102: PetscCall(VecRestoreArrayRead(V, &v));
103: PetscCall(VecRestoreArrayRead(A, &a));
105: PetscCall(MatGetOwnershipRange(P, &i, NULL));
106: PetscCall(MatSetValue(P, i, i, Jac, INSERT_VALUES));
107: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
108: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
109: if (J != P) {
110: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
111: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
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;
129: PetscFunctionBeginUser;
130: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
131: PetscCallMPI(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: PetscCall(PetscOptionsReal("-Cd", "Drag coefficient", "", app.Cd, &app.Cd, NULL));
139: PetscCall(PetscOptionsReal("-Cr", "Restitution coefficient", "", app.Cr, &app.Cr, NULL));
140: PetscCall(PetscOptionsInt("-maxbounces", "Maximum number of bounces", "", app.maxbounces, &app.maxbounces, NULL));
141: PetscOptionsEnd();
143: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
144: /*PetscCall(TSSetSaveTrajectory(ts));*/
145: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
146: PetscCall(TSSetType(ts, TSALPHA2));
148: PetscCall(TSSetMaxTime(ts, PETSC_INFINITY));
149: PetscCall(TSSetTimeStep(ts, 0.1));
150: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
151: PetscCall(TSGetAdapt(ts, &adapt));
152: PetscCall(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: PetscCall(TSSetEventHandler(ts, 2, direction, terminate, Event, PostEvent, &app));
160: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &J));
161: PetscCall(MatSetFromOptions(J));
162: PetscCall(MatSetUp(J));
163: PetscCall(MatCreateVecs(J, NULL, &F));
164: PetscCall(TSSetI2Function(ts, F, I2Function, &app));
165: PetscCall(TSSetI2Jacobian(ts, J, J, I2Jacobian, &app));
166: PetscCall(VecDestroy(&F));
167: PetscCall(MatDestroy(&J));
169: PetscCall(TSGetI2Jacobian(ts, &J, NULL, NULL, NULL));
170: PetscCall(MatCreateVecs(J, &U, NULL));
171: PetscCall(MatCreateVecs(J, &V, NULL));
172: PetscCall(VecGetArray(U, &u));
173: PetscCall(VecGetArray(V, &v));
174: u[0] = 5.0 * rank;
175: v[0] = 20.0;
176: PetscCall(VecRestoreArray(U, &u));
177: PetscCall(VecRestoreArray(V, &v));
179: PetscCall(TS2SetSolution(ts, U, V));
180: PetscCall(TSSetFromOptions(ts));
181: PetscCall(TSSolve(ts, NULL));
183: PetscCall(VecDestroy(&U));
184: PetscCall(VecDestroy(&V));
185: PetscCall(TSDestroy(&ts));
187: PetscCall(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*/