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