Actual source code: ex40.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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:   PetscErrorCode    ierr;
 24:   const PetscScalar *u;

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

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

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

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

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

 71:   f[0] = u[1];
 72:   f[1] = - 9.8;

 74:   VecRestoreArrayRead(U,&u);
 75:   VecRestoreArray(F,&f);
 76:   return(0);
 77: }

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

 90:   VecGetArrayRead(U,&u);

 92:   J[0][0] = 0.0;     J[0][1] = 1.0;
 93:   J[1][0] = 0.0;     J[1][1] = 0.0;
 94:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

 96:   VecRestoreArrayRead(U,&u);

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

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

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

122:   f[0] = udot[0] - u[1];
123:   f[1] = udot[1] + 9.8;

125:   VecRestoreArrayRead(U,&u);
126:   VecRestoreArrayRead(Udot,&udot);
127:   VecRestoreArray(F,&f);
128:   return(0);
129: }

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

142:   VecGetArrayRead(U,&u);
143:   VecGetArrayRead(Udot,&udot);

145:   J[0][0] = a;      J[0][1] = -1.0;
146:   J[1][0] = 0.0;    J[1][1] = a;
147:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

149:   VecRestoreArrayRead(U,&u);
150:   VecRestoreArrayRead(Udot,&udot);

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

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

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:      Initialize program
177:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
179:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
180:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

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

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

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

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

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

245:   /* Set directions and terminate flags for the two events */
246:   direction[0] = -1;            direction[1] = -1;
247:   terminate[0] = PETSC_FALSE;   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 ierr;
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*/