Actual source code: ex40.c

  1: static char help[] = "Serial bouncing ball example to test TS event feature.\n";

  3: /*
  4:   The dynamics of the bouncing ball is described by the ODE
  5:                   u1_t = u2
  6:                   u2_t = -9.8

  8:   There are two events set in this example. The first one checks for the ball hitting the
  9:   ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
 10:   a factor of 0.9. The second event sets a limit on the number of ball bounces.
 11: */

 13: #include <petscts.h>

 15: typedef struct {
 16:   PetscInt maxbounces;
 17:   PetscInt nbounces;
 18: } AppCtx;

 20: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx)
 21: {
 22:   AppCtx            *app = (AppCtx *)ctx;
 23:   const PetscScalar *u;

 26:   /* Event for ball height */
 27:   VecGetArrayRead(U, &u);
 28:   fvalue[0] = u[0];
 29:   /* Event for number of bounces */
 30:   fvalue[1] = app->maxbounces - app->nbounces;
 31:   VecRestoreArrayRead(U, &u);
 32:   return 0;
 33: }

 35: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
 36: {
 37:   AppCtx      *app = (AppCtx *)ctx;
 38:   PetscScalar *u;

 41:   if (event_list[0] == 0) {
 42:     PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds\n", (double)t);
 43:     /* Set new initial conditions with .9 attenuation */
 44:     VecGetArray(U, &u);
 45:     u[0] = 0.0;
 46:     u[1] = -0.9 * u[1];
 47:     VecRestoreArray(U, &u);
 48:   } else if (event_list[0] == 1) {
 49:     PetscPrintf(PETSC_COMM_SELF, "Ball bounced %" PetscInt_FMT " times\n", app->nbounces);
 50:   }
 51:   app->nbounces++;
 52:   return 0;
 53: }

 55: /*
 56:      Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
 57: */
 58: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
 59: {
 60:   PetscScalar       *f;
 61:   const PetscScalar *u;

 64:   /*  The next three lines allow us to access the entries of the vectors directly */
 65:   VecGetArrayRead(U, &u);
 66:   VecGetArray(F, &f);

 68:   f[0] = u[1];
 69:   f[1] = -9.8;

 71:   VecRestoreArrayRead(U, &u);
 72:   VecRestoreArray(F, &f);
 73:   return 0;
 74: }

 76: /*
 77:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
 78: */
 79: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
 80: {
 81:   PetscInt           rowcol[] = {0, 1};
 82:   PetscScalar        J[2][2];
 83:   const PetscScalar *u;

 86:   VecGetArrayRead(U, &u);

 88:   J[0][0] = 0.0;
 89:   J[0][1] = 1.0;
 90:   J[1][0] = 0.0;
 91:   J[1][1] = 0.0;
 92:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);

 94:   VecRestoreArrayRead(U, &u);

 96:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 97:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 98:   if (A != B) {
 99:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
100:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
101:   }
102:   return 0;
103: }

105: /*
106:      Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
107: */
108: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
109: {
110:   PetscScalar       *f;
111:   const PetscScalar *u, *udot;

114:   /*  The next three lines allow us to access the entries of the vectors directly */
115:   VecGetArrayRead(U, &u);
116:   VecGetArrayRead(Udot, &udot);
117:   VecGetArray(F, &f);

119:   f[0] = udot[0] - u[1];
120:   f[1] = udot[1] + 9.8;

122:   VecRestoreArrayRead(U, &u);
123:   VecRestoreArrayRead(Udot, &udot);
124:   VecRestoreArray(F, &f);
125:   return 0;
126: }

128: /*
129:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
130: */
131: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
132: {
133:   PetscInt           rowcol[] = {0, 1};
134:   PetscScalar        J[2][2];
135:   const PetscScalar *u, *udot;

138:   VecGetArrayRead(U, &u);
139:   VecGetArrayRead(Udot, &udot);

141:   J[0][0] = a;
142:   J[0][1] = -1.0;
143:   J[1][0] = 0.0;
144:   J[1][1] = a;
145:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);

147:   VecRestoreArrayRead(U, &u);
148:   VecRestoreArrayRead(Udot, &udot);

150:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
151:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
152:   if (A != B) {
153:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
154:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
155:   }
156:   return 0;
157: }

159: int main(int argc, char **argv)
160: {
161:   TS           ts; /* ODE integrator */
162:   Vec          U;  /* solution will be stored here */
163:   PetscMPIInt  size;
164:   PetscInt     n = 2;
165:   PetscScalar *u;
166:   AppCtx       app;
167:   PetscInt     direction[2];
168:   PetscBool    terminate[2];
169:   PetscBool    rhs_form = PETSC_FALSE, hist = PETSC_TRUE;
170:   TSAdapt      adapt;

172:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173:      Initialize program
174:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176:   PetscInitialize(&argc, &argv, (char *)0, help);
177:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

180:   app.nbounces   = 0;
181:   app.maxbounces = 10;
182:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex40 options", "");
183:   PetscOptionsInt("-maxbounces", "", "", app.maxbounces, &app.maxbounces, NULL);
184:   PetscOptionsBool("-test_adapthistory", "", "", hist, &hist, NULL);
185:   PetscOptionsEnd();

187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188:      Create timestepping solver context
189:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190:   TSCreate(PETSC_COMM_WORLD, &ts);
191:   TSSetType(ts, TSROSW);

193:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194:      Set ODE routines
195:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196:   TSSetProblemType(ts, TS_NONLINEAR);
197:   /* Users are advised against the following branching and code duplication.
198:      For problems without a mass matrix like the one at hand, the RHSFunction
199:      (and companion RHSJacobian) interface is enough to support both explicit
200:      and implicit timesteppers. This tutorial example also deals with the
201:      IFunction/IJacobian interface for demonstration and testing purposes. */
202:   PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL);
203:   if (rhs_form) {
204:     TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
205:     TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL);
206:   } else {
207:     Mat A; /* Jacobian matrix */
208:     MatCreate(PETSC_COMM_WORLD, &A);
209:     MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
210:     MatSetType(A, MATDENSE);
211:     MatSetFromOptions(A);
212:     MatSetUp(A);
213:     TSSetIFunction(ts, NULL, IFunction, NULL);
214:     TSSetIJacobian(ts, A, A, IJacobian, NULL);
215:     MatDestroy(&A);
216:   }

218:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219:      Set initial conditions
220:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221:   VecCreate(PETSC_COMM_WORLD, &U);
222:   VecSetSizes(U, n, PETSC_DETERMINE);
223:   VecSetUp(U);
224:   VecGetArray(U, &u);
225:   u[0] = 0.0;
226:   u[1] = 20.0;
227:   VecRestoreArray(U, &u);
228:   TSSetSolution(ts, U);

230:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231:      Set solver options
232:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233:   if (hist) TSSetSaveTrajectory(ts);
234:   TSSetMaxTime(ts, 30.0);
235:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
236:   TSSetTimeStep(ts, 0.1);
237:   /* The adaptive time step controller could take very large timesteps resulting in
238:      the same event occurring multiple times in the same interval. A maximum step size
239:      limit is enforced here to avoid this issue. */
240:   TSGetAdapt(ts, &adapt);
241:   TSAdaptSetStepLimits(adapt, 0.0, 0.5);

243:   /* Set directions and terminate flags for the two events */
244:   direction[0] = -1;
245:   direction[1] = -1;
246:   terminate[0] = PETSC_FALSE;
247:   terminate[1] = PETSC_TRUE;
248:   TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&app);

250:   TSSetFromOptions(ts);

252:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253:      Run timestepping solver
254:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255:   TSSolve(ts, U);

257:   if (hist) { /* replay following history */
258:     TSTrajectory tj;
259:     PetscReal    tf, t0, dt;

261:     app.nbounces = 0;
262:     TSGetTime(ts, &tf);
263:     TSSetMaxTime(ts, tf);
264:     TSSetStepNumber(ts, 0);
265:     TSRestartStep(ts);
266:     TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
267:     TSSetFromOptions(ts);
268:     TSGetAdapt(ts, &adapt);
269:     TSAdaptSetType(adapt, TSADAPTHISTORY);
270:     TSGetTrajectory(ts, &tj);
271:     TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE);
272:     TSAdaptHistoryGetStep(adapt, 0, &t0, &dt);
273:     /* this example fails with single (or smaller) precision */
274: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
275:     TSAdaptSetType(adapt, TSADAPTBASIC);
276:     TSAdaptSetStepLimits(adapt, 0.0, 0.5);
277:     TSSetFromOptions(ts);
278: #endif
279:     TSSetTime(ts, t0);
280:     TSSetTimeStep(ts, dt);
281:     TSResetTrajectory(ts);
282:     VecGetArray(U, &u);
283:     u[0] = 0.0;
284:     u[1] = 20.0;
285:     VecRestoreArray(U, &u);
286:     TSSolve(ts, U);
287:   }
288:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
290:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291:   VecDestroy(&U);
292:   TSDestroy(&ts);

294:   PetscFinalize();
295:   return 0;
296: }

298: /*TEST

300:     test:
301:       suffix: a
302:       args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
303:       output_file: output/ex40.out

305:     test:
306:       suffix: b
307:       args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
308:       output_file: output/ex40.out

310:     test:
311:       suffix: c
312:       args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
313:       output_file: output/ex40.out

315:     test:
316:       suffix: d
317:       args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
318:       output_file: output/ex40.out

320:     test:
321:       suffix: e
322:       args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
323:       output_file: output/ex40.out

325:     test:
326:       suffix: f
327:       args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
328:       output_file: output/ex40.out

330:     test:
331:       suffix: g
332:       args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
333:       output_file: output/ex40.out

335:     test:
336:       suffix: h
337:       args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
338:       output_file: output/ex40.out

340:     test:
341:       suffix: i
342:       args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
343:       output_file: output/ex40.out

345:     test:
346:       suffix: j
347:       args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
348:       output_file: output/ex40.out

350:     test:
351:       suffix: k
352:       args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir
353:       output_file: output/ex40.out

355:     test:
356:       suffix: l
357:       args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir
358:       args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01
359:       output_file: output/ex40.out

361:     test:
362:       suffix: m
363:       args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false
364:       args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}}

366:     test:
367:       requires: !single
368:       suffix: n
369:       args: -test_adapthistory false
370:       args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view
371:       args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
372:       args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false

374:     test:
375:       requires: !single
376:       suffix: o
377:       args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir
378:       output_file: output/ex40.out
379: TEST*/