1 static char help[] = "Parallel bouncing ball example to test TS event feature.\n";
2
3 /*
4 The dynamics of the bouncing ball is described by the ODE
5 u1_t = u2
6 u2_t = -9.8
7
8 Each processor is assigned one ball.
9
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 */
14
15 #include <petscts.h>
16
EventFunction(TS ts,PetscReal t,Vec U,PetscReal * fvalue,PetscCtx ctx)17 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, PetscCtx ctx)
18 {
19 const PetscScalar *u;
20
21 PetscFunctionBeginUser;
22 /* Event for ball height */
23 PetscCall(VecGetArrayRead(U, &u));
24 fvalue[0] = PetscRealPart(u[0]);
25 PetscCall(VecRestoreArrayRead(U, &u));
26 PetscFunctionReturn(PETSC_SUCCESS);
27 }
28
PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,PetscCtx ctx)29 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx)
30 {
31 PetscScalar *u;
32 PetscMPIInt rank;
33
34 PetscFunctionBeginUser;
35 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
36 if (nevents) {
37 PetscCall(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 PetscCall(VecGetArray(U, &u));
40 u[0] = 1.0 * rank;
41 u[1] = -0.9 * u[1];
42 PetscCall(VecRestoreArray(U, &u));
43 }
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
47 /*
48 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
49 */
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,PetscCtx ctx)50 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx)
51 {
52 PetscScalar *f;
53 const PetscScalar *u;
54
55 PetscFunctionBeginUser;
56 /* The next three lines allow us to access the entries of the vectors directly */
57 PetscCall(VecGetArrayRead(U, &u));
58 PetscCall(VecGetArray(F, &f));
59
60 f[0] = u[1];
61 f[1] = -9.8;
62
63 PetscCall(VecRestoreArrayRead(U, &u));
64 PetscCall(VecRestoreArray(F, &f));
65 PetscFunctionReturn(PETSC_SUCCESS);
66 }
67
68 /*
69 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian.
70 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,PetscCtx ctx)71 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, PetscCtx ctx)
72 {
73 PetscInt rowcol[2], rstart;
74 PetscScalar J[2][2];
75 const PetscScalar *u;
76
77 PetscFunctionBeginUser;
78 PetscCall(VecGetArrayRead(U, &u));
79
80 PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
81 rowcol[0] = rstart;
82 rowcol[1] = rstart + 1;
83
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 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
89
90 PetscCall(VecRestoreArrayRead(U, &u));
91 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
92 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
93 if (A != B) {
94 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
95 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
96 }
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99
100 /*
101 Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
102 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)103 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
104 {
105 PetscScalar *f;
106 const PetscScalar *u, *udot;
107
108 PetscFunctionBeginUser;
109 /* The next three lines allow us to access the entries of the vectors directly */
110 PetscCall(VecGetArrayRead(U, &u));
111 PetscCall(VecGetArrayRead(Udot, &udot));
112 PetscCall(VecGetArray(F, &f));
113
114 f[0] = udot[0] - u[1];
115 f[1] = udot[1] + 9.8;
116
117 PetscCall(VecRestoreArrayRead(U, &u));
118 PetscCall(VecRestoreArrayRead(Udot, &udot));
119 PetscCall(VecRestoreArray(F, &f));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
123 /*
124 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
125 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,PetscCtx ctx)126 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
127 {
128 PetscInt rowcol[2], rstart;
129 PetscScalar J[2][2];
130 const PetscScalar *u, *udot;
131
132 PetscFunctionBeginUser;
133 PetscCall(VecGetArrayRead(U, &u));
134 PetscCall(VecGetArrayRead(Udot, &udot));
135
136 PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
137 rowcol[0] = rstart;
138 rowcol[1] = rstart + 1;
139
140 J[0][0] = a;
141 J[0][1] = -1.0;
142 J[1][0] = 0.0;
143 J[1][1] = a;
144 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
145
146 PetscCall(VecRestoreArrayRead(U, &u));
147 PetscCall(VecRestoreArrayRead(Udot, &udot));
148
149 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
150 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
151 if (A != B) {
152 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
153 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
154 }
155 PetscFunctionReturn(PETSC_SUCCESS);
156 }
157
main(int argc,char ** argv)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;
169
170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 Initialize program
172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173 PetscFunctionBeginUser;
174 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
175 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
176
177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178 Create timestepping solver context
179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
181 PetscCall(TSSetType(ts, TSROSW));
182
183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184 Set ODE routines
185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186 PetscCall(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 PetscCall(PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL));
193 if (rhs_form) {
194 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
195 PetscCall(TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL));
196 } else {
197 PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL));
198 PetscCall(TSSetIJacobian(ts, NULL, NULL, IJacobian, NULL));
199 }
200
201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202 Set initial conditions
203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204 PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
205 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
206 PetscCall(VecSetUp(U));
207 PetscCall(VecGetArray(U, &u));
208 u[0] = 1.0 * rank;
209 u[1] = 20.0;
210 PetscCall(VecRestoreArray(U, &u));
211 PetscCall(TSSetSolution(ts, U));
212
213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214 Set solver options
215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216 PetscCall(TSSetSaveTrajectory(ts));
217 PetscCall(TSSetMaxTime(ts, 30.0));
218 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
219 PetscCall(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 PetscCall(TSGetAdapt(ts, &adapt));
224 PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC));
225 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));
226
227 /* Set direction and terminate flag for the event */
228 PetscCall(TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL));
229
230 PetscCall(TSSetFromOptions(ts));
231
232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233 Run timestepping solver
234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235 PetscCall(TSSolve(ts, U));
236
237 if (hist) { /* replay following history */
238 TSTrajectory tj;
239 PetscReal tf, t0, dt;
240
241 PetscCall(TSGetTime(ts, &tf));
242 PetscCall(TSSetMaxTime(ts, tf));
243 PetscCall(TSSetStepNumber(ts, 0));
244 PetscCall(TSRestartStep(ts));
245 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
246 PetscCall(TSSetFromOptions(ts));
247 PetscCall(TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL));
248 PetscCall(TSGetAdapt(ts, &adapt));
249 PetscCall(TSAdaptSetType(adapt, TSADAPTHISTORY));
250 PetscCall(TSGetTrajectory(ts, &tj));
251 PetscCall(TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE));
252 PetscCall(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 PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC));
256 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));
257 PetscCall(TSSetFromOptions(ts));
258 #endif
259 PetscCall(TSSetTime(ts, t0));
260 PetscCall(TSSetTimeStep(ts, dt));
261 PetscCall(TSResetTrajectory(ts));
262 PetscCall(VecGetArray(U, &u));
263 u[0] = 1.0 * rank;
264 u[1] = 20.0;
265 PetscCall(VecRestoreArray(U, &u));
266 PetscCall(TSSolve(ts, U));
267 }
268 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269 Free work space. All PETSc objects should be destroyed when they are no longer needed.
270 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271 PetscCall(VecDestroy(&U));
272 PetscCall(TSDestroy(&ts));
273
274 PetscCall(PetscFinalize());
275 return 0;
276 }
277
278 /*TEST
279
280 test:
281 suffix: a
282 nsize: 2
283 args: -ts_trajectory_type memory -snes_stol 1e-4
284 filter: sort -b
285
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
291
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
297
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
303
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
309
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
315
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
321
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
328
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
335
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
342
343 TEST*/
344