Actual source code: ex44.c

  1: static char help[] = "Parallel bouncing ball example formulated as a second-order system to test TS event feature.\n";

  3: /*
  4:   The dynamics of the bouncing ball with drag coefficient Cd is described by the ODE

  6:       u_tt = -9.8 - 1/2 Cd (u_t)^2 sign(u_t)

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

 13: #include <petscts.h>

 15: typedef struct {
 16:   PetscReal Cd; /* drag coefficient */
 17:   PetscReal Cr; /* restitution coefficient */
 18:   PetscInt  bounces;
 19:   PetscInt  maxbounces;
 20: } AppCtx;

 22: static PetscErrorCode Event(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx)
 23: {
 24:   AppCtx            *app = (AppCtx *)ctx;
 25:   Vec                V;
 26:   const PetscScalar *u, *v;

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

 41: static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
 42: {
 43:   AppCtx      *app = (AppCtx *)ctx;
 44:   Vec          V;
 45:   PetscScalar *u, *v;
 46:   PetscMPIInt  rank;

 49:   if (!nevents) return 0;
 50:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 51:   if (event_list[0] == 0) {
 52:     PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball hit the ground at t = %5.2f seconds\n", rank, (double)t);
 53:     /* Set new initial conditions with .9 attenuation */
 54:     TS2GetSolution(ts, &U, &V);
 55:     VecGetArray(U, &u);
 56:     VecGetArray(V, &v);
 57:     u[0] = 0.0;
 58:     v[0] = -app->Cr * v[0];
 59:     VecRestoreArray(U, &u);
 60:     VecRestoreArray(V, &v);
 61:     app->bounces++;
 62:   } else if (event_list[0] == 1) {
 63:     PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces);
 64:   }
 65:   return 0;
 66: }

 68: static PetscErrorCode I2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F, void *ctx)
 69: {
 70:   AppCtx            *app = (AppCtx *)ctx;
 71:   const PetscScalar *u, *v, *a;
 72:   PetscScalar        Res, *f;

 75:   VecGetArrayRead(U, &u);
 76:   VecGetArrayRead(V, &v);
 77:   VecGetArrayRead(A, &a);
 78:   Res = a[0] + 9.8 + 0.5 * app->Cd * v[0] * v[0] * PetscSignReal(PetscRealPart(v[0]));
 79:   VecRestoreArrayRead(U, &u);
 80:   VecRestoreArrayRead(V, &v);
 81:   VecRestoreArrayRead(A, &a);

 83:   VecGetArray(F, &f);
 84:   f[0] = Res;
 85:   VecRestoreArray(F, &f);
 86:   return 0;
 87: }

 89: static PetscErrorCode I2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx)
 90: {
 91:   AppCtx            *app = (AppCtx *)ctx;
 92:   const PetscScalar *u, *v, *a;
 93:   PetscInt           i;
 94:   PetscScalar        Jac;

 97:   VecGetArrayRead(U, &u);
 98:   VecGetArrayRead(V, &v);
 99:   VecGetArrayRead(A, &a);
100:   Jac = shiftA + shiftV * app->Cd * v[0];
101:   VecRestoreArrayRead(U, &u);
102:   VecRestoreArrayRead(V, &v);
103:   VecRestoreArrayRead(A, &a);

105:   MatGetOwnershipRange(P, &i, NULL);
106:   MatSetValue(P, i, i, Jac, INSERT_VALUES);
107:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
108:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
109:   if (J != P) {
110:     MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
111:     MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
112:   }
113:   return 0;
114: }

116: int main(int argc, char **argv)
117: {
118:   TS           ts;   /* ODE integrator */
119:   Vec          U, V; /* solution will be stored here */
120:   Vec          F;    /* residual vector */
121:   Mat          J;    /* Jacobian matrix */
122:   PetscMPIInt  rank;
123:   PetscScalar *u, *v;
124:   AppCtx       app;
125:   PetscInt     direction[2];
126:   PetscBool    terminate[2];
127:   TSAdapt      adapt;

130:   PetscInitialize(&argc, &argv, NULL, help);
131:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

133:   app.Cd         = 0.0;
134:   app.Cr         = 0.9;
135:   app.bounces    = 0;
136:   app.maxbounces = 10;
137:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex44 options", "");
138:   PetscOptionsReal("-Cd", "Drag coefficient", "", app.Cd, &app.Cd, NULL);
139:   PetscOptionsReal("-Cr", "Restitution coefficient", "", app.Cr, &app.Cr, NULL);
140:   PetscOptionsInt("-maxbounces", "Maximum number of bounces", "", app.maxbounces, &app.maxbounces, NULL);
141:   PetscOptionsEnd();

143:   TSCreate(PETSC_COMM_WORLD, &ts);
144:   /*TSSetSaveTrajectory(ts);*/
145:   TSSetProblemType(ts, TS_NONLINEAR);
146:   TSSetType(ts, TSALPHA2);

148:   TSSetMaxTime(ts, PETSC_INFINITY);
149:   TSSetTimeStep(ts, 0.1);
150:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
151:   TSGetAdapt(ts, &adapt);
152:   TSAdaptSetStepLimits(adapt, 0.0, 0.5);

154:   direction[0] = -1;
155:   terminate[0] = PETSC_FALSE;
156:   direction[1] = -1;
157:   terminate[1] = PETSC_TRUE;
158:   TSSetEventHandler(ts, 2, direction, terminate, Event, PostEvent, &app);

160:   MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &J);
161:   MatSetFromOptions(J);
162:   MatSetUp(J);
163:   MatCreateVecs(J, NULL, &F);
164:   TSSetI2Function(ts, F, I2Function, &app);
165:   TSSetI2Jacobian(ts, J, J, I2Jacobian, &app);
166:   VecDestroy(&F);
167:   MatDestroy(&J);

169:   TSGetI2Jacobian(ts, &J, NULL, NULL, NULL);
170:   MatCreateVecs(J, &U, NULL);
171:   MatCreateVecs(J, &V, NULL);
172:   VecGetArray(U, &u);
173:   VecGetArray(V, &v);
174:   u[0] = 5.0 * rank;
175:   v[0] = 20.0;
176:   VecRestoreArray(U, &u);
177:   VecRestoreArray(V, &v);

179:   TS2SetSolution(ts, U, V);
180:   TSSetFromOptions(ts);
181:   TSSolve(ts, NULL);

183:   VecDestroy(&U);
184:   VecDestroy(&V);
185:   TSDestroy(&ts);

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

191: /*TEST

193:     test:
194:       suffix: a
195:       args: -ts_alpha_radius {{1.0 0.5}}
196:       output_file: output/ex44.out

198:     test:
199:       suffix: b
200:       args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic
201:       output_file: output/ex44.out

203:     test:
204:       suffix: 2
205:       nsize: 2
206:       args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic
207:       output_file: output/ex44_2.out
208:       filter: sort -b
209:       filter_output: sort -b

211:     test:
212:       requires: !single
213:       args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
214:       args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false

216: TEST*/