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: TSGetDM(ts, &da);
38: DMGetLocalVector(da, &localU);
39: 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);
40: hx = 2.50 / (PetscReal)Mx;
41: sx = 1.0 / (hx * hx);
42: hy = 2.50 / (PetscReal)My;
43: sy = 1.0 / (hy * hy);
45: /*
46: Scatter ghost points to local vector,using the 2-step process
47: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
48: By placing code between these two statements, computations can be
49: done while messages are in transition.
50: */
51: DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
52: DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
54: /*
55: Get pointers to vector data
56: */
57: DMDAVecGetArrayRead(da, localU, &u);
58: DMDAVecGetArray(da, F, &f);
60: /*
61: Get local grid boundaries
62: */
63: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
65: /*
66: Compute function over the locally owned part of the grid
67: */
68: for (j = ys; j < ys + ym; j++) {
69: for (i = xs; i < xs + xm; i++) {
70: uc = u[j][i].u;
71: uxx = (-2.0 * uc + u[j][i - 1].u + u[j][i + 1].u) * sx;
72: uyy = (-2.0 * uc + u[j - 1][i].u + u[j + 1][i].u) * sy;
73: vc = u[j][i].v;
74: vxx = (-2.0 * vc + u[j][i - 1].v + u[j][i + 1].v) * sx;
75: vyy = (-2.0 * vc + u[j - 1][i].v + u[j + 1][i].v) * sy;
76: f[j][i].u = appctx->D1 * (uxx + uyy) - uc * vc * vc + appctx->gamma * (1.0 - uc);
77: f[j][i].v = appctx->D2 * (vxx + vyy) + uc * vc * vc - (appctx->gamma + appctx->kappa) * vc;
78: }
79: }
80: PetscLogFlops(16.0 * xm * ym);
82: /*
83: Restore vectors
84: */
85: DMDAVecRestoreArrayRead(da, localU, &u);
86: DMDAVecRestoreArray(da, F, &f);
87: DMRestoreLocalVector(da, &localU);
88: return 0;
89: }
91: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat BB, void *ctx)
92: {
93: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
94: DM da;
95: PetscInt i, j, Mx, My, xs, ys, xm, ym;
96: PetscReal hx, hy, sx, sy;
97: PetscScalar uc, vc;
98: Field **u;
99: Vec localU;
100: MatStencil stencil[6], rowstencil;
101: PetscScalar entries[6];
103: TSGetDM(ts, &da);
104: DMGetLocalVector(da, &localU);
105: 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);
107: hx = 2.50 / (PetscReal)Mx;
108: sx = 1.0 / (hx * hx);
109: hy = 2.50 / (PetscReal)My;
110: sy = 1.0 / (hy * hy);
112: /*
113: Scatter ghost points to local vector,using the 2-step process
114: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
115: By placing code between these two statements, computations can be
116: done while messages are in transition.
117: */
118: DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
119: DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
121: /*
122: Get pointers to vector data
123: */
124: DMDAVecGetArrayRead(da, localU, &u);
126: /*
127: Get local grid boundaries
128: */
129: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
131: stencil[0].k = 0;
132: stencil[1].k = 0;
133: stencil[2].k = 0;
134: stencil[3].k = 0;
135: stencil[4].k = 0;
136: stencil[5].k = 0;
137: rowstencil.k = 0;
138: /*
139: Compute function over the locally owned part of the grid
140: */
141: for (j = ys; j < ys + ym; j++) {
142: stencil[0].j = j - 1;
143: stencil[1].j = j + 1;
144: stencil[2].j = j;
145: stencil[3].j = j;
146: stencil[4].j = j;
147: stencil[5].j = j;
148: rowstencil.k = 0;
149: rowstencil.j = j;
150: for (i = xs; i < xs + xm; i++) {
151: uc = u[j][i].u;
152: vc = u[j][i].v;
154: /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
155: uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
157: vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
158: vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
159: f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/
161: stencil[0].i = i;
162: stencil[0].c = 0;
163: entries[0] = appctx->D1 * sy;
164: stencil[1].i = i;
165: stencil[1].c = 0;
166: entries[1] = appctx->D1 * sy;
167: stencil[2].i = i - 1;
168: stencil[2].c = 0;
169: entries[2] = appctx->D1 * sx;
170: stencil[3].i = i + 1;
171: stencil[3].c = 0;
172: entries[3] = appctx->D1 * sx;
173: stencil[4].i = i;
174: stencil[4].c = 0;
175: entries[4] = -2.0 * appctx->D1 * (sx + sy) - vc * vc - appctx->gamma;
176: stencil[5].i = i;
177: stencil[5].c = 1;
178: entries[5] = -2.0 * uc * vc;
179: rowstencil.i = i;
180: rowstencil.c = 0;
182: MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
183: if (appctx->aijpc) MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
184: stencil[0].c = 1;
185: entries[0] = appctx->D2 * sy;
186: stencil[1].c = 1;
187: entries[1] = appctx->D2 * sy;
188: stencil[2].c = 1;
189: entries[2] = appctx->D2 * sx;
190: stencil[3].c = 1;
191: entries[3] = appctx->D2 * sx;
192: stencil[4].c = 1;
193: entries[4] = -2.0 * appctx->D2 * (sx + sy) + 2.0 * uc * vc - appctx->gamma - appctx->kappa;
194: stencil[5].c = 0;
195: entries[5] = vc * vc;
196: rowstencil.c = 1;
198: MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
199: if (appctx->aijpc) MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
200: /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
201: }
202: }
204: /*
205: Restore vectors
206: */
207: PetscLogFlops(19.0 * xm * ym);
208: DMDAVecRestoreArrayRead(da, localU, &u);
209: DMRestoreLocalVector(da, &localU);
210: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
211: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
212: MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
213: if (appctx->aijpc) {
214: MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY);
215: MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY);
216: MatSetOption(BB, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
217: }
218: return 0;
219: }
221: /*
222: IFunction - Evaluates implicit nonlinear function, xdot - F(x).
224: Input Parameters:
225: . ts - the TS context
226: . U - input vector
227: . Udot - input vector
228: . ptr - optional user-defined context, as set by TSSetRHSFunction()
230: Output Parameter:
231: . F - function vector
232: */
233: PetscErrorCode IFunction(TS ts, PetscReal ftime, Vec U, Vec Udot, Vec F, void *ptr)
234: {
235: AppCtx *appctx = (AppCtx *)ptr;
236: DM da;
237: PetscInt i, j, Mx, My, xs, ys, xm, ym;
238: PetscReal hx, hy, sx, sy;
239: PetscScalar uc, uxx, uyy, vc, vxx, vyy;
240: Field **u, **f, **udot;
241: Vec localU;
243: TSGetDM(ts, &da);
244: DMGetLocalVector(da, &localU);
245: 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);
246: hx = 2.50 / (PetscReal)Mx;
247: sx = 1.0 / (hx * hx);
248: hy = 2.50 / (PetscReal)My;
249: sy = 1.0 / (hy * hy);
251: /*
252: Scatter ghost points to local vector,using the 2-step process
253: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
254: By placing code between these two statements, computations can be
255: done while messages are in transition.
256: */
257: DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
258: DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
260: /*
261: Get pointers to vector data
262: */
263: DMDAVecGetArrayRead(da, localU, &u);
264: DMDAVecGetArray(da, F, &f);
265: DMDAVecGetArrayRead(da, Udot, &udot);
267: /*
268: Get local grid boundaries
269: */
270: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
272: /*
273: Compute function over the locally owned part of the grid
274: */
275: for (j = ys; j < ys + ym; j++) {
276: for (i = xs; i < xs + xm; i++) {
277: uc = u[j][i].u;
278: uxx = (-2.0 * uc + u[j][i - 1].u + u[j][i + 1].u) * sx;
279: uyy = (-2.0 * uc + u[j - 1][i].u + u[j + 1][i].u) * sy;
280: vc = u[j][i].v;
281: vxx = (-2.0 * vc + u[j][i - 1].v + u[j][i + 1].v) * sx;
282: vyy = (-2.0 * vc + u[j - 1][i].v + u[j + 1][i].v) * sy;
283: f[j][i].u = udot[j][i].u - (appctx->D1 * (uxx + uyy) - uc * vc * vc + appctx->gamma * (1.0 - uc));
284: f[j][i].v = udot[j][i].v - (appctx->D2 * (vxx + vyy) + uc * vc * vc - (appctx->gamma + appctx->kappa) * vc);
285: }
286: }
287: PetscLogFlops(16.0 * xm * ym);
289: /*
290: Restore vectors
291: */
292: DMDAVecRestoreArrayRead(da, localU, &u);
293: DMDAVecRestoreArray(da, F, &f);
294: DMDAVecRestoreArrayRead(da, Udot, &udot);
295: DMRestoreLocalVector(da, &localU);
296: return 0;
297: }
299: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat BB, void *ctx)
300: {
301: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
302: DM da;
303: PetscInt i, j, Mx, My, xs, ys, xm, ym;
304: PetscReal hx, hy, sx, sy;
305: PetscScalar uc, vc;
306: Field **u;
307: Vec localU;
308: MatStencil stencil[6], rowstencil;
309: PetscScalar entries[6];
311: TSGetDM(ts, &da);
312: DMGetLocalVector(da, &localU);
313: 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);
315: hx = 2.50 / (PetscReal)Mx;
316: sx = 1.0 / (hx * hx);
317: hy = 2.50 / (PetscReal)My;
318: sy = 1.0 / (hy * hy);
320: /*
321: Scatter ghost points to local vector,using the 2-step process
322: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
323: By placing code between these two statements, computations can be
324: done while messages are in transition.
325: */
326: DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
327: DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
329: /*
330: Get pointers to vector data
331: */
332: DMDAVecGetArrayRead(da, localU, &u);
334: /*
335: Get local grid boundaries
336: */
337: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
339: stencil[0].k = 0;
340: stencil[1].k = 0;
341: stencil[2].k = 0;
342: stencil[3].k = 0;
343: stencil[4].k = 0;
344: stencil[5].k = 0;
345: rowstencil.k = 0;
346: /*
347: Compute function over the locally owned part of the grid
348: */
349: for (j = ys; j < ys + ym; j++) {
350: stencil[0].j = j - 1;
351: stencil[1].j = j + 1;
352: stencil[2].j = j;
353: stencil[3].j = j;
354: stencil[4].j = j;
355: stencil[5].j = j;
356: rowstencil.k = 0;
357: rowstencil.j = j;
358: for (i = xs; i < xs + xm; i++) {
359: uc = u[j][i].u;
360: vc = u[j][i].v;
362: /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
363: uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
365: vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
366: vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
367: f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/
369: stencil[0].i = i;
370: stencil[0].c = 0;
371: entries[0] = -appctx->D1 * sy;
372: stencil[1].i = i;
373: stencil[1].c = 0;
374: entries[1] = -appctx->D1 * sy;
375: stencil[2].i = i - 1;
376: stencil[2].c = 0;
377: entries[2] = -appctx->D1 * sx;
378: stencil[3].i = i + 1;
379: stencil[3].c = 0;
380: entries[3] = -appctx->D1 * sx;
381: stencil[4].i = i;
382: stencil[4].c = 0;
383: entries[4] = 2.0 * appctx->D1 * (sx + sy) + vc * vc + appctx->gamma + a;
384: stencil[5].i = i;
385: stencil[5].c = 1;
386: entries[5] = 2.0 * uc * vc;
387: rowstencil.i = i;
388: rowstencil.c = 0;
390: MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
391: if (appctx->aijpc) MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
392: stencil[0].c = 1;
393: entries[0] = -appctx->D2 * sy;
394: stencil[1].c = 1;
395: entries[1] = -appctx->D2 * sy;
396: stencil[2].c = 1;
397: entries[2] = -appctx->D2 * sx;
398: stencil[3].c = 1;
399: entries[3] = -appctx->D2 * sx;
400: stencil[4].c = 1;
401: entries[4] = 2.0 * appctx->D2 * (sx + sy) - 2.0 * uc * vc + appctx->gamma + appctx->kappa + a;
402: stencil[5].c = 0;
403: entries[5] = -vc * vc;
404: rowstencil.c = 1;
406: MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
407: if (appctx->aijpc) MatSetValuesStencil(BB, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
408: /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
409: }
410: }
412: /*
413: Restore vectors
414: */
415: PetscLogFlops(19.0 * xm * ym);
416: DMDAVecRestoreArrayRead(da, localU, &u);
417: DMRestoreLocalVector(da, &localU);
418: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
419: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
420: MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
421: if (appctx->aijpc) {
422: MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY);
423: MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY);
424: MatSetOption(BB, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
425: }
426: return 0;
427: }