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 is one event set in this example, which 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. 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: PetscReal Cd; /* drag coefficient */
18: PetscReal Cr; /* restitution coefficient */
19: PetscInt bounces;
20: PetscInt maxbounces;
21: } AppCtx;
23: static PetscErrorCode Event(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx)
24: {
25: Vec V;
26: const PetscScalar *u;
28: PetscFunctionBeginUser;
29: /* Event for ball height */
30: PetscCall(TS2GetSolution(ts, &U, &V));
31: PetscCall(VecGetArrayRead(U, &u));
32: fvalue[0] = PetscRealPart(u[0]);
33: PetscCall(VecRestoreArrayRead(U, &u));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
38: {
39: AppCtx *app = (AppCtx *)ctx;
40: Vec V;
41: PetscScalar *u, *v;
42: PetscMPIInt rank;
43: PetscBool inflag = PETSC_FALSE, outflag;
45: PetscFunctionBeginUser;
46: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
47: if (nevents > 0) {
48: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball hit the ground at t = %5.2f seconds\n", rank, (double)t));
49: /* Set new initial conditions with attenuation Cr */
50: PetscCall(TS2GetSolution(ts, &U, &V));
51: PetscCall(VecGetArray(U, &u));
52: PetscCall(VecGetArray(V, &v));
53: u[0] = 0.0;
54: v[0] = -app->Cr * v[0];
55: PetscCall(VecRestoreArray(U, &u));
56: PetscCall(VecRestoreArray(V, &v));
57: app->bounces++;
58: }
59: if (app->bounces >= app->maxbounces) { // 'app->bounces' may be different on different processes
60: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces));
61: inflag = PETSC_TRUE; // current process requested to terminate
62: }
63: PetscCall(MPIU_Allreduce(&inflag, &outflag, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)ts)));
64: if (outflag) PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER)); // request TS to terminate, sync on all ranks
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(P, MAT_FINAL_ASSEMBLY));
108: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
109: if (J != P) {
110: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
111: PetscCall(MatAssemblyEnd(J, 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[1];
126: PetscBool terminate[1];
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: PetscCall(TSSetEventHandler(ts, 1, direction, terminate, Event, PostEvent, &app)); // each process has one event-function defined
158: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &J));
159: PetscCall(MatSetFromOptions(J));
160: PetscCall(MatSetUp(J));
161: PetscCall(MatCreateVecs(J, NULL, &F));
162: PetscCall(TSSetI2Function(ts, F, I2Function, &app));
163: PetscCall(TSSetI2Jacobian(ts, J, J, I2Jacobian, &app));
164: PetscCall(VecDestroy(&F));
165: PetscCall(MatDestroy(&J));
167: PetscCall(TSGetI2Jacobian(ts, &J, NULL, NULL, NULL));
168: PetscCall(MatCreateVecs(J, &U, NULL));
169: PetscCall(MatCreateVecs(J, &V, NULL));
170: PetscCall(VecGetArray(U, &u));
171: PetscCall(VecGetArray(V, &v));
172: u[0] = 5.0 * rank;
173: v[0] = 20.0;
174: PetscCall(VecRestoreArray(U, &u));
175: PetscCall(VecRestoreArray(V, &v));
177: PetscCall(TS2SetSolution(ts, U, V));
178: PetscCall(TSSetFromOptions(ts));
179: PetscCall(TSSolve(ts, NULL));
181: PetscCall(VecDestroy(&U));
182: PetscCall(VecDestroy(&V));
183: PetscCall(TSDestroy(&ts));
185: PetscCall(PetscFinalize());
186: return 0;
187: }
189: /*TEST
191: test:
192: suffix: a
193: args: -ts_alpha_radius {{1.0 0.5}} -ts_max_time 50
194: output_file: output/ex44.out
196: test:
197: suffix: b
198: args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50
199: output_file: output/ex44.out
201: test:
202: suffix: bmf
203: args: -snes_mf_operator -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50
204: output_file: output/ex44.out
206: test:
207: suffix: 2
208: nsize: 2
209: args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50
210: output_file: output/ex44_2.out
211: filter: sort -b
212: filter_output: sort -b
214: test:
215: requires: !single
216: args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
217: args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false
219: TEST*/