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