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 is one event set in this example, which 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. On reaching the limit on the number of ball bounces,
     11:   the TS run is requested to terminate from the PostEvent() callback.
     12: */

     14: #include <petscts.h>

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

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

     28:   PetscFunctionBeginUser;
     29:   /* Event for ball height */
     30:   PetscCall(TS2GetSolution(ts, &U, &V));
     31:   PetscCall(VecGetArrayRead(U, &u));
     32:   fvalue[0] = PetscRealPart(u[0]);
     33:   PetscCall(VecRestoreArrayRead(U, &u));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
     38: {
     39:   AppCtx      *app = (AppCtx *)ctx;
     40:   Vec          V;
     41:   PetscScalar *u, *v;
     42:   PetscMPIInt  rank;
     43:   PetscBool    inflag = PETSC_FALSE, outflag;

     45:   PetscFunctionBeginUser;
     46:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
     47:   if (nevents > 0) {
     48:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball hit the ground at t = %5.2f seconds\n", rank, (double)t));
     49:     /* Set new initial conditions with attenuation Cr */
     50:     PetscCall(TS2GetSolution(ts, &U, &V));
     51:     PetscCall(VecGetArray(U, &u));
     52:     PetscCall(VecGetArray(V, &v));
     53:     u[0] = 0.0;
     54:     v[0] = -app->Cr * v[0];
     55:     PetscCall(VecRestoreArray(U, &u));
     56:     PetscCall(VecRestoreArray(V, &v));
     57:     app->bounces++;
     58:   }
     59:   if (app->bounces >= app->maxbounces) { // 'app->bounces' may be different on different processes
     60:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces));
     61:     inflag = PETSC_TRUE; // current process requested to terminate
     62:   }
     63:   PetscCallMPI(MPIU_Allreduce(&inflag, &outflag, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)ts)));
     64:   if (outflag) PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER)); // request TS to terminate, sync on all ranks
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     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;

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

     83:   PetscCall(VecGetArray(F, &f));
     84:   f[0] = Res;
     85:   PetscCall(VecRestoreArray(F, &f));
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     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;

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

    105:   PetscCall(MatGetOwnershipRange(P, &i, NULL));
    106:   PetscCall(MatSetValue(P, i, i, Jac, INSERT_VALUES));
    107:   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
    108:   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
    109:   if (J != P) {
    110:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    111:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    112:   }
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    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[1];
    126:   PetscBool    terminate[1];
    127:   TSAdapt      adapt;

    129:   PetscFunctionBeginUser;
    130:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    131:   PetscCallMPI(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:   PetscCall(PetscOptionsReal("-Cd", "Drag coefficient", "", app.Cd, &app.Cd, NULL));
    139:   PetscCall(PetscOptionsReal("-Cr", "Restitution coefficient", "", app.Cr, &app.Cr, NULL));
    140:   PetscCall(PetscOptionsInt("-maxbounces", "Maximum number of bounces", "", app.maxbounces, &app.maxbounces, NULL));
    141:   PetscOptionsEnd();

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

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

    154:   direction[0] = -1;
    155:   terminate[0] = PETSC_FALSE;
    156:   PetscCall(TSSetEventHandler(ts, 1, direction, terminate, Event, PostEvent, &app)); // each process has one event-function defined

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

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

    177:   PetscCall(TS2SetSolution(ts, U, V));
    178:   PetscCall(TSSetFromOptions(ts));
    179:   PetscCall(TSSolve(ts, NULL));

    181:   PetscCall(VecDestroy(&U));
    182:   PetscCall(VecDestroy(&V));
    183:   PetscCall(TSDestroy(&ts));

    185:   PetscCall(PetscFinalize());
    186:   return 0;
    187: }

    189: /*TEST

    191:     test:
    192:       suffix: a
    193:       args: -ts_alpha_radius {{1.0 0.5}} -ts_max_time 50
    194:       output_file: output/ex44.out

    196:     test:
    197:       suffix: b
    198:       args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50
    199:       output_file: output/ex44.out

    201:     test:
    202:       suffix: bmf
    203:       args: -snes_mf_operator -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50
    204:       output_file: output/ex44.out

    206:     test:
    207:       suffix: 2
    208:       nsize: 2
    209:       args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50
    210:       output_file: output/ex44_2.out
    211:       filter: sort -b
    212:       filter_output: sort -b

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

    219: TEST*/