Actual source code: ex36.cxx


  2: /*
  3: Inhomogeneous Laplacian in 3-D. Modeled by the partial differential equation

  5:    -div \rho grad u + \alpha u = f,  0 < x,y,z < 1,

  7: Problem 1: (Default)

  9:   Use the following exact solution with Dirichlet boundary condition

 11:     u = sin(pi*x)*sin(pi*y)*sin(pi*z)

 13:   with Dirichlet boundary conditions

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

 17:   and \rho = 1.0, \alpha = 10.0 uniformly in the domain.

 19:   Use an appropriate forcing function to measure convergence.

 21: Usage:

 23:   Measure convergence rate with uniform refinement with the options: "-problem 1 -error".

 25:     mpiexec -n $NP ./ex36 -problem 1 -error -n 4 -levels 5 -pc_type gamg
 26:     mpiexec -n $NP ./ex36 -problem 1 -error -n 8 -levels 4 -pc_type gamg
 27:     mpiexec -n $NP ./ex36 -problem 1 -error -n 16 -levels 3 -pc_type mg
 28:     mpiexec -n $NP ./ex36 -problem 1 -error -n 32 -levels 2 -pc_type mg
 29:     mpiexec -n $NP ./ex36 -problem 1 -error -n 64 -levels 1 -mg
 30:     mpiexec -n $NP ./ex36 -problem 1 -error -n 128 -levels 0 -mg

 32:   Or with an external mesh file representing [0, 1]^3,

 34:     mpiexec -n $NP ./ex36 -problem 1 -file ./external_mesh.h5m -levels 2 -error -pc_type mg

 36: Problem 2:

 38:   Use the forcing function

 40:      f = e^{-((x-xr)^2+(y-yr)^2+(z-zr)^2)/\nu}

 42:   with Dirichlet boundary conditions

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

 46:   or pure Neumman boundary conditions

 48: Usage:

 50:   Run with different values of \rho and \nu (problem 1) to control diffusion and gaussian source spread. This uses the internal mesh generator implemented in DMMoab.

 52:     mpiexec -n $NP ./ex36 -problem 2 -n 20 -nu 0.02 -rho 0.01
 53:     mpiexec -n $NP ./ex36 -problem 2 -n 40 -nu 0.01 -rho 0.005 -io -ksp_monitor
 54:     mpiexec -n $NP ./ex36 -problem 2 -n 80 -nu 0.01 -rho 0.005 -io -ksp_monitor -pc_type gamg
 55:     mpiexec -n $NP ./ex36 -problem 2 -n 160 -bc neumann -nu 0.005 -rho 0.01 -io
 56:     mpiexec -n $NP ./ex36 -problem 2 -n 320 -bc neumann -nu 0.001 -rho 1 -io

 58:   Or with an external mesh file representing [0, 1]^3,

 60:     mpiexec -n $NP ./ex36 -problem 2 -file ./external_mesh.h5m -levels 1 -pc_type gamg

 62: */

 64: static char help[] = "\
 65:                       Solves a three dimensional inhomogeneous Laplacian equation with a Gaussian source.\n \
 66:                       Usage: ./ex36 -bc dirichlet -nu .01 -n 10\n";

 68: /* PETSc includes */
 69: #include <petscksp.h>
 70: #include <petscdmmoab.h>

 72: typedef enum {
 73:   DIRICHLET,
 74:   NEUMANN
 75: } BCType;

 77: typedef struct {
 78:   PetscInt  problem, dim, n, nlevels;
 79:   PetscReal rho;
 80:   PetscReal bounds[6];
 81:   PetscReal xyzref[3];
 82:   PetscReal nu;
 83:   BCType    bcType;
 84:   char      filename[PETSC_MAX_PATH_LEN];
 85:   PetscBool use_extfile, io, error, usetet, usemg;

 87:   /* Discretization parameters */
 88:   int VPERE;
 89: } UserContext;

 91: static PetscErrorCode ComputeMatrix_MOAB(KSP, Mat, Mat, void *);
 92: static PetscErrorCode ComputeRHS_MOAB(KSP, Vec, void *);
 93: static PetscErrorCode ComputeDiscreteL2Error(KSP ksp, Vec err, UserContext *user);
 94: static PetscErrorCode InitializeOptions(UserContext *user);

 96: int main(int argc, char **argv)
 97: {
 98:   const char *fields[1] = {"T-Variable"};
 99:   DM          dm, dmref, *dmhierarchy;
100:   UserContext user;
101:   PetscInt    k;
102:   KSP         ksp;
103:   PC          pc;
104:   Vec         errv;
105:   Mat         R;
106:   Vec         b, x;

109:   PetscInitialize(&argc, &argv, (char *)0, help);

111:   InitializeOptions(&user);

113:   /* Create the DM object from either a mesh file or from in-memory structured grid */
114:   if (user.use_extfile) {
115:     DMMoabLoadFromFile(PETSC_COMM_WORLD, user.dim, 1, user.filename, "", &dm);
116:   } else {
117:     DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.usetet, NULL, user.n, 1, &dm);
118:   }
119:   DMSetFromOptions(dm);
120:   DMMoabSetFieldNames(dm, 1, fields);

