Actual source code: ex41.c

      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, PetscReal *fvalue, void *ctx)
     18: {
     19:   const PetscScalar *u;

     21:   PetscFunctionBeginUser;
     22:   /* Event for ball height */
     23:   PetscCall(VecGetArrayRead(U, &u));
     24:   fvalue[0] = PetscRealPart(u[0]);
     25:   PetscCall(VecRestoreArrayRead(U, &u));
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
     30: {
     31:   PetscScalar *u;
     32:   PetscMPIInt  rank;

     34:   PetscFunctionBeginUser;
     35:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
     36:   if (nevents) {
     37:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds -> Processor[%d]\n", (double)t, rank));
     38:     /* Set new initial conditions with .9 attenuation */
     39:     PetscCall(VecGetArray(U, &u));
     40:     u[0] = 1.0 * rank;
     41:     u[1] = -0.9 * u[1];
     42:     PetscCall(VecRestoreArray(U, &u));
     43:   }
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: /*
     48:      Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
     49: */
     50: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
     51: {
     52:   PetscScalar       *f;
     53:   const PetscScalar *u;

     55:   PetscFunctionBeginUser;
     56:   /*  The next three lines allow us to access the entries of the vectors directly */
     57:   PetscCall(VecGetArrayRead(U, &u));
     58:   PetscCall(VecGetArray(F, &f));

     60:   f[0] = u[1];
     61:   f[1] = -9.8;

     63:   PetscCall(VecRestoreArrayRead(U, &u));
     64:   PetscCall(VecRestoreArray(F, &f));
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: /*
     69:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian.
     70: */
     71: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
     72: {
     73:   PetscInt           rowcol[2], rstart;
     74:   PetscScalar        J[2][2];
     75:   const PetscScalar *u;

     77:   PetscFunctionBeginUser;
     78:   PetscCall(VecGetArrayRead(U, &u));

     80:   PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
     81:   rowcol[0] = rstart;
     82:   rowcol[1] = rstart + 1;

     84:   J[0][0] = 0.0;
     85:   J[0][1] = 1.0;
     86:   J[1][0] = 0.0;
     87:   J[1][1] = 0.0;
     88:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));

     90:   PetscCall(VecRestoreArrayRead(U, &u));
     91:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     92:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     93:   if (A != B) {
     94:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     95:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     96:   }
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: /*
    101:      Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
    102: */
    103: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
    104: {
    105:   PetscScalar       *f;
    106:   const PetscScalar *u, *udot;

    108:   PetscFunctionBeginUser;
    109:   /*  The next three lines allow us to access the entries of the vectors directly */
    110:   PetscCall(VecGetArrayRead(U, &u));
    111:   PetscCall(VecGetArrayRead(Udot, &udot));
    112:   PetscCall(VecGetArray(F, &f));

    114:   f[0] = udot[0] - u[1];
    115:   f[1] = udot[1] + 9.8;

    117:   PetscCall(VecRestoreArrayRead(U, &u));
    118:   PetscCall(VecRestoreArrayRead(Udot, &udot));
    119:   PetscCall(VecRestoreArray(F, &f));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: /*
    124:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
    125: */
    126: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
    127: {
    128:   PetscInt           rowcol[2], rstart;
    129:   PetscScalar        J[2][2];
    130:   const PetscScalar *u, *udot;

    132:   PetscFunctionBeginUser;
    133:   PetscCall(VecGetArrayRead(U, &u));
    134:   PetscCall(VecGetArrayRead(Udot, &udot));

    136:   PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
    137:   rowcol[0] = rstart;
    138:   rowcol[1] = rstart + 1;

    140:   J[0][0] = a;
    141:   J[0][1] = -1.0;
    142:   J[1][0] = 0.0;
    143:   J[1][1] = a;
    144:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));

    146:   PetscCall(VecRestoreArrayRead(U, &u));
    147:   PetscCall(VecRestoreArrayRead(Udot, &udot));

    149:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    150:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    151:   if (A != B) {
    152:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    153:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    154:   }
    155:   PetscFunctionReturn(PETSC_SUCCESS);
    156: }

    158: int main(int argc, char **argv)
    159: {
    160:   TS           ts; /* ODE integrator */
    161:   Vec          U;  /* solution will be stored here */
    162:   PetscMPIInt  rank;
    163:   PetscInt     n = 2;
    164:   PetscScalar *u;
    165:   PetscInt     direction = -1;
    166:   PetscBool    terminate = PETSC_FALSE;
    167:   PetscBool    rhs_form = PETSC_FALSE, hist = PETSC_TRUE;
    168:   TSAdapt      adapt;

    170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    171:      Initialize program
    172:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    173:   PetscFunctionBeginUser;
    174:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    175:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

    177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    178:      Create timestepping solver context
    179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    180:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    181:   PetscCall(TSSetType(ts, TSROSW));

    183:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    184:      Set ODE routines
    185:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    186:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    187:   /* Users are advised against the following branching and code duplication.
    188:      For problems without a mass matrix like the one at hand, the RHSFunction
    189:      (and companion RHSJacobian) interface is enough to support both explicit
    190:      and implicit timesteppers. This tutorial example also deals with the
    191:      IFunction/IJacobian interface for demonstration and testing purposes. */
    192:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL));
    193:   if (rhs_form) {
    194:     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
    195:     PetscCall(TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL));
    196:   } else {
    197:     PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL));
    198:     PetscCall(TSSetIJacobian(ts, NULL, NULL, IJacobian, NULL));
    199:   }

    201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    202:      Set initial conditions
    203:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    204:   PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
    205:   PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
    206:   PetscCall(VecSetUp(U));
    207:   PetscCall(VecGetArray(U, &u));
    208:   u[0] = 1.0 * rank;
    209:   u[1] = 20.0;
    210:   PetscCall(VecRestoreArray(U, &u));
    211:   PetscCall(TSSetSolution(ts, U));

    213:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    214:      Set solver options
    215:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    216:   PetscCall(TSSetSaveTrajectory(ts));
    217:   PetscCall(TSSetMaxTime(ts, 30.0));
    218:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    219:   PetscCall(TSSetTimeStep(ts, 0.1));
    220:   /* The adaptive time step controller could take very large timesteps resulting in
    221:      the same event occurring multiple times in the same interval. A maximum step size
    222:      limit is enforced here to avoid this issue. */
    223:   PetscCall(TSGetAdapt(ts, &adapt));
    224:   PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC));
    225:   PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));

    227:   /* Set direction and terminate flag for the event */
    228:   PetscCall(TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL));

    230:   PetscCall(TSSetFromOptions(ts));

    232:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    233:      Run timestepping solver
    234:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    235:   PetscCall(TSSolve(ts, U));

    237:   if (hist) { /* replay following history */
    238:     TSTrajectory tj;
    239:     PetscReal    tf, t0, dt;

    241:     PetscCall(TSGetTime(ts, &tf));
    242:     PetscCall(TSSetMaxTime(ts, tf));
    243:     PetscCall(TSSetStepNumber(ts, 0));
    244:     PetscCall(TSRestartStep(ts));
    245:     PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    246:     PetscCall(TSSetFromOptions(ts));
    247:     PetscCall(TSSetEventHandler(ts, 1, &direction, &terminate, EventFunction, PostEventFunction, NULL));
    248:     PetscCall(TSGetAdapt(ts, &adapt));
    249:     PetscCall(TSAdaptSetType(adapt, TSADAPTHISTORY));
    250:     PetscCall(TSGetTrajectory(ts, &tj));
    251:     PetscCall(TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE));
    252:     PetscCall(TSAdaptHistoryGetStep(adapt, 0, &t0, &dt));
    253:     /* this example fails with single (or smaller) precision */
    254: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
    255:     PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC));
    256:     PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));
    257:     PetscCall(TSSetFromOptions(ts));
    258: #endif
    259:     PetscCall(TSSetTime(ts, t0));
    260:     PetscCall(TSSetTimeStep(ts, dt));
    261:     PetscCall(TSResetTrajectory(ts));
    262:     PetscCall(VecGetArray(U, &u));
    263:     u[0] = 1.0 * rank;
    264:     u[1] = 20.0;
    265:     PetscCall(VecRestoreArray(U, &u));
    266:     PetscCall(TSSolve(ts, U));
    267:   }
    268:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    269:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    270:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    271:   PetscCall(VecDestroy(&U));
    272:   PetscCall(TSDestroy(&ts));

    274:   PetscCall(PetscFinalize());
    275:   return 0;
    276: }

    278: /*TEST

    280:    test:
    281:       suffix: a
    282:       nsize: 2
    283:       args: -ts_trajectory_type memory -snes_stol 1e-4
    284:       filter: sort -b

    286:    test:
    287:       suffix: b
    288:       nsize: 2
    289:       args: -ts_trajectory_type memory -ts_type arkimex -snes_stol 1e-4
    290:       filter: sort -b

    292:    test:
    293:       suffix: c
    294:       nsize: 2
    295:       args: -ts_trajectory_type memory -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
    296:       filter: sort -b

    298:    test:
    299:       suffix: d
    300:       nsize: 2
    301:       args: -ts_trajectory_type memory -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
    302:       filter: sort -b

    304:    test:
    305:       suffix: e
    306:       nsize: 2
    307:       args: -ts_trajectory_type memory -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
    308:       filter: sort -b

    310:    test:
    311:       suffix: f
    312:       nsize: 2
    313:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 3bs
    314:       filter: sort -b

    316:    test:
    317:       suffix: g
    318:       nsize: 2
    319:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 5bs
    320:       filter: sort -b

    322:    test:
    323:       suffix: h
    324:       nsize: 2
    325:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 6vr
    326:       filter: sort -b
    327:       output_file: output/ex41_g.out

    329:    test:
    330:       suffix: i
    331:       nsize: 2
    332:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 7vr
    333:       filter: sort -b
    334:       output_file: output/ex41_g.out

    336:    test:
    337:       suffix: j
    338:       nsize: 2
    339:       args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 8vr
    340:       filter: sort -b
    341:       output_file: output/ex41_g.out

    343: TEST*/