Actual source code: reaction_diffusion.c

  1: #include "reaction_diffusion.h"
2: #include <petscdm.h>
3: #include <petscdmda.h>

5: /*F
6:      This example is taken from the book, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
7:       W. Hundsdorf and J.G. Verwer,  Page 21, Pattern Formation with Reaction-Diffusion Equations
8: \begin{eqnarray*}
9:         u_t = D_1 (u_{xx} + u_{yy})  - u*v^2 + \gamma(1 -u)           \\
10:         v_t = D_2 (v_{xx} + v_{yy})  + u*v^2 - (\gamma + \kappa)v
11: \end{eqnarray*}
12:     Unlike in the book this uses periodic boundary conditions instead of Neumann
13:     (since they are easier for finite differences).
14: F*/

16: /*
17:    RHSFunction - Evaluates nonlinear function, F(x).

19:    Input Parameters:
20: .  ts - the TS context
21: .  X - input vector
22: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

24:    Output Parameter:
25: .  F - function vector
26:  */
27: PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec U, Vec F, void *ptr)
28: {
29:   AppCtx     *appctx = (AppCtx *)ptr;
30:   DM          da;
31:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
32:   PetscReal   hx, hy, sx, sy;
33:   PetscScalar uc, uxx, uyy, vc, vxx, vyy;
34:   Field     **u, **f;
35:   Vec         localU;

37:   PetscFunctionBegin;
38:   PetscCall(TSGetDM(ts, &da));
39:   PetscCall(DMGetLocalVector(da, &localU));
40:   PetscCall(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));
41:   hx = 2.50 / (PetscReal)Mx;
42:   sx = 1.0 / (hx * hx);
43:   hy = 2.50 / (PetscReal)My;
44:   sy = 1.0 / (hy * hy);

46:   /*
47:      Scatter ghost points to local vector,using the 2-step process
48:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
49:      By placing code between these two statements, computations can be
50:      done while messages are in transition.
51:   */
52:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
53:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));

55:   /*
56:      Get pointers to vector data
57:   */
59:   PetscCall(DMDAVecGetArray(da, F, &f));

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

66:   /*
67:      Compute function over the locally owned part of the grid
68:   */
69:   for (j = ys; j < ys + ym; j++) {
70:     for (i = xs; i < xs + xm; i++) {
71:       uc        = u[j][i].u;
72:       uxx       = (-2.0 * uc + u[j][i - 1].u + u[j][i + 1].u) * sx;
73:       uyy       = (-2.0 * uc + u[j - 1][i].u + u[j + 1][i].u) * sy;
74:       vc        = u[j][i].v;
75:       vxx       = (-2.0 * vc + u[j][i - 1].v + u[j][i + 1].v) * sx;
76:       vyy       = (-2.0 * vc + u[j - 1][i].v + u[j + 1][i].v) * sy;
77:       f[j][i].u = appctx->D1 * (uxx + uyy) - uc * vc * vc + appctx->gamma * (1.0 - uc);
78:       f[j][i].v = appctx->D2 * (vxx + vyy) + uc * vc * vc - (appctx->gamma + appctx->kappa) * vc;
79:     }
80:   }
81:   PetscCall(PetscLogFlops(16.0 * xm * ym));

83:   /*
84:      Restore vectors
85:   */
87:   PetscCall(DMDAVecRestoreArray(da, F, &f));
88:   PetscCall(DMRestoreLocalVector(da, &localU));
89:   PetscFunctionReturn(PETSC_SUCCESS);
90: }

92: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat BB, void *ctx)
93: {
94:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
95:   DM          da;
96:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
97:   PetscReal   hx, hy, sx, sy;
98:   PetscScalar uc, vc;
99:   Field     **u;
100:   Vec         localU;
101:   MatStencil  stencil[6], rowstencil;
102:   PetscScalar entries[6];

104:   PetscFunctionBegin;
105:   PetscCall(TSGetDM(ts, &da));
106:   PetscCall(DMGetLocalVector(da, &localU));
107:   PetscCall(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));

109:   hx = 2.50 / (PetscReal)Mx;
110:   sx = 1.0 / (hx * hx);
111:   hy = 2.50 / (PetscReal)My;
112:   sy = 1.0 / (hy * hy);

114:   /*
115:      Scatter ghost points to local vector,using the 2-step process
116:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
117:      By placing code between these two statements, computations can be
118:      done while messages are in transition.
119:   */
120:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
121:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));

123:   /*
124:      Get pointers to vector data
125:   */

128:   /*
129:      Get local grid boundaries
130:   */
131:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