122:   /* SetUp the data structures for DMMOAB */
123:   DMSetUp(dm);

125:   DMSetApplicationContext(dm, &user);

127:   KSPCreate(PETSC_COMM_WORLD, &ksp);
128:   KSPSetComputeRHS(ksp, ComputeRHS_MOAB, &user);
129:   KSPSetComputeOperators(ksp, ComputeMatrix_MOAB, &user);

131:   if (user.nlevels) {
132:     KSPGetPC(ksp, &pc);
133:     PetscMalloc(sizeof(DM) * (user.nlevels + 1), &dmhierarchy);
134:     for (k = 0; k <= user.nlevels; k++) dmhierarchy[k] = NULL;

136:     PetscPrintf(PETSC_COMM_WORLD, "Number of mesh hierarchy levels: %d\n", user.nlevels);
137:     DMMoabGenerateHierarchy(dm, user.nlevels, PETSC_NULL);

139:     // coarsest grid = 0
140:     // finest grid = nlevels
141:     dmhierarchy[0]         = dm;
142:     PetscBool usehierarchy = PETSC_FALSE;
143:     if (usehierarchy) {
144:       DMRefineHierarchy(dm, user.nlevels, &dmhierarchy[1]);
145:     } else {
146:       for (k = 1; k <= user.nlevels; k++) DMRefine(dmhierarchy[k - 1], MPI_COMM_NULL, &dmhierarchy[k]);
147:     }
148:     dmref = dmhierarchy[user.nlevels];
149:     PetscObjectReference((PetscObject)dmref);

151:     if (user.usemg) {
152:       PCSetType(pc, PCMG);
153:       PCMGSetLevels(pc, user.nlevels + 1, NULL);
154:       PCMGSetType(pc, PC_MG_MULTIPLICATIVE);
155:       PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH);
156:       PCMGSetCycleType(pc, PC_MG_CYCLE_V);
157:       PCMGSetNumberSmooth(pc, 2);

159:       for (k = 1; k <= user.nlevels; k++) {
160:         DMCreateInterpolation(dmhierarchy[k - 1], dmhierarchy[k], &R, NULL);
161:         PCMGSetInterpolation(pc, k, R);
162:         MatDestroy(&R);
163:       }
164:     }

166:     for (k = 1; k <= user.nlevels; k++) DMDestroy(&dmhierarchy[k]);
167:     PetscFree(dmhierarchy);
168:   } else {
169:     dmref = dm;
170:     PetscObjectReference((PetscObject)dm);
171:   }

