Actual source code: ex32.c
1: /*
2: Laplacian in 2D. Modeled by the partial differential equation
4: div grad u = f, 0 < x,y < 1,
6: with forcing function
8: f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}
10: with pure Neumann boundary conditions
12: The functions are cell-centered
14: This uses multigrid to solve the linear system
16: Contributed by Andrei Draganescu <aidraga@sandia.gov>
18: Note the nice multigrid convergence despite the fact it is only using
19: piecewise constant interpolation/restriction. This is because cell-centered multigrid
20: does not need the same rule:
22: polynomial degree(interpolation) + polynomial degree(restriction) + 2 > degree of PDE
24: that vertex based multigrid needs.
25: */
27: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
29: #include <petscdm.h>
30: #include <petscdmda.h>
31: #include <petscksp.h>
33: extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
34: extern PetscErrorCode ComputeRHS(KSP, Vec, void *);
36: typedef enum {
37: DIRICHLET,
38: NEUMANN
39: } BCType;
41: typedef struct {
42: PetscScalar nu;
43: BCType bcType;
44: } UserContext;
46: int main(int argc, char **argv)
47: {
48: KSP ksp;
49: DM da;
50: UserContext user;
51: const char *bcTypes[2] = {"dirichlet", "neumann"};
52: PetscInt bc;
54: PetscFunctionBeginUser;
55: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
56: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
57: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 2, 2, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da));
58: PetscCall(DMDASetInterpolationType(da, DMDA_Q0));
59: PetscCall(DMSetFromOptions(da));
60: PetscCall(DMSetUp(da));
62: PetscCall(KSPSetDM(ksp, da));
64: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
65: user.nu = 0.1;
66: PetscCall(PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", 0.1, &user.nu, NULL));
67: bc = (PetscInt)NEUMANN;
68: PetscCall(PetscOptionsEList("-bc_type", "Type of boundary condition", "ex29.c", bcTypes, 2, bcTypes[0], &bc, NULL));
69: user.bcType = (BCType)bc;
70: PetscOptionsEnd();
72: PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user));
73: PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user));
74: PetscCall(KSPSetFromOptions(ksp));
75: PetscCall(KSPSetUp(ksp));
76: PetscCall(KSPSolve(ksp, NULL, NULL));
77: PetscCall(KSPDestroy(&ksp));
78: PetscCall(DMDestroy(&da));
79: PetscCall(PetscFinalize());
80: return 0;
81: }
83: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
84: {
85: UserContext *user = (UserContext *)ctx;
86: PetscInt i, j, mx, my, xm, ym, xs, ys;
87: PetscScalar Hx, Hy;
88: PetscScalar **array;
89: DM da;
91: PetscFunctionBeginUser;
92: PetscCall(KSPGetDM(ksp, &da));
93: PetscCall(DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
94: Hx = 1.0 / (PetscReal)mx;
95: Hy = 1.0 / (PetscReal)my;
96: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
97: PetscCall(DMDAVecGetArray(da, b, &array));
98: for (j = ys; j < ys + ym; j++) {
99: for (i = xs; i < xs + xm; i++) array[j][i] = PetscExpScalar(-(((PetscReal)i + 0.5) * Hx) * (((PetscReal)i + 0.5) * Hx) / user->nu) * PetscExpScalar(-(((PetscReal)j + 0.5) * Hy) * (((PetscReal)j + 0.5) * Hy) / user->nu) * Hx * Hy;
100: }
101: PetscCall(DMDAVecRestoreArray(da, b, &array));
102: PetscCall(VecAssemblyBegin(b));
103: PetscCall(VecAssemblyEnd(b));
105: /* force right-hand side to be consistent for singular matrix */
106: /* note this is really a hack, normally the model would provide you with a consistent right handside */
107: if (user->bcType == NEUMANN) {
108: MatNullSpace nullspace;
110: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
111: PetscCall(MatNullSpaceRemove(nullspace, b));
112: PetscCall(MatNullSpaceDestroy(&nullspace));
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
118: {
119: UserContext *user = (UserContext *)ctx;
120: PetscInt i, j, mx, my, xm, ym, xs, ys, num, numi, numj;
121: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
122: MatStencil row, col[5];
123: DM da;
125: PetscFunctionBeginUser;
126: PetscCall(KSPGetDM(ksp, &da));
127: PetscCall(DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
128: Hx = 1.0 / (PetscReal)mx;
129: Hy = 1.0 / (PetscReal)my;
130: HxdHy = Hx / Hy;
131: HydHx = Hy / Hx;
132: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
133: for (j = ys; j < ys + ym; j++) {
134: for (i = xs; i < xs + xm; i++) {
135: row.i = i;
136: row.j = j;
137: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
138: if (user->bcType == DIRICHLET) {
139: v[0] = 2.0 * (HxdHy + HydHx);
140: PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
141: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dirichlet boundary conditions not supported !");
142: } else if (user->bcType == NEUMANN) {
143: num = 0;
144: numi = 0;
145: numj = 0;
146: if (j != 0) {
147: v[num] = -HxdHy;
148: col[num].i = i;
149: col[num].j = j - 1;
150: num++;
151: numj++;
152: }
153: if (i != 0) {
154: v[num] = -HydHx;
155: col[num].i = i - 1;
156: col[num].j = j;
157: num++;
158: numi++;
159: }
160: if (i != mx - 1) {
161: v[num] = -HydHx;
162: col[num].i = i + 1;
163: col[num].j = j;
164: num++;
165: numi++;
166: }
167: if (j != my - 1) {
168: v[num] = -HxdHy;
169: col[num].i = i;
170: col[num].j = j + 1;
171: num++;
172: numj++;
173: }
174: v[num] = (PetscReal)numj * HxdHy + (PetscReal)numi * HydHx;
175: col[num].i = i;
176: col[num].j = j;
177: num++;
178: PetscCall(MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES));
179: }
180: } else {
181: v[0] = -HxdHy;
182: col[0].i = i;
183: col[0].j = j - 1;
184: v[1] = -HydHx;
185: col[1].i = i - 1;
186: col[1].j = j;
187: v[2] = 2.0 * (HxdHy + HydHx);
188: col[2].i = i;
189: col[2].j = j;
190: v[3] = -HydHx;
191: col[3].i = i + 1;
192: col[3].j = j;
193: v[4] = -HxdHy;
194: col[4].i = i;
195: col[4].j = j + 1;
196: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES));
197: }
198: }
199: }
200: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
201: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
202: if (user->bcType == NEUMANN) {
203: MatNullSpace nullspace;
205: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
206: PetscCall(MatSetNullSpace(J, nullspace));
207: PetscCall(MatNullSpaceDestroy(&nullspace));
208: }
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*TEST
214: test:
215: suffix: 2
216: requires: !single
217: args: -pc_type mg -pc_mg_levels 5 -ksp_monitor_true_residual -ksp_rtol 1.e-10 -ksp_type cg -mg_levels_pc_type sor -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 2 -mg_coarse_pc_type svd -da_refine 5
219: test:
220: suffix: 3
221: requires: !single
222: args: -pc_type mg -pc_mg_levels 2 -ksp_monitor_true_residual -ksp_rtol 1.e-10 -ksp_type cg -mg_levels_pc_type sor -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 2 -mg_coarse_pc_type svd -da_refine 5
224: test:
225: suffix: 4
226: requires: !single
227: args: -pc_type mg -pc_mg_levels 2 -ksp_monitor_true_residual -ksp_rtol 1.e-10 -ksp_type cg -mg_levels_pc_type sor -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 2 -mg_coarse_pc_type svd -da_refine 4
229: TEST*/