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*/