Actual source code: ex20.c

      1: static char help[] = "Solves the van der Pol equation.\n\
      2: Input parameters include:\n";

      4: /* ------------------------------------------------------------------------

      6:    This program solves the van der Pol DAE ODE equivalent
      7:        y' = z                 (1)
      8:        z' = \mu ((1-y^2)z-y)
      9:    on the domain 0 <= x <= 1, with the boundary conditions
     10:        y(0) = 2, y'(0) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2),
     11:    and
     12:        \mu = 10^6 ( y'(0) ~ -0.6666665432100101).
     13:    This is a nonlinear equation. The well prepared initial condition gives errors that are not dominated by the first few steps of the method when \mu is large.

     15:    Notes:
     16:    This code demonstrates the TS solver interface to an ODE -- RHSFunction for explicit form and IFunction for implicit form.

     18:   ------------------------------------------------------------------------- */

     20: #include <petscts.h>

     22: typedef struct _n_User *User;
     23: struct _n_User {
     24:   PetscReal mu;
     25:   PetscReal next_output;
     26: };

     28: /*
     29:    User-defined routines
     30: */
     31: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ctx)
     32: {
     33:   User               user = (User)ctx;
     34:   PetscScalar       *f;
     35:   const PetscScalar *x;

     37:   PetscFunctionBeginUser;
     38:   PetscCall(VecGetArrayRead(X, &x));
     39:   PetscCall(VecGetArray(F, &f));
     40:   f[0] = x[1];
     41:   f[1] = user->mu * (1. - x[0] * x[0]) * x[1] - x[0];
     42:   PetscCall(VecRestoreArrayRead(X, &x));
     43:   PetscCall(VecRestoreArray(F, &f));
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
     48: {
     49:   User               user = (User)ctx;
     50:   const PetscScalar *x, *xdot;
     51:   PetscScalar       *f;

     53:   PetscFunctionBeginUser;
     54:   PetscCall(VecGetArrayRead(X, &x));
     55:   PetscCall(VecGetArrayRead(Xdot, &xdot));
     56:   PetscCall(VecGetArray(F, &f));
     57:   f[0] = xdot[0] - x[1];
     58:   f[1] = xdot[1] - user->mu * ((1.0 - x[0] * x[0]) * x[1] - x[0]);
     59:   PetscCall(VecRestoreArrayRead(X, &x));
     60:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
     61:   PetscCall(VecRestoreArray(F, &f));
     62:   PetscFunctionReturn(PETSC_SUCCESS);
     63: }

     65: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
     66: {
     67:   User               user     = (User)ctx;
     68:   PetscInt           rowcol[] = {0, 1};
     69:   const PetscScalar *x;
     70:   PetscScalar        J[2][2];

     72:   PetscFunctionBeginUser;
     73:   PetscCall(VecGetArrayRead(X, &x));
     74:   J[0][0] = a;
     75:   J[0][1] = -1.0;
     76:   J[1][0] = user->mu * (2.0 * x[0] * x[1] + 1.0);
     77:   J[1][1] = a - user->mu * (1.0 - x[0] * x[0]);
     78:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
     79:   PetscCall(VecRestoreArrayRead(X, &x));

     81:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     82:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     83:   if (A != B) {
     84:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     85:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     86:   }
     87:   PetscFunctionReturn(PETSC_SUCCESS);
     88: }

     90: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
     91: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx)
     92: {
     93:   const PetscScalar *x;
     94:   PetscReal          tfinal, dt;
     95:   User               user = (User)ctx;
     96:   Vec                interpolatedX;

     98:   PetscFunctionBeginUser;
     99:   PetscCall(TSGetTimeStep(ts, &dt));
    100:   PetscCall(TSGetMaxTime(ts, &tfinal));

    102:   while (user->next_output <= t && user->next_output <= tfinal) {
    103:     PetscCall(VecDuplicate(X, &interpolatedX));
    104:     PetscCall(TSInterpolate(ts, user->next_output, interpolatedX));
    105:     PetscCall(VecGetArrayRead(interpolatedX, &x));
    106:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
    107:     PetscCall(VecRestoreArrayRead(interpolatedX, &x));
    108:     PetscCall(VecDestroy(&interpolatedX));
    109:     user->next_output += 0.1;
    110:   }
    111:   PetscFunctionReturn(PETSC_SUCCESS);
    112: }

    114: int main(int argc, char **argv)
    115: {
    116:   TS             ts; /* nonlinear solver */
    117:   Vec            x;  /* solution, residual vectors */
    118:   Mat            A;  /* Jacobian matrix */
    119:   PetscInt       steps;
    120:   PetscReal      ftime   = 0.5;
    121:   PetscBool      monitor = PETSC_FALSE, implicitform = PETSC_TRUE;
    122:   PetscScalar   *x_ptr;
    123:   PetscMPIInt    size;
    124:   struct _n_User user;

    126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    127:      Initialize program
    128:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    129:   PetscFunctionBeginUser;
    130:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    131:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    132:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

    134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    135:     Set runtime options
    136:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    137:   user.next_output = 0.0;
    138:   user.mu          = 1.0e3;
    139:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
    140:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &implicitform, NULL));
    141:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Physical parameters", NULL);
    142:   PetscCall(PetscOptionsReal("-mu", "Stiffness parameter", "<1.0e6>", user.mu, &user.mu, NULL));
    143:   PetscOptionsEnd();

    145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    146:     Create necessary matrix and vectors, solve same ODE on every process
    147:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    148:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    149:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
    150:   PetscCall(MatSetFromOptions(A));
    151:   PetscCall(MatSetUp(A));

    153:   PetscCall(MatCreateVecs(A, &x, NULL));

    155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    156:      Create timestepping solver context
    157:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    158:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    159:   if (implicitform) {
    160:     PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
    161:     PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user));
    162:     PetscCall(TSSetType(ts, TSBEULER));
    163:   } else {
    164:     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
    165:     PetscCall(TSSetType(ts, TSRK));
    166:   }
    167:   PetscCall(TSSetMaxTime(ts, ftime));
    168:   PetscCall(TSSetTimeStep(ts, 0.001));
    169:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    170:   if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));

    172:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    173:      Set initial conditions
    174:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    175:   PetscCall(VecGetArray(x, &x_ptr));
    176:   x_ptr[0] = 2.0;
    177:   x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
    178:   PetscCall(VecRestoreArray(x, &x_ptr));

    180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    181:      Set runtime options
    182:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    183:   PetscCall(TSSetFromOptions(ts));

    185:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    186:      Solve nonlinear system
    187:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    188:   PetscCall(TSSolve(ts, x));
    189:   PetscCall(TSGetSolveTime(ts, &ftime));
    190:   PetscCall(TSGetStepNumber(ts, &steps));
    191:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "steps %" PetscInt_FMT ", ftime %g\n", steps, (double)ftime));
    192:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));

    194:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    195:      Free work space.  All PETSc objects should be destroyed when they
    196:      are no longer needed.
    197:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    198:   PetscCall(MatDestroy(&A));
    199:   PetscCall(VecDestroy(&x));
    200:   PetscCall(TSDestroy(&ts));

    202:   PetscCall(PetscFinalize());
    203:   return 0;
    204: }

    206: /*TEST

    208:     test:
    209:       requires: !single
    210:       args: -mu 1e6

    212:     test:
    213:       requires: !single
    214:       suffix: 2
    215:       args: -implicitform false -ts_type rk -ts_rk_type 5dp -ts_adapt_type dsp

    217:     test:
    218:       requires: !single
    219:       suffix: 3
    220:       args: -implicitform false -ts_type rk -ts_rk_type 5dp -ts_adapt_type dsp -ts_adapt_dsp_filter H0312

    222: TEST*/