2: static char help[] = "Solves 1D wave equation using multigrid.\n\n";

  4: #include <petscdm.h>
  5: #include <petscdmda.h>
  6: #include <petscksp.h>

  8: extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
  9: extern PetscErrorCode ComputeRHS(KSP, Vec, void *);
 10: extern PetscErrorCode ComputeInitialSolution(DM, Vec);

 12: int main(int argc, char **argv)
 13: {
 14:   PetscInt i;
 15:   KSP      ksp;
 16:   DM       da;
 17:   Vec      x;

 20:   PetscInitialize(&argc, &argv, (char *)0, help);
 21:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 22:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 3, 2, 1, 0, &da);
 23:   DMSetFromOptions(da);
 24:   DMSetUp(da);
 25:   KSPSetDM(ksp, da);
 26:   KSPSetComputeRHS(ksp, ComputeRHS, NULL);
 27:   KSPSetComputeOperators(ksp, ComputeMatrix, NULL);

 29:   KSPSetFromOptions(ksp);
 30:   DMCreateGlobalVector(da, &x);
 31:   ComputeInitialSolution(da, x);
 32:   DMSetApplicationContext(da, x);
 33:   KSPSetUp(ksp);
 34:   VecView(x, PETSC_VIEWER_DRAW_WORLD);
 35:   for (i = 0; i < 10; i++) {
 36:     KSPSolve(ksp, NULL, x);
 37:     VecView(x, PETSC_VIEWER_DRAW_WORLD);
 38:   }
 39:   VecDestroy(&x);
 40:   KSPDestroy(&ksp);
 41:   DMDestroy(&da);
 42:   PetscFinalize();
 43:   return 0;
 44: }

 46: PetscErrorCode ComputeInitialSolution(DM da, Vec x)
 47: {
 48:   PetscInt    mx, col[2], xs, xm, i;
 49:   PetscScalar Hx, val[2];

 52:   DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 53:   Hx = 2.0 * PETSC_PI / (PetscReal)(mx);
 54:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);

 56:   for (i = xs; i < xs + xm; i++) {
 57:     col[0] = 2 * i;
 58:     col[1] = 2 * i + 1;
 59:     val[0] = val[1] = PetscSinScalar(((PetscScalar)i) * Hx);
 60:     VecSetValues(x, 2, col, val, INSERT_VALUES);
 61:   }
 62:   VecAssemblyBegin(x);
 63:   VecAssemblyEnd(x);
 64:   return 0;
 65: }

 67: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
 68: {
 69:   PetscInt    mx;
 70:   PetscScalar h;
 71:   Vec         x;
 72:   DM          da;

 75:   KSPGetDM(ksp, &da);
 76:   DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 77:   DMGetApplicationContext(da, &x);
 78:   h = 2.0 * PETSC_PI / ((mx));
 79:   VecCopy(x, b);
 80:   VecScale(b, h);
 81:   return 0;
 82: }

 84: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
 85: {
 86:   PetscInt    i, mx, xm, xs;
 87:   PetscScalar v[7], Hx;
 88:   MatStencil  row, col[7];
 89:   PetscScalar lambda;
 90:   DM          da;

 93:   KSPGetDM(ksp, &da);
 94:   PetscArrayzero(col, 7);
 95:   DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 96:   Hx = 2.0 * PETSC_PI / (PetscReal)(mx);
 97:   DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
 98:   lambda = 2.0 * Hx;
 99:   for (i = xs; i < xs + xm; i++) {
100:     row.i    = i;
101:     row.j    = 0;
102:     row.k    = 0;
103:     row.c    = 0;
104:     v[0]     = Hx;
105:     col[0].i = i;
106:     col[0].c = 0;
107:     v[1]     = lambda;
108:     col[1].i = i - 1;
109:     col[1].c = 1;
110:     v[2]     = -lambda;
111:     col[2].i = i + 1;
112:     col[2].c = 1;
113:     MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES);

115:     row.i    = i;
116:     row.j    = 0;
117:     row.k    = 0;
118:     row.c    = 1;
119:     v[0]     = lambda;
120:     col[0].i = i - 1;
121:     col[0].c = 0;
122:     v[1]     = Hx;
123:     col[1].i = i;
124:     col[1].c = 1;
125:     v[2]     = -lambda;
126:     col[2].i = i + 1;
127:     col[2].c = 0;
128:     MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES);
129:   }
130:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
131:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
132:   MatView(jac, PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
133:   return 0;
134: }

136: /*TEST

138:    test:
139:       args: -ksp_monitor_short -pc_type mg -pc_mg_type full -ksp_type fgmres -da_refine 2 -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 1 -mg_levels_pc_type ilu

141: TEST*/