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, (char *)0, help));
56: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
57: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 12, 12, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da));
58: PetscCall(DMSetFromOptions(da));
59: PetscCall(DMSetUp(da));
60: PetscCall(DMDASetInterpolationType(da, DMDA_Q0));
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(KSPSolve(ksp, NULL, NULL));
76: PetscCall(KSPDestroy(&ksp));
77: PetscCall(DMDestroy(&da));
78: PetscCall(PetscFinalize());
79: return 0;
80: }
82: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
83: {
84: UserContext *user = (UserContext *)ctx;
85: PetscInt i, j, mx, my, xm, ym, xs, ys;
86: PetscScalar Hx, Hy;
87: PetscScalar **array;
88: DM da;
90: PetscFunctionBeginUser;
91: PetscCall(KSPGetDM(ksp, &da));
92: PetscCall(DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
93: Hx = 1.0 / (PetscReal)(mx);
94: Hy = 1.0 / (PetscReal)(my);
95: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
96: PetscCall(DMDAVecGetArray(da, b, &array));
97: for (j = ys; j < ys + ym; j++) {
98: 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;
99: }
100: PetscCall(DMDAVecRestoreArray(da, b, &array));
101: PetscCall(VecAssemblyBegin(b));
102: PetscCall(VecAssemblyEnd(b));
104: /* force right-hand side to be consistent for singular matrix */
105: /* note this is really a hack, normally the model would provide you with a consistent right handside */
106: if (user->bcType == NEUMANN) {
107: MatNullSpace nullspace;
109: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
110: PetscCall(MatNullSpaceRemove(nullspace, b));
111: PetscCall(MatNullSpaceDestroy(&nullspace));
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
117: {
118: UserContext *user = (UserContext *)ctx;
119: PetscInt i, j, mx, my, xm, ym, xs, ys, num, numi, numj;
120: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
121: MatStencil row, col[5];
122: DM da;
124: PetscFunctionBeginUser;
125: PetscCall(KSPGetDM(ksp, &da));
126: PetscCall(DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
127: Hx = 1.0 / (PetscReal)(mx);
128: Hy = 1.0 / (PetscReal)(my);
129: HxdHy = Hx / Hy;
130: HydHx = Hy / Hx;
131: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
132: for (j = ys; j < ys + ym; j++) {
133: for (i = xs; i < xs + xm; i++) {
134: row.i = i;
135: row.j = j;
136: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
137: if (user->bcType == DIRICHLET) {
138: v[0] = 2.0 * (HxdHy + HydHx);
139: PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
140: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dirichlet boundary conditions not supported !");
141: } else if (user->bcType == NEUMANN) {
142: num = 0;
143: numi = 0;
144: numj = 0;
145: if (j != 0) {
146: v[num] = -HxdHy;
147: col[num].i = i;
148: col[num].j = j - 1;
149: num++;
150: numj++;
151: }
152: if (i != 0) {
153: v[num] = -HydHx;
154: col[num].i = i - 1;
155: col[num].j = j;
156: num++;
157: numi++;
158: }
159: if (i != mx - 1) {
160: v[num] = -HydHx;
161: col[num].i = i + 1;
162: col[num].j = j;
163: num++;
164: numi++;
165: }
166: if (j != my - 1) {
167: v[num] = -HxdHy;
168: col[num].i = i;
169: col[num].j = j + 1;
170: num++;
171: numj++;
172: }
173: v[num] = (PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx;
174: col[num].i = i;
175: col[num].j = j;
176: num++;
177: PetscCall(MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES));
178: }
179: } else {
180: v[0] = -HxdHy;
181: col[0].i = i;
182: col[0].j = j - 1;
183: v[1] = -HydHx;
184: col[1].i = i - 1;
185: col[1].j = j;
186: v[2] = 2.0 * (HxdHy + HydHx);
187: col[2].i = i;
188: col[2].j = j;
189: v[3] = -HydHx;
190: col[3].i = i + 1;
191: col[3].j = j;
192: v[4] = -HxdHy;
193: col[4].i = i;
194: col[4].j = j + 1;
195: PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES));
196: }
197: }
198: }
199: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
200: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
201: if (user->bcType == NEUMANN) {
202: MatNullSpace nullspace;
204: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
205: PetscCall(MatSetNullSpace(J, nullspace));
206: PetscCall(MatNullSpaceDestroy(&nullspace));
207: }
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*TEST
213: test:
214: args: -pc_type mg -pc_mg_type full -ksp_type fgmres -ksp_monitor_short -pc_mg_levels 3 -mg_coarse_pc_factor_shift_type nonzero
216: TEST*/