Actual source code: ex25.c
1: /*
2: Partial differential equation
4: d (1 + e*sine(2*pi*k*x)) d u = 1, 0 < x < 1,
5: -- ---
6: dx dx
7: with boundary conditions
9: u = 0 for x = 0, x = 1
11: This uses multigrid to solve the linear system
13: */
15: static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.\n\n";
17: #include <petscdm.h>
18: #include <petscdmda.h>
19: #include <petscksp.h>
21: static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
22: static PetscErrorCode ComputeRHS(KSP, Vec, void *);
24: typedef struct {
25: PetscInt k;
26: PetscScalar e;
27: } AppCtx;
29: int main(int argc, char **argv)
30: {
31: KSP ksp;
32: DM da;
33: AppCtx user;
34: Mat A;
35: Vec b, b2;
36: Vec x;
37: PetscReal nrm;
39: PetscFunctionBeginUser;
40: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
41: user.k = 1;
42: user.e = .99;
43: PetscCall(PetscOptionsGetInt(NULL, 0, "-k", &user.k, 0));
44: PetscCall(PetscOptionsGetScalar(NULL, 0, "-e", &user.e, 0));
46: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
47: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 128, 1, 1, 0, &da));
48: PetscCall(DMSetFromOptions(da));
49: PetscCall(DMSetUp(da));
50: PetscCall(KSPSetDM(ksp, da));
51: PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user));
52: PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user));
53: PetscCall(KSPSetFromOptions(ksp));
54: PetscCall(KSPSolve(ksp, NULL, NULL));
56: PetscCall(KSPGetOperators(ksp, &A, NULL));
57: PetscCall(KSPGetSolution(ksp, &x));
58: PetscCall(KSPGetRhs(ksp, &b));
59: PetscCall(VecDuplicate(b, &b2));
60: PetscCall(MatMult(A, x, b2));
61: PetscCall(VecAXPY(b2, -1.0, b));
62: PetscCall(VecNorm(b2, NORM_MAX, &nrm));
63: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)nrm));
65: PetscCall(VecDestroy(&b2));
66: PetscCall(KSPDestroy(&ksp));
67: PetscCall(DMDestroy(&da));
68: PetscCall(PetscFinalize());
69: return 0;
70: }
72: static PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
73: {
74: PetscInt mx, idx[2];
75: PetscScalar h, v[2];
76: DM da;
78: PetscFunctionBeginUser;
79: PetscCall(KSPGetDM(ksp, &da));
80: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
81: h = 1.0 / (mx - 1);
82: PetscCall(VecSet(b, h));
83: idx[0] = 0;
84: idx[1] = mx - 1;
85: v[0] = v[1] = 0.0;
86: PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES));
87: PetscCall(VecAssemblyBegin(b));
88: PetscCall(VecAssemblyEnd(b));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
93: {
94: AppCtx *user = (AppCtx *)ctx;
95: PetscInt i, mx, xm, xs;
96: PetscScalar v[3], h, xlow, xhigh;
97: MatStencil row, col[3];
98: DM da;
100: PetscFunctionBeginUser;
101: PetscCall(KSPGetDM(ksp, &da));
102: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
103: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
104: h = 1.0 / (mx - 1);
106: for (i = xs; i < xs + xm; i++) {
107: row.i = i;
108: if (i == 0 || i == mx - 1) {
109: v[0] = 2.0 / h;
110: PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
111: } else {
112: xlow = h * (PetscReal)i - .5 * h;
113: xhigh = xlow + h;
114: v[0] = (-1.0 - user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xlow)) / h;
115: col[0].i = i - 1;
116: 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;
117: col[1].i = row.i;
118: v[2] = (-1.0 - user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xhigh)) / h;
119: col[2].i = i + 1;
120: PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES));
121: }
122: }
123: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
124: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*TEST
130: test:
131: 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
132: requires: !single
134: test:
135: suffix: 2
136: nsize: 2
137: 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
138: requires: !single
140: TEST*/