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;

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

111:   PetscCall(InitializeOptions(&user));

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

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

125:   PetscCall(DMSetApplicationContext(dm, &user));

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

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

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

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

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

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

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

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

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

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

187:   if (user.io) {
188:     /* Write out the solution along with the mesh */
189:     PetscCall(DMMoabSetGlobalFieldVector(dmref, x));
190: #ifdef MOAB_HAVE_HDF5
191:     PetscCall(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:     PetscCall(DMMoabOutput(dmref, "ex36.vtk", ""));
199: #endif
200:   }

202:   /* Cleanup objects */
203:   PetscCall(KSPDestroy(&ksp));
204:   PetscCall(DMDestroy(&dmref));
205:   PetscCall(DMDestroy(&dm));
206:   PetscCall(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:   PetscFunctionBegin;
271:   PetscCall(KSPGetDM(ksp, &dm));

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

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

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

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

288:     /* Get connectivity information: */
289:     PetscCall(DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect));
290:     PetscCheck(nconn == 4 || nconn == 8, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only HEX8/TET4 element bases are supported in the current example. n(Connectivity)=%" PetscInt_FMT ".", nconn);

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

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

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

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

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

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

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

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

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

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

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

341:   /* Restore vectors */
342:   PetscCall(VecAssemblyBegin(b));
343:   PetscCall(VecAssemblyEnd(b));
344:   PetscCall(PetscFree3(phi, phypts, jxw));
345:   PetscCall(PetscQuadratureDestroy(&quadratureObj));
346:   PetscFunctionReturn(PETSC_SUCCESS);
347: }

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

365:   PetscFunctionBeginUser;
366:   PetscCall(KSPGetDM(ksp, &dm));

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

375:   PetscCall(DMMoabFEMCreateQuadratureDefault(user->dim, user->VPERE, &quadratureObj));
376:   PetscCall(PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL));
377:   PetscCall(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));

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

383:     /* Get connectivity information: */
384:     PetscCall(DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect));
385:     PetscCheck(nconn == 4 || nconn == 8, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only HEX8/TET4 element bases are supported in the current example. n(Connectivity)=%" PetscInt_FMT ".", nconn);

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

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

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

397:     PetscCall(PetscArrayzero(array, nconn * nconn));

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

406:       const PetscInt offset = q * nconn;

408:       for (i = 0; i < nconn; ++i) {
409:         for (j = 0; j < nconn; ++j) {
410:           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]));
411:         }
412:       }
413:     }

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

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

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

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

439:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
440:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

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

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

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

465:   PetscFunctionBegin;
466:   PetscCall(KSPGetDM(ksp, &dm));

468:   /* get the solution vector */
469:   PetscCall(KSPGetSolution(ksp, &sol));

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

481:   PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));

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

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

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

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

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

506:   /* Restore vectors */
507:   PetscCall(VecRestoreArrayRead(sol, &x));
508:   if (err) PetscCall(VecRestoreArray(err, &e));
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }

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

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

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

556:   if (user->problem < 1 || user->problem > 2) user->problem = 1;
557:   user->bcType = (BCType)bc;
558:   user->VPERE  = (user->usetet ? 4 : 8);
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: /*TEST

564:    build:
565:       requires: moab

567:    test:
568:       args: -levels 1 -nu .01 -n 4 -mg -ksp_converged_reason

570:    test:
571:       suffix: 2
572:       nsize: 2
573:       requires: hdf5
574:       args: -levels 2 -nu .01 -n 2 -mg -ksp_converged_reason

576: TEST*/