Actual source code: ex29.c


  2: /*
  3: Added at the request of Marc Garbey.

  5: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation

  7:    -div \rho grad u = f,  0 < x,y < 1,

  9: with forcing function

 11:    f = e^{-x^2/\nu} e^{-y^2/\nu}

 13: with Dirichlet boundary conditions

 15:    u = f(x,y) for x = 0, x = 1, y = 0, y = 1

 17: or pure Neumman boundary conditions

 19: This uses multigrid to solve the linear system
 20: */

 22: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

 24: #include <petscdm.h>
 25: #include <petscdmda.h>
 26: #include <petscksp.h>

 28: extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
 29: extern PetscErrorCode ComputeRHS(KSP, Vec, void *);

 31: typedef enum {
 32:   DIRICHLET,
 33:   NEUMANN
 34: } BCType;

 36: typedef struct {
 37:   PetscReal rho;
 38:   PetscReal nu;
 39:   BCType    bcType;
 40: } UserContext;

 42: int main(int argc, char **argv)
 43: {
 44:   KSP         ksp;
 45:   DM          da;
 46:   UserContext user;
 47:   const char *bcTypes[2] = {"dirichlet", "neumann"};
 48:   PetscInt    bc;
 49:   Vec         b, x;
 50:   PetscBool   testsolver = PETSC_FALSE;

 53:   PetscInitialize(&argc, &argv, (char *)0, help);
 54:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 55:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 3, 3, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da);
 56:   DMSetFromOptions(da);
 57:   DMSetUp(da);
 58:   DMDASetUniformCoordinates(da, 0, 1, 0, 1, 0, 0);
 59:   DMDASetFieldName(da, 0, "Pressure");

 61:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
 62:   user.rho = 1.0;
 63:   PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, NULL);
 64:   user.nu = 0.1;
 65:   PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, NULL);
 66:   bc = (PetscInt)DIRICHLET;
 67:   PetscOptionsEList("-bc_type", "Type of boundary condition", "ex29.c", bcTypes, 2, bcTypes[0], &bc, NULL);
 68:   user.bcType = (BCType)bc;
 69:   PetscOptionsBool("-testsolver", "Run solver multiple times, useful for performance studies of solver", "ex29.c", testsolver, &testsolver, NULL);
 70:   PetscOptionsEnd();

 72:   KSPSetComputeRHS(ksp, ComputeRHS, &user);
 73:   KSPSetComputeOperators(ksp, ComputeMatrix, &user);
 74:   KSPSetDM(ksp, da);
 75:   KSPSetFromOptions(ksp);
 76:   KSPSetUp(ksp);
 77:   KSPSolve(ksp, NULL, NULL);

 79:   if (testsolver) {
 80:     KSPGetSolution(ksp, &x);
 81:     KSPGetRhs(ksp, &b);
 82:     KSPSetDMActive(ksp, PETSC_FALSE);
 83:     KSPSolve(ksp, b, x);
 84:     {
 85: #if defined(PETSC_USE_LOG)
 86:       PetscLogStage stage;
 87: #endif
 88:       PetscInt i, n = 20;

 90:       PetscLogStageRegister("Solve only", &stage);
 91:       PetscLogStagePush(stage);
 92:       for (i = 0; i < n; i++) KSPSolve(ksp, b, x);
 93:       PetscLogStagePop();
 94:     }
 95:   }

 97:   DMDestroy(&da);
 98:   KSPDestroy(&ksp);
 99:   PetscFinalize();
100:   return 0;
101: }

103: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
104: {
105:   UserContext  *user = (UserContext *)ctx;
106:   PetscInt      i, j, mx, my, xm, ym, xs, ys;
107:   PetscScalar   Hx, Hy;
108:   PetscScalar **array;
109:   DM            da;

112:   KSPGetDM(ksp, &da);
113:   DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
114:   Hx = 1.0 / (PetscReal)(mx - 1);
115:   Hy = 1.0 / (PetscReal)(my - 1);
116:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
117:   DMDAVecGetArray(da, b, &array);
118:   for (j = ys; j < ys + ym; j++) {
119:     for (i = xs; i < xs + xm; i++) array[j][i] = PetscExpScalar(-((PetscReal)i * Hx) * ((PetscReal)i * Hx) / user->nu) * PetscExpScalar(-((PetscReal)j * Hy) * ((PetscReal)j * Hy) / user->nu) * Hx * Hy;
120:   }
121:   DMDAVecRestoreArray(da, b, &array);
122:   VecAssemblyBegin(b);
123:   VecAssemblyEnd(b);

125:   /* force right hand side to be consistent for singular matrix */
126:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
127:   if (user->bcType == NEUMANN) {
128:     MatNullSpace nullspace;

130:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
131:     MatNullSpaceRemove(nullspace, b);
132:     MatNullSpaceDestroy(&nullspace);
133:   }
134:   return 0;
135: }

137: PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
138: {
140:   if ((i > mx / 3.0) && (i < 2.0 * mx / 3.0) && (j > my / 3.0) && (j < 2.0 * my / 3.0)) {
141:     *rho = centerRho;
142:   } else {
143:     *rho = 1.0;
144:   }
145:   return 0;
146: }