173:   KSPSetDM(ksp, dmref);
174:   KSPSetFromOptions(ksp);

176:   /* Perform the actual solve */
177:   KSPSolve(ksp, NULL, NULL);
178:   KSPGetSolution(ksp, &x);
179:   KSPGetRhs(ksp, &b);

181:   if (user.error) {
182:     VecDuplicate(b, &errv);
183:     ComputeDiscreteL2Error(ksp, errv, &user);
184:     VecDestroy(&errv);
185:   }

187:   if (user.io) {
188:     /* Write out the solution along with the mesh */
189:     DMMoabSetGlobalFieldVector(dmref, x);
190: #ifdef MOAB_HAVE_HDF5
191:     DMMoabOutput(dmref, "ex36.h5m", "");
192: #else
193:     /* MOAB does not support true parallel writers that aren't HDF5 based
194:        And so if you are using VTK as the output format in parallel,
195:        the data could be jumbled due to the order in which the processors
196:        write out their parts of the mesh and solution tags
197:     */
198:     DMMoabOutput(dmref, "ex36.vtk", "");
199: #endif
200:   }

202:   /* Cleanup objects */
203:   KSPDestroy(&ksp);
204:   DMDestroy(&dmref);
205:   DMDestroy(&dm);
206:   PetscFinalize();
207:   return 0;
208: }

210: PetscReal ComputeDiffusionCoefficient(PetscReal coords[3], UserContext *user)
211: {
212:   if (user->problem == 2) {
213:     if ((coords[0] > 1.0 / 3.0) && (coords[0] < 2.0 / 3.0) && (coords[1] > 1.0 / 3.0) && (coords[1] < 2.0 / 3.0) && (coords[2] > 1.0 / 3.0) && (coords[2] < 2.0 / 3.0)) return user->rho;
214:     else return 1.0;
215:   } else return 1.0; /* problem = 1 */
216: }

218: PetscReal ComputeReactionCoefficient(PetscReal coords[3], UserContext *user)
219: {
220:   if (user->problem == 2) {
221:     if ((coords[0] > 1.0 / 3.0) && (coords[0] < 2.0 / 3.0) && (coords[1] > 1.0 / 3.0) && (coords[1] < 2.0 / 3.0) && (coords[2] > 1.0 / 3.0) && (coords[2] < 2.0 / 3.0)) return 10.0;
222:     else return 0.0;
223:   } else return 5.0; /* problem = 1 */
224: }

226: double ExactSolution(PetscReal coords[3], UserContext *user)
227: {
228:   if (user->problem == 2) {
229:     const PetscScalar xx = (coords[0] - user->xyzref[0]) * (coords[0] - user->xyzref[0]);
230:     const PetscScalar yy = (coords[1] - user->xyzref[1]) * (coords[1] - user->xyzref[1]);
231:     const PetscScalar zz = (coords[2] - user->xyzref[2]) * (coords[2] - user->xyzref[2]);
232:     return PetscExpScalar(-(xx + yy + zz) / user->nu);
233:   } else return sin(PETSC_PI * coords[0]) * sin(PETSC_PI * coords[1]) * sin(PETSC_PI * coords[2]);
234: }

236: PetscReal exact_solution(PetscReal x, PetscReal y, PetscReal z)
237: {
238:   PetscReal coords[3] = {x, y, z};
239:   return ExactSolution(coords, 0);
240: }

242: double ForcingFunction(PetscReal coords[3], UserContext *user)
243: {
244:   const PetscReal exact = ExactSolution(coords, user);
245:   if (user->problem == 2) {
246:     const PetscReal duxyz = ((coords[0] - user->xyzref[0]) + (coords[1] - user->xyzref[1]) + (coords[2] - user->xyzref[2]));
247:     return (4.0 / user->nu * duxyz * duxyz - 6.0) * exact / user->nu;
248:   } else {
249:     const PetscReal reac = ComputeReactionCoefficient(coords, user);
250:     return (3.0 * PETSC_PI * PETSC_PI + reac) * exact;
251:   }
252: }

