Actual source code: ex5.c

      2: static char help[] = "Demonstrates Pattern Formation with Reaction-Diffusion Equations.\n";

This example is taken from the book, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
W. Hundsdorf and J.G. Verwer, Page 21, Pattern Formation with Reaction-Diffusion Equations
\begin{eqnarray*}
u_t = D_1 (u_{xx} + u_{yy}) - u*v^2 + \gamma(1 -u) \\
v_t = D_2 (v_{xx} + v_{yy}) + u*v^2 - (\gamma + \kappa)v
\end{eqnarray*}
Unlike in the book this uses periodic boundary conditions instead of Neumann
(since they are easier for finite differences).

     15: /*
     16:       Helpful runtime monitor options:
     17:            -ts_monitor_draw_solution
     18:            -draw_save -draw_save_movie

     20:       Helpful runtime linear solver options:
     21:            -pc_type mg -pc_mg_galerkin pmat -da_refine 1 -snes_monitor -ksp_monitor -ts_view  (note that these Jacobians are so well-conditioned multigrid may not be the best solver)

     23:       Point your browser to localhost:8080 to monitor the simulation
     24:            ./ex5  -ts_view_pre saws  -stack_view saws -draw_save -draw_save_single_file -x_virtual -ts_monitor_draw_solution -saws_root .

     26: */

     28: /*

     30:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
     31:    Include "petscts.h" so that we can use SNES numerical (ODE) integrators.  Note that this
     32:    file automatically includes:
     33:      petscsys.h       - base PETSc routines   petscvec.h  - vectors
     34:      petscmat.h - matrices                    petscis.h   - index sets
     35:      petscksp.h - Krylov subspace methods     petscpc.h   - preconditioners
     36:      petscviewer.h - viewers                  petscsnes.h - nonlinear solvers
     37: */
     38: #include <petscdm.h>
     39: #include <petscdmda.h>
     40: #include <petscts.h>

     42: /*  Simple C struct that allows us to access the two velocity (x and y directions) values easily in the code */
     43: typedef struct {
     44:   PetscScalar u,v;
     45: } Field;

     47: /*  Data structure to store the model parameters */
     48: typedef struct {
     49:   PetscReal D1,D2,gamma,kappa;
     50: } AppCtx;

     52: /*
     53:    User-defined routines
     54: */
     55: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*),InitialConditions(DM,Vec);
     56: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);

     58: int main(int argc,char **argv)
     59: {
     60:   TS             ts;                  /* ODE integrator */
     61:   Vec            x;                   /* solution */
     63:   DM             da;
     64:   AppCtx         appctx;

     66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     67:      Initialize program
     68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     69:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     71:   appctx.D1    = 8.0e-5;
     72:   appctx.D2    = 4.0e-5;
     73:   appctx.gamma = .024;
     74:   appctx.kappa = .06;

     76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     77:      Create distributed array (DMDA) to manage parallel grid and vectors
     78:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     79:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,65,65,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
     80:   DMSetFromOptions(da);
     81:   DMSetUp(da);
     82:   DMDASetFieldName(da,0,"u");
     83:   DMDASetFieldName(da,1,"v");

     85:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     86:      Create global vector from DMDA; this will be used to store the solution
     87:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     88:   DMCreateGlobalVector(da,&x);

     90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     91:      Create timestepping solver context
     92:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     93:   TSCreate(PETSC_COMM_WORLD,&ts);
     94:   TSSetType(ts,TSARKIMEX);
     95:   TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
     96:   TSSetDM(ts,da);
     97:   TSSetProblemType(ts,TS_NONLINEAR);
     98:   TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);
     99:   TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);

    101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    102:      Set initial conditions
    103:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    104:   InitialConditions(da,x);
    105:   TSSetSolution(ts,x);

    107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    108:      Set solver options
    109:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    110:   TSSetMaxTime(ts,2000.0);
    111:   TSSetTimeStep(ts,.0001);
    112:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
    113:   TSSetFromOptions(ts);

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:      Solve ODE system
    117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    118:   TSSolve(ts,x);

    120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121:      Free work space.
    122:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    123:   VecDestroy(&x);
    124:   TSDestroy(&ts);
    125:   DMDestroy(&da);

    127:   PetscFinalize();
    128:   return ierr;
    129: }
    130: /* ------------------------------------------------------------------- */
    131: /*
    132:    RHSFunction - Evaluates nonlinear function, that defines the right
    133:      hand side of the ODE

    135:    Input Parameters:
    136: .  ts - the TS context
    137: .  time - current time
    138: .  U - input vector
    139: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

    141:    Output Parameter:
    142: .  F - function vector
    143:  */
    144: PetscErrorCode RHSFunction(TS ts,PetscReal time,Vec U,Vec F,void *ptr)
    145: {
    146:   AppCtx         *appctx = (AppCtx*)ptr;
    147:   DM             da;
    149:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
    150:   PetscReal      hx,hy,sx,sy;
    151:   PetscScalar    uc,uxx,uyy,vc,vxx,vyy;
    152:   Field          **u,**f;
    153:   Vec            localU;

    156:   TSGetDM(ts,&da);
    157:   /* Get local (ghosted) work vector */
    158:   DMGetLocalVector(da,&localU);
    159:   /* Get information about mesh needed for discretization */
    160:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    162:   hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx);
    163:   hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy);

    165:   /*
    166:      Scatter ghost points to local vector, using the 2-step process
    167:   */
    168:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
    169:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

    171:   /*
    172:      Get pointers to actual vector data
    173:   */
    174:   DMDAVecGetArrayRead(da,localU,&u);
    175:   DMDAVecGetArray(da,F,&f);

    177:   /*
    178:      Get local grid boundaries; this is the region that this process owns and must operate on
    179:   */
    180:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

    182:   /*
    183:      Compute function over the locally owned part of the grid with standard finite differences
    184:   */
    185:   for (j=ys; j<ys+ym; j++) {
    186:     for (i=xs; i<xs+xm; i++) {
    187:       uc        = u[j][i].u;
    188:       uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
    189:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
    190:       vc        = u[j][i].v;
    191:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
    192:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
    193:       f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);
    194:       f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc;
    195:     }
    196:   }
    197:   PetscLogFlops(16.0*xm*ym);

    199:   /*
    200:      Restore access to vectors and return no longer needed work vector
    201:   */
    202:   DMDAVecRestoreArrayRead(da,localU,&u);
    203:   DMDAVecRestoreArray(da,F,&f);
    204:   DMRestoreLocalVector(da,&localU);
    205:   return(0);
    206: }

    208: /* ------------------------------------------------------------------- */
    209: PetscErrorCode InitialConditions(DM da,Vec U)
    210: {
    212:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
    213:   Field          **u;
    214:   PetscReal      hx,hy,x,y;

    217:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    219:   hx = 2.5/(PetscReal)(Mx);
    220:   hy = 2.5/(PetscReal)(My);

    222:   /*
    223:      Get pointers to actual vector data
    224:   */
    225:   DMDAVecGetArray(da,U,&u);

    227:   /*
    228:      Get local grid boundaries
    229:   */
    230:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

    232:   /*
    233:      Compute function over the locally owned part of the grid
    234:   */
    235:   for (j=ys; j<ys+ym; j++) {
    236:     y = j*hy;
    237:     for (i=xs; i<xs+xm; i++) {
    238:       x = i*hx;
    239:       if (PetscGTE(x,1.0) && PetscLTE(x,1.5) && PetscGTE(y,1.0) && PetscLTE(y,1.5)) u[j][i].v = PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0)/4.0;
    240:       else u[j][i].v = 0.0;

    242:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
    243:     }
    244:   }

    246:   /*
    247:      Restore access to vector
    248:   */
    249:   DMDAVecRestoreArray(da,U,&u);
    250:   return(0);
    251: }

    253: /*
    254:    RHSJacobian - Evaluates the Jacobian of the right hand side
    255:      function of the ODE.

    257:    Input Parameters:
    258: .  ts - the TS context
    259: .  time - current time
    260: .  U - input vector
    261: .  ptr - optional user-defined context, as set by TSSetRHSJacobian()

    263:    Output Parameter:
    264: .  A - the Jacobian
    265: .  BB - optional additional matrix where an approximation to the Jacobian
    266:         may be stored from which the preconditioner is constructed
    267:  */
    268: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
    269: {
    270:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
    271:   DM             da;
    273:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
    274:   PetscReal      hx,hy,sx,sy;
    275:   PetscScalar    uc,vc;
    276:   Field          **u;
    277:   Vec            localU;
    278:   MatStencil     stencil[6],rowstencil;
    279:   PetscScalar    entries[6];

    282:   TSGetDM(ts,&da);
    283:   DMGetLocalVector(da,&localU);
    284:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    286:   hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx);
    287:   hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy);

    289:   /*
    290:      Scatter ghost points to local vector,using the 2-step process
    291:   */
    292:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
    293:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

    295:   /*
    296:      Get pointers to vector data
    297:   */
    298:   DMDAVecGetArrayRead(da,localU,&u);

    300:   /*
    301:      Get local grid boundaries
    302:   */
    303:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

    305:   stencil[0].k = 0;
    306:   stencil[1].k = 0;
    307:   stencil[2].k = 0;
    308:   stencil[3].k = 0;
    309:   stencil[4].k = 0;
    310:   stencil[5].k = 0;
    311:   rowstencil.k = 0;
    312:   /*
    313:      Compute function over the locally owned part of the grid
    314:   */
    315:   for (j=ys; j<ys+ym; j++) {

    317:     stencil[0].j = j-1;
    318:     stencil[1].j = j+1;
    319:     stencil[2].j = j;
    320:     stencil[3].j = j;
    321:     stencil[4].j = j;
    322:     stencil[5].j = j;
    323:     rowstencil.k = 0; rowstencil.j = j;
    324:     for (i=xs; i<xs+xm; i++) {
    325:       uc = u[j][i].u;
    326:       vc = u[j][i].v;

    328:       /*      uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
    329:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;

    331:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
    332:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
    333:        f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/

    335:       stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy;
    336:       stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy;
    337:       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx;
    338:       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx;
    339:       stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma;
    340:       stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc;
    341:       rowstencil.i = i; rowstencil.c = 0;

    343:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);

    345:       stencil[0].c = 1; entries[0] = appctx->D2*sy;
    346:       stencil[1].c = 1; entries[1] = appctx->D2*sy;
    347:       stencil[2].c = 1; entries[2] = appctx->D2*sx;
    348:       stencil[3].c = 1; entries[3] = appctx->D2*sx;
    349:       stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa;
    350:       stencil[5].c = 0; entries[5] = vc*vc;
    351:       rowstencil.c = 1;

    353:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
    354:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
    355:     }
    356:   }

    358:   /*
    359:      Restore vectors
    360:   */
    361:   PetscLogFlops(19.0*xm*ym);
    362:   DMDAVecRestoreArrayRead(da,localU,&u);
    363:   DMRestoreLocalVector(da,&localU);
    364:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    365:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    366:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
    367:   return(0);
    368: }

    371: /*TEST

    373:    test:
    374:       args: -ts_view  -ts_monitor -ts_max_time 500
    375:       requires: double
    376:       timeoutfactor: 3

    378:    test:
    379:       suffix: 2
    380:       args: -ts_view  -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution
    381:       requires: x double
    382:       output_file: output/ex5_1.out
    383:       timeoutfactor: 3

    385: TEST*/