148: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
149: {
150:   UserContext *user = (UserContext *)ctx;
151:   PetscReal    centerRho;
152:   PetscInt     i, j, mx, my, xm, ym, xs, ys;
153:   PetscScalar  v[5];
154:   PetscReal    Hx, Hy, HydHx, HxdHy, rho;
155:   MatStencil   row, col[5];
156:   DM           da;
157:   PetscBool    check_matis = PETSC_FALSE;

160:   KSPGetDM(ksp, &da);
161:   centerRho = user->rho;
162:   DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
163:   Hx    = 1.0 / (PetscReal)(mx - 1);
164:   Hy    = 1.0 / (PetscReal)(my - 1);
165:   HxdHy = Hx / Hy;
166:   HydHx = Hy / Hx;
167:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
168:   for (j = ys; j < ys + ym; j++) {
169:     for (i = xs; i < xs + xm; i++) {
170:       row.i = i;
171:       row.j = j;
172:       ComputeRho(i, j, mx, my, centerRho, &rho);
173:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
174:         if (user->bcType == DIRICHLET) {
175:           v[0] = 2.0 * rho * (HxdHy + HydHx);
176:           MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES);
177:         } else if (user->bcType == NEUMANN) {
178:           PetscInt numx = 0, numy = 0, num = 0;
179:           if (j != 0) {
180:             v[num]     = -rho * HxdHy;
181:             col[num].i = i;
182:             col[num].j = j - 1;
183:             numy++;
184:             num++;
185:           }
186:           if (i != 0) {
187:             v[num]     = -rho * HydHx;
188:             col[num].i = i - 1;
189:             col[num].j = j;
190:             numx++;
191:             num++;
192:           }
193:           if (i != mx - 1) {
194:             v[num]     = -rho * HydHx;
195:             col[num].i = i + 1;
196:             col[num].j = j;
197:             numx++;
198:             num++;
199:           }
200:           if (j != my - 1) {
201:             v[num]     = -rho * HxdHy;
202:             col[num].i = i;
203:             col[num].j = j + 1;
204:             numy++;
205:             num++;
206:           }
207:           v[num]     = numx * rho * HydHx + numy * rho * HxdHy;
208:           col[num].i = i;
209:           col[num].j = j;
210:           num++;
211:           MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES);
212:         }
213:       } else {
214:         v[0]     = -rho * HxdHy;
215:         col[0].i = i;
216:         col[0].j = j - 1;
217:         v[1]     = -rho * HydHx;
218:         col[1].i = i - 1;
219:         col[1].j = j;
220:         v[2]     = 2.0 * rho * (HxdHy + HydHx);
221:         col[2].i = i;
222:         col[2].j = j;
223:         v[3]     = -rho * HydHx;
224:         col[3].i = i + 1;
225:         col[3].j = j;
226:         v[4]     = -rho * HxdHy;
227:         col[4].i = i;
228:         col[4].j = j + 1;
229:         MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES);
230:       }
231:     }
232:   }
233:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
234:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
235:   MatViewFromOptions(jac, NULL, "-view_mat");
236:   PetscOptionsGetBool(NULL, NULL, "-check_matis", &check_matis, NULL);
237:   if (check_matis) {
238:     void (*f)(void);
239:     Mat       J2;
240:     MatType   jtype;
241:     PetscReal nrm;

243:     MatGetType(jac, &jtype);
244:     MatConvert(jac, MATIS, MAT_INITIAL_MATRIX, &J2);
245:     MatViewFromOptions(J2, NULL, "-view_conv");
246:     MatConvert(J2, jtype, MAT_INPLACE_MATRIX, &J2);
247:     MatGetOperation(jac, MATOP_VIEW, &f);
248:     MatSetOperation(J2, MATOP_VIEW, f);
249:     MatSetDM(J2, da);
250:     MatViewFromOptions(J2, NULL, "-view_conv_assembled");
251:     MatAXPY(J2, -1., jac, DIFFERENT_NONZERO_PATTERN);
252:     MatNorm(J2, NORM_FROBENIUS, &nrm);
253:     PetscPrintf(PETSC_COMM_WORLD, "Error MATIS %g\n", (double)nrm);
254:     MatViewFromOptions(J2, NULL, "-view_conv_err");
255:     MatDestroy(&J2);
256:   }
257:   if (user->bcType == NEUMANN) {
258:     MatNullSpace nullspace;

260:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
261:     MatSetNullSpace(J, nullspace);
262:     MatNullSpaceDestroy(&nullspace);
263:   }
264:   return 0;
265: }

267: /*TEST

269:    test:
270:       args: -pc_type mg -pc_mg_type full -ksp_type fgmres -ksp_monitor_short -da_refine 8 -ksp_rtol 1.e-3

272:    test:
273:       suffix: 2
274:       args: -bc_type neumann -pc_type mg -pc_mg_type full -ksp_type fgmres -ksp_monitor_short -da_refine 8 -mg_coarse_pc_factor_shift_type nonzero
275:       requires: !single

277:    test:
278:       suffix: telescope
279:       nsize: 4
280:       args: -ksp_monitor_short -da_grid_x 257 -da_grid_y 257 -pc_type mg -pc_mg_galerkin pmat -pc_mg_levels 4 -ksp_type richardson -mg_levels_ksp_type chebyshev -mg_levels_pc_type jacobi -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_ignore_kspcomputeoperators -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_galerkin pmat -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_levels_ksp_type chebyshev -mg_coarse_telescope_mg_levels_pc_type jacobi -mg_coarse_pc_telescope_reduction_factor 4

282:    test:
283:       suffix: 3
284:       args: -ksp_view -da_refine 2 -pc_type mg -pc_mg_distinct_smoothup -mg_levels_up_pc_type jacobi

286:    test:
287:       suffix: 4
288:       args: -ksp_view -da_refine 2 -pc_type mg -pc_mg_distinct_smoothup -mg_levels_up_ksp_max_it 3 -mg_levels_ksp_max_it 4

290:    test:
291:       suffix: 5
292:       nsize: 2
293:       requires: hypre !complex
294:       args: -pc_type mg  -da_refine 2 -ksp_monitor  -matptap_via hypre -pc_mg_galerkin both

296:    test:
297:       suffix: 6
298:       args: -pc_type svd -pc_svd_monitor ::all

300: TEST*/