Actual source code: ex25.c

  1: /*
  2:  Partial differential equation

  4:    d  (1 + e*sine(2*pi*k*x)) d u = 1, 0 < x < 1,
  5:    --                        ---
  6:    dx                        dx
  7: with boundary conditions

  9:    u = 0 for x = 0, x = 1

 11:    This uses multigrid to solve the linear system

 13: */

 15: static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.\n\n";

 17: #include <petscdm.h>
 18: #include <petscdmda.h>
 19: #include <petscksp.h>

 21: static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
 22: static PetscErrorCode ComputeRHS(KSP, Vec, void *);

 24: typedef struct {
 25:   PetscInt    k;
 26:   PetscScalar e;
 27: } AppCtx;

 29: int main(int argc, char **argv)
 30: {
 31:   KSP       ksp;
 32:   DM        da;
 33:   AppCtx    user;
 34:   Mat       A;
 35:   Vec       b, b2;
 36:   Vec       x;
 37:   PetscReal nrm;

 39:   PetscFunctionBeginUser;
 40:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 41:   user.k = 1;
 42:   user.e = .99;
 43:   PetscCall(PetscOptionsGetInt(NULL, 0, "-k", &user.k, 0));
 44:   PetscCall(PetscOptionsGetScalar(NULL, 0, "-e", &user.e, 0));

 46:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 47:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 128, 1, 1, 0, &da));
 48:   PetscCall(DMSetFromOptions(da));
 49:   PetscCall(DMSetUp(da));
 50:   PetscCall(KSPSetDM(ksp, da));
 51:   PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user));
 52:   PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user));
 53:   PetscCall(KSPSetFromOptions(ksp));
 54:   PetscCall(KSPSolve(ksp, NULL, NULL));

 56:   PetscCall(KSPGetOperators(ksp, &A, NULL));
 57:   PetscCall(KSPGetSolution(ksp, &x));
 58:   PetscCall(KSPGetRhs(ksp, &b));
 59:   PetscCall(VecDuplicate(b, &b2));
 60:   PetscCall(MatMult(A, x, b2));
 61:   PetscCall(VecAXPY(b2, -1.0, b));
 62:   PetscCall(VecNorm(b2, NORM_MAX, &nrm));
 63:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)nrm));

 65:   PetscCall(VecDestroy(&b2));
 66:   PetscCall(KSPDestroy(&ksp));
 67:   PetscCall(DMDestroy(&da));
 68:   PetscCall(PetscFinalize());
 69:   return 0;
 70: }

 72: static PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
 73: {
 74:   PetscInt    mx, idx[2];
 75:   PetscScalar h, v[2];
 76:   DM          da;

 78:   PetscFunctionBeginUser;
 79:   PetscCall(KSPGetDM(ksp, &da));
 80:   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
 81:   h = 1.0 / (mx - 1);
 82:   PetscCall(VecSet(b, h));
 83:   idx[0] = 0;
 84:   idx[1] = mx - 1;
 85:   v[0] = v[1] = 0.0;
 86:   PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES));
 87:   PetscCall(VecAssemblyBegin(b));
 88:   PetscCall(VecAssemblyEnd(b));
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
 93: {
 94:   AppCtx     *user = (AppCtx *)ctx;
 95:   PetscInt    i, mx, xm, xs;
 96:   PetscScalar v[3], h, xlow, xhigh;
 97:   MatStencil  row, col[3];
 98:   DM          da;

100:   PetscFunctionBeginUser;
101:   PetscCall(KSPGetDM(ksp, &da));
102:   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
103:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
104:   h = 1.0 / (mx - 1);

106:   for (i = xs; i < xs + xm; i++) {
107:     row.i = i;
108:     if (i == 0 || i == mx - 1) {
109:       v[0] = 2.0 / h;
110:       PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
111:     } else {
112:       xlow     = h * (PetscReal)i - .5 * h;
113:       xhigh    = xlow + h;
114:       v[0]     = (-1.0 - user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xlow)) / h;
115:       col[0].i = i - 1;
116:       v[1]     = (2.0 + user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xlow) + user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xhigh)) / h;
117:       col[1].i = row.i;
118:       v[2]     = (-1.0 - user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xhigh)) / h;
119:       col[2].i = i + 1;
120:       PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES));
121:     }
122:   }
123:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
124:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: /*TEST

130:    test:
131:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -ksp_view -pc_mg_type full
132:       requires: !single

134:    test:
135:       suffix: 2
136:       nsize: 2
137:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -ksp_view -pc_mg_type full
138:       requires: !single

140: TEST*/