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