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;
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*/