Actual source code: ex41.c

petsc-master 2017-02-25
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,"Processor [%d]: Ball hit the ground at t = %5.2f seconds\n",rank,(double)t);
 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
 51: */
 52: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
 53: {
 54:   PetscErrorCode    ierr;
 55:   PetscScalar       *f;
 56:   const PetscScalar *u,*udot;

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

 64:   f[0] = udot[0] - u[1];
 65:   f[1] = udot[1] + 9.8;

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

 73: /*
 74:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
 75: */
 76: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
 77: {
 78:   PetscErrorCode    ierr;
 79:   PetscInt          rowcol[2];
 80:   const PetscScalar *u,*udot;
 81:   PetscScalar       J[2][2];
 82:   PetscInt          rstart;

 85:   VecGetArrayRead(U,&u);
 86:   VecGetArrayRead(Udot,&udot);

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

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

 95:   VecRestoreArrayRead(U,&u);
 96:   VecRestoreArrayRead(Udot,&udot);

 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: int main(int argc,char **argv)
108: {
109:   TS             ts;            /* ODE integrator */
110:   Vec            U;             /* solution will be stored here */
111:   Mat            A;             /* Jacobian matrix */
113:   PetscMPIInt    rank;
114:   PetscInt       n = 2;
115:   PetscScalar    *u;
116:   PetscInt       direction=-1;
117:   PetscBool      terminate=PETSC_FALSE;
118:   TSAdapt        adapt;

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Initialize program
122:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
124:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:     Create necessary matrix and vectors
128:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129:   MatCreate(PETSC_COMM_WORLD,&A);
130:   MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
131:   MatSetType(A,MATDENSE);
132:   MatSetFromOptions(A);
133:   MatSetUp(A);

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

137:   VecGetArray(U,&u);
138:   u[0] = 1.0*rank;
139:   u[1] = 20.0;
140:   VecRestoreArray(U,&u);

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:      Create timestepping solver context
144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145:   TSCreate(PETSC_COMM_WORLD,&ts);
146:   TSSetSaveTrajectory(ts);
147:   TSSetProblemType(ts,TS_NONLINEAR);
148:   TSSetType(ts,TSROSW);
149:   TSSetIFunction(ts,NULL,(TSIFunction) IFunction,NULL);
150:   TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,NULL);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:      Set initial conditions
154:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155:   TSSetSolution(ts,U);

157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158:      Set solver options
159:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160:   TSSetDuration(ts,1000,30.0);
161:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
162:   TSSetInitialTimeStep(ts,0.0,0.1);

164:   TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);

166:   /* The adapative time step controller could take very large timesteps resulting in
167:      the same event occuring multiple times in the same interval. A maximum step size
168:      limit is enforced here to avoid this issue.
169:   */
170:   TSGetAdapt(ts,&adapt);
171:   TSAdaptSetStepLimits(adapt,0.0,0.5);

173:   TSSetFromOptions(ts);

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:      Run timestepping solver
177:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178:   TSSolve(ts,U);

180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
182:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183:   MatDestroy(&A);
184:   VecDestroy(&U);
185:   TSDestroy(&ts);

187:   PetscFinalize();
188:   return ierr;
189: }