Actual source code: ex40.c

petsc-master 2017-02-19
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
 59: */
 60: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
 61: {
 62:   PetscErrorCode    ierr;
 63:   PetscScalar       *f;
 64:   const PetscScalar *u,*udot;

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

 72:   f[0] = udot[0] - u[1];
 73:   f[1] = udot[1] + 9.8;

 75:   VecRestoreArrayRead(U,&u);
 76:   VecRestoreArrayRead(Udot,&udot);
 77:   VecRestoreArray(F,&f);
 78:   return(0);
 79: }

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

 92:   VecGetArrayRead(U,&u);
 93:   VecGetArrayRead(Udot,&udot);

 95:   J[0][0] = a;                       J[0][1] = -1;
 96:   J[1][0] = 0.0;                     J[1][1] = a;
 97:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

 99:   VecRestoreArrayRead(U,&u);
100:   VecRestoreArrayRead(Udot,&udot);

102:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
103:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
104:   if (A != B) {
105:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
106:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
107:   }
108:   return(0);
109: }

111: int main(int argc,char **argv)
112: {
113:   TS             ts;            /* ODE integrator */
114:   Vec            U;             /* solution will be stored here */
115:   Mat            A;             /* Jacobian matrix */
117:   PetscMPIInt    size;
118:   PetscInt       n = 2;
119:   PetscScalar    *u;
120:   AppCtx         app;
121:   PetscInt       direction[2];
122:   PetscBool      terminate[2];
123:   TSAdapt        adapt;

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:      Initialize program
127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
129:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
130:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

132:   app.nbounces = 0;
133:   app.maxbounces = 10;
134:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");
135:   PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
136:   PetscOptionsEnd();

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:     Create necessary matrix and vectors
140:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:   MatCreate(PETSC_COMM_WORLD,&A);
142:   MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
143:   MatSetType(A,MATDENSE);
144:   MatSetFromOptions(A);
145:   MatSetUp(A);

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

149:   VecGetArray(U,&u);
150:   u[0] = 0.0;
151:   u[1] = 20.0;
152:   VecRestoreArray(U,&u);

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:      Create timestepping solver context
156:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157:   TSCreate(PETSC_COMM_WORLD,&ts);
158:   TSSetSaveTrajectory(ts);
159:   TSSetProblemType(ts,TS_NONLINEAR);
160:   TSSetType(ts,TSROSW);
161:   TSSetIFunction(ts,NULL,(TSIFunction) IFunction,NULL);
162:   TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,NULL);

164:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165:      Set initial conditions
166:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167:   TSSetSolution(ts,U);

169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:      Set solver options
171:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172:   TSSetDuration(ts,1000,30.0);
173:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
174:   TSSetInitialTimeStep(ts,0.0,0.1);

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

181:   /* The adapative time step controller could take very large timesteps resulting in
182:      the same event occuring multiple times in the same interval. A maximum step size
183:      limit is enforced here to avoid this issue.
184:   */
185:   TSGetAdapt(ts,&adapt);
186:   TSAdaptSetStepLimits(adapt,0.0,0.5);

188:   TSSetFromOptions(ts);

190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191:      Run timestepping solver
192:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193:   TSSolve(ts,U);

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198:   MatDestroy(&A);
199:   VecDestroy(&U);
200:   TSDestroy(&ts);

202:   PetscFinalize();
203:   return ierr;
204: }