254: PetscErrorCode ComputeRHS_MOAB(KSP ksp, Vec b, void *ptr)
255: {
256:   UserContext              *user = (UserContext *)ptr;
257:   DM                        dm;
258:   PetscInt                  dof_indices[8], nc, npoints;
259:   PetscBool                 dbdry[8];
260:   PetscReal                 vpos[8 * 3];
261:   PetscInt                  i, q, nconn;
262:   const moab::EntityHandle *connect;
263:   const moab::Range        *elocal;
264:   moab::Interface          *mbImpl;
265:   PetscScalar               localv[8];
266:   PetscReal                *phi, *phypts, *jxw;
267:   PetscBool                 elem_on_boundary;
268:   PetscQuadrature           quadratureObj;

270:   KSPGetDM(ksp, &dm);

272:   /* reset the RHS */
273:   VecSet(b, 0.0);

275:   DMMoabFEMCreateQuadratureDefault(user->dim, user->VPERE, &quadratureObj);
276:   PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL);
277:   PetscMalloc3(user->VPERE * npoints, &phi, npoints * 3, &phypts, npoints, &jxw);

279:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
280:   DMMoabGetInterface(dm, &mbImpl);
281:   DMMoabGetLocalElements(dm, &elocal);

283:   /* loop over local elements */
284:   for (moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
285:     const moab::EntityHandle ehandle = *iter;

287:     /* Get connectivity information: */
288:     DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect);

291:     /* get the coordinates of the element vertices */
292:     DMMoabGetVertexCoordinates(dm, nconn, connect, vpos);

294:     /* get the local DoF numbers to appropriately set the element contribution in the operator */
295:     DMMoabGetFieldDofsLocal(dm, nconn, connect, 0, dof_indices);

297:     /* compute the quadrature points transformed to the physical space and then
298:        compute the basis functions to compute local operators */
299:     DMMoabFEMComputeBasis(user->dim, nconn, vpos, quadratureObj, phypts, jxw, phi, NULL);

301:     PetscArrayzero(localv, nconn);
302:     /* Compute function over the locally owned part of the grid */
303:     for (q = 0; q < npoints; ++q) {
304:       const double   ff     = ForcingFunction(&phypts[3 * q], user);
305:       const PetscInt offset = q * nconn;

307:       for (i = 0; i < nconn; ++i) localv[i] += jxw[q] * phi[offset + i] * ff;
308:     }

310:     /* check if element is on the boundary */
311:     DMMoabIsEntityOnBoundary(dm, ehandle, &elem_on_boundary);

313:     /* apply dirichlet boundary conditions */
314:     if (elem_on_boundary && user->bcType == DIRICHLET) {
315:       /* get the list of nodes on boundary so that we can enforce dirichlet conditions strongly */
316:       DMMoabCheckBoundaryVertices(dm, nconn, connect, dbdry);

318:       for (i = 0; i < nconn; ++i) {
319:         if (dbdry[i]) { /* dirichlet node */
320:           /* think about strongly imposing dirichlet */
321:           localv[i] = ForcingFunction(&vpos[3 * i], user);
322:         }
323:       }
324:     }

326:     /* set the values directly into appropriate locations. Can alternately use VecSetValues */
327:     VecSetValuesLocal(b, nconn, dof_indices, localv, ADD_VALUES);
328:   }

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

335:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
336:     MatNullSpaceRemove(nullspace, b);
337:     MatNullSpaceDestroy(&nullspace);
338:   }

340:   /* Restore vectors */
341:   VecAssemblyBegin(b);
342:   VecAssemblyEnd(b);
343:   PetscFree3(phi, phypts, jxw);
344:   PetscQuadratureDestroy(&quadratureObj);
345:   return 0;
346: }

