Actual source code: ex25.c


  2: /*
  3:  Partial differential equation

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

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

 12:    This uses multigrid to solve the linear system

 14: */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

129: /*TEST

131:    test:
132:       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
133:       requires: !single

135:    test:
136:       suffix: 2
137:       nsize: 2
138:       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
139:       requires: !single

141: TEST*/