Actual source code: ex41.c

  1: static char help[] = "Parallel 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:   Each processor is assigned one ball.

 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: */

 15: #include <petscts.h>

 17: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx)
 18: {
 19:   const PetscScalar *u;

 22:   /* Event for ball height */
 23:   VecGetArrayRead(U, &u);
 24:   fvalue[0] = u[0];
 25:   VecRestoreArrayRead(U, &u);
 26:   return 0;
 27: }

 29: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
 30: {
 31:   PetscScalar *u;
 32:   PetscMPIInt  rank;

 35:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 36:   if (nevents) {
 37:     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:     VecGetArray(U, &u);
 40:     u[0] = 1.0 * rank;
 41:     u[1] = -0.9 * u[1];
 42:     VecRestoreArray(U, &u);
 43:   }
 44:   return 0;
 45: }

 47: /*
 48:      Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
 49: */
 50: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
 51: {
 52:   PetscScalar       *f;
 53:   const PetscScalar *u;

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

 60:   f[0] = u[1];
 61:   f[1] = -9.8;

 63:   VecRestoreArrayRead(U, &u);
 64:   VecRestoreArray(F, &f);
 65:   return 0;
 66: }

 68: /*
 69:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian.
 70: */
 71: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
 72: {
 73:   PetscInt           rowcol[2], rstart;
 74:   PetscScalar        J[2][2];
 75:   const PetscScalar *u;

 78:   VecGetArrayRead(U, &u);

 80:   MatGetOwnershipRange(B, &rstart, NULL);
 81:   rowcol[0] = rstart;
 82:   rowcol[1] = rstart + 1;

 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:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);

 90:   VecRestoreArrayRead(U, &u);
 91:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 92:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 93:   if (A != B) {
 94:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 95:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 96:   }
 97:   return 0;
 98: }

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

109:   /*  The next three lines allow us to access the entries of the vectors directly */
110:   VecGetArrayRead(U, &u);
111:   VecGetArrayRead(Udot, &udot);
112:   VecGetArray(F, &f);

114:   f[0] = udot[0] - u[1];
115:   f[1] = udot[1] + 9.8;

117:   VecRestoreArrayRead(U, &u);
118:   VecRestoreArrayRead(Udot, &udot);
119:   VecRestoreArray(F, &f);
120:   return 0;
121: }

123: /*
124:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
125: */
126: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
127: {
128:   PetscInt           rowcol[2], rstart;
129:   PetscScalar        J[2][2];
130:   const PetscScalar *u, *udot;

133:   VecGetArrayRead(U, &u);
134:   VecGetArrayRead(Udot, &udot);

136:   MatGetOwnershipRange(B, &rstart, NULL);
137:   rowcol[0] = rstart;
138:   rowcol[1] = rstart + 1;

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

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

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

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;

170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171:      Initialize program
172:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174:   PetscInitialize(&argc, &argv, (char *)0, help);
175:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Create timestepping solver context
179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180:   TSCreate(PETSC_COMM_WORLD, &ts);
181:   TSSetType(ts, TSROSW);

183:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184:      Set ODE routines
185:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186:   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:   PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL);
193:   if (rhs_form) {
194:     TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
195:     TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL);
196:   } else {
197:     TSSetIFunction(ts, NULL, IFunction, NULL);
198:     TSSetIJacobian(ts, NULL, NULL, IJacobian, NULL);
199:   }

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:      Set initial conditions
203:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204:   VecCreate(PETSC_COMM_WORLD, &U);
205:   VecSetSizes(U, n, PETSC_DETERMINE);
206:   VecSetUp(U);
207:   VecGetArray(U, &u);
208:   u[0] = 1.0 * rank;
209:   u[1] = 20.0;
210:   VecRestoreArray(U, &u);
211:   TSSetSolution(ts, U);

213:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214:      Set solver options
215:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216:   TSSetSaveTrajectory(ts);
217:   TSSetMaxTime(ts, 30.0);
218:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
219:   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:   TSGetAdapt(ts, &adapt);
224:   TSAdaptSetType(adapt, TSADAPTBASIC);
225:   TSAdaptSetStepLimits(adapt, 0.0, 0.5);

227:   /* Set direction and terminate flag for the event */
228:   TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL);

230:   TSSetFromOptions(ts);

232:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233:      Run timestepping solver
234:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235:   TSSolve(ts, U);

237:   if (hist) { /* replay following history */
238:     TSTrajectory tj;
239:     PetscReal    tf, t0, dt;

241:     TSGetTime(ts, &tf);
242:     TSSetMaxTime(ts, tf);
243:     TSSetStepNumber(ts, 0);
244:     TSRestartStep(ts);
245:     TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
246:     TSSetFromOptions(ts);
247:     TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL);
248:     TSGetAdapt(ts, &adapt);
249:     TSAdaptSetType(adapt, TSADAPTHISTORY);
250:     TSGetTrajectory(ts, &tj);
251:     TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE);
252:     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:     TSAdaptSetType(adapt, TSADAPTBASIC);
256:     TSAdaptSetStepLimits(adapt, 0.0, 0.5);
257:     TSSetFromOptions(ts);
258: #endif
259:     TSSetTime(ts, t0);
260:     TSSetTimeStep(ts, dt);
261:     TSResetTrajectory(ts);
262:     VecGetArray(U, &u);
263:     u[0] = 1.0 * rank;
264:     u[1] = 20.0;
265:     VecRestoreArray(U, &u);
266:     TSSolve(ts, U);
267:   }
268:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
270:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271:   VecDestroy(&U);
272:   TSDestroy(&ts);

274:   PetscFinalize();
275:   return 0;
276: }

278: /*TEST

280:    test:
281:       suffix: a
282:       nsize: 2
283:       args: -ts_trajectory_type memory -snes_stol 1e-4
284:       filter: sort -b

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

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

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

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

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

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

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

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

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

343: TEST*/