Actual source code: ct_vdp_imex.c

      1: /*
      2:  * ex_vdp.c
      3:  *
      4:  *  Created on: Jun 1, 2012
      5:  *      Author: Hong Zhang
      6:  */
      7: static char help[] = "Solves the van der Pol equation. \n Input parameters include:\n";

      9: /*
     10:  * This program solves the van der Pol equation
     11:  * y' = z                               (1)
     12:  * z' = (((1-y^2)*z-y)/eps              (2)
     13:  * on the domain 0<=x<=0.5, with the initial conditions
     14:  * y(0) = 2,
     15:  * z(0) = -2/3 + 10/81*eps - 292/2187*eps^2-1814/19683*eps^3
     16:  * IMEX schemes are applied to the split equation
     17:  * [y'] = [z]  + [0                 ]
     18:  * [z']   [0]    [(((1-y^2)*z-y)/eps]
     19:  *
     20:  * F(x)= [z]
     21:  *       [0]
     22:  *
     23:  * G(x)= [y'] -   [0                 ]
     24:  *       [z']     [(((1-y^2)*z-y)/eps]
     25:  *
     26:  * JG(x) =  G_x + a G_xdot
     27:  */

     29: #include <petscdmda.h>
     30: #include <petscts.h>

     32: typedef struct _User *User;
     33: struct _User {
     34:   PetscReal mu; /*stiffness control coefficient: epsilon*/
     35: };

     37: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
     38: static PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     39: static PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);

     41: int main(int argc, char **argv)
     42: {
     43:   TS           ts;
     44:   Vec          x; /* solution vector */
     45:   Mat          A; /* Jacobian */
     46:   PetscInt     steps, mx, eimex_rowcol[2], two;
     47:   PetscScalar *x_ptr;
     48:   PetscReal    ftime, dt, norm;
     49:   Vec          ref;
     50:   struct _User user; /* user-defined work context */
     51:   PetscViewer  viewer;

     53:   PetscFunctionBeginUser;
     54:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     55:   /* Initialize user application context */
     56:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "van der Pol options", "");
     57:   user.mu = 1e0;
     58:   PetscCall(PetscOptionsReal("-eps", "Stiffness controller", "", user.mu, &user.mu, NULL));
     59:   PetscOptionsEnd();

     61:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     62:    Set runtime options
     63:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     64:   /*
     65:    PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL));
     66:    */

     68:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     69:    Create necessary matrix and vectors, solve same ODE on every process
     70:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     71:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     72:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
     73:   PetscCall(MatSetFromOptions(A));
     74:   PetscCall(MatSetUp(A));
     75:   PetscCall(MatCreateVecs(A, &x, NULL));

     77:   PetscCall(MatCreateVecs(A, &ref, NULL));
     78:   PetscCall(VecGetArray(ref, &x_ptr));
     79:   /*
     80:    * [0,1], mu=10^-3
     81:    */
     82:   /*
     83:    x_ptr[0] = -1.8881254106283;
     84:    x_ptr[1] =  0.7359074233370;*/

     86:   /*
     87:    * [0,0.5],mu=10^-3
     88:    */
     89:   /*
     90:    x_ptr[0] = 1.596980778659137;
     91:    x_ptr[1] = -1.029103015879544;
     92:    */
     93:   /*
     94:    * [0,0.5],mu=1
     95:    */
     96:   x_ptr[0] = 1.619084329683235;
     97:   x_ptr[1] = -0.803530465176385;

     99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    100:    Create timestepping solver context
    101:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    102:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    103:   PetscCall(TSSetType(ts, TSEIMEX));
    104:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
    105:   PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
    106:   PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user));

    108:   dt    = 0.00001;
    109:   ftime = 1.1;
    110:   PetscCall(TSSetTimeStep(ts, dt));
    111:   PetscCall(TSSetMaxTime(ts, ftime));
    112:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    114:    Set initial conditions
    115:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    116:   PetscCall(VecGetArray(x, &x_ptr));
    117:   x_ptr[0] = 2.;
    118:   x_ptr[1] = -2. / 3. + 10. / 81. * (user.mu) - 292. / 2187. * (user.mu) * (user.mu) - 1814. / 19683. * (user.mu) * (user.mu) * (user.mu);
    119:   PetscCall(TSSetSolution(ts, x));
    120:   PetscCall(VecGetSize(x, &mx));

    122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    123:    Set runtime options
    124:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    125:   PetscCall(TSSetFromOptions(ts));

    127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    128:    Solve nonlinear system
    129:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    130:   PetscCall(TSSolve(ts, x));
    131:   PetscCall(TSGetTime(ts, &ftime));
    132:   PetscCall(TSGetStepNumber(ts, &steps));

    134:   PetscCall(VecAXPY(x, -1.0, ref));
    135:   PetscCall(VecNorm(x, NORM_2, &norm));
    136:   PetscCall(TSGetTimeStep(ts, &dt));

    138:   eimex_rowcol[0] = 0;
    139:   eimex_rowcol[1] = 0;
    140:   two             = 2;
    141:   PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-ts_eimex_row_col", eimex_rowcol, &two, NULL));
    142:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "order %11s %18s %37s\n", "dt", "norm", "final solution components 0 and 1"));
    143:   PetscCall(VecGetArray(x, &x_ptr));
    144:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(%" PetscInt_FMT ",%" PetscInt_FMT ") %10.8f %18.15f %18.15f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm, (double)PetscRealPart(x_ptr[0]), (double)PetscRealPart(x_ptr[1])));
    145:   PetscCall(VecRestoreArray(x, &x_ptr));

    147:   /* Write line in convergence log */
    148:   PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
    149:   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
    150:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_APPEND));
    151:   PetscCall(PetscViewerFileSetName(viewer, "eimex_nonstiff_vdp.txt"));
    152:   PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %10.8f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm));
    153:   PetscCall(PetscViewerDestroy(&viewer));

    155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    156:    Free work space.
    157:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    158:   PetscCall(MatDestroy(&A));
    159:   PetscCall(VecDestroy(&x));
    160:   PetscCall(VecDestroy(&ref));
    161:   PetscCall(TSDestroy(&ts));
    162:   PetscCall(PetscFinalize());
    163:   return 0;
    164: }

    166: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    167: {
    168:   PetscScalar       *f;
    169:   const PetscScalar *x;

    171:   PetscFunctionBegin;
    172:   PetscCall(VecGetArrayRead(X, &x));
    173:   PetscCall(VecGetArray(F, &f));
    174:   f[0] = x[1];
    175:   f[1] = 0.0;
    176:   PetscCall(VecRestoreArrayRead(X, &x));
    177:   PetscCall(VecRestoreArray(F, &f));
    178:   PetscFunctionReturn(PETSC_SUCCESS);
    179: }

    181: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
    182: {
    183:   User               user = (User)ptr;
    184:   PetscScalar       *f;
    185:   const PetscScalar *x, *xdot;

    187:   PetscFunctionBegin;
    188:   PetscCall(VecGetArrayRead(X, &x));
    189:   PetscCall(VecGetArrayRead(Xdot, &xdot));
    190:   PetscCall(VecGetArray(F, &f));
    191:   f[0] = xdot[0];
    192:   f[1] = xdot[1] - ((1. - x[0] * x[0]) * x[1] - x[0]) / user->mu;
    193:   PetscCall(VecRestoreArrayRead(X, &x));
    194:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
    195:   PetscCall(VecRestoreArray(F, &f));
    196:   PetscFunctionReturn(PETSC_SUCCESS);
    197: }

    199: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ptr)
    200: {
    201:   User               user     = (User)ptr;
    202:   PetscReal          mu       = user->mu;
    203:   PetscInt           rowcol[] = {0, 1};
    204:   PetscScalar        J[2][2];
    205:   const PetscScalar *x;

    207:   PetscFunctionBegin;
    208:   PetscCall(VecGetArrayRead(X, &x));
    209:   J[0][0] = a;
    210:   J[0][1] = 0;
    211:   J[1][0] = (2. * x[0] * x[1] + 1.) / mu;
    212:   J[1][1] = a - (1. - x[0] * x[0]) / mu;
    213:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
    214:   PetscCall(VecRestoreArrayRead(X, &x));

    216:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    217:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    218:   if (A != B) {
    219:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    220:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    221:   }
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

    225: /*TEST

    227:    test:
    228:      args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_row_col 3,3 -ts_monitor_lg_solution
    229:      requires: x

    231:    test:
    232:      suffix: adapt
    233:      args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_lg_solution
    234:      requires: x

    236:    test:
    237:      suffix: loop
    238:      args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt {{0.005 0.001 0.0005}separate output} -ts_max_steps {{100 500 1000}separate output} -ts_eimex_row_col {{1,1 2,1 3,1 2,2 3,2 3,3}separate output}
    239:      requires: x

    241:  TEST*/