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: }