1 /*
2 Partial differential equation
3
4 d (1 + e*sine(2*pi*k*x)) d u = 1, 0 < x < 1,
5 -- ---
6 dx dx
7 with boundary conditions
8
9 u = 0 for x = 0, x = 1
10
11 This uses multigrid to solve the linear system
12
13 */
14
15 static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.\n\n";
16
17 #include <petscdm.h>
18 #include <petscdmda.h>
19 #include <petscksp.h>
20
21 static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
22 static PetscErrorCode ComputeRHS(KSP, Vec, void *);
23
24 typedef struct {
25 PetscInt k;
26 PetscScalar e;
27 } AppCtx;
28
main(int argc,char ** argv)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;
38
39 PetscFunctionBeginUser;
40 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
45
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));
55
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));
64
65 PetscCall(VecDestroy(&b2));
66 PetscCall(KSPDestroy(&ksp));
67 PetscCall(DMDestroy(&da));
68 PetscCall(PetscFinalize());
69 return 0;
70 }
71
ComputeRHS(KSP ksp,Vec b,PetscCtx ctx)72 static PetscErrorCode ComputeRHS(KSP ksp, Vec b, PetscCtx ctx)
73 {
74 PetscInt mx, idx[2];
75 PetscScalar h, v[2];
76 DM da;
77
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 }
91
ComputeMatrix(KSP ksp,Mat J,Mat jac,PetscCtx ctx)92 static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, PetscCtx 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;
99
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);
105
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 }
127
128 /*TEST
129
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
133
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
139
140 TEST*/
141