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: */
58: PetscCall(DMDAVecGetArrayRead(da, localU, &u));
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: */
86: PetscCall(DMDAVecRestoreArrayRead(da, localU, &u));
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: */
126: PetscCall(DMDAVecGetArrayRead(da, localU, &u));
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));
210: PetscCall(DMDAVecRestoreArrayRead(da, localU, &u));
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: */
266: PetscCall(DMDAVecGetArrayRead(da, localU, &u));
267: PetscCall(DMDAVecGetArray(da, F, &f));
268: PetscCall(DMDAVecGetArrayRead(da, Udot, &udot));
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: */
295: PetscCall(DMDAVecRestoreArrayRead(da, localU, &u));
296: PetscCall(DMDAVecRestoreArray(da, F, &f));
297: PetscCall(DMDAVecRestoreArrayRead(da, Udot, &udot));
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: */
336: PetscCall(DMDAVecGetArrayRead(da, localU, &u));
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));
420: PetscCall(DMDAVecRestoreArrayRead(da, localU, &u));
421: PetscCall(DMRestoreLocalVector(da, &localU));
422: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
423: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
424: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
425: if (appctx->aijpc) {
426: PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
427: PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
428: PetscCall(MatSetOption(BB, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
429: }
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }