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