Actual source code: ex70.c

      1: static char help[] = "Poiseuille flow problem. Viscous, laminar flow in a 2D channel with parabolic velocity\n\
      2:                       profile and linear pressure drop, exact solution of the 2D Stokes\n";

      4: /*
      5:      M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
      6:    author : Christiaan M. Klaij

      8:    Poiseuille flow problem.

     10:    Viscous, laminar flow in a 2D channel with parabolic velocity
     11:    profile and linear pressure drop, exact solution of the 2D Stokes
     12:    equations.

     14:    Discretized with the cell-centered finite-volume method on a
     15:    Cartesian grid with co-located variables. Variables ordered as
     16:    [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with
     17:    PCFIELDSPLIT. Lower factorization is used to mimic the Semi-Implicit
     18:    Method for Pressure Linked Equations (SIMPLE) used as preconditioner
     19:    instead of solver.

     21:    Disclaimer: does not contain the pressure-weighed interpolation
     22:    method needed to suppress spurious pressure modes in real-life
     23:    problems.

     25:    Usage:
     26:      mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none

     28:      Runs with PCFIELDSPLIT on 32x48 grid, no PC for the Schur
     29:      complement because A11 is zero. FGMRES is needed because
     30:      PCFIELDSPLIT is a variable preconditioner.

     32:      mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc

     34:      Same as above but with user defined PC for the true Schur
     35:      complement. PC based on the SIMPLE-type approximation (inverse of
     36:      A00 approximated by inverse of its diagonal).

     38:      mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_ksp

     40:      Replace the true Schur complement with a user defined Schur
     41:      complement based on the SIMPLE-type approximation. Same matrix is
     42:      used as PC.

     44:      mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi

     46:      Out-of-the-box SIMPLE-type preconditioning. The major advantage
     47:      is that the user neither needs to provide the approximation of
     48:      the Schur complement, nor the corresponding preconditioner.
     49: */

     51: #include <petscksp.h>

     53: typedef struct {
     54:   PetscBool userPC, userKSP, matsymmetric; /* user defined preconditioner and matrix for the Schur complement */
     55:   PetscInt  nx, ny;                        /* nb of cells in x- and y-direction */
     56:   PetscReal hx, hy;                        /* mesh size in x- and y-direction */
     57:   Mat       A;                             /* block matrix */
     58:   Mat       subA[4];                       /* the four blocks */
     59:   Mat       myS;                           /* the approximation of the Schur complement */
     60:   Vec       x, b, y;                       /* solution, rhs and temporary vector */
     61:   IS        isg[2];                        /* index sets of split "0" and "1" */
     62: } Stokes;

     64: PetscErrorCode StokesSetupMatBlock00(Stokes *); /* setup the block Q */
     65: PetscErrorCode StokesSetupMatBlock01(Stokes *); /* setup the block G */
     66: PetscErrorCode StokesSetupMatBlock10(Stokes *); /* setup the block D (equal to the transpose of G) */
     67: PetscErrorCode StokesSetupMatBlock11(Stokes *); /* setup the block C (equal to zero) */

     69: PetscErrorCode StokesGetPosition(Stokes *, PetscInt, PetscInt *, PetscInt *); /* row number j*nx+i corresponds to position (i,j) in grid */

     71: PetscErrorCode StokesStencilLaplacian(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Laplacian operator */
     72: PetscErrorCode StokesStencilGradientX(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Gradient operator (x-component) */
     73: PetscErrorCode StokesStencilGradientY(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Gradient operator (y-component) */

     75: PetscErrorCode StokesRhs(Stokes *);                                        /* rhs vector */
     76: PetscErrorCode StokesRhsMomX(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of velocity (x-component) */
     77: PetscErrorCode StokesRhsMomY(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of velocity (y-component) */
     78: PetscErrorCode StokesRhsMass(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of pressure */

     80: PetscErrorCode StokesSetupApproxSchur(Stokes *); /* approximation of the Schur complement */

     82: PetscErrorCode StokesExactSolution(Stokes *); /* exact solution vector */
     83: PetscErrorCode StokesWriteSolution(Stokes *); /* write solution to file */

     85: /* exact solution for the velocity (x-component, y-component is zero) */
     86: PetscScalar StokesExactVelocityX(const PetscScalar y)
     87: {
     88:   return 4.0 * y * (1.0 - y);
     89: }

     91: /* exact solution for the pressure */
     92: PetscScalar StokesExactPressure(const PetscScalar x)
     93: {
     94:   return 8.0 * (2.0 - x);
     95: }

     97: PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp)
     98: {
     99:   KSP     *subksp;
    100:   PC       pc;
    101:   PetscInt n = 1;

    103:   PetscFunctionBeginUser;
    104:   PetscCall(KSPGetPC(ksp, &pc));
    105:   PetscCall(PCFieldSplitSetIS(pc, "0", s->isg[0]));
    106:   PetscCall(PCFieldSplitSetIS(pc, "1", s->isg[1]));
    107:   if (s->userPC) PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS));
    108:   if (s->userKSP) {
    109:     PetscCall(PCSetUp(pc));
    110:     PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
    111:     PetscCall(KSPSetOperators(subksp[1], s->myS, s->myS));
    112:     PetscCall(PetscFree(subksp));
    113:   }
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: PetscErrorCode StokesWriteSolution(Stokes *s)
    118: {
    119:   PetscMPIInt        size;
    120:   PetscInt           n, i, j;
    121:   const PetscScalar *array;

    123:   PetscFunctionBeginUser;
    124:   /* write data (*warning* only works sequential) */
    125:   PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
    126:   if (size == 1) {
    127:     PetscViewer viewer;
    128:     PetscCall(VecGetArrayRead(s->x, &array));
    129:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer));
    130:     PetscCall(PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n"));
    131:     for (j = 0; j < s->ny; j++) {
    132:       for (i = 0; i < s->nx; i++) {
    133:         n = j * s->nx + i;
    134:         PetscCall(PetscViewerASCIIPrintf(viewer, "%.12g %.12g %.12g %.12g %.12g\n", (double)(i * s->hx + s->hx / 2), (double)(j * s->hy + s->hy / 2), (double)PetscRealPart(array[n]), (double)PetscRealPart(array[n + s->nx * s->ny]), (double)PetscRealPart(array[n + 2 * s->nx * s->ny])));
    135:       }
    136:     }
    137:     PetscCall(VecRestoreArrayRead(s->x, &array));
    138:     PetscCall(PetscViewerDestroy(&viewer));
    139:   }
    140:   PetscFunctionReturn(PETSC_SUCCESS);
    141: }

    143: PetscErrorCode StokesSetupIndexSets(Stokes *s)
    144: {
    145:   PetscFunctionBeginUser;
    146:   /* the two index sets */
    147:   PetscCall(MatNestGetISs(s->A, s->isg, NULL));
    148:   PetscFunctionReturn(PETSC_SUCCESS);
    149: }

    151: PetscErrorCode StokesSetupVectors(Stokes *s)
    152: {
    153:   PetscFunctionBeginUser;
    154:   /* solution vector x */
    155:   PetscCall(VecCreate(PETSC_COMM_WORLD, &s->x));
    156:   PetscCall(VecSetSizes(s->x, PETSC_DECIDE, 3 * s->nx * s->ny));
    157:   PetscCall(VecSetType(s->x, VECMPI));

    159:   /* exact solution y */
    160:   PetscCall(VecDuplicate(s->x, &s->y));
    161:   PetscCall(StokesExactSolution(s));

    163:   /* rhs vector b */
    164:   PetscCall(VecDuplicate(s->x, &s->b));
    165:   PetscCall(StokesRhs(s));
    166:   PetscFunctionReturn(PETSC_SUCCESS);
    167: }

    169: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
    170: {
    171:   PetscInt n;

    173:   PetscFunctionBeginUser;
    174:   /* cell number n=j*nx+i has position (i,j) in grid */
    175:   n  = row % (s->nx * s->ny);
    176:   *i = n % s->nx;
    177:   *j = (n - (*i)) / s->nx;
    178:   PetscFunctionReturn(PETSC_SUCCESS);
    179: }

    181: PetscErrorCode StokesExactSolution(Stokes *s)
    182: {
    183:   PetscInt    row, start, end, i, j;
    184:   PetscScalar val;
    185:   Vec         y0, y1;

    187:   PetscFunctionBeginUser;
    188:   /* velocity part */
    189:   PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
    190:   PetscCall(VecGetOwnershipRange(y0, &start, &end));
    191:   for (row = start; row < end; row++) {
    192:     PetscCall(StokesGetPosition(s, row, &i, &j));
    193:     if (row < s->nx * s->ny) {
    194:       val = StokesExactVelocityX(j * s->hy + s->hy / 2);
    195:     } else {
    196:       val = 0;
    197:     }
    198:     PetscCall(VecSetValue(y0, row, val, INSERT_VALUES));
    199:   }
    200:   PetscCall(VecAssemblyBegin(y0));
    201:   PetscCall(VecAssemblyEnd(y0));
    202:   PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));

    204:   /* pressure part */
    205:   PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
    206:   PetscCall(VecGetOwnershipRange(y1, &start, &end));
    207:   for (row = start; row < end; row++) {
    208:     PetscCall(StokesGetPosition(s, row, &i, &j));
    209:     val = StokesExactPressure(i * s->hx + s->hx / 2);
    210:     PetscCall(VecSetValue(y1, row, val, INSERT_VALUES));
    211:   }
    212:   PetscCall(VecAssemblyBegin(y1));
    213:   PetscCall(VecAssemblyEnd(y1));
    214:   PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: PetscErrorCode StokesRhs(Stokes *s)
    219: {
    220:   PetscInt    row, start, end, i, j;
    221:   PetscScalar val;
    222:   Vec         b0, b1;

    224:   PetscFunctionBeginUser;
    225:   /* velocity part */
    226:   PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
    227:   PetscCall(VecGetOwnershipRange(b0, &start, &end));
    228:   for (row = start; row < end; row++) {
    229:     PetscCall(StokesGetPosition(s, row, &i, &j));
    230:     if (row < s->nx * s->ny) {
    231:       PetscCall(StokesRhsMomX(s, i, j, &val));
    232:     } else {
    233:       PetscCall(StokesRhsMomY(s, i, j, &val));
    234:     }
    235:     PetscCall(VecSetValue(b0, row, val, INSERT_VALUES));
    236:   }
    237:   PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));

    239:   /* pressure part */
    240:   PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
    241:   PetscCall(VecGetOwnershipRange(b1, &start, &end));
    242:   for (row = start; row < end; row++) {
    243:     PetscCall(StokesGetPosition(s, row, &i, &j));
    244:     PetscCall(StokesRhsMass(s, i, j, &val));
    245:     if (s->matsymmetric) val = -1.0 * val;
    246:     PetscCall(VecSetValue(b1, row, val, INSERT_VALUES));
    247:   }
    248:   PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
    253: {
    254:   PetscInt    row, start, end, sz, i, j;
    255:   PetscInt    cols[5];
    256:   PetscScalar vals[5];

    258:   PetscFunctionBeginUser;
    259:   /* A[0] is 2N-by-2N */
    260:   PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0]));
    261:   PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_"));
    262:   PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny));
    263:   PetscCall(MatSetType(s->subA[0], MATMPIAIJ));
    264:   PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL));
    265:   PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end));

    267:   for (row = start; row < end; row++) {
    268:     PetscCall(StokesGetPosition(s, row, &i, &j));
    269:     /* first part: rows 0 to (nx*ny-1) */
    270:     PetscCall(StokesStencilLaplacian(s, i, j, &sz, cols, vals));
    271:     /* second part: rows (nx*ny) to (2*nx*ny-1) */
    272:     if (row >= s->nx * s->ny) {
    273:       for (i = 0; i < sz; i++) cols[i] += s->nx * s->ny;
    274:     }
    275:     for (i = 0; i < sz; i++) vals[i] = -1.0 * vals[i]; /* dynamic viscosity coef mu=-1 */
    276:     PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES));
    277:   }
    278:   PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY));
    279:   PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY));
    280:   PetscFunctionReturn(PETSC_SUCCESS);
    281: }

    283: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
    284: {
    285:   PetscInt    row, start, end, sz, i, j;
    286:   PetscInt    cols[5];
    287:   PetscScalar vals[5];

    289:   PetscFunctionBeginUser;
    290:   /* A[1] is 2N-by-N */
    291:   PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1]));
    292:   PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_"));
    293:   PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny));
    294:   PetscCall(MatSetType(s->subA[1], MATMPIAIJ));
    295:   PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL));
    296:   PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end));

    298:   PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));

    300:   for (row = start; row < end; row++) {
    301:     PetscCall(StokesGetPosition(s, row, &i, &j));
    302:     /* first part: rows 0 to (nx*ny-1) */
    303:     if (row < s->nx * s->ny) {
    304:       PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals));
    305:     } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
    306:       PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals));
    307:     }
    308:     PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES));
    309:   }
    310:   PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY));
    311:   PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
    316: {
    317:   PetscFunctionBeginUser;
    318:   /* A[2] is minus transpose of A[1] */
    319:   PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]));
    320:   if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0));
    321:   PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_"));
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }

    325: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
    326: {
    327:   PetscFunctionBeginUser;
    328:   /* A[3] is N-by-N null matrix */
    329:   PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3]));
    330:   PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_"));
    331:   PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny));
    332:   PetscCall(MatSetType(s->subA[3], MATMPIAIJ));
    333:   PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL));
    334:   PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY));
    335:   PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY));
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }

    339: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
    340: {
    341:   Vec diag;

    343:   PetscFunctionBeginUser;
    344:   /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */
    345:   /* note: A11 is zero */
    346:   /* note: in real life this matrix would be build directly, */
    347:   /* i.e. without MatMatMult */

    349:   /* inverse of diagonal of A00 */
    350:   PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
    351:   PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny));
    352:   PetscCall(VecSetType(diag, VECMPI));
    353:   PetscCall(MatGetDiagonal(s->subA[0], diag));
    354:   PetscCall(VecReciprocal(diag));

    356:   /* compute: - A10 inv(DIAGFORM(A00)) A01 */
    357:   PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */
    358:   PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DETERMINE, &s->myS));
    359:   PetscCall(MatScale(s->myS, -1.0));

    361:   /* restore A10 */
    362:   PetscCall(MatGetDiagonal(s->subA[0], diag));
    363:   PetscCall(MatDiagonalScale(s->subA[1], diag, NULL));
    364:   PetscCall(VecDestroy(&diag));
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }

    368: PetscErrorCode StokesSetupMatrix(Stokes *s)
    369: {
    370:   PetscFunctionBeginUser;
    371:   PetscCall(StokesSetupMatBlock00(s));
    372:   PetscCall(StokesSetupMatBlock01(s));
    373:   PetscCall(StokesSetupMatBlock10(s));
    374:   PetscCall(StokesSetupMatBlock11(s));
    375:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A));
    376:   PetscCall(StokesSetupApproxSchur(s));
    377:   PetscFunctionReturn(PETSC_SUCCESS);
    378: }

    380: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
    381: {
    382:   PetscInt    p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx;
    383:   PetscScalar ae = s->hy / s->hx, aeb = 0;
    384:   PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2);
    385:   PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2);
    386:   PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2);

    388:   PetscFunctionBeginUser;
    389:   if (i == 0 && j == 0) { /* south-west corner */
    390:     *sz     = 3;
    391:     cols[0] = p;
    392:     vals[0] = -(ae + awb + asb + an);
    393:     cols[1] = e;
    394:     vals[1] = ae;
    395:     cols[2] = n;
    396:     vals[2] = an;
    397:   } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
    398:     *sz     = 3;
    399:     cols[0] = s2;
    400:     vals[0] = as;
    401:     cols[1] = p;
    402:     vals[1] = -(ae + awb + as + anb);
    403:     cols[2] = e;
    404:     vals[2] = ae;
    405:   } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
    406:     *sz     = 3;
    407:     cols[0] = w;
    408:     vals[0] = aw;
    409:     cols[1] = p;
    410:     vals[1] = -(aeb + aw + asb + an);
    411:     cols[2] = n;
    412:     vals[2] = an;
    413:   } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
    414:     *sz     = 3;
    415:     cols[0] = s2;
    416:     vals[0] = as;
    417:     cols[1] = w;
    418:     vals[1] = aw;
    419:     cols[2] = p;
    420:     vals[2] = -(aeb + aw + as + anb);
    421:   } else if (i == 0) { /* west boundary */
    422:     *sz     = 4;
    423:     cols[0] = s2;
    424:     vals[0] = as;
    425:     cols[1] = p;
    426:     vals[1] = -(ae + awb + as + an);
    427:     cols[2] = e;
    428:     vals[2] = ae;
    429:     cols[3] = n;
    430:     vals[3] = an;
    431:   } else if (i == s->nx - 1) { /* east boundary */
    432:     *sz     = 4;
    433:     cols[0] = s2;
    434:     vals[0] = as;
    435:     cols[1] = w;
    436:     vals[1] = aw;
    437:     cols[2] = p;
    438:     vals[2] = -(aeb + aw + as + an);
    439:     cols[3] = n;
    440:     vals[3] = an;
    441:   } else if (j == 0) { /* south boundary */
    442:     *sz     = 4;
    443:     cols[0] = w;
    444:     vals[0] = aw;
    445:     cols[1] = p;
    446:     vals[1] = -(ae + aw + asb + an);
    447:     cols[2] = e;
    448:     vals[2] = ae;
    449:     cols[3] = n;
    450:     vals[3] = an;
    451:   } else if (j == s->ny - 1) { /* north boundary */
    452:     *sz     = 4;
    453:     cols[0] = s2;
    454:     vals[0] = as;
    455:     cols[1] = w;
    456:     vals[1] = aw;
    457:     cols[2] = p;
    458:     vals[2] = -(ae + aw + as + anb);
    459:     cols[3] = e;
    460:     vals[3] = ae;
    461:   } else { /* interior */
    462:     *sz     = 5;
    463:     cols[0] = s2;
    464:     vals[0] = as;
    465:     cols[1] = w;
    466:     vals[1] = aw;
    467:     cols[2] = p;
    468:     vals[2] = -(ae + aw + as + an);
    469:     cols[3] = e;
    470:     vals[3] = ae;
    471:     cols[4] = n;
    472:     vals[4] = an;
    473:   }
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
    478: {
    479:   PetscInt    p = j * s->nx + i, w = p - 1, e = p + 1;
    480:   PetscScalar ae = s->hy / 2, aeb = s->hy;
    481:   PetscScalar aw = -s->hy / 2, awb = 0;

    483:   PetscFunctionBeginUser;
    484:   if (i == 0 && j == 0) { /* south-west corner */
    485:     *sz     = 2;
    486:     cols[0] = p;
    487:     vals[0] = -(ae + awb);
    488:     cols[1] = e;
    489:     vals[1] = ae;
    490:   } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
    491:     *sz     = 2;
    492:     cols[0] = p;
    493:     vals[0] = -(ae + awb);
    494:     cols[1] = e;
    495:     vals[1] = ae;
    496:   } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
    497:     *sz     = 2;
    498:     cols[0] = w;
    499:     vals[0] = aw;
    500:     cols[1] = p;
    501:     vals[1] = -(aeb + aw);
    502:   } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
    503:     *sz     = 2;
    504:     cols[0] = w;
    505:     vals[0] = aw;
    506:     cols[1] = p;
    507:     vals[1] = -(aeb + aw);
    508:   } else if (i == 0) { /* west boundary */
    509:     *sz     = 2;
    510:     cols[0] = p;
    511:     vals[0] = -(ae + awb);
    512:     cols[1] = e;
    513:     vals[1] = ae;
    514:   } else if (i == s->nx - 1) { /* east boundary */
    515:     *sz     = 2;
    516:     cols[0] = w;
    517:     vals[0] = aw;
    518:     cols[1] = p;
    519:     vals[1] = -(aeb + aw);
    520:   } else if (j == 0) { /* south boundary */
    521:     *sz     = 3;
    522:     cols[0] = w;
    523:     vals[0] = aw;
    524:     cols[1] = p;
    525:     vals[1] = -(ae + aw);
    526:     cols[2] = e;
    527:     vals[2] = ae;
    528:   } else if (j == s->ny - 1) { /* north boundary */
    529:     *sz     = 3;
    530:     cols[0] = w;
    531:     vals[0] = aw;
    532:     cols[1] = p;
    533:     vals[1] = -(ae + aw);
    534:     cols[2] = e;
    535:     vals[2] = ae;
    536:   } else { /* interior */
    537:     *sz     = 3;
    538:     cols[0] = w;
    539:     vals[0] = aw;
    540:     cols[1] = p;
    541:     vals[1] = -(ae + aw);
    542:     cols[2] = e;
    543:     vals[2] = ae;
    544:   }
    545:   PetscFunctionReturn(PETSC_SUCCESS);
    546: }

    548: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
    549: {
    550:   PetscInt    p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx;
    551:   PetscScalar as = -s->hx / 2, asb = 0;
    552:   PetscScalar an = s->hx / 2, anb = 0;

    554:   PetscFunctionBeginUser;
    555:   if (i == 0 && j == 0) { /* south-west corner */
    556:     *sz     = 2;
    557:     cols[0] = p;
    558:     vals[0] = -(asb + an);
    559:     cols[1] = n;
    560:     vals[1] = an;
    561:   } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
    562:     *sz     = 2;
    563:     cols[0] = s2;
    564:     vals[0] = as;
    565:     cols[1] = p;
    566:     vals[1] = -(as + anb);
    567:   } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
    568:     *sz     = 2;
    569:     cols[0] = p;
    570:     vals[0] = -(asb + an);
    571:     cols[1] = n;
    572:     vals[1] = an;
    573:   } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
    574:     *sz     = 2;
    575:     cols[0] = s2;
    576:     vals[0] = as;
    577:     cols[1] = p;
    578:     vals[1] = -(as + anb);
    579:   } else if (i == 0) { /* west boundary */
    580:     *sz     = 3;
    581:     cols[0] = s2;
    582:     vals[0] = as;
    583:     cols[1] = p;
    584:     vals[1] = -(as + an);
    585:     cols[2] = n;
    586:     vals[2] = an;
    587:   } else if (i == s->nx - 1) { /* east boundary */
    588:     *sz     = 3;
    589:     cols[0] = s2;
    590:     vals[0] = as;
    591:     cols[1] = p;
    592:     vals[1] = -(as + an);
    593:     cols[2] = n;
    594:     vals[2] = an;
    595:   } else if (j == 0) { /* south boundary */
    596:     *sz     = 2;
    597:     cols[0] = p;
    598:     vals[0] = -(asb + an);
    599:     cols[1] = n;
    600:     vals[1] = an;
    601:   } else if (j == s->ny - 1) { /* north boundary */
    602:     *sz     = 2;
    603:     cols[0] = s2;
    604:     vals[0] = as;
    605:     cols[1] = p;
    606:     vals[1] = -(as + anb);
    607:   } else { /* interior */
    608:     *sz     = 3;
    609:     cols[0] = s2;
    610:     vals[0] = as;
    611:     cols[1] = p;
    612:     vals[1] = -(as + an);
    613:     cols[2] = n;
    614:     vals[2] = an;
    615:   }
    616:   PetscFunctionReturn(PETSC_SUCCESS);
    617: }

    619: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
    620: {
    621:   PetscScalar y   = j * s->hy + s->hy / 2;
    622:   PetscScalar awb = s->hy / (s->hx / 2);

    624:   PetscFunctionBeginUser;
    625:   if (i == 0) { /* west boundary */
    626:     *val = awb * StokesExactVelocityX(y);
    627:   } else {
    628:     *val = 0.0;
    629:   }
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }

    633: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
    634: {
    635:   PetscFunctionBeginUser;
    636:   *val = 0.0;
    637:   PetscFunctionReturn(PETSC_SUCCESS);
    638: }

    640: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
    641: {
    642:   PetscScalar y   = j * s->hy + s->hy / 2;
    643:   PetscScalar aeb = s->hy;

    645:   PetscFunctionBeginUser;
    646:   if (i == 0) { /* west boundary */
    647:     *val = aeb * StokesExactVelocityX(y);
    648:   } else {
    649:     *val = 0.0;
    650:   }
    651:   PetscFunctionReturn(PETSC_SUCCESS);
    652: }

    654: PetscErrorCode StokesCalcResidual(Stokes *s)
    655: {
    656:   PetscReal val;
    657:   Vec       b0, b1;

    659:   PetscFunctionBeginUser;
    660:   /* residual Ax-b (*warning* overwrites b) */
    661:   PetscCall(VecScale(s->b, -1.0));
    662:   PetscCall(MatMultAdd(s->A, s->x, s->b, s->b));

    664:   /* residual velocity */
    665:   PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
    666:   PetscCall(VecNorm(b0, NORM_2, &val));
    667:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val));
    668:   PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));

    670:   /* residual pressure */
    671:   PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
    672:   PetscCall(VecNorm(b1, NORM_2, &val));
    673:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val));
    674:   PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));

    676:   /* total residual */
    677:   PetscCall(VecNorm(s->b, NORM_2, &val));
    678:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val));
    679:   PetscFunctionReturn(PETSC_SUCCESS);
    680: }

    682: PetscErrorCode StokesCalcError(Stokes *s)
    683: {
    684:   PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny);
    685:   PetscReal   val;
    686:   Vec         y0, y1;

    688:   PetscFunctionBeginUser;
    689:   /* error y-x */
    690:   PetscCall(VecAXPY(s->y, -1.0, s->x));

    692:   /* error in velocity */
    693:   PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
    694:   PetscCall(VecNorm(y0, NORM_2, &val));
    695:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale))));
    696:   PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));

    698:   /* error in pressure */
    699:   PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
    700:   PetscCall(VecNorm(y1, NORM_2, &val));
    701:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale))));
    702:   PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));

    704:   /* total error */
    705:   PetscCall(VecNorm(s->y, NORM_2, &val));
    706:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart(val / scale)));
    707:   PetscFunctionReturn(PETSC_SUCCESS);
    708: }

    710: int main(int argc, char **argv)
    711: {
    712:   Stokes s;
    713:   KSP    ksp;

    715:   PetscFunctionBeginUser;
    716:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    717:   s.nx = 4;
    718:   s.ny = 6;
    719:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL));
    720:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL));
    721:   s.hx           = 2.0 / s.nx;
    722:   s.hy           = 1.0 / s.ny;
    723:   s.matsymmetric = PETSC_FALSE;
    724:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL));
    725:   s.userPC = s.userKSP = PETSC_FALSE;
    726:   PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC));
    727:   PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP));

    729:   PetscCall(StokesSetupMatrix(&s));
    730:   PetscCall(StokesSetupIndexSets(&s));
    731:   PetscCall(StokesSetupVectors(&s));

    733:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    734:   PetscCall(KSPSetOperators(ksp, s.A, s.A));
    735:   PetscCall(KSPSetFromOptions(ksp));
    736:   PetscCall(StokesSetupPC(&s, ksp));
    737:   PetscCall(KSPSolve(ksp, s.b, s.x));

    739:   /* don't trust, verify! */
    740:   PetscCall(StokesCalcResidual(&s));
    741:   PetscCall(StokesCalcError(&s));
    742:   PetscCall(StokesWriteSolution(&s));

    744:   PetscCall(KSPDestroy(&ksp));
    745:   PetscCall(MatDestroy(&s.subA[0]));
    746:   PetscCall(MatDestroy(&s.subA[1]));
    747:   PetscCall(MatDestroy(&s.subA[2]));
    748:   PetscCall(MatDestroy(&s.subA[3]));
    749:   PetscCall(MatDestroy(&s.A));
    750:   PetscCall(VecDestroy(&s.x));
    751:   PetscCall(VecDestroy(&s.b));
    752:   PetscCall(VecDestroy(&s.y));
    753:   PetscCall(MatDestroy(&s.myS));
    754:   PetscCall(PetscFinalize());
    755:   return 0;
    756: }

    758: /*TEST

    760:    test:
    761:       nsize: 2
    762:       args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none

    764:    test:
    765:       suffix: 2
    766:       nsize: 2
    767:       args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc

    769:    test:
    770:       suffix: 3
    771:       nsize: 2
    772:       args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc

    774:    test:
    775:       suffix: 4
    776:       nsize: 2
    777:       args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi

    779:    test:
    780:       suffix: 4_pcksp
    781:       nsize: 2
    782:       args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi

    784:    test:
    785:       suffix: 5
    786:       nsize: 2
    787:       args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type gkb -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-10

    789:    test:
    790:       suffix: 6
    791:       nsize: 2
    792:       args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type gkb -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-10

    794:    test:
    795:       suffix: 7
    796:       nsize: 2
    797:       args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type gkb -pc_fieldsplit_gkb_tol 1e-4 -pc_fieldsplit_gkb_nu 5 -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-6

    799: TEST*/