Actual source code: ex32.c

      1: /*
      2: Laplacian in 2D. Modeled by the partial differential equation

      4:    div  grad u = f,  0 < x,y < 1,

      6: with forcing function

      8:    f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}

     10: with pure Neumann boundary conditions

     12: The functions are cell-centered

     14: This uses multigrid to solve the linear system

     16:        Contributed by Andrei Draganescu <aidraga@sandia.gov>

     18: Note the nice multigrid convergence despite the fact it is only using
     19: piecewise constant interpolation/restriction. This is because cell-centered multigrid
     20: does not need the same rule:

     22:     polynomial degree(interpolation) + polynomial degree(restriction) + 2 > degree of PDE

     24: that vertex based multigrid needs.
     25: */

     27: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

     29: #include <petscdm.h>
     30: #include <petscdmda.h>
     31: #include <petscksp.h>

     33: extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
     34: extern PetscErrorCode ComputeRHS(KSP, Vec, void *);

     36: typedef enum {
     37:   DIRICHLET,
     38:   NEUMANN
     39: } BCType;

     41: typedef struct {
     42:   PetscScalar nu;
     43:   BCType      bcType;
     44: } UserContext;

     46: int main(int argc, char **argv)
     47: {
     48:   KSP         ksp;
     49:   DM          da;
     50:   UserContext user;
     51:   const char *bcTypes[2] = {"dirichlet", "neumann"};
     52:   PetscInt    bc;

     54:   PetscFunctionBeginUser;
     55:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     56:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
     57:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 2, 2, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da));
     58:   PetscCall(DMDASetInterpolationType(da, DMDA_Q0));
     59:   PetscCall(DMSetFromOptions(da));
     60:   PetscCall(DMSetUp(da));

     62:   PetscCall(KSPSetDM(ksp, da));

     64:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
     65:   user.nu = 0.1;
     66:   PetscCall(PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", 0.1, &user.nu, NULL));
     67:   bc = (PetscInt)NEUMANN;
     68:   PetscCall(PetscOptionsEList("-bc_type", "Type of boundary condition", "ex29.c", bcTypes, 2, bcTypes[0], &bc, NULL));
     69:   user.bcType = (BCType)bc;
     70:   PetscOptionsEnd();

     72:   PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user));
     73:   PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user));
     74:   PetscCall(KSPSetFromOptions(ksp));
     75:   PetscCall(KSPSetUp(ksp));
     76:   PetscCall(KSPSolve(ksp, NULL, NULL));
     77:   PetscCall(KSPDestroy(&ksp));
     78:   PetscCall(DMDestroy(&da));
     79:   PetscCall(PetscFinalize());
     80:   return 0;
     81: }

     83: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
     84: {
     85:   UserContext  *user = (UserContext *)ctx;
     86:   PetscInt      i, j, mx, my, xm, ym, xs, ys;
     87:   PetscScalar   Hx, Hy;
     88:   PetscScalar **array;
     89:   DM            da;

     91:   PetscFunctionBeginUser;
     92:   PetscCall(KSPGetDM(ksp, &da));
     93:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
     94:   Hx = 1.0 / (PetscReal)mx;
     95:   Hy = 1.0 / (PetscReal)my;
     96:   PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
     97:   PetscCall(DMDAVecGetArray(da, b, &array));
     98:   for (j = ys; j < ys + ym; j++) {
     99:     for (i = xs; i < xs + xm; i++) array[j][i] = PetscExpScalar(-(((PetscReal)i + 0.5) * Hx) * (((PetscReal)i + 0.5) * Hx) / user->nu) * PetscExpScalar(-(((PetscReal)j + 0.5) * Hy) * (((PetscReal)j + 0.5) * Hy) / user->nu) * Hx * Hy;
    100:   }
    101:   PetscCall(DMDAVecRestoreArray(da, b, &array));
    102:   PetscCall(VecAssemblyBegin(b));
    103:   PetscCall(VecAssemblyEnd(b));

    105:   /* force right-hand side to be consistent for singular matrix */
    106:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
    107:   if (user->bcType == NEUMANN) {
    108:     MatNullSpace nullspace;

    110:     PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
    111:     PetscCall(MatNullSpaceRemove(nullspace, b));
    112:     PetscCall(MatNullSpaceDestroy(&nullspace));
    113:   }
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
    118: {
    119:   UserContext *user = (UserContext *)ctx;
    120:   PetscInt     i, j, mx, my, xm, ym, xs, ys, num, numi, numj;
    121:   PetscScalar  v[5], Hx, Hy, HydHx, HxdHy;
    122:   MatStencil   row, col[5];
    123:   DM           da;

    125:   PetscFunctionBeginUser;
    126:   PetscCall(KSPGetDM(ksp, &da));
    127:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    128:   Hx    = 1.0 / (PetscReal)mx;
    129:   Hy    = 1.0 / (PetscReal)my;
    130:   HxdHy = Hx / Hy;
    131:   HydHx = Hy / Hx;
    132:   PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
    133:   for (j = ys; j < ys + ym; j++) {
    134:     for (i = xs; i < xs + xm; i++) {
    135:       row.i = i;
    136:       row.j = j;
    137:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
    138:         if (user->bcType == DIRICHLET) {
    139:           v[0] = 2.0 * (HxdHy + HydHx);
    140:           PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
    141:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dirichlet boundary conditions not supported !");
    142:         } else if (user->bcType == NEUMANN) {
    143:           num  = 0;
    144:           numi = 0;
    145:           numj = 0;
    146:           if (j != 0) {
    147:             v[num]     = -HxdHy;
    148:             col[num].i = i;
    149:             col[num].j = j - 1;
    150:             num++;
    151:             numj++;
    152:           }
    153:           if (i != 0) {
    154:             v[num]     = -HydHx;
    155:             col[num].i = i - 1;
    156:             col[num].j = j;
    157:             num++;
    158:             numi++;
    159:           }
    160:           if (i != mx - 1) {
    161:             v[num]     = -HydHx;
    162:             col[num].i = i + 1;
    163:             col[num].j = j;
    164:             num++;
    165:             numi++;
    166:           }
    167:           if (j != my - 1) {
    168:             v[num]     = -HxdHy;
    169:             col[num].i = i;
    170:             col[num].j = j + 1;
    171:             num++;
    172:             numj++;
    173:           }
    174:           v[num]     = (PetscReal)numj * HxdHy + (PetscReal)numi * HydHx;
    175:           col[num].i = i;
    176:           col[num].j = j;
    177:           num++;
    178:           PetscCall(MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES));
    179:         }
    180:       } else {
    181:         v[0]     = -HxdHy;
    182:         col[0].i = i;
    183:         col[0].j = j - 1;
    184:         v[1]     = -HydHx;
    185:         col[1].i = i - 1;
    186:         col[1].j = j;
    187:         v[2]     = 2.0 * (HxdHy + HydHx);
    188:         col[2].i = i;
    189:         col[2].j = j;
    190:         v[3]     = -HydHx;
    191:         col[3].i = i + 1;
    192:         col[3].j = j;
    193:         v[4]     = -HxdHy;
    194:         col[4].i = i;
    195:         col[4].j = j + 1;
    196:         PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES));
    197:       }
    198:     }
    199:   }
    200:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    201:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    202:   if (user->bcType == NEUMANN) {
    203:     MatNullSpace nullspace;

    205:     PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
    206:     PetscCall(MatSetNullSpace(J, nullspace));
    207:     PetscCall(MatNullSpaceDestroy(&nullspace));
    208:   }
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: /*TEST

    214:    test:
    215:      suffix: 2
    216:      requires: !single
    217:      args: -pc_type mg -pc_mg_levels 5 -ksp_monitor_true_residual -ksp_rtol 1.e-10 -ksp_type cg -mg_levels_pc_type sor -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 2 -mg_coarse_pc_type svd -da_refine 5

    219:    test:
    220:      suffix: 3
    221:      requires: !single
    222:      args: -pc_type mg -pc_mg_levels 2 -ksp_monitor_true_residual -ksp_rtol 1.e-10 -ksp_type cg -mg_levels_pc_type sor -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 2 -mg_coarse_pc_type svd -da_refine 5

    224:    test:
    225:      suffix: 4
    226:      requires: !single
    227:      args: -pc_type mg -pc_mg_levels 2 -ksp_monitor_true_residual -ksp_rtol 1.e-10 -ksp_type cg -mg_levels_pc_type sor -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 2 -mg_coarse_pc_type svd -da_refine 4

    229: TEST*/