Actual source code: ex5.c
petsc-3.12.2 2019-11-22
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*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 ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.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*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*/