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*/