348: PetscErrorCode ComputeMatrix_MOAB(KSP ksp, Mat J, Mat jac, void *ctx)
349: {
350:   UserContext              *user = (UserContext *)ctx;
351:   DM                        dm;
352:   PetscInt                  i, j, q, nconn, nglobale, nglobalv, nc, npoints, hlevel;
353:   PetscInt                  dof_indices[8];
354:   PetscReal                 vpos[8 * 3], rho, alpha;
355:   PetscBool                 dbdry[8];
356:   const moab::EntityHandle *connect;
357:   const moab::Range        *elocal;
358:   moab::Interface          *mbImpl;
359:   PetscBool                 elem_on_boundary;
360:   PetscScalar               array[8 * 8];
361:   PetscReal                *phi, *dphi[3], *phypts, *jxw;
362:   PetscQuadrature           quadratureObj;

365:   KSPGetDM(ksp, &dm);

367:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
368:   DMMoabGetInterface(dm, &mbImpl);
369:   DMMoabGetLocalElements(dm, &elocal);
370:   DMMoabGetSize(dm, &nglobale, &nglobalv);
371:   DMMoabGetHierarchyLevel(dm, &hlevel);
372:   PetscPrintf(PETSC_COMM_WORLD, "ComputeMatrix: Level = %d, N(elements) = %d, N(vertices) = %d \n", hlevel, nglobale, nglobalv);

374:   DMMoabFEMCreateQuadratureDefault(user->dim, user->VPERE, &quadratureObj);
375:   PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL);
376:   PetscMalloc6(user->VPERE * npoints, &phi, user->VPERE * npoints, &dphi[0], user->VPERE * npoints, &dphi[1], user->VPERE * npoints, &dphi[2], npoints * 3, &phypts, npoints, &jxw);

378:   /* loop over local elements */
379:   for (moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
380:     const moab::EntityHandle ehandle = *iter;

382:     /* Get connectivity information: */
383:     DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect);

386:     /* get the coordinates of the element vertices */
387:     DMMoabGetVertexCoordinates(dm, nconn, connect, vpos);

389:     /* get the local DoF numbers to appropriately set the element contribution in the operator */
390:     DMMoabGetFieldDofsLocal(dm, nconn, connect, 0, dof_indices);

392:     /* compute the quadrature points transformed to the physical space and
393:        compute the basis functions and the derivatives wrt x, y and z directions */
394:     DMMoabFEMComputeBasis(user->dim, nconn, vpos, quadratureObj, phypts, jxw, phi, dphi);

396:     PetscArrayzero(array, nconn * nconn);

398:     /* Compute function over the locally owned part of the grid */
399:     for (q = 0; q < npoints; ++q) {
400:       /* compute the inhomogeneous diffusion coefficient at the quadrature point
401:           -- for large spatial variations, embed this property evaluation inside quadrature loop */
402:       rho   = ComputeDiffusionCoefficient(&phypts[q * 3], user);
403:       alpha = ComputeReactionCoefficient(&phypts[q * 3], user);

405:       const PetscInt offset = q * nconn;

407:       for (i = 0; i < nconn; ++i) {
408:         for (j = 0; j < nconn; ++j) {
409:           array[i * nconn + j] += jxw[q] * (rho * (dphi[0][offset + i] * dphi[0][offset + j] + dphi[1][offset + i] * dphi[1][offset + j] + dphi[2][offset + i] * dphi[2][offset + j]) + alpha * (phi[offset + i] * phi[offset + j]));
410:         }
411:       }
412:     }

414:     /* check if element is on the boundary */
415:     DMMoabIsEntityOnBoundary(dm, ehandle, &elem_on_boundary);

