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 "reaction_diffusion.h"
 39: #include <petscdm.h>
 40: #include <petscdmda.h>

 42: /* ------------------------------------------------------------------- */
 43: PetscErrorCode InitialConditions(DM da,Vec U)
 44: {
 45:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
 46:   Field          **u;
 47:   PetscReal      hx,hy,x,y;

 49:   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);

 51:   hx = 2.5/(PetscReal)(Mx);
 52:   hy = 2.5/(PetscReal)(My);

 54:   /*
 55:      Get pointers to actual vector data
 56:   */
 57:   DMDAVecGetArray(da,U,&u);

 59:   /*
 60:      Get local grid boundaries
 61:   */
 62:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

 64:   /*
 65:      Compute function over the locally owned part of the grid
 66:   */
 67:   for (j=ys; j<ys+ym; j++) {
 68:     y = j*hy;
 69:     for (i=xs; i<xs+xm; i++) {
 70:       x = i*hx;
 71:       if (PetscApproximateGTE(x,1.0) && PetscApproximateLTE(x,1.5) && PetscApproximateGTE(y,1.0) && PetscApproximateLTE(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;
 72:       else u[j][i].v = 0.0;

 74:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
 75:     }
 76:   }

 78:   /*
 79:      Restore access to vector
 80:   */
 81:   DMDAVecRestoreArray(da,U,&u);
 82:   return 0;
 83: }

 85: int main(int argc,char **argv)
 86: {
 87:   TS     ts;                    /* ODE integrator */
 88:   Vec    x;                     /* solution */
 89:   DM     da;
 90:   AppCtx appctx;

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Initialize program
 94:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   PetscInitialize(&argc,&argv,(char*)0,help);
 97:   appctx.D1    = 8.0e-5;
 98:   appctx.D2    = 4.0e-5;
 99:   appctx.gamma = .024;
100:   appctx.kappa = .06;
101:   appctx.aijpc = PETSC_FALSE;

103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104:      Create distributed array (DMDA) to manage parallel grid and vectors
105:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,65,65,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
107:   DMSetFromOptions(da);
108:   DMSetUp(da);
109:   DMDASetFieldName(da,0,"u");
110:   DMDASetFieldName(da,1,"v");

112:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113:      Create global vector from DMDA; this will be used to store the solution
114:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115:   DMCreateGlobalVector(da,&x);

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:      Create timestepping solver context
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   TSCreate(PETSC_COMM_WORLD,&ts);
121:   TSSetType(ts,TSARKIMEX);
122:   TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
123:   TSSetDM(ts,da);
124:   TSSetProblemType(ts,TS_NONLINEAR);
125:   TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);
126:   TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:      Set initial conditions
130:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131:   InitialConditions(da,x);
132:   TSSetSolution(ts,x);

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Set solver options
136:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   TSSetMaxTime(ts,2000.0);
138:   TSSetTimeStep(ts,.0001);
139:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
140:   TSSetFromOptions(ts);

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:      Solve ODE system
144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145:   TSSolve(ts,x);

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:      Free work space.
149:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150:   VecDestroy(&x);
151:   TSDestroy(&ts);
152:   DMDestroy(&da);

154:   PetscFinalize();
155:   return 0;
156: }

158: /*TEST

160:    build:
161:      depends: reaction_diffusion.c

163:    test:
164:       args: -ts_view  -ts_monitor -ts_max_time 500
165:       requires: double
166:       timeoutfactor: 3

168:    test:
169:       suffix: 2
170:       args: -ts_view  -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution
171:       requires: x double
172:       output_file: output/ex5_1.out
173:       timeoutfactor: 3

175: TEST*/