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:  * Processors:1
     11:  */

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

     33: #include <petscdmda.h>
     34: #include <petscts.h>

     36: typedef struct _User *User;
     37: struct _User {
     38:   PetscReal mu;  /*stiffness control coefficient: epsilon*/
     39: };

     41: static PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
     42: static PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
     43: static PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

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

     59:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
     60:   /* Initialize user application context */
     61:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"van der Pol options","");
     62:   user.mu      = 1e0;
     63:   PetscOptionsReal("-eps","Stiffness controller","",user.mu,&user.mu,NULL);
     64:   PetscOptionsEnd();

     66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     67:    Set runtime options
     68:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     69:   /*
     70:    PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
     71:    */

     73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     74:    Create necessary matrix and vectors, solve same ODE on every process
     75:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     76:   MatCreate(PETSC_COMM_WORLD,&A);
     77:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
     78:   MatSetFromOptions(A);
     79:   MatSetUp(A);
     80:   MatCreateVecs(A,&x,NULL);

     82:   MatCreateVecs(A,&ref,NULL);
     83:   VecGetArray(ref,&x_ptr);
     84:   /*
     85:    * [0,1], mu=10^-3
     86:    */
     87:   /*
     88:    x_ptr[0] = -1.8881254106283;
     89:    x_ptr[1] =  0.7359074233370;*/

     91:   /*
     92:    * [0,0.5],mu=10^-3
     93:    */
     94:   /*
     95:    x_ptr[0] = 1.596980778659137;
     96:    x_ptr[1] = -1.029103015879544;
     97:    */
     98:   /*
     99:    * [0,0.5],mu=1
    100:    */
    101:   x_ptr[0] = 1.619084329683235;
    102:   x_ptr[1] = -0.803530465176385;

    104:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    105:    Create timestepping solver context
    106:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    107:   TSCreate(PETSC_COMM_WORLD,&ts);
    108:   TSSetType(ts,TSEIMEX);
    109:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
    110:   TSSetIFunction(ts,NULL,IFunction,&user);
    111:   TSSetIJacobian(ts,A,A,IJacobian,&user);

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

    128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    129:    Set runtime options
    130:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    131:   TSSetFromOptions(ts);

    133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    134:    Solve nonlinear system
    135:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    136:   TSSolve(ts,x);
    137:   TSGetTime(ts,&ftime);
    138:   TSGetStepNumber(ts,&steps);

    140:   VecAXPY(x,-1.0,ref);
    141:   VecNorm(x,NORM_2,&norm);
    142:   TSGetTimeStep(ts,&dt);

    144:   eimex_rowcol[0] = 0; eimex_rowcol[1] = 0; two = 2;
    145:   PetscOptionsGetIntArray(NULL,NULL,"-ts_eimex_row_col",eimex_rowcol,&two,NULL);
    146:   PetscPrintf(PETSC_COMM_WORLD,"order %11s %18s %37s\n","dt","norm","final solution components 0 and 1");
    147:   VecGetArray(x,&x_ptr);
    148:   PetscPrintf(PETSC_COMM_WORLD,"(%D,%D) %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]));
    149:   VecRestoreArray(x,&x_ptr);

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

    159:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    160:    Free work space.
    161:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    162:   MatDestroy(&A);
    163:   VecDestroy(&x);
    164:   VecDestroy(&ref);
    165:   TSDestroy(&ts);
    166:   PetscFinalize();
    167:   return ierr;
    168: }

    171: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
    172: {
    173:   PetscErrorCode    ierr;
    174:   PetscScalar       *f;
    175:   const PetscScalar *x;

    178:   VecGetArrayRead(X,&x);
    179:   VecGetArray(F,&f);
    180:   f[0] = x[1];
    181:   f[1] = 0.0;
    182:   VecRestoreArrayRead(X,&x);
    183:   VecRestoreArray(F,&f);
    184:   return(0);
    185: }

    187: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
    188: {
    189:   User              user = (User)ptr;
    190:   PetscScalar       *f;
    191:   const PetscScalar *x,*xdot;
    192:   PetscErrorCode    ierr;

    195:   VecGetArrayRead(X,&x);
    196:   VecGetArrayRead(Xdot,&xdot);
    197:   VecGetArray(F,&f);
    198:   f[0] = xdot[0];
    199:   f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu;
    200:   VecRestoreArrayRead(X,&x);
    201:   VecRestoreArrayRead(Xdot,&xdot);
    202:   VecRestoreArray(F,&f);
    203:   return(0);
    204: }

    206: static PetscErrorCode IJacobian(TS  ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ptr)
    207: {
    208:   PetscErrorCode    ierr;
    209:   User              user = (User)ptr;
    210:   PetscReal         mu = user->mu;
    211:   PetscInt          rowcol[] = {0,1};
    212:   PetscScalar       J[2][2];
    213:   const PetscScalar *x;

    216:   VecGetArrayRead(X,&x);
    217:   J[0][0] = a;
    218:   J[0][1] = 0;
    219:   J[1][0] = (2.*x[0]*x[1]+1.)/mu;
    220:   J[1][1] = a - (1. - x[0]*x[0])/mu;
    221:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
    222:   VecRestoreArrayRead(X,&x);

    224:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    225:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    226:   if (A != B) {
    227:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    228:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    229:   }
    230:   return(0);
    231: }

    233: /*TEST

    235:    test:
    236:      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
    237:      requires: x

    239:    test:
    240:      suffix: adapt
    241:      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
    242:      requires: x

    244:    test:
    245:      suffix: loop
    246:      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}
    247:      requires: x

    249:  TEST*/