133:   stencil[0].k = 0;
134:   stencil[1].k = 0;
135:   stencil[2].k = 0;
136:   stencil[3].k = 0;
137:   stencil[4].k = 0;
138:   stencil[5].k = 0;
139:   rowstencil.k = 0;
140:   /*
141:      Compute function over the locally owned part of the grid
142:   */
143:   for (j = ys; j < ys + ym; j++) {
144:     stencil[0].j = j - 1;
145:     stencil[1].j = j + 1;
146:     stencil[2].j = j;
147:     stencil[3].j = j;
148:     stencil[4].j = j;
149:     stencil[5].j = j;
150:     rowstencil.k = 0;
151:     rowstencil.j = j;
152:     for (i = xs; i < xs + xm; i++) {
153:       uc = u[j][i].u;
154:       vc = u[j][i].v;

156:       /*      uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
157:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;

159:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
160:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
161:        f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/

163:       stencil[0].i = i;
164:       stencil[0].c = 0;
165:       entries[0]   = appctx->D1 * sy;
166:       stencil[1].i = i;
167:       stencil[1].c = 0;
168:       entries[1]   = appctx->D1 * sy;
169:       stencil[2].i = i - 1;
170:       stencil[2].c = 0;
171:       entries[2]   = appctx->D1 * sx;
172:       stencil[3].i = i + 1;
173:       stencil[3].c = 0;
174:       entries[3]   = appctx->D1 * sx;
175:       stencil[4].i = i;
176:       stencil[4].c = 0;
177:       entries[4]   = -2.0 * appctx->D1 * (sx + sy) - vc * vc - appctx->gamma;
178:       stencil[5].i = i;
179:       stencil[5].c = 1;
180:       entries[5]   = -2.0 * uc * vc;
181:       rowstencil.i = i;
182:       rowstencil.c = 0;

184:       PetscCall(MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
185:       if (appctx->aijpc) PetscCall(MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
186:       stencil[0].c = 1;
187:       entries[0]   = appctx->D2 * sy;
188:       stencil[1].c = 1;
189:       entries[1]   = appctx->D2 * sy;
190:       stencil[2].c = 1;
191:       entries[2]   = appctx->D2 * sx;
192:       stencil[3].c = 1;
193:       entries[3]   = appctx->D2 * sx;
194:       stencil[4].c = 1;
195:       entries[4]   = -2.0 * appctx->D2 * (sx + sy) + 2.0 * uc * vc - appctx->gamma - appctx->kappa;
196:       stencil[5].c = 0;
197:       entries[5]   = vc * vc;
198:       rowstencil.c = 1;

200:       PetscCall(MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
201:       if (appctx->aijpc) PetscCall(MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
202:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
203:     }
204:   }

206:   /*
207:      Restore vectors
208:   */
209:   PetscCall(PetscLogFlops(19.0 * xm * ym));
211:   PetscCall(DMRestoreLocalVector(da, &localU));
212:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
213:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
214:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
215:   if (appctx->aijpc) {
216:     PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
217:     PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
218:     PetscCall(MatSetOption(BB, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
219:   }
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: /*
224:    IFunction - Evaluates implicit nonlinear function, xdot - F(x).

226:    Input Parameters:
227: .  ts - the TS context
228: .  U - input vector
229: .  Udot - input vector
230: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

232:    Output Parameter:
233: .  F - function vector
234:  */
235: PetscErrorCode IFunction(TS ts, PetscReal ftime, Vec U, Vec Udot, Vec F, void *ptr)
236: {
237:   AppCtx     *appctx = (AppCtx *)ptr;
238:   DM          da;
239:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
240:   PetscReal   hx, hy, sx, sy;
241:   PetscScalar uc, uxx, uyy, vc, vxx, vyy;
242:   Field     **u, **f, **udot;
243:   Vec         localU;

245:   PetscFunctionBegin;
246:   PetscCall(TSGetDM(ts, &da));
247:   PetscCall(DMGetLocalVector(da, &localU));
248:   PetscCall(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));
249:   hx = 2.50 / (PetscReal)Mx;
250:   sx = 1.0 / (hx * hx);
251:   hy = 2.50 / (PetscReal)My;
252:   sy = 1.0 / (hy * hy);

254:   /*
255:      Scatter ghost points to local vector,using the 2-step process
256:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
257:      By placing code between these two statements, computations can be
258:      done while messages are in transition.
259:   */
260:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
261:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));

263:   /*
264:      Get pointers to vector data
265:   */
267:   PetscCall(DMDAVecGetArray(da, F, &f));

270:   /*
271:      Get local grid boundaries
272:   */
273:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

275:   /*
276:      Compute function over the locally owned part of the grid
277:   */
278:   for (j = ys; j < ys + ym; j++) {
279:     for (i = xs; i < xs + xm; i++) {
280:       uc        = u[j][i].u;
281:       uxx       = (-2.0 * uc + u[j][i - 1].u + u[j][i + 1].u) * sx;
282:       uyy       = (-2.0 * uc + u[j - 1][i].u + u[j + 1][i].u) * sy;
283:       vc        = u[j][i].v;
284:       vxx       = (-2.0 * vc + u[j][i - 1].v + u[j][i + 1].v) * sx;
285:       vyy       = (-2.0 * vc + u[j - 1][i].v + u[j + 1][i].v) * sy;
286:       f[j][i].u = udot[j][i].u - (appctx->D1 * (uxx + uyy) - uc * vc * vc + appctx->gamma * (1.0 - uc));
287:       f[j][i].v = udot[j][i].v - (appctx->D2 * (vxx + vyy) + uc * vc * vc - (appctx->gamma + appctx->kappa) * vc);
288:     }
289:   }
290:   PetscCall(PetscLogFlops(16.0 * xm * ym));

292:   /*
293:      Restore vectors
294:   */
296:   PetscCall(DMDAVecRestoreArray(da, F, &f));
298:   PetscCall(DMRestoreLocalVector(da, &localU));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat BB, void *ctx)
303: {
304:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
305:   DM          da;
306:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
307:   PetscReal   hx, hy, sx, sy;
308:   PetscScalar uc, vc;
309:   Field     **u;
310:   Vec         localU;
311:   MatStencil  stencil[6], rowstencil;
312:   PetscScalar entries[6];

314:   PetscFunctionBegin;
315:   PetscCall(TSGetDM(ts, &da));
316:   PetscCall(DMGetLocalVector(da, &localU));
317:   PetscCall(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));

319:   hx = 2.50 / (PetscReal)Mx;
320:   sx = 1.0 / (hx * hx);
321:   hy = 2.50 / (PetscReal)My;
322:   sy = 1.0 / (hy * hy);

324:   /*
325:      Scatter ghost points to local vector,using the 2-step process
326:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
327:      By placing code between these two statements, computations can be
328:      done while messages are in transition.
329:   */
330:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
331:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));

