Actual source code: ex40.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  2: static char help[] = "Serial bouncing ball example to test TS event feature.\n";

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

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

 14: #include <petscts.h>

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

 23: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
 24: {
 25:   AppCtx            *app=(AppCtx*)ctx;
 26:   PetscErrorCode    ierr;
 27:   const PetscScalar *u;

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

 41: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
 42: {
 43:   AppCtx         *app=(AppCtx*)ctx;
 45:   PetscScalar    *u;

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

 65: /*
 66:      Defines the ODE passed to the ODE solver
 67: */
 68: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
 69: {
 70:   PetscErrorCode    ierr;
 71:   PetscScalar       *f;
 72:   const PetscScalar *u,*udot;

 75:   /*  The next three lines allow us to access the entries of the vectors directly */
 76:   VecGetArrayRead(U,&u);
 77:   VecGetArrayRead(Udot,&udot);
 78:   VecGetArray(F,&f);

 80:   f[0] = udot[0] - u[1];
 81:   f[1] = udot[1] + 9.8;

 83:   VecRestoreArrayRead(U,&u);
 84:   VecRestoreArrayRead(Udot,&udot);
 85:   VecRestoreArray(F,&f);
 86:   return(0);
 87: }

 91: /*
 92:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
 93: */
 94: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
 95: {
 96:   PetscErrorCode    ierr;
 97:   PetscInt          rowcol[] = {0,1};
 98:   PetscScalar       J[2][2];
 99:   const PetscScalar *u,*udot;

102:   VecGetArrayRead(U,&u);
103:   VecGetArrayRead(Udot,&udot);

105:   J[0][0] = a;                       J[0][1] = -1;
106:   J[1][0] = 0.0;                     J[1][1] = a;
107:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

109:   VecRestoreArrayRead(U,&u);
110:   VecRestoreArrayRead(Udot,&udot);

112:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
113:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
114:   if (A != B) {
115:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
116:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
117:   }
118:   return(0);
119: }

123: int main(int argc,char **argv)
124: {
125:   TS             ts;            /* ODE integrator */
126:   Vec            U;             /* solution will be stored here */
127:   Mat            A;             /* Jacobian matrix */
129:   PetscMPIInt    size;
130:   PetscInt       n = 2;
131:   PetscScalar    *u;
132:   AppCtx         app;
133:   PetscInt       direction[2];
134:   PetscBool      terminate[2];
135:   TSAdapt        adapt;

137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:      Initialize program
139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140:   PetscInitialize(&argc,&argv,(char*)0,help);
141:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
142:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

144:   app.nbounces = 0;
145:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");
146:   {
147:     app.maxbounces = 10;
148:     PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
149:   }
150:   PetscOptionsEnd();

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:     Create necessary matrix and vectors
154:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155:   MatCreate(PETSC_COMM_WORLD,&A);
156:   MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
157:   MatSetType(A,MATDENSE);
158:   MatSetFromOptions(A);
159:   MatSetUp(A);

161:   MatCreateVecs(A,&U,NULL);

163:   VecGetArray(U,&u);
164:   u[0] = 0.0;
165:   u[1] = 20.0;
166:   VecRestoreArray(U,&u);

168:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169:      Create timestepping solver context
170:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171:   TSCreate(PETSC_COMM_WORLD,&ts);
172:   TSSetSaveTrajectory(ts);
173:   TSSetProblemType(ts,TS_NONLINEAR);
174:   TSSetType(ts,TSROSW);
175:   TSSetIFunction(ts,NULL,(TSIFunction) IFunction,NULL);
176:   TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,NULL);

178:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179:      Set initial conditions
180:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181:   TSSetSolution(ts,U);

183:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184:      Set solver options
185:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186:   TSSetDuration(ts,1000,30.0);
187:   TSSetInitialTimeStep(ts,0.0,0.1);
188:   TSSetFromOptions(ts);
189: 
190:   /* Set directions and terminate flags for the two events */
191:   direction[0] = -1;            direction[1] = -1;
192:   terminate[0] = PETSC_FALSE;   terminate[1] = PETSC_TRUE;
193:   TSSetEventMonitor(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);

195:   TSGetAdapt(ts,&adapt);
196:   /* The adapative time step controller could take very large timesteps resulting in 
197:      the same event occuring multiple times in the same interval. A max. step 
198:      limit is enforced here to avoid this issue. 
199:   */
200:   TSAdaptSetStepLimits(adapt,0.0,0.5);
201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:      Run timestepping solver
203:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204:     TSSolve(ts,U);

206:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
208:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209: 
210:   MatDestroy(&A);
211:   VecDestroy(&U);
212:   TSDestroy(&ts);
213: 
214:   PetscFinalize();
215:   return(0);
216: }