417:     /* apply dirichlet boundary conditions */
418:     if (elem_on_boundary && user->bcType == DIRICHLET) {
419:       /* get the list of nodes on boundary so that we can enforce dirichlet conditions strongly */
420:       DMMoabCheckBoundaryVertices(dm, nconn, connect, dbdry);

422:       for (i = 0; i < nconn; ++i) {
423:         if (dbdry[i]) { /* dirichlet node */
424:           /* think about strongly imposing dirichlet */
425:           for (j = 0; j < nconn; ++j) {
426:             /* TODO: symmetrize the system - need the RHS */
427:             array[i * nconn + j] = 0.0;
428:           }
429:           array[i * nconn + i] = 1.0;
430:         }
431:       }
432:     }

434:     /* set the values directly into appropriate locations. Can alternately use VecSetValues */
435:     MatSetValuesLocal(jac, nconn, dof_indices, nconn, dof_indices, array, ADD_VALUES);
436:   }

438:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
439:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

441:   if (user->bcType == NEUMANN && false) {
442:     MatNullSpace nullspace;

444:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
445:     MatSetNullSpace(jac, nullspace);
446:     MatNullSpaceDestroy(&nullspace);
447:   }
448:   PetscFree6(phi, dphi[0], dphi[1], dphi[2], phypts, jxw);
449:   PetscQuadratureDestroy(&quadratureObj);
450:   return 0;
451: }

453: PetscErrorCode ComputeDiscreteL2Error(KSP ksp, Vec err, UserContext *user)
454: {
455:   DM                 dm;
456:   Vec                sol;
457:   PetscScalar        vpos[3];
458:   const PetscScalar *x;
459:   PetscScalar       *e;
460:   PetscReal          l2err = 0.0, linferr = 0.0, global_l2, global_linf;
461:   PetscInt           dof_index, N;
462:   const moab::Range *ownedvtx;

464:   KSPGetDM(ksp, &dm);

466:   /* get the solution vector */
467:   KSPGetSolution(ksp, &sol);

469:   /* Get the internal reference to the vector arrays */
470:   VecGetArrayRead(sol, &x);
471:   VecGetSize(sol, &N);
472:   if (err) {
473:     /* reset the error vector */
474:     VecSet(err, 0.0);
475:     /* get array reference */
476:     VecGetArray(err, &e);
477:   }

479:   DMMoabGetLocalVertices(dm, &ownedvtx, NULL);

481:   /* Compute function over the locally owned part of the grid */
482:   for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
483:     const moab::EntityHandle vhandle = *iter;
484:     DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index);

486:     /* compute the mid-point of the element and use a 1-point lumped quadrature */
487:     DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos);

489:     /* compute the discrete L2 error against the exact solution */
490:     const PetscScalar lerr = (ExactSolution(vpos, user) - x[dof_index]);
491:     l2err += lerr * lerr;
492:     if (linferr < fabs(lerr)) linferr = fabs(lerr);

494:     if (err) {
495:       /* set the discrete L2 error against the exact solution */
496:       e[dof_index] = lerr;
497:     }
498:   }

