Actual source code: ex28.c
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;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
21: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
22: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 3, 2, 1, 0, &da));
23: PetscCall(DMSetFromOptions(da));
24: PetscCall(DMSetUp(da));
25: PetscCall(KSPSetDM(ksp, da));
26: PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, NULL));
27: PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, NULL));
29: PetscCall(KSPSetFromOptions(ksp));
30: PetscCall(DMCreateGlobalVector(da, &x));
31: PetscCall(ComputeInitialSolution(da, x));
32: PetscCall(DMSetApplicationContext(da, x));
33: PetscCall(KSPSetUp(ksp));
34: PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
35: for (i = 0; i < 10; i++) {
36: PetscCall(KSPSolve(ksp, NULL, x));
37: PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
38: }
39: PetscCall(VecDestroy(&x));
40: PetscCall(KSPDestroy(&ksp));
41: PetscCall(DMDestroy(&da));
42: PetscCall(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];
51: PetscFunctionBeginUser;
52: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
53: Hx = 2.0 * PETSC_PI / (PetscReal)(mx);
54: PetscCall(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: PetscCall(VecSetValues(x, 2, col, val, INSERT_VALUES));
61: }
62: PetscCall(VecAssemblyBegin(x));
63: PetscCall(VecAssemblyEnd(x));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
68: {
69: PetscInt mx;
70: PetscScalar h;
71: Vec x;
72: DM da;
74: PetscFunctionBeginUser;
75: PetscCall(KSPGetDM(ksp, &da));
76: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
77: PetscCall(DMGetApplicationContext(da, &x));
78: h = 2.0 * PETSC_PI / ((mx));
79: PetscCall(VecCopy(x, b));
80: PetscCall(VecScale(b, h));
81: PetscFunctionReturn(PETSC_SUCCESS);
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;
92: PetscFunctionBeginUser;
93: PetscCall(KSPGetDM(ksp, &da));
94: PetscCall(PetscArrayzero(col, 7));
95: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
96: Hx = 2.0 * PETSC_PI / (PetscReal)(mx);
97: PetscCall(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: PetscCall(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: PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES));
129: }
130: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
132: PetscCall(MatView(jac, PETSC_VIEWER_BINARY_(PETSC_COMM_SELF)));
133: PetscFunctionReturn(PETSC_SUCCESS);
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*/