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;

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

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

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

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

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

 65:   f[0] = u[1];
 66:   f[1] = - 9.8;

 68:   VecRestoreArrayRead(U,&u);
 69:   VecRestoreArray(F,&f);
 70:   return 0;
 71: }

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

 82:   VecGetArrayRead(U,&u);

 84:   J[0][0] = 0.0;     J[0][1] = 1.0;
 85:   J[1][0] = 0.0;     J[1][1] = 0.0;
 86:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

 88:   VecRestoreArrayRead(U,&u);

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

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

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

112:   f[0] = udot[0] - u[1];
113:   f[1] = udot[1] + 9.8;

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

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

130:   VecGetArrayRead(U,&u);
131:   VecGetArrayRead(Udot,&udot);

133:   J[0][0] = a;      J[0][1] = -1.0;
134:   J[1][0] = 0.0;    J[1][1] = a;
135:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

137:   VecRestoreArrayRead(U,&u);
138:   VecRestoreArrayRead(Udot,&udot);

140:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
141:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
142:   if (A != B) {
143:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
144:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
145:   }
146:   return 0;
147: }

149: int main(int argc,char **argv)
150: {
151:   TS             ts;            /* ODE integrator */
152:   Vec            U;             /* solution will be stored here */
154:   PetscMPIInt    size;
155:   PetscInt       n = 2;
156:   PetscScalar    *u;
157:   AppCtx         app;
158:   PetscInt       direction[2];
159:   PetscBool      terminate[2];
160:   PetscBool      rhs_form=PETSC_FALSE,hist=PETSC_TRUE;
161:   TSAdapt        adapt;

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Initialize program
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166:   PetscInitialize(&argc,&argv,(char*)0,help);
167:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

170:   app.nbounces = 0;
171:   app.maxbounces = 10;
172:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex40 options","");
173:   PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
174:   PetscOptionsBool("-test_adapthistory","","",hist,&hist,NULL);
175:   PetscOptionsEnd();

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:     Mat A; /* Jacobian matrix */
198:     MatCreate(PETSC_COMM_WORLD,&A);
199:     MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
200:     MatSetType(A,MATDENSE);
201:     MatSetFromOptions(A);
202:     MatSetUp(A);
203:     TSSetIFunction(ts,NULL,IFunction,NULL);
204:     TSSetIJacobian(ts,A,A,IJacobian,NULL);
205:     MatDestroy(&A);
206:   }

208:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209:      Set initial conditions
210:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
211:   VecCreate(PETSC_COMM_WORLD,&U);
212:   VecSetSizes(U,n,PETSC_DETERMINE);
213:   VecSetUp(U);
214:   VecGetArray(U,&u);
215:   u[0] = 0.0;
216:   u[1] = 20.0;
217:   VecRestoreArray(U,&u);
218:   TSSetSolution(ts,U);

220:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221:      Set solver options
222:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
223:   if (hist) TSSetSaveTrajectory(ts);
224:   TSSetMaxTime(ts,30.0);
225:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
226:   TSSetTimeStep(ts,0.1);
227:   /* The adaptive time step controller could take very large timesteps resulting in
228:      the same event occurring multiple times in the same interval. A maximum step size
229:      limit is enforced here to avoid this issue. */
230:   TSGetAdapt(ts,&adapt);
231:   TSAdaptSetStepLimits(adapt,0.0,0.5);

233:   /* Set directions and terminate flags for the two events */
234:   direction[0] = -1;            direction[1] = -1;
235:   terminate[0] = PETSC_FALSE;   terminate[1] = PETSC_TRUE;
236:   TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);

238:   TSSetFromOptions(ts);

240:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241:      Run timestepping solver
242:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243:   TSSolve(ts,U);

245:   if (hist) { /* replay following history */
246:     TSTrajectory tj;
247:     PetscReal    tf,t0,dt;

249:     app.nbounces = 0;
250:     TSGetTime(ts,&tf);
251:     TSSetMaxTime(ts,tf);
252:     TSSetStepNumber(ts,0);
253:     TSRestartStep(ts);
254:     TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
255:     TSSetFromOptions(ts);
256:     TSGetAdapt(ts,&adapt);
257:     TSAdaptSetType(adapt,TSADAPTHISTORY);
258:     TSGetTrajectory(ts,&tj);
259:     TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE);
260:     TSAdaptHistoryGetStep(adapt,0,&t0,&dt);
261:     /* this example fails with single (or smaller) precision */
262: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL__FP16)
263:     TSAdaptSetType(adapt,TSADAPTBASIC);
264:     TSAdaptSetStepLimits(adapt,0.0,0.5);
265:     TSSetFromOptions(ts);
266: #endif
267:     TSSetTime(ts,t0);
268:     TSSetTimeStep(ts,dt);
269:     TSResetTrajectory(ts);
270:     VecGetArray(U,&u);
271:     u[0] = 0.0;
272:     u[1] = 20.0;
273:     VecRestoreArray(U,&u);
274:     TSSolve(ts,U);
275:   }
276:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
278:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279:   VecDestroy(&U);
280:   TSDestroy(&ts);

282:   PetscFinalize();
283:   return 0;
284: }

286: /*TEST

288:     test:
289:       suffix: a
290:       args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
291:       output_file: output/ex40.out

293:     test:
294:       suffix: b
295:       args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
296:       output_file: output/ex40.out

298:     test:
299:       suffix: c
300:       args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
301:       output_file: output/ex40.out

303:     test:
304:       suffix: d
305:       args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
306:       output_file: output/ex40.out

308:     test:
309:       suffix: e
310:       args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
311:       output_file: output/ex40.out

313:     test:
314:       suffix: f
315:       args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
316:       output_file: output/ex40.out

318:     test:
319:       suffix: g
320:       args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
321:       output_file: output/ex40.out

323:     test:
324:       suffix: h
325:       args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
326:       output_file: output/ex40.out

328:     test:
329:       suffix: i
330:       args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
331:       output_file: output/ex40.out

333:     test:
334:       suffix: j
335:       args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
336:       output_file: output/ex40.out

338:     test:
339:       suffix: k
340:       args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir
341:       output_file: output/ex40.out

343:     test:
344:       suffix: l
345:       args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir
346:       args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01
347:       output_file: output/ex40.out

349:     test:
350:       suffix: m
351:       args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false
352:       args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}}

354:     test:
355:       requires: !single
356:       suffix: n
357:       args: -test_adapthistory false
358:       args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view
359:       args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
360:       args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false

362:     test:
363:       requires: !single
364:       suffix: o
365:       args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir
366:       output_file: output/ex40.out
367: TEST*/