Actual source code: ex12.c
petsc-3.3-p7 2013-05-11
1: static char help[] = "Bratu nonlinear PDE in 2d.\n\
2: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
3: domain, using a parallel unstructured mesh (DMMESH) to discretize it.\n\
4: The command line options include:\n\
5: -lambda <parameter>, where <parameter> indicates the problem's nonlinearity\n\
6: problem SFI: <parameter> = Bratu parameter (0 <= lambda <= 6.81)\n\n";
8: /*T
9: Concepts: SNES^parallel Bratu example
10: Concepts: DMMESH^using unstructured grids;
11: Processors: n
12: T*/
14: /* ------------------------------------------------------------------------
16: Solid Fuel Ignition (SFI) problem. This problem is modeled by
17: the partial differential equation
19: -Laplacian u - lambda*exp(u) = f(x,y), 0 < x,y < 1,
21: with boundary conditions
23: u = 0 for x = 0, x = 1, y = 0, y = 1.
25: A finite element approximation with the usual P_1 linear basis
26: is used to discretize the boundary value problem to obtain a nonlinear
27: system of equations.
29: Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options]
30: e.g.,
31: ./ex5 -draw_pause -1
32: mpiexec -n 2 ./ex5 -log_summary
34: ------------------------------------------------------------------------- */
36: /*
37: Include "petscdmmesh.h" so that we can use unstructured meshes (DMMESHs).
38: Include "petscsnes.h" so that we can use SNES solvers. Note that this
39: file automatically includes:
40: petscsys.h - base PETSc routines
41: petscvec.h - vectors petscmat.h - matrices
42: petscis.h - index sets petscksp.h - Krylov subspace methods
43: petscviewer.h - viewers petscpc.h - preconditioners
44: */
45: #include <petscdmmesh.h>
46: #include <petscsnes.h>
48: typedef enum {RUN_FULL, RUN_TEST, RUN_MESH} RunType;
49: typedef enum {NEUMANN, DIRICHLET} BCType;
51: /*------------------------------------------------------------------------------
52: This code can be generated using config/PETSc/FEM.py
54: import PETSc.FEM
55: from FIAT.reference_element import default_simplex
56: from FIAT.lagrange import Lagrange
58: generator = PETSc.FEM.QuadratureGenerator()
59: generator.setup()
60: dim = 2
61: order = 1
62: elements = [Lagrange(default_simplex(dim), order))]
63: generator.run(elements, filename)
64: -----------------------------------------------------------------------------*/
65: #include <stdlib.h>
67: #define NUM_QUADRATURE_POINTS_0 1
69: /* Quadrature points
70: - (x1,y1,x2,y2,...) */
71: static PetscReal points_0[1] = {0.0};
73: /* Quadrature weights
74: - (v1,v2,...) */
75: static PetscReal weights_0[1] = {2.0};
77: #define NUM_BASIS_FUNCTIONS_0 2
79: /* Nodal basis function evaluations
80: - basis function is fastest varying, then point */
81: static PetscReal Basis_0[2] = {
82: 0.5,
83: 0.5};
85: /* Nodal basis function derivative evaluations,
86: - derivative direction fastest varying, then basis function, then point */
87: static PetscReal BasisDerivatives_0[2] = {
88: -0.5,
89: 0.5};
91: #define NUM_DUAL_POINTS_0 2
93: /* Dual points
94: - (x1,y1,x2,y2,...) */
95: static PetscReal dualPoints_0[2] = {
96: -1.0,
97: 1.0};
103: PetscReal IntegrateDualBasis_gen_0(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
104: {
105: PetscReal refCoords[1];
106: PetscReal coords[1];
108: switch(dualIndex) {
109: case 0:
110: refCoords[0] = -1.0;
111: break;
112: case 1:
113: refCoords[0] = 1.0;
114: break;
115: default:
116: printf("dualIndex: %d\n", dualIndex);
117: throw ALE::Exception("Bad dual index");
118: }
119: for(int d = 0; d < 1; d++)
120: {
121: coords[d] = v0[d];
122: for(int e = 0; e < 1; e++)
123: {
124: coords[d] += J[d * 1 + e] * (refCoords[e] + 1.0);
125: }
126: }
127: return (*func)(coords);
128: }
134: PetscReal IntegrateBdDualBasis_gen_0(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
135: {
136: PetscReal refCoords[1];
137: PetscReal coords[2];
139: switch(dualIndex) {
140: case 0:
141: refCoords[0] = -1.0;
142: break;
143: case 1:
144: refCoords[0] = 1.0;
145: break;
146: default:
147: printf("dualIndex: %d\n", dualIndex);
148: throw ALE::Exception("Bad dual index");
149: }
150: for(int d = 0; d < 2; d++)
151: {
152: coords[d] = v0[d];
153: for(int e = 0; e < 1; e++)
154: {
155: coords[d] += J[d * 2 + e] * (refCoords[e] + 1.0);
156: }
157: }
158: return (*func)(coords);
159: }
161: #define NUM_QUADRATURE_POINTS_1 1
163: /* Quadrature points
164: - (x1,y1,x2,y2,...) */
165: static PetscReal points_1[2] = {
166: -0.333333333333,
167: -0.333333333333};
169: /* Quadrature weights
170: - (v1,v2,...) */
171: static PetscReal weights_1[1] = {2.0};
173: #define NUM_BASIS_FUNCTIONS_1 3
175: /* Nodal basis function evaluations
176: - basis function is fastest varying, then point */
177: static PetscReal Basis_1[3] = {
178: 0.333333333333,
179: 0.333333333333,
180: 0.333333333333};
182: /* Nodal basis function derivative evaluations,
183: - derivative direction fastest varying, then basis function, then point */
184: static PetscReal BasisDerivatives_1[6] = {
185: -0.5,
186: -0.5,
187: 0.5,
188: 0.0,
189: 0.0,
190: 0.5};
192: #define NUM_DUAL_POINTS_1 3
194: /* Dual points
195: - (x1,y1,x2,y2,...) */
196: static PetscReal dualPoints_1[6] = {
197: -1.0,
198: -1.0,
199: 1.0,
200: -1.0,
201: -1.0,
202: 1.0};
208: PetscReal IntegrateDualBasis_gen_1(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
209: {
210: PetscReal refCoords[2];
211: PetscReal coords[2];
213: switch(dualIndex) {
214: case 0:
215: refCoords[0] = -1.0;
216: refCoords[1] = -1.0;
217: break;
218: case 1:
219: refCoords[0] = 1.0;
220: refCoords[1] = -1.0;
221: break;
222: case 2:
223: refCoords[0] = -1.0;
224: refCoords[1] = 1.0;
225: break;
226: default:
227: printf("dualIndex: %d\n", dualIndex);
228: throw ALE::Exception("Bad dual index");
229: }
230: for(int d = 0; d < 2; d++)
231: {
232: coords[d] = v0[d];
233: for(int e = 0; e < 2; e++)
234: {
235: coords[d] += J[d * 2 + e] * (refCoords[e] + 1.0);
236: }
237: }
238: return (*func)(coords);
239: }
245: PetscReal IntegrateBdDualBasis_gen_1(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
246: {
247: PetscReal refCoords[2];
248: PetscReal coords[3];
250: switch(dualIndex) {
251: case 0:
252: refCoords[0] = -1.0;
253: refCoords[1] = -1.0;
254: break;
255: case 1:
256: refCoords[0] = 1.0;
257: refCoords[1] = -1.0;
258: break;
259: case 2:
260: refCoords[0] = -1.0;
261: refCoords[1] = 1.0;
262: break;
263: default:
264: printf("dualIndex: %d\n", dualIndex);
265: throw ALE::Exception("Bad dual index");
266: }
267: for(int d = 0; d < 3; d++)
268: {
269: coords[d] = v0[d];
270: for(int e = 0; e < 2; e++)
271: {
272: coords[d] += J[d * 3 + e] * (refCoords[e] + 1.0);
273: }
274: }
275: return (*func)(coords);
276: }
279: #define NUM_QUADRATURE_POINTS_2 1
281: /* Quadrature points
282: - (x1,y1,x2,y2,...) */
283: static PetscReal points_2[3] = {
284: -0.5,
285: -0.5,
286: -0.5};
288: /* Quadrature weights
289: - (v1,v2,...) */
290: static PetscReal weights_2[1] = {1.33333333333};
292: #define NUM_BASIS_FUNCTIONS_2 4
294: /* Nodal basis function evaluations
295: - basis function is fastest varying, then point */
296: static PetscReal Basis_2[4] = {
297: 0.25,
298: 0.25,
299: 0.25,
300: 0.25};
302: /* Nodal basis function derivative evaluations,
303: - derivative direction fastest varying, then basis function, then point */
304: static PetscReal BasisDerivatives_2[12] = {
305: -0.5,
306: -0.5,
307: -0.5,
308: 0.5,
309: 0.0,
310: 0.0,
311: 0.0,
312: 0.5,
313: 0.0,
314: 0.0,
315: 0.0,
316: 0.5};
318: #define NUM_DUAL_POINTS_2 4
320: /* Dual points
321: - (x1,y1,x2,y2,...) */
322: static PetscReal dualPoints_2[12] = {
323: -1.0,
324: -1.0,
325: -1.0,
326: 1.0,
327: -1.0,
328: -1.0,
329: -1.0,
330: 1.0,
331: -1.0,
332: -1.0,
333: -1.0,
334: 1.0};
340: PetscReal IntegrateDualBasis_gen_2(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
341: {
342: PetscReal refCoords[3];
343: PetscReal coords[3];
345: switch(dualIndex) {
346: case 0:
347: refCoords[0] = -1.0;
348: refCoords[1] = -1.0;
349: refCoords[2] = -1.0;
350: break;
351: case 1:
352: refCoords[0] = 1.0;
353: refCoords[1] = -1.0;
354: refCoords[2] = -1.0;
355: break;
356: case 2:
357: refCoords[0] = -1.0;
358: refCoords[1] = 1.0;
359: refCoords[2] = -1.0;
360: break;
361: case 3:
362: refCoords[0] = -1.0;
363: refCoords[1] = -1.0;
364: refCoords[2] = 1.0;
365: break;
366: default:
367: printf("dualIndex: %d\n", dualIndex);
368: throw ALE::Exception("Bad dual index");
369: }
370: for(int d = 0; d < 3; d++)
371: {
372: coords[d] = v0[d];
373: for(int e = 0; e < 3; e++)
374: {
375: coords[d] += J[d * 3 + e] * (refCoords[e] + 1.0);
376: }
377: }
378: return (*func)(coords);
379: }
385: PetscReal IntegrateBdDualBasis_gen_2(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
386: {
387: PetscReal refCoords[3];
388: PetscReal coords[4];
390: switch(dualIndex) {
391: case 0:
392: refCoords[0] = -1.0;
393: refCoords[1] = -1.0;
394: refCoords[2] = -1.0;
395: break;
396: case 1:
397: refCoords[0] = 1.0;
398: refCoords[1] = -1.0;
399: refCoords[2] = -1.0;
400: break;
401: case 2:
402: refCoords[0] = -1.0;
403: refCoords[1] = 1.0;
404: refCoords[2] = -1.0;
405: break;
406: case 3:
407: refCoords[0] = -1.0;
408: refCoords[1] = -1.0;
409: refCoords[2] = 1.0;
410: break;
411: default:
412: printf("dualIndex: %d\n", dualIndex);
413: throw ALE::Exception("Bad dual index");
414: }
415: for(int d = 0; d < 4; d++)
416: {
417: coords[d] = v0[d];
418: for(int e = 0; e < 3; e++)
419: {
420: coords[d] += J[d * 4 + e] * (refCoords[e] + 1.0);
421: }
422: }
423: return (*func)(coords);
424: }
426: /*------------------------------------------------------------------------------
427: end of generated code
428: -----------------------------------------------------------------------------*/
430: /*
431: User-defined application context - contains data needed by the
432: application-provided call-back routines, FormJacobianLocal() and
433: FormFunctionLocal().
434: */
435: typedef struct {
436: DM dm; /* The unstructured mesh data structure */
437: PetscInt debug; /* The debugging level */
438: PetscMPIInt rank; /* The process rank */
439: PetscMPIInt numProcs; /* The number of processes */
440: RunType run; /* The run type */
441: PetscInt dim; /* The topological mesh dimension */
442: PetscBool interpolate; /* Generate intermediate mesh elements */
443: PetscReal refinementLimit; /* The largest allowable cell volume */
444: char partitioner[2048]; /* The graph partitioner */
445: /* Element quadrature */
446: PetscQuadrature q;
447: /* Problem specific parameters */
448: BCType bcType; /* The type of boundary conditions */
449: PetscReal lambda; /* The Bratu problem parameter */
450: PetscScalar (*rhsFunc)(const PetscReal []); /* The rhs function f(x,y,z) */
451: PetscScalar (*exactFunc)(const PetscReal []); /* The exact solution function u(x,y,z) */
452: Vec exactSol; /* The discrete exact solution */
453: Vec error; /* The discrete cell-wise error */
454: } AppCtx;
456: /*
457: User-defined routines
458: */
459: extern PetscErrorCode FormInitialGuess(Vec X, PetscScalar (*guessFunc)(const PetscReal []), InsertMode mode, AppCtx *user);
460: extern PetscErrorCode FormFunctionLocal(DM dm, Vec X, Vec F, AppCtx *user);
461: extern PetscErrorCode FormJacobianLocal(DM dm, Vec X, Mat J, AppCtx *user);
463: PetscReal lambda = 0.0;
464: PetscScalar guess(const PetscReal coords[]) {
465: PetscScalar scale = lambda/(lambda+1.0);
466: return scale*(0.5 - fabs(coords[0]-0.5))*(0.5 - fabs(coords[1]-0.5));
467: }
469: PetscScalar zero(const PetscReal coords[]) {
470: return 0.0;
471: }
473: PetscScalar constant(const PetscReal x[]) {
474: return -4.0;
475: };
477: PetscScalar nonlinear_2d(const PetscReal x[]) {
478: return -4.0 - lambda*PetscExpScalar(x[0]*x[0] + x[1]*x[1]);
479: };
481: PetscScalar linear_2d(const PetscReal x[]) {
482: return -6.0*(x[0] - 0.5) - 6.0*(x[1] - 0.5);
483: };
485: PetscScalar quadratic_2d(const PetscReal x[]) {
486: return x[0]*x[0] + x[1]*x[1];
487: };
489: PetscScalar cubic_2d(const PetscReal x[]) {
490: return x[0]*x[0]*x[0] - 1.5*x[0]*x[0] + x[1]*x[1]*x[1] - 1.5*x[1]*x[1] + 0.5;
491: };
493: PetscScalar nonlinear_3d(const PetscReal x[]) {
494: return -4.0 - lambda*PetscExpScalar((2.0/3.0)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]));
495: };
497: PetscScalar linear_3d(const PetscReal x[]) {
498: return -6.0*(x[0] - 0.5) - 6.0*(x[1] - 0.5) - 6.0*(x[2] - 0.5);
499: };
501: PetscScalar quadratic_3d(const PetscReal x[]) {
502: return (2.0/3.0)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
503: };
505: PetscScalar cubic_3d(const PetscReal x[]) {
506: return x[0]*x[0]*x[0] - 1.5*x[0]*x[0] + x[1]*x[1]*x[1] - 1.5*x[1]*x[1] + x[2]*x[2]*x[2] - 1.5*x[2]*x[2] + 0.75;
507: };
511: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {
512: const char *runTypes[3] = {"full", "test", "mesh"};
513: const char *bcTypes[2] = {"neumann", "dirichlet"};
514: PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0.0;
515: PetscInt run, bc;
519: options->debug = 0;
520: options->run = RUN_FULL;
521: options->dim = 2;
522: options->interpolate = PETSC_FALSE;
523: options->refinementLimit = 0.0;
524: options->bcType = DIRICHLET;
525: options->lambda = 6.0;
526: options->rhsFunc = zero;
528: MPI_Comm_size(comm, &options->numProcs);
529: MPI_Comm_rank(comm, &options->rank);
530: PetscOptionsBegin(comm, "", "Bratu Problem Options", "DMMESH");
531: PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, PETSC_NULL);
532: run = options->run;
533: PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 3, runTypes[options->run], &run, PETSC_NULL);
534: options->run = (RunType) run;
535: PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, PETSC_NULL);
536: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, PETSC_NULL);
537: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, PETSC_NULL);
538: PetscStrcpy(options->partitioner, "chaco");
539: PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, PETSC_NULL);
540: bc = options->bcType;
541: PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,2,bcTypes[options->bcType],&bc,PETSC_NULL);
542: options->bcType = (BCType) bc;
543: PetscOptionsReal("-lambda", "The parameter controlling nonlinearity", "ex12.c", options->lambda, &options->lambda, PETSC_NULL);
544: if (options->lambda >= bratu_lambda_max || options->lambda < bratu_lambda_min) {
545: SETERRQ3(PETSC_COMM_WORLD, 1, "Lambda, %g, is out of range, [%g, %g)", options->lambda, bratu_lambda_min, bratu_lambda_max);
546: }
547: PetscOptionsEnd();
548: lambda = options->lambda;
549: return(0);
550: };
554: PetscErrorCode SetupQuadrature(AppCtx *user) {
556: switch(user->dim) {
557: case 1:
558: user->q.numQuadPoints = NUM_QUADRATURE_POINTS_0;
559: user->q.quadPoints = points_0;
560: user->q.quadWeights = weights_0;
561: user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
562: user->q.basis = Basis_0;
563: user->q.basisDer = BasisDerivatives_0;
564: break;
565: case 2:
566: user->q.numQuadPoints = NUM_QUADRATURE_POINTS_1;
567: user->q.quadPoints = points_1;
568: user->q.quadWeights = weights_1;
569: user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_1;
570: user->q.basis = Basis_1;
571: user->q.basisDer = BasisDerivatives_1;
572: break;
573: case 3:
574: user->q.numQuadPoints = NUM_QUADRATURE_POINTS_2;
575: user->q.quadPoints = points_2;
576: user->q.quadWeights = weights_2;
577: user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_2;
578: user->q.basis = Basis_2;
579: user->q.basisDer = BasisDerivatives_2;
580: break;
581: default:
582: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
583: }
584: return(0);
585: }
589: PetscErrorCode SetupSection(AppCtx *user) {
590: PetscSection section;
591: /* These can be generated using config/PETSc/FEM.py */
592: PetscInt numDof_0[2] = {1, 0};
593: PetscInt numDof_1[3] = {1, 0, 0};
594: PetscInt numDof_2[4] = {1, 0, 0, 0};
595: PetscInt *numDof = PETSC_NULL;
596: PetscInt numBC = 0;
597: PetscInt bcField[1] = {0};
598: IS bcPoints[1] = {PETSC_NULL};
602: switch(user->dim) {
603: case 1:
604: numDof = numDof_0;
605: break;
606: case 2:
607: numDof = numDof_1;
608: break;
609: case 3:
610: numDof = numDof_2;
611: break;
612: default:
613: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid spatial dimension %d", user->dim);
614: }
615: if (user->bcType == DIRICHLET) {
616: numBC = 1;
617: DMMeshGetStratumIS(user->dm, "marker", 1, &bcPoints[0]);
618: }
619: DMMeshCreateSection(user->dm, user->dim, 1, PETSC_NULL, numDof, numBC, bcField, bcPoints, §ion);
620: DMMeshSetDefaultSection(user->dm, section);
621: if (user->bcType == DIRICHLET) {
622: ISDestroy(&bcPoints[0]);
623: }
624: return(0);
625: }
629: PetscErrorCode SetupExactSolution(AppCtx *user) {
631: switch(user->dim) {
632: case 2:
633: if (user->bcType == DIRICHLET) {
634: if (user->lambda > 0.0) {
635: user->rhsFunc = nonlinear_2d;
636: user->exactFunc = quadratic_2d;
637: } else {
638: user->rhsFunc = constant;
639: user->exactFunc = quadratic_2d;
640: }
641: } else {
642: user->rhsFunc = linear_2d;
643: user->exactFunc = cubic_2d;
644: }
645: break;
646: case 3:
647: if (user->bcType == DIRICHLET) {
648: if (user->lambda > 0.0) {
649: user->rhsFunc = nonlinear_3d;
650: user->exactFunc = quadratic_3d;
651: } else {
652: user->rhsFunc = constant;
653: user->exactFunc = quadratic_3d;
654: }
655: } else {
656: user->rhsFunc = linear_3d;
657: user->exactFunc = cubic_3d;
658: }
659: break;
660: default:
661: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
662: }
663: return(0);
664: }
668: PetscErrorCode ComputeError(Vec X, PetscReal *error, AppCtx *user) {
669: PetscScalar (*exactFunc)(const PetscReal []) = user->exactFunc;
670: const PetscInt debug = user->debug;
671: const PetscInt dim = user->dim;
672: const PetscInt numQuadPoints = user->q.numQuadPoints;
673: const PetscReal *quadPoints = user->q.quadPoints;
674: const PetscReal *quadWeights = user->q.quadWeights;
675: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
676: const PetscReal *basis = user->q.basis;
677: Vec localX;
678: PetscReal *coords, *v0, *J, *invJ, detJ;
679: PetscReal localError;
680: PetscInt cStart, cEnd;
681: PetscErrorCode ierr;
684: DMGetLocalVector(user->dm, &localX);
685: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
686: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
687: PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
688: DMMeshGetHeightStratum(user->dm, 0, &cStart, &cEnd);
689: for(PetscInt c = cStart; c < cEnd; ++c) {
690: const PetscScalar *x;
691: PetscReal elemError = 0.0;
693: DMMeshComputeCellGeometry(user->dm, c, v0, J, invJ, &detJ);
694: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
695: DMMeshVecGetClosure(user->dm, localX, c, &x);
696: if (debug) {DMPrintCellVector(c, "Solution", numBasisFuncs, x);}
697: for(int q = 0; q < numQuadPoints; ++q) {
698: for(int d = 0; d < dim; d++) {
699: coords[d] = v0[d];
700: for(int e = 0; e < dim; e++) {
701: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
702: }
703: }
704: const PetscScalar funcVal = (*exactFunc)(coords);
705: PetscReal interpolant = 0.0;
706: for(int f = 0; f < numBasisFuncs; ++f) {
707: interpolant += x[f]*basis[q*numBasisFuncs+f];
708: }
709: elemError += PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ;
710: }
711: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d error %g\n", c, elemError);}
712: localError += elemError;
713: }
714: PetscFree4(coords,v0,J,invJ);
715: DMRestoreLocalVector(user->dm, &localX);
716: MPI_Allreduce(&localError, error, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);
717: *error = PetscSqrtReal(*error);
718: return(0);
719: }
723: int main(int argc, char **argv)
724: {
725: SNES snes; /* nonlinear solver */
726: Vec u,r; /* solution, residual vectors */
727: Mat A,J; /* Jacobian matrix */
728: AppCtx user; /* user-defined work context */
729: PetscInt its; /* iterations for convergence */
730: PetscReal error; /* L_2 error in the solution */
733: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
734: Initialize program
735: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
736: PetscInitialize(&argc, &argv, PETSC_NULL, help);
738: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
739: Initialize problem parameters
740: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
741: ProcessOptions(PETSC_COMM_WORLD, &user);
743: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
744: Create nonlinear solver context
745: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
746: SNESCreate(PETSC_COMM_WORLD, &snes);
748: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
749: Create unstructured mesh (DMMESH) to manage parallel grid and vectors
750: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
751: DMMeshCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.interpolate, &user.dm);
752: {
753: DM refinedMesh = PETSC_NULL;
754: DM distributedMesh = PETSC_NULL;
756: /* Refine mesh using a volume constraint */
757: DMMeshRefine(user.dm, user.refinementLimit, user.interpolate, &refinedMesh);
758: if (refinedMesh) {
759: DMDestroy(&user.dm);
760: user.dm = refinedMesh;
761: }
762: #if 0
763: {
764: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
765: DMMeshGetMesh(user.dm, oldMesh);
766: oldMesh->setDebug(1);
767: }
768: #endif
769: /* Distribute mesh over processes */
770: DMMeshDistribute(user.dm, user.partitioner, &distributedMesh);
771: if (distributedMesh) {
772: DMDestroy(&user.dm);
773: user.dm = distributedMesh;
774: }
775: /* Mark boundary cells for higher order element calculations */
776: if (user.bcType == DIRICHLET) {
777: DMMeshMarkBoundaryCells(user.dm, "marker", 1, 2);
778: }
779: }
780: DMSetFromOptions(user.dm);
781: SNESSetDM(snes, user.dm);
783: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
784: Setup dof layout.
786: For a DMDA, this is automatic given the number of dof at each vertex.
787: However, for a DMMesh, we need to specify this.
788: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
789: SetupExactSolution(&user);
790: SetupQuadrature(&user);
791: SetupSection(&user);
793: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
794: Extract global vectors from DM; then duplicate for remaining
795: vectors that are the same types
796: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
797: DMCreateGlobalVector(user.dm, &u);
798: VecDuplicate(u, &r);
800: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
801: Create matrix data structure; set Jacobian evaluation routine
803: Set Jacobian matrix data structure and default Jacobian evaluation
804: routine. User can override with:
805: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
806: (unless user explicitly sets preconditioner)
807: -snes_mf_operator : form preconditioning matrix as set by the user,
808: but use matrix-free approx for Jacobian-vector
809: products within Newton-Krylov method
810: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
811: /* J can be type of MATAIJ, MATBAIJ or MATSBAIJ */
812: DMCreateMatrix(user.dm, MATAIJ, &J);
813: A = J;
814: SNESSetJacobian(snes, A, J, SNESDMMeshComputeJacobian, &user);
816: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
817: Set local function evaluation routine
818: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
819: DMMeshSetLocalFunction(user.dm, (DMMeshLocalFunction1) FormFunctionLocal);
820: DMMeshSetLocalJacobian(user.dm, (DMMeshLocalJacobian1) FormJacobianLocal);
821: SNESSetFunction(snes, r, SNESDMMeshComputeFunction, &user);
823: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
824: Customize nonlinear solver; set runtime options
825: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
826: SNESSetFromOptions(snes);
828: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
829: Setup boundary conditions
830: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
831: FormInitialGuess(u, user.exactFunc, INSERT_ALL_VALUES, &user);
832: if (user.run == RUN_FULL) {
833: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
834: Evaluate initial guess
835: Note: The user should initialize the vector u, with the initial guess
836: for the nonlinear solver prior to calling SNESSolve(). In particular,
837: to employ an initial guess of zero, the user should explicitly set
838: this vector to zero by calling VecSet().
839: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
840: FormInitialGuess(u, guess, INSERT_VALUES, &user);
841: if (user.debug) {
842: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
843: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
844: }
845: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
846: Solve nonlinear system
847: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
848: SNESSolve(snes, PETSC_NULL, u);
849: SNESGetIterationNumber(snes, &its);
850: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
851: ComputeError(u, &error, &user);
852: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
853: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
854: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
855: } else {
856: PetscReal res;
858: /* Check discretization error */
859: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
860: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
861: ComputeError(u, &error, &user);
862: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
863: /* Check residual */
864: SNESDMMeshComputeFunction(snes, u, r, &user);
865: VecNorm(r, NORM_2, &res);
866: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
867: }
869: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
870: Output results
871: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
872: if (0) {
873: PetscViewer viewer;
875: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
876: /*PetscViewerSetType(viewer, PETSCVIEWERDRAW);
877: PetscViewerDrawSetInfo(viewer, PETSC_NULL, "Solution", PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE); */
878: PetscViewerSetType(viewer, PETSCVIEWERASCII);
879: PetscViewerFileSetName(viewer, "ex12_sol.vtk");
880: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
881: DMView(user.dm, viewer);
882: VecView(u, viewer);
883: PetscViewerDestroy(&viewer);
884: }
886: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
887: Free work space. All PETSc objects should be destroyed when they
888: are no longer needed.
889: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
890: if (A != J) {
891: MatDestroy(&A);
892: }
893: MatDestroy(&J);
894: VecDestroy(&u);
895: VecDestroy(&r);
896: SNESDestroy(&snes);
897: DMDestroy(&user.dm);
898: PetscFinalize();
899: return 0;
900: }
904: /*
905: FormInitialGuess - Forms initial approximation.
907: Input Parameters:
908: + user - user-defined application context
909: - guessFunc - The coordinate function to use for the guess
911: Output Parameter:
912: . X - vector
913: */
914: PetscErrorCode FormInitialGuess(Vec X, PetscScalar (*guessFunc)(const PetscReal []), InsertMode mode, AppCtx *user)
915: {
916: Vec localX, coordinates;
917: PetscSection section, cSection;
918: PetscInt vStart, vEnd;
922: DMGetLocalVector(user->dm, &localX);
923: DMMeshGetDepthStratum(user->dm, 0, &vStart, &vEnd);
924: DMMeshGetDefaultSection(user->dm, §ion);
925: DMMeshGetCoordinateSection(user->dm, &cSection);
926: DMMeshGetCoordinateVec(user->dm, &coordinates);
927: for(PetscInt v = vStart; v < vEnd; ++v) {
928: PetscScalar values[1];
929: PetscScalar *coords;
931: VecGetValuesSection(coordinates, cSection, v, &coords);
932: values[0] = (*guessFunc)(coords);
933: VecSetValuesSection(localX, section, v, values, mode);
934: }
935: VecDestroy(&coordinates);
936: PetscSectionDestroy(&cSection);
938: DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X);
939: DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X);
940: DMRestoreLocalVector(user->dm, &localX);
941: #if 0
942: /* This is necessary for higher order elements */
943: MeshGetSectionReal(mesh, "exactSolution", &this->_options.exactSol.section);
944: const Obj<PETSC_MESH_TYPE::real_section_type>& s = this->_mesh->getRealSection("exactSolution");
945: this->_mesh->setupField(s);
946: const Obj<PETSC_MESH_TYPE::label_sequence>& cells = this->_mesh->heightStratum(0);
947: const Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = this->_mesh->getRealSection("coordinates");
948: const int localDof = this->_mesh->sizeWithBC(s, *cells->begin());
949: PETSC_MESH_TYPE::real_section_type::value_type *values = new PETSC_MESH_TYPE::real_section_type::value_type[localDof];
950: PetscReal *v0 = new PetscReal[dim()];
951: PetscReal *J = new PetscReal[dim()*dim()];
952: PetscReal detJ;
953: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV((int) pow(this->_mesh->getSieve()->getMaxConeSize(), this->_mesh->depth())+1, true);
955: for(PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
956: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), *c_iter, pV);
957: const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
958: const int oSize = pV.getSize();
959: int v = 0;
961: this->_mesh->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
962: for(int cl = 0; cl < oSize; ++cl) {
963: const int pointDim = s->getFiberDimension(oPoints[cl]);
965: if (pointDim) {
966: for(int d = 0; d < pointDim; ++d, ++v) {
967: values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
968: }
969: }
970: }
971: this->_mesh->update(s, *c_iter, values);
972: pV.clear();
973: }
974: delete [] values;
975: delete [] v0;
976: delete [] J;
977: #endif
978: return(0);
979: }
983: /*
984: FormFunctionLocal - Evaluates nonlinear function, F(x).
985: */
986: PetscErrorCode FormFunctionLocal(DM dm, Vec X, Vec F, AppCtx *user)
987: {
988: PetscScalar (*rhsFunc)(const PetscReal []) = user->rhsFunc;
989: const PetscInt debug = user->debug;
990: const PetscInt dim = user->dim;
991: const PetscInt numQuadPoints = user->q.numQuadPoints;
992: const PetscReal *quadPoints = user->q.quadPoints;
993: const PetscReal *quadWeights = user->q.quadWeights;
994: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
995: const PetscReal *basis = user->q.basis;
996: const PetscReal *basisDer = user->q.basisDer;
997: const PetscReal lambda = user->lambda;
998: PetscReal *coords, *v0, *J, *invJ, detJ;
999: PetscScalar *realSpaceDer, *fieldGrad, *elemVec;
1000: PetscInt cStart, cEnd;
1001: PetscErrorCode ierr;
1004: VecSet(F, 0.0);
1005: PetscMalloc3(dim,PetscScalar,&realSpaceDer,dim,PetscScalar,&fieldGrad,numBasisFuncs,PetscScalar,&elemVec);
1006: PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
1007: DMMeshGetHeightStratum(user->dm, 0, &cStart, &cEnd);
1008: for(PetscInt c = cStart; c < cEnd; ++c) {
1009: const PetscScalar *x;
1011: PetscMemzero(elemVec, numBasisFuncs * sizeof(PetscScalar));
1012: DMMeshComputeCellGeometry(user->dm, c, v0, J, invJ, &detJ);
1013: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1014: DMMeshVecGetClosure(user->dm, X, c, &x);
1015: if (debug) {DMPrintCellVector(c, "Solution", numBasisFuncs, x);}
1017: for(int q = 0; q < numQuadPoints; ++q) {
1018: PetscScalar fieldVal = 0.0;
1020: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
1021: for(int d = 0; d < dim; ++d) {
1022: fieldGrad[d] = 0.0;
1023: coords[d] = v0[d];
1024: for(int e = 0; e < dim; ++e) {
1025: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1026: }
1027: if (debug) {PetscPrintf(PETSC_COMM_SELF, " coords[%d] %g\n", d, coords[d]);}
1028: }
1029: for(int f = 0; f < numBasisFuncs; ++f) {
1030: fieldVal += x[f]*basis[q*numBasisFuncs+f];
1032: for(int d = 0; d < dim; ++d) {
1033: realSpaceDer[d] = 0.0;
1034: for(int e = 0; e < dim; ++e) {
1035: realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
1036: }
1037: fieldGrad[d] += realSpaceDer[d]*x[f];
1038: }
1039: }
1040: if (debug) {
1041: for(int d = 0; d < dim; ++d) {
1042: PetscPrintf(PETSC_COMM_SELF, " fieldGrad[%d] %g\n", d, fieldGrad[d]);
1043: }
1044: }
1045: const PetscScalar funcVal = (*rhsFunc)(coords);
1046: for(int f = 0; f < numBasisFuncs; ++f) {
1047: /* Constant term: -f(x) */
1048: elemVec[f] -= basis[q*numBasisFuncs+f]*funcVal*quadWeights[q]*detJ;
1049: /* Linear term: -\Delta u */
1050: PetscScalar product = 0.0;
1051: for(int d = 0; d < dim; ++d) {
1052: realSpaceDer[d] = 0.0;
1053: for(int e = 0; e < dim; ++e) {
1054: realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
1055: }
1056: product += realSpaceDer[d]*fieldGrad[d];
1057: }
1058: elemVec[f] += product*quadWeights[q]*detJ;
1059: /* Nonlinear term: -\lambda e^{u} */
1060: elemVec[f] -= basis[q*numBasisFuncs+f]*lambda*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
1061: }
1062: }
1063: if (debug) {DMPrintCellVector(c, "Residual", numBasisFuncs, elemVec);}
1064: DMMeshVecSetClosure(user->dm, F, c, elemVec, ADD_VALUES);
1065: }
1066: PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
1067: PetscFree3(realSpaceDer,fieldGrad,elemVec);
1068: PetscFree4(coords,v0,J,invJ);
1070: PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");
1071: for(int p = 0; p < user->numProcs; ++p) {
1072: if (p == user->rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
1073: PetscBarrier((PetscObject) user->dm);
1074: }
1075: return(0);
1076: }
1080: /*
1081: FormJacobianLocal - Evaluates Jacobian matrix.
1082: */
1083: PetscErrorCode FormJacobianLocal(DM dm, Vec X, Mat Jac, AppCtx *user)
1084: {
1085: const PetscInt debug = user->debug;
1086: const PetscInt dim = user->dim;
1087: const PetscInt numQuadPoints = user->q.numQuadPoints;
1088: const PetscReal *quadWeights = user->q.quadWeights;
1089: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
1090: const PetscReal *basis = user->q.basis;
1091: const PetscReal *basisDer = user->q.basisDer;
1092: const PetscReal lambda = user->lambda;
1093: PetscReal *v0, *J, *invJ, detJ;
1094: PetscScalar *realSpaceTestDer, *realSpaceBasisDer, *elemMat;
1095: PetscInt cStart, cEnd;
1096: PetscErrorCode ierr;
1099: MatZeroEntries(Jac);
1100: PetscMalloc3(dim,PetscScalar,&realSpaceTestDer,dim,PetscScalar,&realSpaceBasisDer,numBasisFuncs*numBasisFuncs,PetscScalar,&elemMat);
1101: PetscMalloc3(dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
1102: DMMeshGetHeightStratum(user->dm, 0, &cStart, &cEnd);
1103: for(PetscInt c = cStart; c < cEnd; ++c) {
1104: const PetscScalar *x;
1106: PetscMemzero(elemMat, numBasisFuncs*numBasisFuncs * sizeof(PetscScalar));
1107: DMMeshComputeCellGeometry(user->dm, c, v0, J, invJ, &detJ);
1108: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1109: DMMeshVecGetClosure(user->dm, X, c, &x);
1111: for(int q = 0; q < numQuadPoints; ++q) {
1112: PetscScalar fieldVal = 0.0;
1114: for(int f = 0; f < numBasisFuncs; ++f) {
1115: fieldVal += x[f]*basis[q*numBasisFuncs+f];
1116: }
1117: for(int f = 0; f < numBasisFuncs; ++f) {
1118: for(int d = 0; d < dim; ++d) {
1119: realSpaceTestDer[d] = 0.0;
1120: for(int e = 0; e < dim; ++e) {
1121: realSpaceTestDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
1122: }
1123: }
1124: for(int g = 0; g < numBasisFuncs; ++g) {
1125: for(int d = 0; d < dim; ++d) {
1126: realSpaceBasisDer[d] = 0.0;
1127: for(int e = 0; e < dim; ++e) {
1128: realSpaceBasisDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+g)*dim+e];
1129: }
1130: }
1131: /* Linear term: -\Delta u */
1132: PetscScalar product = 0.0;
1133: for(int d = 0; d < dim; ++d) product += realSpaceTestDer[d]*realSpaceBasisDer[d];
1134: elemMat[f*numBasisFuncs+g] += product*quadWeights[q]*detJ;
1135: /* Nonlinear term: -\lambda e^{u} */
1136: elemMat[f*numBasisFuncs+g] -= basis[q*numBasisFuncs+f]*basis[q*numBasisFuncs+g]*lambda*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
1137: }
1138: }
1139: }
1140: if (debug) {DMPrintCellMatrix(c, "Jacobian", numBasisFuncs, numBasisFuncs, elemMat);}
1141: DMMeshMatSetClosure(user->dm, Jac, c, elemMat, ADD_VALUES);
1142: }
1143: PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
1144: PetscFree3(realSpaceTestDer,realSpaceBasisDer,elemMat);
1145: PetscFree3(v0,J,invJ);
1147: /* Assemble matrix, using the 2-step process:
1148: MatAssemblyBegin(), MatAssemblyEnd(). */
1149: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
1150: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
1151: /* Tell the matrix we will never add a new nonzero location to the
1152: matrix. If we do, it will generate an error. */
1153: MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1154: if (user->bcType == NEUMANN) {
1155: MatNullSpace nullsp;
1157: MatNullSpaceCreate(((PetscObject) Jac)->comm, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
1158: MatSetNullSpace(Jac, nullsp);
1159: }
1160: return(0);
1161: }