Actual source code: ex41.c

petsc-master 2019-06-15
Report Typos and Errors
  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:   PetscErrorCode    ierr;
 20:   const PetscScalar *u;

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

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

 37:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 38:   if (nevents) {
 39:     PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %5.2f seconds -> Processor[%d]\n",(double)t,rank);
 40:     /* Set new initial conditions with .9 attenuation */
 41:     VecGetArray(U,&u);
 42:     u[0] =  1.0*rank;
 43:     u[1] = -0.9*u[1];
 44:     VecRestoreArray(U,&u);
 45:   }
 46:   return(0);
 47: }

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

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

 63:   f[0] = u[1];
 64:   f[1] = - 9.8;

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

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

 82:   VecGetArrayRead(U,&u);

 84:   MatGetOwnershipRange(A,&rstart,NULL);
 85:   rowcol[0] = rstart; rowcol[1] = rstart+1;

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

 91:   VecRestoreArrayRead(U,&u);

 93:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 94:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 95:   if (A != B) {
 96:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 97:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 98:   }
 99:   return(0);
100: }

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

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

117:   f[0] = udot[0] - u[1];
118:   f[1] = udot[1] + 9.8;

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

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

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

140:   MatGetOwnershipRange(A,&rstart,NULL);
141:   rowcol[0] = rstart; rowcol[1] = rstart+1;

143:   J[0][0] = a;        J[0][1] = -1.0;
144:   J[1][0] = 0.0;      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 */
164:   PetscMPIInt    rank;
165:   PetscInt       n = 2;
166:   PetscScalar    *u;
167:   PetscInt       direction=-1;
168:   PetscBool      terminate=PETSC_FALSE;
169:   PetscBool      rhs_form=PETSC_FALSE,hist=PETSC_TRUE;
170:   TSAdapt        adapt;

172:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173:      Initialize program
174:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
176:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

184:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185:      Set ODE routines
186:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187:   TSSetProblemType(ts,TS_NONLINEAR);
188:   /* Users are advised against the following branching and code duplication.
189:      For problems without a mass matrix like the one at hand, the RHSFunction
190:      (and companion RHSJacobian) interface is enough to support both explicit
191:      and implicit timesteppers. This tutorial example also deals with the
192:      IFunction/IJacobian interface for demonstration and testing purposes. */
193:   PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL);
194:   if (rhs_form) {
195:     TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
196:     TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL);
197:   } else {
198:     TSSetIFunction(ts,NULL,IFunction,NULL);
199:     TSSetIJacobian(ts,NULL,NULL,IJacobian,NULL);
200:   }

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

214:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215:      Set solver options
216:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217:   TSSetSaveTrajectory(ts);
218:   TSSetMaxTime(ts,30.0);
219:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
220:   TSSetTimeStep(ts,0.1);
221:   /* The adapative time step controller could take very large timesteps resulting in
222:      the same event occuring multiple times in the same interval. A maximum step size
223:      limit is enforced here to avoid this issue. */
224:   TSGetAdapt(ts,&adapt);
225:   TSAdaptSetType(adapt,TSADAPTBASIC);
226:   TSAdaptSetStepLimits(adapt,0.0,0.5);

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

231:   TSSetFromOptions(ts);

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

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

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

275:   PetscFinalize();
276:   return ierr;
277: }

279: /*TEST

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

287:    test:
288:       suffix: b
289:       nsize: 2
290:       args: -ts_trajectory_type memory -ts_type arkimex -snes_stol 1e-4
291:       filter: sort -b

293:    test:
294:       suffix: c
295:       nsize: 2
296:       args: -ts_trajectory_type memory -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
297:       filter: sort -b

299:    test:
300:       suffix: d
301:       nsize: 2
302:       args: -ts_trajectory_type memory -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
303:       filter: sort -b

305:    test:
306:       suffix: e
307:       nsize: 2
308:       args: -ts_trajectory_type memory -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
309:       filter: sort -b

311:    test:
312:       suffix: f
313:       nsize: 2
314:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 3bs
315:       filter: sort -b

317:    test:
318:       suffix: g
319:       nsize: 2
320:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 5bs
321:       filter: sort -b

323:    test:
324:       suffix: h
325:       nsize: 2
326:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 6vr
327:       filter: sort -b
328:       output_file: output/ex41_g.out

330:    test:
331:       suffix: i
332:       nsize: 2
333:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 7vr
334:       filter: sort -b
335:       output_file: output/ex41_g.out

337:    test:
338:       suffix: j
339:       nsize: 2
340:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 8vr
341:       filter: sort -b
342:       output_file: output/ex41_g.out

344: TEST*/