Actual source code: allen_cahn.c

      1: static char help[] ="Solves the time dependent Allen-Cahn equation with IMEX methods";

      3: /*
      4:  * allen_cahn.c
      5:  *
      6:  *  Created on: Jun 8, 2012
      7:  *      Author: Hong Zhang
      8:  */

     10: #include <petscts.h>

     12: /*
     13:  * application context
     14:  */
     15: typedef struct {
     16:   PetscReal   param;        /* parameter */
     17:   PetscReal   xleft,xright;  /* range in x-direction */
     18:   PetscInt    mx;           /* Discretization in x-direction */
     19: }AppCtx;

     22: static PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
     23: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
     24: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void *ctx);
     25: static PetscErrorCode FormInitialSolution(TS,Vec,void*);

     28: int main(int argc, char **argv)
     29: {
     30:   TS                ts;
     31:   Vec               x; /*solution vector*/
     32:   Mat               A; /*Jacobian*/
     33:   PetscInt          steps,mx;
     34:   PetscErrorCode    ierr;
     35:   PetscReal         ftime;
     36:   AppCtx            user;       /* user-defined work context */

     38:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;

     40:   /* Initialize user application context */
     41:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Allen-Cahn equation","");
     42:   user.param = 9e-4;
     43:   user.xleft = -1.;
     44:   user.xright = 2.;
     45:   user.mx = 400;
     46:   PetscOptionsReal("-eps","parameter","",user.param,&user.param,NULL);
     47:   PetscOptionsEnd();

     49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     50:    Set runtime options
     51:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     52:   /*
     53:    * PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
     54:    */
     55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     56:    Create necessary matrix and vectors, solve same ODE on every process
     57:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     58:   MatCreate(PETSC_COMM_WORLD,&A);
     59:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,user.mx,user.mx);
     60:   MatSetFromOptions(A);
     61:   MatSetUp(A);
     62:   MatCreateVecs(A,&x,NULL);

     64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     65:    Create time stepping solver context
     66:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     67:   TSCreate(PETSC_COMM_WORLD,&ts);
     68:   TSSetType(ts,TSEIMEX);
     69:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
     70:   TSSetIFunction(ts,NULL,FormIFunction,&user);
     71:   TSSetIJacobian(ts,A,A,FormIJacobian,&user);
     72:   ftime = 22;
     73:   TSSetMaxTime(ts,ftime);
     74:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

     76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     77:    Set initial conditions
     78:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     79:   FormInitialSolution(ts,x,&user);
     80:   TSSetSolution(ts,x);
     81:   VecGetSize(x,&mx);

     83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     84:    Set runtime options
     85:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     86:   TSSetFromOptions(ts);

     88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     89:    Solve nonlinear system
     90:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     91:   TSSolve(ts,x);
     92:   TSGetTime(ts,&ftime);
     93:   TSGetStepNumber(ts,&steps);
     94:   PetscPrintf(PETSC_COMM_WORLD,"eps %g, steps %D, ftime %g\n",(double)user.param,steps,(double)ftime);
     95:   /*   VecView(x,PETSC_VIEWER_STDOUT_WORLD);*/

     97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     98:    Free work space.
     99:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    100:   MatDestroy(&A);
    101:   VecDestroy(&x);
    102:   TSDestroy(&ts);
    103:   PetscFinalize();
    104:   return ierr;
    105: }

    108: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
    109: {
    110:   PetscErrorCode    ierr;
    111:   AppCtx            *user = (AppCtx*)ptr;
    112:   PetscScalar       *f;
    113:   const PetscScalar *x;
    114:   PetscInt          i,mx;
    115:   PetscReal         hx,eps;

    118:   mx = user->mx;
    119:   eps = user->param;
    120:   hx = (user->xright-user->xleft)/(mx-1);
    121:   VecGetArrayRead(X,&x);
    122:   VecGetArray(F,&f);
    123:   f[0] = 2.*eps*(x[1]-x[0])/(hx*hx); /*boundary*/
    124:   for (i=1;i<mx-1;i++) {
    125:     f[i]= eps*(x[i+1]-2.*x[i]+x[i-1])/(hx*hx);
    126:   }
    127:   f[mx-1] = 2.*eps*(x[mx-2]- x[mx-1])/(hx*hx); /*boundary*/
    128:   VecRestoreArrayRead(X,&x);
    129:   VecRestoreArray(F,&f);
    130:   return(0);
    131: }

    134: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
    135: {
    136:   PetscErrorCode    ierr;
    137:   AppCtx            *user = (AppCtx*)ptr;
    138:   PetscScalar       *f;
    139:   const PetscScalar *x,*xdot;
    140:   PetscInt          i,mx;

    143:   mx = user->mx;
    144:   VecGetArrayRead(X,&x);
    145:   VecGetArrayRead(Xdot,&xdot);
    146:   VecGetArray(F,&f);

    148:   for (i=0;i<mx;i++) {
    149:     f[i]= xdot[i] - x[i]*(1.-x[i]*x[i]);
    150:   }

    152:   VecRestoreArrayRead(X,&x);
    153:   VecRestoreArrayRead(Xdot,&xdot);
    154:   VecRestoreArray(F,&f);
    155:   return(0);
    156: }

    159: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U, Vec Udot, PetscReal a, Mat J,Mat Jpre,void *ctx)
    160: {
    161:   PetscErrorCode    ierr;
    162:   AppCtx            *user = (AppCtx *)ctx;
    163:   PetscScalar       v;
    164:   const PetscScalar *x;
    165:   PetscInt          i,col;

    168:   VecGetArrayRead(U,&x);
    169:   for (i=0; i < user->mx; i++) {
    170:     v = a - 1. + 3.*x[i]*x[i];
    171:     col = i;
    172:     MatSetValues(J,1,&i,1,&col,&v,INSERT_VALUES);
    173:   }
    174:   VecRestoreArrayRead(U,&x);

    176:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    177:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    178:   if (J != Jpre) {
    179:     MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
    180:     MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
    181:   }
    182:   /*  MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/
    183:   return(0);
    184: }

    187: static PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ctx)
    188: {
    189:   AppCtx *user = (AppCtx*)ctx;
    190:   PetscInt       i;
    191:   PetscScalar    *x;
    192:   PetscReal      hx,x_map;

    196:   hx = (user->xright-user->xleft)/(PetscReal)(user->mx-1);
    197:   VecGetArray(U,&x);
    198:   for (i=0;i<user->mx;i++) {
    199:     x_map = user->xleft + i*hx;
    200:     if (x_map >= 0.7065) {
    201:       x[i] = PetscTanhReal((x_map-0.8)/(2.*PetscSqrtReal(user->param)));
    202:     } else if (x_map >= 0.4865) {
    203:       x[i] = PetscTanhReal((0.613-x_map)/(2.*PetscSqrtReal(user->param)));
    204:     } else if (x_map >= 0.28) {
    205:       x[i] = PetscTanhReal((x_map-0.36)/(2.*PetscSqrtReal(user->param)));
    206:     } else if (x_map >= -0.7) {
    207:       x[i] = PetscTanhReal((0.2-x_map)/(2.*PetscSqrtReal(user->param)));
    208:     } else if (x_map >= -1) {
    209:       x[i] = PetscTanhReal((x_map+0.9)/(2.*PetscSqrtReal(user->param)));
    210:     }
    211:   }
    212:   VecRestoreArray(U,&x);
    213:   return(0);
    214: }

    216: /*TEST

    218:      test:
    219:        args:  -ts_rtol 1e-04 -ts_dt 0.025 -pc_type lu -ksp_error_if_not_converged TRUE  -ts_type eimex -ts_adapt_type none -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_draw_solution
    220:        requires: x
    221:        timeoutfactor: 3

    223: TEST*/