Actual source code: ex51.c

      1: static char help[] = "Small ODE to test TS accuracy.\n";

      3: /*
      4:   The ODE
      5:                   u1_t = cos(t),
      6:                   u2_t = sin(u2)
      7:   with analytical solution
      8:                   u1(t) = sin(t),
      9:                   u2(t) = 2 * atan(exp(t) * tan(0.5))
     10:   is used to test the accuracy of TS schemes.
     11: */

     13: #include <petscts.h>

     15: /*
     16:      Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
     17: */
     18: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *s)
     19: {
     20:   PetscScalar       *f;
     21:   const PetscScalar *u;

     23:   PetscFunctionBeginUser;
     24:   PetscCall(VecGetArrayRead(U, &u));
     25:   PetscCall(VecGetArray(F, &f));

     27:   f[0] = PetscCosReal(t);
     28:   f[1] = PetscSinScalar(u[1]);

     30:   PetscCall(VecRestoreArrayRead(U, &u));
     31:   PetscCall(VecRestoreArray(F, &f));
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: /*
     36:      Defines the exact solution.
     37: */
     38: static PetscErrorCode ExactSolution(PetscReal t, Vec U)
     39: {
     40:   PetscScalar *u;

     42:   PetscFunctionBeginUser;
     43:   PetscCall(VecGetArray(U, &u));

     45:   u[0] = PetscSinReal(t);
     46:   u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5));

     48:   PetscCall(VecRestoreArray(U, &u));
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: int main(int argc, char **argv)
     53: {
     54:   TS           ts;  /* ODE integrator */
     55:   Vec          U;   /* numerical solution will be stored here */
     56:   Vec          Uex; /* analytical (exact) solution will be stored here */
     57:   PetscMPIInt  size;
     58:   PetscInt     n = 2;
     59:   PetscScalar *u;
     60:   PetscReal    t, final_time = 1.0, dt = 0.25;
     61:   PetscReal    error;
     62:   TSAdapt      adapt;

     64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     65:      Initialize program
     66:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     67:   PetscFunctionBeginUser;
     68:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     69:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     70:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

     72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     73:      Create timestepping solver context
     74:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     75:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     76:   PetscCall(TSSetType(ts, TSROSW));

     78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     79:      Set ODE routines
     80:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     81:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
     82:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));

     84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     85:      Set initial conditions
     86:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     87:   PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
     88:   PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
     89:   PetscCall(VecSetUp(U));
     90:   PetscCall(VecGetArray(U, &u));
     91:   u[0] = 0.0;
     92:   u[1] = 1.0;
     93:   PetscCall(VecRestoreArray(U, &u));
     94:   PetscCall(TSSetSolution(ts, U));

     96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     97:      Set solver options
     98:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     99:   PetscCall(TSSetSaveTrajectory(ts));
    100:   PetscCall(TSSetMaxTime(ts, final_time));
    101:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    102:   PetscCall(TSSetTimeStep(ts, dt));
    103:   /* The adaptive time step controller is forced to take constant time steps. */
    104:   PetscCall(TSGetAdapt(ts, &adapt));
    105:   PetscCall(TSAdaptSetType(adapt, TSADAPTNONE));

    107:   PetscCall(TSSetFromOptions(ts));

    109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    110:      Run timestepping solver and compute error
    111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    112:   PetscCall(TSSolve(ts, U));
    113:   PetscCall(TSGetTime(ts, &t));

    115:   if (PetscAbsReal(t - final_time) > 100 * PETSC_MACHINE_EPSILON) {
    116:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference of %g between the prescribed final time %g and the actual final time.\n", (double)(final_time - t), (double)final_time));
    117:   }
    118:   PetscCall(VecDuplicate(U, &Uex));
    119:   PetscCall(ExactSolution(t, Uex));

    121:   PetscCall(VecAYPX(Uex, -1.0, U));
    122:   PetscCall(VecNorm(Uex, NORM_2, &error));
    123:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error at final time: %.2E\n", (double)error));

    125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    126:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    128:   PetscCall(VecDestroy(&U));
    129:   PetscCall(VecDestroy(&Uex));
    130:   PetscCall(TSDestroy(&ts));

    132:   PetscCall(PetscFinalize());
    133:   return 0;
    134: }

    136: /*TEST

    138:     test:
    139:       suffix: 3bs
    140:       args: -ts_type rk -ts_rk_type 3bs
    141:       requires: !single

    143:     test:
    144:       suffix: 5bs
    145:       args: -ts_type rk -ts_rk_type 5bs
    146:       requires: !single

    148:     test:
    149:       suffix: 5dp
    150:       args: -ts_type rk -ts_rk_type 5dp
    151:       requires: !single

    153:     test:
    154:       suffix: 6vr
    155:       args: -ts_type rk -ts_rk_type 6vr
    156:       requires: !single

    158:     test:
    159:       suffix: 7vr
    160:       args: -ts_type rk -ts_rk_type 7vr
    161:       requires: !single

    163:     test:
    164:       suffix: 8vr
    165:       args: -ts_type rk -ts_rk_type 8vr
    166:       requires: !single

    168: TEST*/