1
2 static char help[] = "Demonstrates Pattern Formation with Reaction-Diffusion Equations.\n";
3
4 /*F
5 This example is taken from the book, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
6 W. Hundsdorf and J.G. Verwer, Page 21, Pattern Formation with Reaction-Diffusion Equations
7 \begin{eqnarray*}
8 u_t = D_1 (u_{xx} + u_{yy}) - u*v^2 + \gamma(1 -u) \\
9 v_t = D_2 (v_{xx} + v_{yy}) + u*v^2 - (\gamma + \kappa)v
10 \end{eqnarray*}
11 Unlike in the book this uses periodic boundary conditions instead of Neumann
12 (since they are easier for finite differences).
13 F*/
14
15 /*
16 Helpful runtime monitor options:
17 -ts_monitor_draw_solution
18 -draw_save -draw_save_movie
19
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)
22
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 .
25
26 */
27
28 /*
29
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>
41
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;
46
47 /* Data structure to store the model parameters */
48 typedef struct {
49 PetscReal D1,D2,gamma,kappa;
50 } AppCtx;
51
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*);
57
main(int argc,char ** argv)58 int main(int argc,char **argv)
59 {
60 TS ts; /* ODE integrator */
61 Vec x; /* solution */
62 PetscErrorCode ierr;
63 DM da;
64 AppCtx appctx;
65
66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67 Initialize program
68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
70 PetscFunctionBeginUser;
71 appctx.D1 = 8.0e-5;
72 appctx.D2 = 4.0e-5;
73 appctx.gamma = .024;
74 appctx.kappa = .06;
75
76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77 Create distributed array (DMDA) to manage parallel grid and vectors
78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,65,65,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);
80 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
81 ierr = DMSetUp(da);CHKERRQ(ierr);
82 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
83 ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr);
84
85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 Create global vector from DMDA; this will be used to store the solution
87 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
89
90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 Create timestepping solver context
92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
94 ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);
95 ierr = TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);CHKERRQ(ierr);
96 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
97 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
98 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);CHKERRQ(ierr);
99 ierr = TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr);
100
101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102 Set initial conditions
103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104 ierr = InitialConditions(da,x);CHKERRQ(ierr);
105 ierr = TSSetSolution(ts,x);CHKERRQ(ierr);
106
107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108 Set solver options
109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110 ierr = TSSetMaxTime(ts,2000.0);CHKERRQ(ierr);
111 ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr);
112 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
113 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
114
115 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116 Solve ODE system
117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118 ierr = TSSolve(ts,x);CHKERRQ(ierr);
119
120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121 Free work space.
122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123 ierr = VecDestroy(&x);CHKERRQ(ierr);
124 ierr = TSDestroy(&ts);CHKERRQ(ierr);
125 ierr = DMDestroy(&da);CHKERRQ(ierr);
126
127 ierr = PetscFinalize();
128 return ierr;
129 }
130 /* ------------------------------------------------------------------- */
131 /*
132 RHSFunction - Evaluates nonlinear function, that defines the right
133 hand side of the ODE
134
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()
140
141 Output Parameter:
142 . F - function vector
143 */
RHSFunction(TS ts,PetscReal time,Vec U,Vec F,void * ptr)144 PetscErrorCode RHSFunction(TS ts,PetscReal time,Vec U,Vec F,void *ptr)
145 {
146 AppCtx *appctx = (AppCtx*)ptr;
147 DM da;
148 PetscErrorCode ierr;
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;
154
155 PetscFunctionBegin;
156 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
157 /* Get local (ghosted) work vector */
158 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
159 /* Get information about mesh needed for discretization */
160 ierr = 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);CHKERRQ(ierr);
161
162 hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx);
163 hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy);
164
165 /*
166 Scatter ghost points to local vector, using the 2-step process
167 */
168 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
169 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
170
171 /*
172 Get pointers to actual vector data
173 */
174 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
175 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
176
177 /*
178 Get local grid boundaries; this is the region that this process owns and must operate on
179 */
180 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
181
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 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
198
199 /*
200 Restore access to vectors and return no longer needed work vector
201 */
202 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
203 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
204 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
205 PetscFunctionReturn(0);
206 }
207
208 /* ------------------------------------------------------------------- */
InitialConditions(DM da,Vec U)209 PetscErrorCode InitialConditions(DM da,Vec U)
210 {
211 PetscErrorCode ierr;
212 PetscInt i,j,xs,ys,xm,ym,Mx,My;
213 Field **u;
214 PetscReal hx,hy,x,y;
215
216 PetscFunctionBegin;
217 ierr = 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);CHKERRQ(ierr);
218
219 hx = 2.5/(PetscReal)(Mx);
220 hy = 2.5/(PetscReal)(My);
221
222 /*
223 Get pointers to actual vector data
224 */
225 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
226
227 /*
228 Get local grid boundaries
229 */
230 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
231
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;
241
242 u[j][i].u = 1.0 - 2.0*u[j][i].v;
243 }
244 }
245
246 /*
247 Restore access to vector
248 */
249 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
250 PetscFunctionReturn(0);
251 }
252
253 /*
254 RHSJacobian - Evaluates the Jacobian of the right hand side
255 function of the ODE.
256
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()
262
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 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void * ctx)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;
272 PetscErrorCode ierr;
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];
280
281 PetscFunctionBegin;
282 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
283 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
284 ierr = 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);CHKERRQ(ierr);
285
286 hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx);
287 hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy);
288
289 /*
290 Scatter ghost points to local vector,using the 2-step process
291 */
292 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
293 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
294
295 /*
296 Get pointers to vector data
297 */
298 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
299
300 /*
301 Get local grid boundaries
302 */
303 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
304
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++) {
316
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;
327
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;
330
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);*/
334
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;
342
343 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
344
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;
352
353 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
354 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
355 }
356 }
357
358 /*
359 Restore vectors
360 */
361 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr);
362 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
363 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
364 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
365 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
366 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
367 PetscFunctionReturn(0);
368 }
369
370
371 /*TEST
372
373 test:
374 args: -ts_view -ts_monitor -ts_max_time 500
375 requires: double
376 timeoutfactor: 3
377
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
384
385 TEST*/
386