Actual source code: ex41.c

petsc-master 2016-07-30
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>

 19: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
 20: {
 21:   PetscErrorCode    ierr;
 22:   const PetscScalar *u;

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

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

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

 55: /*
 56:      Defines the ODE passed to the ODE solver
 57: */
 58: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
 59: {
 60:   PetscErrorCode    ierr;
 61:   PetscScalar       *f;
 62:   const PetscScalar *u,*udot;

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

 70:   f[0] = udot[0] - u[1];
 71:   f[1] = udot[1] + 9.8;

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

 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[2];
 88:   const PetscScalar *u,*udot;
 89:   PetscScalar       J[2][2];
 90:   PetscInt          rstart;

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

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

 99:   J[0][0] = a;                       J[0][1] = -1;
100:   J[1][0] = 0.0;                     J[1][1] = a;
101:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

103:   VecRestoreArrayRead(U,&u);
104:   VecRestoreArrayRead(Udot,&udot);

106:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108:   if (A != B) {
109:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
110:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
111:   }
112:   return(0);
113: }

117: int main(int argc,char **argv)
118: {
119:   TS             ts;            /* ODE integrator */
120:   Vec            U;             /* solution will be stored here */
121:   Mat            A;             /* Jacobian matrix */
123:   PetscMPIInt    rank;
124:   PetscInt       n = 2;
125:   PetscScalar    *u;
126:   PetscInt       direction=-1;
127:   PetscBool      terminate=PETSC_FALSE;
128:   TSAdapt        adapt;

130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:      Initialize program
132:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
134:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

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

147:   VecGetArray(U,&u);
148:   u[0] = 1.0*rank;
149:   u[1] = 20.0;
150:   VecRestoreArray(U,&u);

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

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

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

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

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

183:   TSSetFromOptions(ts);

185:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186:      Run timestepping solver
187:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188:   TSSolve(ts,U);

190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
192:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193:   MatDestroy(&A);
194:   VecDestroy(&U);
195:   TSDestroy(&ts);

197:   PetscFinalize();
198:   return ierr;
199: }