Actual source code: ex36.cxx
1: /*
2: Inhomogeneous Laplacian in 3-D. Modeled by the partial differential equation
4: -div \rho grad u + \alpha u = f, 0 < x,y,z < 1,
6: Problem 1: (Default)
8: Use the following exact solution with Dirichlet boundary condition
10: u = sin(pi*x)*sin(pi*y)*sin(pi*z)
12: with Dirichlet boundary conditions
14: u = f(x,y,z) for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
16: and \rho = 1.0, \alpha = 10.0 uniformly in the domain.
18: Use an appropriate forcing function to measure convergence.
20: Usage:
22: Measure convergence rate with uniform refinement with the options: "-problem 1 -error".
24: mpiexec -n $NP ./ex36 -problem 1 -error -n 4 -levels 5 -pc_type gamg
25: mpiexec -n $NP ./ex36 -problem 1 -error -n 8 -levels 4 -pc_type gamg
26: mpiexec -n $NP ./ex36 -problem 1 -error -n 16 -levels 3 -pc_type mg
27: mpiexec -n $NP ./ex36 -problem 1 -error -n 32 -levels 2 -pc_type mg
28: mpiexec -n $NP ./ex36 -problem 1 -error -n 64 -levels 1 -mg
29: mpiexec -n $NP ./ex36 -problem 1 -error -n 128 -levels 0 -mg
31: Or with an external mesh file representing [0, 1]^3,
33: mpiexec -n $NP ./ex36 -problem 1 -file ./external_mesh.h5m -levels 2 -error -pc_type mg
35: Problem 2:
37: Use the forcing function
39: f = e^{-((x-xr)^2+(y-yr)^2+(z-zr)^2)/\nu}
41: with Dirichlet boundary conditions
43: u = f(x,y,z) for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
45: or pure Neumman boundary conditions
47: Usage:
49: 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.
51: mpiexec -n $NP ./ex36 -problem 2 -n 20 -nu 0.02 -rho 0.01
52: mpiexec -n $NP ./ex36 -problem 2 -n 40 -nu 0.01 -rho 0.005 -io -ksp_monitor
53: mpiexec -n $NP ./ex36 -problem 2 -n 80 -nu 0.01 -rho 0.005 -io -ksp_monitor -pc_type gamg
54: mpiexec -n $NP ./ex36 -problem 2 -n 160 -bc neumann -nu 0.005 -rho 0.01 -io
55: mpiexec -n $NP ./ex36 -problem 2 -n 320 -bc neumann -nu 0.001 -rho 1 -io
57: Or with an external mesh file representing [0, 1]^3,
59: mpiexec -n $NP ./ex36 -problem 2 -file ./external_mesh.h5m -levels 1 -pc_type gamg
61: */
63: static char help[] = "\
64: Solves a three dimensional inhomogeneous Laplacian equation with a Gaussian source.\n \
65: Usage: ./ex36 -bc dirichlet -nu .01 -n 10\n";
67: /* PETSc includes */
68: #include <petscksp.h>
69: #include <petscdmmoab.h>
71: typedef enum {
72: DIRICHLET,
73: NEUMANN
74: } BCType;
76: typedef struct {
77: PetscInt problem, dim, n, nlevels;
78: PetscReal rho;
79: PetscReal bounds[6];
80: PetscReal xyzref[3];
81: PetscReal nu;
82: BCType bcType;
83: char filename[PETSC_MAX_PATH_LEN];
84: PetscBool use_extfile, io, error, usetet, usemg;
86: /* Discretization parameters */
87: int VPERE;
88: } UserContext;
90: static PetscErrorCode ComputeMatrix_MOAB(KSP, Mat, Mat, void *);
91: static PetscErrorCode ComputeRHS_MOAB(KSP, Vec, void *);
92: static PetscErrorCode ComputeDiscreteL2Error(KSP ksp, Vec err, UserContext *user);
93: static PetscErrorCode InitializeOptions(UserContext *user);
95: int main(int argc, char **argv)
96: {
97: const char *fields[1] = {"T-Variable"};
98: DM dm, dmref, *dmhierarchy;
99: UserContext user;
100: PetscInt k;
101: KSP ksp;
102: PC pc;
103: Vec errv;
104: Mat R;
105: Vec b, x;
107: PetscFunctionBeginUser;
108: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
110: PetscCall(InitializeOptions(&user));
112: /* Create the DM object from either a mesh file or from in-memory structured grid */
113: if (user.use_extfile) {
114: PetscCall(DMMoabLoadFromFile(PETSC_COMM_WORLD, user.dim, 1, user.filename, "", &dm));
115: } else {
116: PetscCall(DMMoabCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.usetet, NULL, user.n, 1, &dm));
117: }
118: PetscCall(DMSetFromOptions(dm));
119: PetscCall(DMMoabSetFieldNames(dm, 1, fields));
121: /* SetUp the data structures for DMMOAB */
122: PetscCall(DMSetUp(dm));
124: PetscCall(DMSetApplicationContext(dm, &user));
126: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
127: PetscCall(KSPSetComputeRHS(ksp, ComputeRHS_MOAB, &user));
128: PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix_MOAB, &user));
130: if (user.nlevels) {
131: PetscCall(KSPGetPC(ksp, &pc));
132: PetscCall(PetscMalloc(sizeof(DM) * (user.nlevels + 1), &dmhierarchy));
133: for (k = 0; k <= user.nlevels; k++) dmhierarchy[k] = NULL;
135: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of mesh hierarchy levels: %d\n", user.nlevels));
136: PetscCall(DMMoabGenerateHierarchy(dm, user.nlevels, NULL));
138: // coarsest grid = 0
139: // finest grid = nlevels
140: dmhierarchy[0] = dm;
141: PetscBool usehierarchy = PETSC_FALSE;
142: if (usehierarchy) {
143: PetscCall(DMRefineHierarchy(dm, user.nlevels, &dmhierarchy[1]));
144: } else {
145: for (k = 1; k <= user.nlevels; k++) PetscCall(DMRefine(dmhierarchy[k - 1], MPI_COMM_NULL, &dmhierarchy[k]));
146: }
147: dmref = dmhierarchy[user.nlevels];
148: PetscCall(PetscObjectReference((PetscObject)dmref));
150: if (user.usemg) {
151: PetscCall(PCSetType(pc, PCMG));
152: PetscCall(PCMGSetLevels(pc, user.nlevels + 1, NULL));
153: PetscCall(PCMGSetType(pc, PC_MG_MULTIPLICATIVE));
154: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
155: PetscCall(PCMGSetCycleType(pc, PC_MG_CYCLE_V));
156: PetscCall(PCMGSetNumberSmooth(pc, 2));
158: for (k = 1; k <= user.nlevels; k++) {
159: PetscCall(DMCreateInterpolation(dmhierarchy[k - 1], dmhierarchy[k], &R, NULL));
160: PetscCall(PCMGSetInterpolation(pc, k, R));
161: PetscCall(MatDestroy(&R));
162: }
163: }
165: for (k = 1; k <= user.nlevels; k++) PetscCall(DMDestroy(&dmhierarchy[k]));
166: PetscCall(PetscFree(dmhierarchy));
167: } else {
168: dmref = dm;
169: PetscCall(PetscObjectReference((PetscObject)dm));
170: }
172: PetscCall(KSPSetDM(ksp, dmref));
173: PetscCall(KSPSetFromOptions(ksp));
175: /* Perform the actual solve */
176: PetscCall(KSPSolve(ksp, NULL, NULL));
177: PetscCall(KSPGetSolution(ksp, &x));
178: PetscCall(KSPGetRhs(ksp, &b));
180: if (user.error) {
181: PetscCall(VecDuplicate(b, &errv));
182: PetscCall(ComputeDiscreteL2Error(ksp, errv, &user));
183: PetscCall(VecDestroy(&errv));
184: }
186: if (user.io) {
187: /* Write out the solution along with the mesh */
188: PetscCall(DMMoabSetGlobalFieldVector(dmref, x));
189: #ifdef MOAB_HAVE_HDF5
190: PetscCall(DMMoabOutput(dmref, "ex36.h5m", ""));
191: #else
192: /* MOAB does not support true parallel writers that aren't HDF5 based
193: And so if you are using VTK as the output format in parallel,
194: the data could be jumbled due to the order in which the processors
195: write out their parts of the mesh and solution tags
196: */
197: PetscCall(DMMoabOutput(dmref, "ex36.vtk", ""));
198: #endif
199: }
201: /* Cleanup objects */
202: PetscCall(KSPDestroy(&ksp));
203: PetscCall(DMDestroy(&dmref));
204: PetscCall(DMDestroy(&dm));
205: PetscCall(PetscFinalize());
206: return 0;
207: }
209: PetscReal ComputeDiffusionCoefficient(PetscReal coords[3], UserContext *user)
210: {
211: if (user->problem == 2) {
212: 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;
213: else return 1.0;
214: } else return 1.0; /* problem = 1 */
215: }
217: PetscReal ComputeReactionCoefficient(PetscReal coords[3], UserContext *user)
218: {
219: if (user->problem == 2) {
220: 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;
221: else return 0.0;
222: } else return 5.0; /* problem = 1 */
223: }
225: double ExactSolution(PetscReal coords[3], UserContext *user)
226: {
227: if (user->problem == 2) {
228: const PetscScalar xx = (coords[0] - user->xyzref[0]) * (coords[0] - user->xyzref[0]);
229: const PetscScalar yy = (coords[1] - user->xyzref[1]) * (coords[1] - user->xyzref[1]);
230: const PetscScalar zz = (coords[2] - user->xyzref[2]) * (coords[2] - user->xyzref[2]);
231: return PetscExpScalar(-(xx + yy + zz) / user->nu);
232: } else return sin(PETSC_PI * coords[0]) * sin(PETSC_PI * coords[1]) * sin(PETSC_PI * coords[2]);
233: }
235: PetscReal exact_solution(PetscReal x, PetscReal y, PetscReal z)
236: {
237: PetscReal coords[3] = {x, y, z};
238: return ExactSolution(coords, 0);
239: }
241: double ForcingFunction(PetscReal coords[3], UserContext *user)
242: {
243: const PetscReal exact = ExactSolution(coords, user);
244: if (user->problem == 2) {
245: const PetscReal duxyz = ((coords[0] - user->xyzref[0]) + (coords[1] - user->xyzref[1]) + (coords[2] - user->xyzref[2]));
246: return (4.0 / user->nu * duxyz * duxyz - 6.0) * exact / user->nu;
247: } else {
248: const PetscReal reac = ComputeReactionCoefficient(coords, user);
249: return (3.0 * PETSC_PI * PETSC_PI + reac) * exact;
250: }
251: }
253: PetscErrorCode ComputeRHS_MOAB(KSP ksp, Vec b, void *ptr)
254: {
255: UserContext *user = (UserContext *)ptr;
256: DM dm;
257: PetscInt dof_indices[8], nc, npoints;
258: PetscBool dbdry[8];
259: PetscReal vpos[8 * 3];
260: PetscInt i, q, nconn;
261: const moab::EntityHandle *connect;
262: const moab::Range *elocal;
263: moab::Interface *mbImpl;
264: PetscScalar localv[8];
265: PetscReal *phi, *phypts, *jxw;
266: PetscBool elem_on_boundary;
267: PetscQuadrature quadratureObj;
269: PetscFunctionBegin;
270: PetscCall(KSPGetDM(ksp, &dm));
272: /* reset the RHS */
273: PetscCall(VecSet(b, 0.0));
275: PetscCall(DMMoabFEMCreateQuadratureDefault(user->dim, user->VPERE, &quadratureObj));
276: PetscCall(PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL));
277: PetscCall(PetscMalloc3(user->VPERE * npoints, &phi, npoints * 3, &phypts, npoints, &jxw));
279: /* get the essential MOAB mesh related quantities needed for FEM assembly */
280: PetscCall(DMMoabGetInterface(dm, &mbImpl));
281: PetscCall(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: PetscCall(DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect));
289: 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);
291: /* get the coordinates of the element vertices */
292: PetscCall(DMMoabGetVertexCoordinates(dm, nconn, connect, vpos));
294: /* get the local DoF numbers to appropriately set the element contribution in the operator */
295: PetscCall(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: PetscCall(DMMoabFEMComputeBasis(user->dim, nconn, vpos, quadratureObj, phypts, jxw, phi, NULL));
301: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
336: PetscCall(MatNullSpaceRemove(nullspace, b));
337: PetscCall(MatNullSpaceDestroy(&nullspace));
338: }
340: /* Restore vectors */
341: PetscCall(VecAssemblyBegin(b));
342: PetscCall(VecAssemblyEnd(b));
343: PetscCall(PetscFree3(phi, phypts, jxw));
344: PetscCall(PetscQuadratureDestroy(&quadratureObj));
345: PetscFunctionReturn(PETSC_SUCCESS);
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;
364: PetscFunctionBeginUser;
365: PetscCall(KSPGetDM(ksp, &dm));
367: /* get the essential MOAB mesh related quantities needed for FEM assembly */
368: PetscCall(DMMoabGetInterface(dm, &mbImpl));
369: PetscCall(DMMoabGetLocalElements(dm, &elocal));
370: PetscCall(DMMoabGetSize(dm, &nglobale, &nglobalv));
371: PetscCall(DMMoabGetHierarchyLevel(dm, &hlevel));
372: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ComputeMatrix: Level = %d, N(elements) = %d, N(vertices) = %d \n", hlevel, nglobale, nglobalv));
374: PetscCall(DMMoabFEMCreateQuadratureDefault(user->dim, user->VPERE, &quadratureObj));
375: PetscCall(PetscQuadratureGetData(quadratureObj, NULL, &nc, &npoints, NULL, NULL));
376: 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));
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: PetscCall(DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect));
384: 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);
386: /* get the coordinates of the element vertices */
387: PetscCall(DMMoabGetVertexCoordinates(dm, nconn, connect, vpos));
389: /* get the local DoF numbers to appropriately set the element contribution in the operator */
390: PetscCall(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: PetscCall(DMMoabFEMComputeBasis(user->dim, nconn, vpos, quadratureObj, phypts, jxw, phi, dphi));
396: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValuesLocal(jac, nconn, dof_indices, nconn, dof_indices, array, ADD_VALUES));
436: }
438: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
439: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
441: if (user->bcType == NEUMANN && false) {
442: MatNullSpace nullspace;
444: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
445: PetscCall(MatSetNullSpace(jac, nullspace));
446: PetscCall(MatNullSpaceDestroy(&nullspace));
447: }
448: PetscCall(PetscFree6(phi, dphi[0], dphi[1], dphi[2], phypts, jxw));
449: PetscCall(PetscQuadratureDestroy(&quadratureObj));
450: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
465: PetscCall(KSPGetDM(ksp, &dm));
467: /* get the solution vector */
468: PetscCall(KSPGetSolution(ksp, &sol));
470: /* Get the internal reference to the vector arrays */
471: PetscCall(VecGetArrayRead(sol, &x));
472: PetscCall(VecGetSize(sol, &N));
473: if (err) {
474: /* reset the error vector */
475: PetscCall(VecSet(err, 0.0));
476: /* get array reference */
477: PetscCall(VecGetArray(err, &e));
478: }
480: PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));
482: /* Compute function over the locally owned part of the grid */
483: for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
484: const moab::EntityHandle vhandle = *iter;
485: PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index));
487: /* compute the mid-point of the element and use a 1-point lumped quadrature */
488: PetscCall(DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos));
490: /* compute the discrete L2 error against the exact solution */
491: const PetscScalar lerr = (ExactSolution(vpos, user) - x[dof_index]);
492: l2err += lerr * lerr;
493: if (linferr < fabs(lerr)) linferr = fabs(lerr);
495: if (err) {
496: /* set the discrete L2 error against the exact solution */
497: e[dof_index] = lerr;
498: }
499: }
501: PetscCallMPI(MPI_Allreduce(&l2err, &global_l2, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD));
502: PetscCallMPI(MPI_Allreduce(&linferr, &global_linf, 1, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD));
503: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Computed Errors: L_2 = %f, L_inf = %f\n", sqrt(global_l2 / N), global_linf));
505: /* Restore vectors */
506: PetscCall(VecRestoreArrayRead(sol, &x));
507: if (err) PetscCall(VecRestoreArray(err, &e));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: PetscErrorCode InitializeOptions(UserContext *user)
512: {
513: const char *bcTypes[2] = {"dirichlet", "neumann"};
514: PetscInt bc;
516: PetscFunctionBegin;
517: /* set default parameters */
518: user->dim = 3; /* 3-Dimensional problem */
519: user->problem = 1;
520: user->n = 2;
521: user->nlevels = 2;
522: user->rho = 0.1;
523: user->bounds[0] = user->bounds[2] = user->bounds[4] = 0.0;
524: user->bounds[1] = user->bounds[3] = user->bounds[5] = 1.0;
525: user->xyzref[0] = user->bounds[1] / 2;
526: user->xyzref[1] = user->bounds[3] / 2;
527: user->xyzref[2] = user->bounds[5] / 2;
528: user->nu = 0.05;
529: user->usemg = PETSC_FALSE;
530: user->io = PETSC_FALSE;
531: user->usetet = PETSC_FALSE;
532: user->error = PETSC_FALSE;
533: bc = (PetscInt)DIRICHLET;
535: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "ex36.cxx");
536: PetscCall(PetscOptionsInt("-problem", "The type of problem being solved (controls forcing function)", "ex36.cxx", user->problem, &user->problem, NULL));
537: PetscCall(PetscOptionsInt("-n", "The elements in each direction", "ex36.cxx", user->n, &user->n, NULL));
538: PetscCall(PetscOptionsInt("-levels", "Number of levels in the multigrid hierarchy", "ex36.cxx", user->nlevels, &user->nlevels, NULL));
539: PetscCall(PetscOptionsReal("-rho", "The conductivity", "ex36.cxx", user->rho, &user->rho, NULL));
540: PetscCall(PetscOptionsReal("-x", "The domain size in x-direction", "ex36.cxx", user->bounds[1], &user->bounds[1], NULL));
541: PetscCall(PetscOptionsReal("-y", "The domain size in y-direction", "ex36.cxx", user->bounds[3], &user->bounds[3], NULL));
542: PetscCall(PetscOptionsReal("-z", "The domain size in y-direction", "ex36.cxx", user->bounds[5], &user->bounds[5], NULL));
543: PetscCall(PetscOptionsReal("-xref", "The x-coordinate of Gaussian center (for -problem 1)", "ex36.cxx", user->xyzref[0], &user->xyzref[0], NULL));
544: PetscCall(PetscOptionsReal("-yref", "The y-coordinate of Gaussian center (for -problem 1)", "ex36.cxx", user->xyzref[1], &user->xyzref[1], NULL));
545: PetscCall(PetscOptionsReal("-zref", "The y-coordinate of Gaussian center (for -problem 1)", "ex36.cxx", user->xyzref[2], &user->xyzref[2], NULL));
546: PetscCall(PetscOptionsReal("-nu", "The width of the Gaussian source (for -problem 1)", "ex36.cxx", user->nu, &user->nu, NULL));
547: PetscCall(PetscOptionsBool("-mg", "Use multigrid preconditioner", "ex36.cxx", user->usemg, &user->usemg, NULL));
548: PetscCall(PetscOptionsBool("-io", "Write out the solution and mesh data", "ex36.cxx", user->io, &user->io, NULL));
549: PetscCall(PetscOptionsBool("-tet", "Use tetrahedra to discretize the domain", "ex36.cxx", user->usetet, &user->usetet, NULL));
550: PetscCall(PetscOptionsBool("-error", "Compute the discrete L_2 and L_inf errors of the solution", "ex36.cxx", user->error, &user->error, NULL));
551: PetscCall(PetscOptionsEList("-bc", "Type of boundary condition", "ex36.cxx", bcTypes, 2, bcTypes[0], &bc, NULL));
552: PetscCall(PetscOptionsString("-file", "The mesh file for the problem", "ex36.cxx", "", user->filename, sizeof(user->filename), &user->use_extfile));
553: PetscOptionsEnd();
555: if (user->problem < 1 || user->problem > 2) user->problem = 1;
556: user->bcType = (BCType)bc;
557: user->VPERE = (user->usetet ? 4 : 8);
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /*TEST
563: build:
564: requires: moab !complex
566: test:
567: args: -levels 1 -nu .01 -n 4 -mg -ksp_converged_reason
569: test:
570: suffix: 2
571: nsize: 2
572: requires: hdf5
573: args: -levels 2 -nu .01 -n 2 -mg -ksp_converged_reason
575: TEST*/