Actual source code: ex32.c
2: /*
3: Laplacian in 2D. Modeled by the partial differential equation
5: div grad u = f, 0 < x,y < 1,
7: with forcing function
9: f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}
11: with pure Neumann boundary conditions
13: The functions are cell-centered
15: This uses multigrid to solve the linear system
17: Contributed by Andrei Draganescu <aidraga@sandia.gov>
19: Note the nice multigrid convergence despite the fact it is only using
20: piecewise constant interpolation/restriction. This is because cell-centered multigrid
21: does not need the same rule:
23: polynomial degree(interpolation) + polynomial degree(restriction) + 2 > degree of PDE
25: that vertex based multigrid needs.
26: */
28: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
30: #include <petscdm.h>
31: #include <petscdmda.h>
32: #include <petscksp.h>
34: extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
35: extern PetscErrorCode ComputeRHS(KSP, Vec, void *);
37: typedef enum {
38: DIRICHLET,
39: NEUMANN
40: } BCType;
42: typedef struct {
43: PetscScalar nu;
44: BCType bcType;
45: } UserContext;
47: int main(int argc, char **argv)
48: {
49: KSP ksp;
50: DM da;
51: UserContext user;
52: const char *bcTypes[2] = {"dirichlet", "neumann"};
53: PetscInt bc;
56: PetscInitialize(&argc, &argv, (char *)0, help);
57: KSPCreate(PETSC_COMM_WORLD, &ksp);
58: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 12, 12, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da);
59: DMSetFromOptions(da);
60: DMSetUp(da);
61: DMDASetInterpolationType(da, DMDA_Q0);
63: KSPSetDM(ksp, da);
65: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
66: user.nu = 0.1;
67: PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", 0.1, &user.nu, NULL);
68: bc = (PetscInt)NEUMANN;
69: PetscOptionsEList("-bc_type", "Type of boundary condition", "ex29.c", bcTypes, 2, bcTypes[0], &bc, NULL);
70: user.bcType = (BCType)bc;
71: PetscOptionsEnd();
73: KSPSetComputeRHS(ksp, ComputeRHS, &user);
74: KSPSetComputeOperators(ksp, ComputeMatrix, &user);
75: KSPSetFromOptions(ksp);
76: KSPSolve(ksp, NULL, NULL);
77: KSPDestroy(&ksp);
78: DMDestroy(&da);
79: 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;
92: KSPGetDM(ksp, &da);
93: 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: DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
97: 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: DMDAVecRestoreArray(da, b, &array);
102: VecAssemblyBegin(b);
103: 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: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
111: MatNullSpaceRemove(nullspace, b);
112: MatNullSpaceDestroy(&nullspace);
113: }
114: return 0;
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;
126: KSPGetDM(ksp, &da);
127: 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: 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: 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: 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: MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES);
197: }
198: }
199: }
200: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
201: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
202: if (user->bcType == NEUMANN) {
203: MatNullSpace nullspace;
205: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
206: MatSetNullSpace(J, nullspace);
207: MatNullSpaceDestroy(&nullspace);
208: }
209: return 0;
210: }
212: /*TEST
214: test:
215: 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
217: TEST*/