333:   /*
334:      Get pointers to vector data
335:   */

338:   /*
339:      Get local grid boundaries
340:   */
341:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

343:   stencil[0].k = 0;
344:   stencil[1].k = 0;
345:   stencil[2].k = 0;
346:   stencil[3].k = 0;
347:   stencil[4].k = 0;
348:   stencil[5].k = 0;
349:   rowstencil.k = 0;
350:   /*
351:      Compute function over the locally owned part of the grid
352:   */
353:   for (j = ys; j < ys + ym; j++) {
354:     stencil[0].j = j - 1;
355:     stencil[1].j = j + 1;
356:     stencil[2].j = j;
357:     stencil[3].j = j;
358:     stencil[4].j = j;
359:     stencil[5].j = j;
360:     rowstencil.k = 0;
361:     rowstencil.j = j;
362:     for (i = xs; i < xs + xm; i++) {
363:       uc = u[j][i].u;
364:       vc = u[j][i].v;

366:       /*      uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
367:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;

369:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
370:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
371:        f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/

373:       stencil[0].i = i;
374:       stencil[0].c = 0;
375:       entries[0]   = -appctx->D1 * sy;
376:       stencil[1].i = i;
377:       stencil[1].c = 0;
378:       entries[1]   = -appctx->D1 * sy;
379:       stencil[2].i = i - 1;
380:       stencil[2].c = 0;
381:       entries[2]   = -appctx->D1 * sx;
382:       stencil[3].i = i + 1;
383:       stencil[3].c = 0;
384:       entries[3]   = -appctx->D1 * sx;
385:       stencil[4].i = i;
386:       stencil[4].c = 0;
387:       entries[4]   = 2.0 * appctx->D1 * (sx + sy) + vc * vc + appctx->gamma + a;
388:       stencil[5].i = i;
389:       stencil[5].c = 1;
390:       entries[5]   = 2.0 * uc * vc;
391:       rowstencil.i = i;
392:       rowstencil.c = 0;

394:       PetscCall(MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
395:       if (appctx->aijpc) PetscCall(MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
396:       stencil[0].c = 1;
397:       entries[0]   = -appctx->D2 * sy;
398:       stencil[1].c = 1;
399:       entries[1]   = -appctx->D2 * sy;
400:       stencil[2].c = 1;
401:       entries[2]   = -appctx->D2 * sx;
402:       stencil[3].c = 1;
403:       entries[3]   = -appctx->D2 * sx;
404:       stencil[4].c = 1;
405:       entries[4]   = 2.0 * appctx->D2 * (sx + sy) - 2.0 * uc * vc + appctx->gamma + appctx->kappa + a;
406:       stencil[5].c = 0;
407:       entries[5]   = -vc * vc;
408:       rowstencil.c = 1;

410:       PetscCall(MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
411:       if (appctx->aijpc) PetscCall(MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES));
412:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
413:     }
414:   }

416:   /*
417:      Restore vectors
418:   */
419:   PetscCall(PetscLogFlops(19.0 * xm * ym));