Actual source code: ex40.c

      1: static char help[] = "Serial 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:   There is one event set in this example, which checks for the ball hitting the
      9:   ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
     10:   a factor of 0.9. 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:   PetscInt maxbounces;
     18:   PetscInt nbounces;
     19: } AppCtx;

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

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

     33: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
     34: {
     35:   AppCtx      *app = (AppCtx *)ctx;
     36:   PetscScalar *u;

     38:   PetscFunctionBeginUser;
     39:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds\n", (double)t));
     40:   /* Set new initial conditions with .9 attenuation */
     41:   PetscCall(VecGetArray(U, &u));
     42:   u[0] = 0.0;
     43:   u[1] = -0.9 * u[1];
     44:   PetscCall(VecRestoreArray(U, &u));
     45:   app->nbounces++;
     46:   if (app->nbounces >= app->maxbounces) {
     47:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball bounced %" PetscInt_FMT " times\n", app->nbounces));
     48:     PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER)); // request TS to terminate; since the program is serial, no need to sync this call
     49:   }
     50:   PetscFunctionReturn(PETSC_SUCCESS);
     51: }

     53: /*
     54:      Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
     55: */
     56: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
     57: {
     58:   PetscScalar       *f;
     59:   const PetscScalar *u;

     61:   PetscFunctionBeginUser;
     62:   /*  The following lines allow us to access the entries of the vectors directly */
     63:   PetscCall(VecGetArrayRead(U, &u));
     64:   PetscCall(VecGetArray(F, &f));

     66:   f[0] = u[1];
     67:   f[1] = -9.8;

     69:   PetscCall(VecRestoreArrayRead(U, &u));
     70:   PetscCall(VecRestoreArray(F, &f));
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /*
     75:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
     76: */
     77: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
     78: {
     79:   PetscInt           rowcol[] = {0, 1};
     80:   PetscScalar        J[2][2];
     81:   const PetscScalar *u;

     83:   PetscFunctionBeginUser;
     84:   PetscCall(VecGetArrayRead(U, &u));

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

     92:   PetscCall(VecRestoreArrayRead(U, &u));

     94:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     95:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     96:   if (A != B) {
     97:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     98:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     99:   }
    100:   PetscFunctionReturn(PETSC_SUCCESS);
    101: }

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

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

    117:   f[0] = udot[0] - u[1];
    118:   f[1] = udot[1] + 9.8;

    120:   PetscCall(VecRestoreArrayRead(U, &u));
    121:   PetscCall(VecRestoreArrayRead(Udot, &udot));
    122:   PetscCall(VecRestoreArray(F, &f));
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

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

    135:   PetscFunctionBeginUser;
    136:   PetscCall(VecGetArrayRead(U, &u));
    137:   PetscCall(VecGetArrayRead(Udot, &udot));

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

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

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

    157: int main(int argc, char **argv)
    158: {
    159:   TS           ts; /* ODE integrator */
    160:   Vec          U;  /* solution will be stored here */
    161:   PetscMPIInt  size;
    162:   PetscInt     n = 2;
    163:   PetscScalar *u;
    164:   AppCtx       app;
    165:   PetscInt     direction[1];
    166:   PetscBool    terminate[1];
    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_size(PETSC_COMM_WORLD, &size));
    176:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    178:   app.nbounces   = 0;
    179:   app.maxbounces = 10;
    180:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex40 options", "");
    181:   PetscCall(PetscOptionsInt("-maxbounces", "", "", app.maxbounces, &app.maxbounces, NULL));
    182:   PetscCall(PetscOptionsBool("-test_adapthistory", "", "", hist, &hist, NULL));
    183:   PetscOptionsEnd();

    185:   Mat A; /* Jacobian matrix */
    186:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    187:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
    188:   PetscCall(MatSetType(A, MATDENSE));
    189:   PetscCall(MatSetFromOptions(A));
    190:   PetscCall(MatSetUp(A));
    191:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    192:      Create timestepping solver context
    193:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    194:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    195:   PetscCall(TSSetType(ts, TSROSW));

    197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    198:      Set ODE routines
    199:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    200:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    201:   /* Users are advised against the following branching and code duplication.
    202:      For problems without a mass matrix like the one at hand, the RHSFunction
    203:      (and companion RHSJacobian) interface is enough to support both explicit
    204:      and implicit timesteppers. This tutorial example also deals with the
    205:      IFunction/IJacobian interface for demonstration and testing purposes. */
    206:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL));
    207:   if (rhs_form) {
    208:     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
    209:     PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
    210:   } else {
    211:     PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL));
    212:     PetscCall(TSSetIJacobian(ts, A, A, IJacobian, NULL));
    213:   }

    215:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    216:      Set initial conditions
    217:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    218:   PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
    219:   PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
    220:   PetscCall(VecSetUp(U));
    221:   PetscCall(VecGetArray(U, &u));
    222:   u[0] = 0.0;
    223:   u[1] = 20.0;
    224:   PetscCall(VecRestoreArray(U, &u));
    225:   PetscCall(TSSetSolution(ts, U));

    227:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    228:      Set solver options
    229:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    230:   if (hist) PetscCall(TSSetSaveTrajectory(ts));
    231:   PetscCall(TSSetMaxTime(ts, 30.0));
    232:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    233:   PetscCall(TSSetTimeStep(ts, 0.1));
    234:   /* The adaptive time step controller could take very large timesteps
    235:      jumping over the next event zero-crossing point. A maximum step size
    236:      limit is enforced here to avoid this issue. */
    237:   PetscCall(TSGetAdapt(ts, &adapt));
    238:   PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));

    240:   /* Set direction and terminate flag for the event */
    241:   direction[0] = -1;
    242:   terminate[0] = PETSC_FALSE;
    243:   PetscCall(TSSetEventHandler(ts, 1, direction, terminate, EventFunction, PostEventFunction, (void *)&app));

    245:   PetscCall(TSSetFromOptions(ts));

    247:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    248:      Run timestepping solver
    249:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    250:   PetscCall(TSSolve(ts, U));

    252:   if (hist) { /* replay following history */
    253:     TSTrajectory tj;
    254:     PetscReal    tf, t0, dt;

    256:     app.nbounces = 0;
    257:     PetscCall(TSGetTime(ts, &tf));
    258:     PetscCall(TSSetMaxTime(ts, tf));
    259:     PetscCall(TSSetStepNumber(ts, 0));
    260:     PetscCall(TSRestartStep(ts));
    261:     PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    262:     PetscCall(TSSetFromOptions(ts));
    263:     PetscCall(TSGetAdapt(ts, &adapt));
    264:     PetscCall(TSAdaptSetType(adapt, TSADAPTHISTORY));
    265:     PetscCall(TSGetTrajectory(ts, &tj));
    266:     PetscCall(TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE));
    267:     PetscCall(TSAdaptHistoryGetStep(adapt, 0, &t0, &dt));
    268:     /* this example fails with single (or smaller) precision */
    269: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
    270:     /*
    271:        In the first TSSolve() the final time 'tf' is the event location found after a few event handler iterations.
    272:        If 'tf' is set as the max time for the second run, the TS solver may approach this point by
    273:        slightly different steps, resulting in a slightly different solution and fvalue[] at 'tf',
    274:        so that the event may not be triggered at 'tf' anymore. Fix: apply safety factor 1.05
    275:     */
    276:     PetscCall(TSSetMaxTime(ts, tf * 1.05));
    277:     PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC));
    278:     PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5));
    279:     PetscCall(TSSetFromOptions(ts));
    280: #endif
    281:     PetscCall(TSSetTime(ts, t0));
    282:     PetscCall(TSSetTimeStep(ts, dt));
    283:     PetscCall(TSResetTrajectory(ts));
    284:     PetscCall(VecGetArray(U, &u));
    285:     u[0] = 0.0;
    286:     u[1] = 20.0;
    287:     PetscCall(VecRestoreArray(U, &u));
    288:     PetscCall(TSSolve(ts, U));
    289:   }
    290:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    291:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    292:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    293:   PetscCall(MatDestroy(&A));
    294:   PetscCall(VecDestroy(&U));
    295:   PetscCall(TSDestroy(&ts));

    297:   PetscCall(PetscFinalize());
    298:   return 0;
    299: }

    301: /*TEST

    303:     test:
    304:       suffix: a
    305:       args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
    306:       output_file: output/ex40.out

    308:     test:
    309:       suffix: b
    310:       args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
    311:       output_file: output/ex40.out

    313:     test:
    314:       suffix: c
    315:       args: -snes_mf_operator -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
    316:       output_file: output/ex40.out

    318:     test:
    319:       suffix: cr
    320:       args: -rhs-form -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_cr_dir
    321:       output_file: output/ex40.out

    323:     test:
    324:       suffix: crmf
    325:       args: -rhs-form -snes_mf_operator -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_crmf_dir
    326:       output_file: output/ex40.out

    328:     test:
    329:       suffix: d
    330:       args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
    331:       output_file: output/ex40.out

    333:     test:
    334:       suffix: e
    335:       args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
    336:       output_file: output/ex40.out

    338:     test:
    339:       suffix: f
    340:       args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
    341:       output_file: output/ex40.out

    343:     test:
    344:       suffix: g
    345:       args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
    346:       output_file: output/ex40.out

    348:     test:
    349:       suffix: h
    350:       args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
    351:       output_file: output/ex40.out

    353:     test:
    354:       suffix: i
    355:       args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
    356:       output_file: output/ex40.out

    358:     test:
    359:       suffix: j
    360:       args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
    361:       output_file: output/ex40.out

    363:     test:
    364:       suffix: k
    365:       args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir
    366:       output_file: output/ex40.out

    368:     test:
    369:       suffix: l
    370:       args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir
    371:       args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01
    372:       output_file: output/ex40.out

    374:     test:
    375:       suffix: m
    376:       args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false
    377:       args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}}

    379:     test:
    380:       requires: !single
    381:       suffix: n
    382:       args: -test_adapthistory false
    383:       args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view
    384:       args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
    385:       args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false

    387:     test:
    388:       requires: !single
    389:       suffix: o
    390:       args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir
    391:       output_file: output/ex40.out
    392: TEST*/