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

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

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

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

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

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

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

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

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

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

 73:   PetscFunctionBeginUser;
 74:   PetscCall(KSPGetDM(ksp, &da));
 75:   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
 76:   PetscCall(DMGetApplicationContext(da, &x));
 77:   h = 2.0 * PETSC_PI / mx;
 78:   PetscCall(VecCopy(x, b));
 79:   PetscCall(VecScale(b, h));
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

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

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

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

135: /*TEST

137:    test:
138:       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

140: TEST*/