Actual source code: ex16.c

      2: static char help[] = "Solves the van der Pol equation and demonstrate IMEX.\n\
      3: Input parameters include:\n\
      4:       -mu : stiffness parameter\n\n";

      6: /*
      7:    Concepts: TS^time-dependent nonlinear problems
      8:    Concepts: TS^van der Pol equation
      9:    Processors: 1
     10: */
     11: /* ------------------------------------------------------------------------

     13:    This program solves the van der Pol equation
     14:        y'' - \mu ((1-y^2)*y' - y) = 0        (1)
     15:    on the domain 0 <= x <= 1, with the boundary conditions
     16:        y(0) = 2, y'(0) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2),
     17:    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.

     19:    Notes:
     20:    This code demonstrates the TS solver interface to two variants of
     21:    linear problems, u_t = f(u,t), namely turning (1) into a system of
     22:    first order differential equations,

     24:    [ y' ] = [          z            ]
     25:    [ z' ]   [ \mu ((1 - y^2) z - y) ]

     27:    which then we can write as a vector equation

     29:    [ u_1' ] = [             u_2           ]  (2)
     30:    [ u_2' ]   [ \mu (1 - u_1^2) u_2 - u_1 ]

     32:    which is now in the desired form of u_t = f(u,t). One way that we
     33:    can split f(u,t) in (2) is to split by component,

     35:    [ u_1' ] = [ u_2 ] + [            0                ]
     36:    [ u_2' ]   [  0  ]   [ \mu ((1 - u_1^2) u_2 - u_1) ]

     38:    where

     40:    [ G(u,t) ] = [ u_2 ]
     41:                 [  0  ]

     43:    and

     45:    [ F(u',u,t) ] = [ u_1' ] - [            0                ]
     46:                    [ u_2' ]   [ \mu ((1 - u_1^2) u_2 - u_1) ]

     48:    Using the definition of the Jacobian of F (from the PETSc user manual),
     49:    in the equation F(u',u,t) = G(u,t),

     51:               dF   dF
     52:    J(F) = a * -- - --
     53:               du'  du

     55:    where d is the partial derivative. In this example,

     57:    dF   [ 1 ; 0 ]
     58:    -- = [       ]
     59:    du'  [ 0 ; 1 ]

     61:    dF   [       0             ;         0        ]
     62:    -- = [                                        ]
     63:    du   [ -\mu (2*u_1*u_2 + 1);  \mu (1 - u_1^2) ]

     65:    Hence,

     67:           [      a             ;          0          ]
     68:    J(F) = [                                          ]
     69:           [ \mu (2*u_1*u_2 + 1); a - \mu (1 - u_1^2) ]

     71:   ------------------------------------------------------------------------- */

     73: #include <petscts.h>

     75: typedef struct _n_User *User;
     76: struct _n_User {
     77:   PetscReal mu;
     78:   PetscBool imex;
     79:   PetscReal next_output;
     80: };

     82: /*
     83: *  User-defined routines
     84: */
     85: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx)
     86: {
     87:   PetscErrorCode    ierr;
     88:   User              user = (User)ctx;
     89:   PetscScalar       *f;
     90:   const PetscScalar *x;

     93:   VecGetArrayRead(X,&x);
     94:   VecGetArray(F,&f);
     95:   f[0] = (user->imex ? x[1] : 0);
     96:   f[1] = 0.0;
     97:   VecRestoreArrayRead(X,&x);
     98:   VecRestoreArray(F,&f);
     99:   return(0);
    100: }

    102: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
    103: {
    104:   PetscErrorCode    ierr;
    105:   User              user = (User)ctx;
    106:   const PetscScalar *x,*xdot;
    107:   PetscScalar       *f;

    110:   VecGetArrayRead(X,&x);
    111:   VecGetArrayRead(Xdot,&xdot);
    112:   VecGetArray(F,&f);
    113:   f[0] = xdot[0] + (user->imex ? 0 : x[1]);
    114:   f[1] = xdot[1] - user->mu*((1. - x[0]*x[0])*x[1] - x[0]);
    115:   VecRestoreArrayRead(X,&x);
    116:   VecRestoreArrayRead(Xdot,&xdot);
    117:   VecRestoreArray(F,&f);
    118:   return(0);
    119: }

    121: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
    122: {
    123:   PetscErrorCode    ierr;
    124:   User              user     = (User)ctx;
    125:   PetscReal         mu       = user->mu;
    126:   PetscInt          rowcol[] = {0,1};
    127:   const PetscScalar *x;
    128:   PetscScalar       J[2][2];

    131:   VecGetArrayRead(X,&x);
    132:   J[0][0] = a;                    J[0][1] = (user->imex ? 0 : 1.);
    133:   J[1][0] = mu*(2.*x[0]*x[1]+1.);   J[1][1] = a - mu*(1. - x[0]*x[0]);
    134:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
    135:   VecRestoreArrayRead(X,&x);

    137:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    138:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    139:   if (A != B) {
    140:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    141:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    142:   }
    143:   return(0);
    144: }

    146: static PetscErrorCode RegisterMyARK2(void)
    147: {

    151:   {
    152:     const PetscReal
    153:       A[3][3] = {{0,0,0},
    154:                  {0.41421356237309504880,0,0},
    155:                  {0.75,0.25,0}},
    156:       At[3][3] = {{0,0,0},
    157:                   {0.12132034355964257320,0.29289321881345247560,0},
    158:                   {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}},
    159:       *bembedt = NULL,*bembed = NULL;
    160:     TSARKIMEXRegister("myark2",2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembed,0,NULL,NULL);
    161:   }
    162:   return(0);
    163: }

    165: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
    166: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
    167: {
    168:   PetscErrorCode    ierr;
    169:   const PetscScalar *x;
    170:   PetscReal         tfinal, dt;
    171:   User              user = (User)ctx;
    172:   Vec               interpolatedX;

    175:   TSGetTimeStep(ts,&dt);
    176:   TSGetMaxTime(ts,&tfinal);

    178:   while (user->next_output <= t && user->next_output <= tfinal) {
    179:     VecDuplicate(X,&interpolatedX);
    180:     TSInterpolate(ts,user->next_output,interpolatedX);
    181:     VecGetArrayRead(interpolatedX,&x);
    182:     PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",user->next_output,step,t,dt,(double)PetscRealPart(x[0]),(double)PetscRealPart(x[1]));
    183:     VecRestoreArrayRead(interpolatedX,&x);
    184:     VecDestroy(&interpolatedX);

    186:     user->next_output += 0.1;
    187:   }
    188:   return(0);
    189: }

    191: int main(int argc,char **argv)
    192: {
    193:   TS             ts;            /* nonlinear solver */
    194:   Vec            x;             /* solution, residual vectors */
    195:   Mat            A;             /* Jacobian matrix */
    196:   PetscInt       steps;
    197:   PetscReal      ftime = 0.5;
    198:   PetscBool      monitor = PETSC_FALSE;
    199:   PetscScalar    *x_ptr;
    200:   PetscMPIInt    size;
    201:   struct _n_User user;

    204:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    205:      Initialize program
    206:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    207:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
    208:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    209:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");

    211:   RegisterMyARK2();

    213:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    214:     Set runtime options
    215:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    216:   user.mu          = 1000.0;
    217:   user.imex        = PETSC_TRUE;
    218:   user.next_output = 0.0;

    220:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
    221:   PetscOptionsGetBool(NULL,NULL,"-imex",&user.imex,NULL);
    222:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);

    224:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    225:     Create necessary matrix and vectors, solve same ODE on every process
    226:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    227:   MatCreate(PETSC_COMM_WORLD,&A);
    228:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
    229:   MatSetFromOptions(A);
    230:   MatSetUp(A);
    231:   MatCreateVecs(A,&x,NULL);

    233:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    234:      Create timestepping solver context
    235:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    236:   TSCreate(PETSC_COMM_WORLD,&ts);
    237:   TSSetType(ts,TSBEULER);
    238:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
    239:   TSSetIFunction(ts,NULL,IFunction,&user);
    240:   TSSetIJacobian(ts,A,A,IJacobian,&user);
    241:   TSSetMaxTime(ts,ftime);
    242:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
    243:   if (monitor) {
    244:     TSMonitorSet(ts,Monitor,&user,NULL);
    245:   }

    247:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    248:      Set initial conditions
    249:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    250:   VecGetArray(x,&x_ptr);
    251:   x_ptr[0] = 2.0;
    252:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
    253:   VecRestoreArray(x,&x_ptr);
    254:   TSSetTimeStep(ts,0.01);

    256:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    257:      Set runtime options
    258:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    259:   TSSetFromOptions(ts);

    261:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    262:      Solve nonlinear system
    263:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    264:   TSSolve(ts,x);
    265:   TSGetSolveTime(ts,&ftime);
    266:   TSGetStepNumber(ts,&steps);
    267:   PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,steps,(double)ftime);
    268:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

    270:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    271:      Free work space.  All PETSc objects should be destroyed when they
    272:      are no longer needed.
    273:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    274:   MatDestroy(&A);
    275:   VecDestroy(&x);
    276:   TSDestroy(&ts);

    278:   PetscFinalize();
    279:   return ierr;
    280: }

    282: /*TEST

    284:     test:
    285:       args: -ts_type arkimex -ts_arkimex_type myark2 -ts_adapt_type none
    286:       requires: !single

    288: TEST*/