1 static char help[] = "Parallel bouncing ball example formulated as a second-order system to test TS event feature.\n";
2
3 /*
4 The dynamics of the bouncing ball with drag coefficient Cd is described by the ODE
5
6 u_tt = -9.8 - 1/2 Cd (u_t)^2 sign(u_t)
7
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 */
13
14 #include <petscts.h>
15
16 typedef struct {
17 PetscReal Cd; /* drag coefficient */
18 PetscReal Cr; /* restitution coefficient */
19 PetscInt bounces;
20 PetscInt maxbounces;
21 } AppCtx;
22
Event(TS ts,PetscReal t,Vec U,PetscReal * fvalue,PetscCtx ctx)23 static PetscErrorCode Event(TS ts, PetscReal t, Vec U, PetscReal *fvalue, PetscCtx ctx)
24 {
25 Vec V;
26 const PetscScalar *u;
27
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 }
36
PostEvent(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,PetscCtx ctx)37 static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx)
38 {
39 AppCtx *app = (AppCtx *)ctx;
40 Vec V;
41 PetscScalar *u, *v;
42 PetscMPIInt rank;
43 PetscBool inflag = PETSC_FALSE, outflag;
44
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 PetscCallMPI(MPIU_Allreduce(&inflag, &outflag, 1, MPI_C_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 }
67
I2Function(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec F,PetscCtx ctx)68 static PetscErrorCode I2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F, PetscCtx ctx)
69 {
70 AppCtx *app = (AppCtx *)ctx;
71 const PetscScalar *u, *v, *a;
72 PetscScalar Res, *f;
73
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));
82
83 PetscCall(VecGetArray(F, &f));
84 f[0] = Res;
85 PetscCall(VecRestoreArray(F, &f));
86 PetscFunctionReturn(PETSC_SUCCESS);
87 }
88
I2Jacobian(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P,PetscCtx ctx)89 static PetscErrorCode I2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, PetscCtx ctx)
90 {
91 AppCtx *app = (AppCtx *)ctx;
92 const PetscScalar *u, *v, *a;
93 PetscInt i;
94 PetscScalar Jac;
95
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));
104
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 }
115
main(int argc,char ** argv)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;
128
129 PetscFunctionBeginUser;
130 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
131 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
132
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();
142
143 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
144 /*PetscCall(TSSetSaveTrajectory(ts));*/
145 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
146 PetscCall(TSSetType(ts, TSALPHA2));
147
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));
153
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
157
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));
166
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));
176
177 PetscCall(TS2SetSolution(ts, U, V));
178 PetscCall(TSSetFromOptions(ts));
179 PetscCall(TSSolve(ts, NULL));
180
181 PetscCall(VecDestroy(&U));
182 PetscCall(VecDestroy(&V));
183 PetscCall(TSDestroy(&ts));
184
185 PetscCall(PetscFinalize());
186 return 0;
187 }
188
189 /*TEST
190
191 test:
192 suffix: a
193 args: -ts_alpha_radius {{1.0 0.5}} -ts_max_time 50
194 output_file: output/ex44.out
195
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
200
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
205
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
213
214 test:
215 requires: !single
216 args: -ts_time_step 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
217 args: -ts_max_steps 1 -ts_max_step_rejections {{0 1 2}separate_output} -ts_error_if_step_fails false
218
219 TEST*/
220