500:   MPI_Allreduce(&l2err, &global_l2, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
501:   MPI_Allreduce(&linferr, &global_linf, 1, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
502:   PetscPrintf(PETSC_COMM_WORLD, "Computed Errors: L_2 = %f, L_inf = %f\n", sqrt(global_l2 / N), global_linf);

504:   /* Restore vectors */
505:   VecRestoreArrayRead(sol, &x);
506:   if (err) VecRestoreArray(err, &e);
507:   return 0;
508: }

510: PetscErrorCode InitializeOptions(UserContext *user)
511: {
512:   const char *bcTypes[2] = {"dirichlet", "neumann"};
513:   PetscInt    bc;

515:   /* set default parameters */
516:   user->dim       = 3; /* 3-Dimensional problem */
517:   user->problem   = 1;
518:   user->n         = 2;
519:   user->nlevels   = 2;
520:   user->rho       = 0.1;
521:   user->bounds[0] = user->bounds[2] = user->bounds[4] = 0.0;
522:   user->bounds[1] = user->bounds[3] = user->bounds[5] = 1.0;
523:   user->xyzref[0]                                     = user->bounds[1] / 2;
524:   user->xyzref[1]                                     = user->bounds[3] / 2;
525:   user->xyzref[2]                                     = user->bounds[5] / 2;
526:   user->nu                                            = 0.05;
527:   user->usemg                                         = PETSC_FALSE;
528:   user->io                                            = PETSC_FALSE;
529:   user->usetet                                        = PETSC_FALSE;
530:   user->error                                         = PETSC_FALSE;
531:   bc                                                  = (PetscInt)DIRICHLET;

533:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "ex36.cxx");
534:   PetscOptionsInt("-problem", "The type of problem being solved (controls forcing function)", "ex36.cxx", user->problem, &user->problem, NULL);
535:   PetscOptionsInt("-n", "The elements in each direction", "ex36.cxx", user->n, &user->n, NULL);
536:   PetscOptionsInt("-levels", "Number of levels in the multigrid hierarchy", "ex36.cxx", user->nlevels, &user->nlevels, NULL);
537:   PetscOptionsReal("-rho", "The conductivity", "ex36.cxx", user->rho, &user->rho, NULL);
538:   PetscOptionsReal("-x", "The domain size in x-direction", "ex36.cxx", user->bounds[1], &user->bounds[1], NULL);
539:   PetscOptionsReal("-y", "The domain size in y-direction", "ex36.cxx", user->bounds[3], &user->bounds[3], NULL);
540:   PetscOptionsReal("-z", "The domain size in y-direction", "ex36.cxx", user->bounds[5], &user->bounds[5], NULL);
541:   PetscOptionsReal("-xref", "The x-coordinate of Gaussian center (for -problem 1)", "ex36.cxx", user->xyzref[0], &user->xyzref[0], NULL);
542:   PetscOptionsReal("-yref", "The y-coordinate of Gaussian center (for -problem 1)", "ex36.cxx", user->xyzref[1], &user->xyzref[1], NULL);
543:   PetscOptionsReal("-zref", "The y-coordinate of Gaussian center (for -problem 1)", "ex36.cxx", user->xyzref[2], &user->xyzref[2], NULL);
544:   PetscOptionsReal("-nu", "The width of the Gaussian source (for -problem 1)", "ex36.cxx", user->nu, &user->nu, NULL);
545:   PetscOptionsBool("-mg", "Use multigrid preconditioner", "ex36.cxx", user->usemg, &user->usemg, NULL);
546:   PetscOptionsBool("-io", "Write out the solution and mesh data", "ex36.cxx", user->io, &user->io, NULL);
547:   PetscOptionsBool("-tet", "Use tetrahedra to discretize the domain", "ex36.cxx", user->usetet, &user->usetet, NULL);
548:   PetscOptionsBool("-error", "Compute the discrete L_2 and L_inf errors of the solution", "ex36.cxx", user->error, &user->error, NULL);
549:   PetscOptionsEList("-bc", "Type of boundary condition", "ex36.cxx", bcTypes, 2, bcTypes[0], &bc, NULL);
550:   PetscOptionsString("-file", "The mesh file for the problem", "ex36.cxx", "", user->filename, sizeof(user->filename), &user->use_extfile);
551:   PetscOptionsEnd();

553:   if (user->problem < 1 || user->problem > 2) user->problem = 1;
554:   user->bcType = (BCType)bc;
555:   user->VPERE  = (user->usetet ? 4 : 8);
556:   return 0;
557: }

559: /*TEST

561:    build:
562:       requires: moab

564:    test:
565:       args: -levels 1 -nu .01 -n 4 -mg -ksp_converged_reason

567:    test:
568:       suffix: 2
569:       nsize: 2
570:       requires: hdf5
571:       args: -levels 2 -nu .01 -n 2 -mg -ksp_converged_reason

573: TEST*/