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;
40: PetscFunctionBeginUser;
41: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
42: user.k = 1;
43: user.e = .99;
44: PetscCall(PetscOptionsGetInt(NULL, 0, "-k", &user.k, 0));
45: PetscCall(PetscOptionsGetScalar(NULL, 0, "-e", &user.e, 0));
47: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
48: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 128, 1, 1, 0, &da));
49: PetscCall(DMSetFromOptions(da));
50: PetscCall(DMSetUp(da));
51: PetscCall(KSPSetDM(ksp, da));
52: PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user));
53: PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user));
54: PetscCall(KSPSetFromOptions(ksp));
55: PetscCall(KSPSolve(ksp, NULL, NULL));
57: PetscCall(KSPGetOperators(ksp, &A, NULL));
58: PetscCall(KSPGetSolution(ksp, &x));
59: PetscCall(KSPGetRhs(ksp, &b));
60: PetscCall(VecDuplicate(b, &b2));
61: PetscCall(MatMult(A, x, b2));
62: PetscCall(VecAXPY(b2, -1.0, b));
63: PetscCall(VecNorm(b2, NORM_MAX, &nrm));
64: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)nrm));
66: PetscCall(VecDestroy(&b2));
67: PetscCall(KSPDestroy(&ksp));
68: PetscCall(DMDestroy(&da));
69: PetscCall(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;
79: PetscFunctionBeginUser;
80: PetscCall(KSPGetDM(ksp, &da));
81: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
82: h = 1.0 / ((mx - 1));
83: PetscCall(VecSet(b, h));
84: idx[0] = 0;
85: idx[1] = mx - 1;
86: v[0] = v[1] = 0.0;
87: PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES));
88: PetscCall(VecAssemblyBegin(b));
89: PetscCall(VecAssemblyEnd(b));
90: PetscFunctionReturn(PETSC_SUCCESS);
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;
101: PetscFunctionBeginUser;
102: PetscCall(KSPGetDM(ksp, &da));
103: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
104: PetscCall(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: PetscCall(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: PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES));
122: }
123: }
124: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
125: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
126: PetscFunctionReturn(PETSC_SUCCESS);
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*/