Actual source code: ex12.c
petsc-3.4.5 2014-06-29
1: static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements.\n\
2: We solve the Poisson problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: #include <petscdmplex.h>
6: #include <petscsnes.h>
7: #if defined(PETSC_HAVE_EXODUSII)
8: #include <exodusII.h>
9: #endif
11: /*------------------------------------------------------------------------------
12: This code can be generated using 'bin/pythonscripts/PetscGenerateFEMQuadrature.py dim order dim 1 laplacian dim order dim 1 boundary src/snes/examples/tutorials/ex12.h'
13: -----------------------------------------------------------------------------*/
14: #include "ex12.h"
15: #include "ex12_bd.h"
17: typedef enum {NEUMANN, DIRICHLET} BCType;
18: typedef enum {RUN_FULL, RUN_TEST} RunType;
20: typedef struct {
21: DM dm; /* REQUIRED in order to use SNES evaluation functions */
22: PetscFEM fem; /* REQUIRED to use DMPlexComputeResidualFEM() */
23: PetscInt debug; /* The debugging level */
24: PetscMPIInt rank; /* The process rank */
25: PetscMPIInt numProcs; /* The number of processes */
26: RunType runType; /* Whether to run tests, or solve the full problem */
27: PetscBool jacobianMF; /* Whether to calculate the Jacobian action on the fly */
28: PetscLogEvent createMeshEvent;
29: PetscBool showInitial, showSolution;
30: /* Domain and mesh definition */
31: PetscInt dim; /* The topological mesh dimension */
32: char filename[2048]; /* The optional ExodusII file */
33: PetscBool interpolate; /* Generate intermediate mesh elements */
34: PetscReal refinementLimit; /* The largest allowable cell volume */
35: char partitioner[2048]; /* The graph partitioner */
36: /* GPU partitioning */
37: PetscInt numBatches; /* The number of cell batches per kernel */
38: PetscInt numBlocks; /* The number of concurrent blocks per kernel */
39: /* Element quadrature */
40: PetscQuadrature q[NUM_FIELDS];
41: PetscQuadrature qbd[NUM_FIELDS];
42: /* Problem definition */
43: void (*f0Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]); /* f0_u(x,y,z), and f0_p(x,y,z) */
44: void (*f1Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]); /* f1_u(x,y,z), and f1_p(x,y,z) */
45: void (*g0Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[]); /* g0_uu(x,y,z), g0_up(x,y,z), g0_pu(x,y,z), and g0_pp(x,y,z) */
46: void (*g1Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[]); /* g1_uu(x,y,z), g1_up(x,y,z), g1_pu(x,y,z), and g1_pp(x,y,z) */
47: void (*g2Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[]); /* g2_uu(x,y,z), g2_up(x,y,z), g2_pu(x,y,z), and g2_pp(x,y,z) */
48: void (*g3Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[]); /* g3_uu(x,y,z), g3_up(x,y,z), g3_pu(x,y,z), and g3_pp(x,y,z) */
49: void (*exactFuncs[NUM_BASIS_COMPONENTS_TOTAL])(const PetscReal x[], PetscScalar *u); /* The exact solution function u(x,y,z), v(x,y,z), and p(x,y,z) */
50: void (*f0BdFuncs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]); /* f0_u(x,y,z), and f0_p(x,y,z) */
51: void (*f1BdFuncs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]); /* f1_u(x,y,z), and f1_p(x,y,z) */
52: BCType bcType;
53: } AppCtx;
55: void zero(const PetscReal coords[], PetscScalar *u)
56: {
57: *u = 0.0;
58: }
60: /*
61: In 2D for Dirichlet conditions, we use exact solution:
63: u = x^2 + y^2
64: f = 4
66: so that
68: -\Delta u + f = -4 + 4 = 0
70: For Neumann conditions, we have
72: \nabla u \cdot -\hat y |_{y=0} = -(2y)|_{y=0} = 0 (bottom)
73: \nabla u \cdot \hat y |_{y=1} = (2y)|_{y=1} = 2 (top)
74: \nabla u \cdot -\hat x |_{x=0} = -(2x)|_{x=0} = 0 (left)
75: \nabla u \cdot \hat x |_{x=1} = (2x)|_{x=1} = 2 (right)
77: Which we can express as
79: \nabla u \cdot \hat n|_\Gamma = 2 (x + y)
80: */
81: void quadratic_u_2d(const PetscReal x[], PetscScalar *u)
82: {
83: *u = x[0]*x[0] + x[1]*x[1];
84: }
86: void f0_u(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
87: {
88: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
89: PetscInt comp;
91: for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 4.0;
92: }
94: void f0_bd_u(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f0[])
95: {
96: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
97: PetscInt comp;
98: PetscScalar val = 0.0;
100: if ((fabs(x[0] - 1.0) < 1.0e-9) || (fabs(x[1] - 1.0) < 1.0e-9)) {val = -2.0;}
101: for (comp = 0; comp < Ncomp; ++comp) f0[comp] = val;
102: }
104: void f0_bd_zero(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f0[])
105: {
106: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
107: PetscInt comp;
109: for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 0.0;
110: }
112: void f1_bd_zero(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f1[])
113: {
114: const PetscInt Ncomp = SPATIAL_DIM_0*NUM_BASIS_COMPONENTS_0;
115: PetscInt comp;
117: for (comp = 0; comp < Ncomp; ++comp) f1[comp] = 0.0;
118: }
120: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
121: void f1_u(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
122: {
123: const PetscInt dim = SPATIAL_DIM_0;
124: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
125: PetscInt comp, d;
127: for (comp = 0; comp < Ncomp; ++comp) {
128: for (d = 0; d < dim; ++d) {
129: f1[comp*dim+d] = gradU[comp*dim+d];
130: }
131: }
132: }
134: /* < \nabla v, \nabla u + {\nabla u}^T >
135: This just gives \nabla u, give the perdiagonal for the transpose */
136: void g3_uu(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[])
137: {
138: const PetscInt dim = SPATIAL_DIM_0;
139: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
140: PetscInt compI, d;
142: for (compI = 0; compI < Ncomp; ++compI) {
143: for (d = 0; d < dim; ++d) {
144: g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
145: }
146: }
147: }
149: /*
150: In 3D we use exact solution:
152: u = x^2 + y^2 + z^2
153: f = 6
155: so that
157: -\Delta u + f = -6 + 6 = 0
158: */
159: void quadratic_u_3d(const PetscReal x[], PetscScalar *u)
160: {
161: *u = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
162: }
166: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
167: {
168: const char *bcTypes[2] = {"neumann", "dirichlet"};
169: const char *runTypes[2] = {"full", "test"};
170: PetscInt bc, run;
171: PetscBool flg;
175: options->debug = 0;
176: options->runType = RUN_FULL;
177: options->dim = 2;
178: options->interpolate = PETSC_FALSE;
179: options->refinementLimit = 0.0;
180: options->bcType = DIRICHLET;
181: options->numBatches = 1;
182: options->numBlocks = 1;
183: options->jacobianMF = PETSC_FALSE;
184: options->showInitial = PETSC_FALSE;
185: options->showSolution = PETSC_TRUE;
187: options->fem.quad = (PetscQuadrature*) &options->q;
188: options->fem.f0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f0Funcs;
189: options->fem.f1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f1Funcs;
190: options->fem.g0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g0Funcs;
191: options->fem.g1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g1Funcs;
192: options->fem.g2Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g2Funcs;
193: options->fem.g3Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g3Funcs;
194: options->fem.bcFuncs = (PetscScalar (**)(const PetscReal[])) &options->exactFuncs;
195: options->fem.quadBd = (PetscQuadrature*) &options->qbd;
196: options->fem.f0BdFuncs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[])) &options->f0BdFuncs;
197: options->fem.f1BdFuncs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[])) &options->f1BdFuncs;
199: MPI_Comm_size(comm, &options->numProcs);
200: MPI_Comm_rank(comm, &options->rank);
201: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
202: PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, NULL);
203: run = options->runType;
204: PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 2, runTypes[options->runType], &run, NULL);
206: options->runType = (RunType) run;
208: PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
209: PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
210: #if !defined(PETSC_HAVE_EXODUSII)
211: if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
212: #endif
213: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
214: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
215: PetscStrcpy(options->partitioner, "chaco");
216: PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
217: bc = options->bcType;
218: PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);
220: options->bcType = (BCType) bc;
222: PetscOptionsInt("-gpu_batches", "The number of cell batches per kernel", "ex12.c", options->numBatches, &options->numBatches, NULL);
223: PetscOptionsInt("-gpu_blocks", "The number of concurrent blocks per kernel", "ex12.c", options->numBlocks, &options->numBlocks, NULL);
224: PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
225: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
226: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
227: PetscOptionsEnd();
229: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
230: return(0);
231: }
235: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
236: {
237: Vec lv;
238: PetscInt p;
239: PetscMPIInt rank, numProcs;
243: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
244: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
245: DMGetLocalVector(dm, &lv);
246: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
247: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
248: PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
249: for (p = 0; p < numProcs; ++p) {
250: if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
251: PetscBarrier((PetscObject) dm);
252: }
253: DMRestoreLocalVector(dm, &lv);
254: return(0);
255: }
259: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
260: {
261: PetscInt dim = user->dim;
262: const char *filename = user->filename;
263: PetscBool interpolate = user->interpolate;
264: PetscReal refinementLimit = user->refinementLimit;
265: const char *partitioner = user->partitioner;
266: size_t len;
270: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
271: PetscStrlen(filename, &len);
272: if (!len) {
273: DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
274: } else {
275: #if defined(PETSC_HAVE_EXODUSII)
276: int CPU_word_size = 0, IO_word_size = 0, exoid;
277: float version;
278: PetscMPIInt rank;
280: MPI_Comm_rank(comm, &rank);
281: if (!rank) {
282: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
283: if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
284: } else exoid = -1; /* Not used */
285: DMPlexCreateExodus(comm, exoid, interpolate, dm);
286: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
287: if (!rank) {ex_close(exoid);}
288: /* Must have boundary marker for Dirichlet conditions */
289: #endif
290: }
291: {
292: DM refinedMesh = NULL;
293: DM distributedMesh = NULL;
295: /* Refine mesh using a volume constraint */
296: DMPlexSetRefinementLimit(*dm, refinementLimit);
297: DMRefine(*dm, comm, &refinedMesh);
298: if (refinedMesh) {
299: DMDestroy(dm);
300: *dm = refinedMesh;
301: }
302: /* Distribute mesh over processes */
303: DMPlexDistribute(*dm, partitioner, 0, &distributedMesh);
304: if (distributedMesh) {
305: DMDestroy(dm);
306: *dm = distributedMesh;
307: }
308: }
309: DMSetFromOptions(*dm);
310: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
311: user->dm = *dm;
312: return(0);
313: }
317: PetscErrorCode SetupQuadrature(AppCtx *user)
318: {
320: user->fem.quad[0].numQuadPoints = NUM_QUADRATURE_POINTS_0;
321: user->fem.quad[0].quadPoints = points_0;
322: user->fem.quad[0].quadWeights = weights_0;
323: user->fem.quad[0].numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
324: user->fem.quad[0].numComponents = NUM_BASIS_COMPONENTS_0;
325: user->fem.quad[0].basis = Basis_0;
326: user->fem.quad[0].basisDer = BasisDerivatives_0;
327: user->fem.quadBd[0].numQuadPoints = NUM_QUADRATURE_POINTS_0_BD;
328: user->fem.quadBd[0].quadPoints = points_0_BD;
329: user->fem.quadBd[0].quadWeights = weights_0_BD;
330: user->fem.quadBd[0].numBasisFuncs = NUM_BASIS_FUNCTIONS_0_BD;
331: user->fem.quadBd[0].numComponents = NUM_BASIS_COMPONENTS_0_BD;
332: user->fem.quadBd[0].basis = Basis_0_BD;
333: user->fem.quadBd[0].basisDer = BasisDerivatives_0_BD;
334: return(0);
335: }
339: /*
340: There is a problem here with uninterpolated meshes. The index in numDof[] is not dimension in this case,
341: but sieve depth.
342: */
343: PetscErrorCode SetupSection(DM dm, AppCtx *user)
344: {
345: PetscSection section;
346: const PetscInt numFields = NUM_FIELDS;
347: PetscInt dim = user->dim;
348: const char *bdLabel = user->bcType == NEUMANN ? "boundary" : "marker";
349: PetscInt numBC = 0;
350: PetscInt numComp[NUM_FIELDS] = {NUM_BASIS_COMPONENTS_0};
351: PetscInt bcFields[1] = {0};
352: IS bcPoints[1] = {NULL};
353: PetscInt numDof[NUM_FIELDS*(SPATIAL_DIM_0+1)];
354: PetscInt f, d;
355: PetscBool has;
359: if (dim != SPATIAL_DIM_0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_0);
360: for (d = 0; d <= dim; ++d) {
361: numDof[0*(dim+1)+d] = numDof_0[d];
362: }
363: for (f = 0; f < numFields; ++f) {
364: for (d = 1; d < dim; ++d) {
365: if ((numDof[f*(dim+1)+d] > 0) && !user->interpolate) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
366: }
367: }
368: DMPlexHasLabel(dm, bdLabel, &has);
369: if (!has) {
370: DMLabel label;
372: DMPlexCreateLabel(dm, bdLabel);
373: DMPlexGetLabel(dm, bdLabel, &label);
374: DMPlexMarkBoundaryFaces(dm, label);
375: if (user->bcType == DIRICHLET) {
376: DMPlexLabelComplete(dm, label);
377: }
378: }
379: if (user->bcType == DIRICHLET) {
380: numBC = 1;
381: DMPlexGetStratumIS(dm, bdLabel, 1, &bcPoints[0]);
382: }
383: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, §ion);
384: PetscSectionSetFieldName(section, 0, "potential");
385: DMSetDefaultSection(dm, section);
386: if (user->bcType == DIRICHLET) {
387: ISDestroy(&bcPoints[0]);
388: }
389: return(0);
390: }
394: PetscErrorCode SetupExactSolution(DM dm, AppCtx *user)
395: {
396: PetscFEM *fem = &user->fem;
400: fem->f0Funcs[0] = f0_u;
401: fem->f1Funcs[0] = f1_u;
402: fem->g0Funcs[0] = NULL;
403: fem->g1Funcs[0] = NULL;
404: fem->g2Funcs[0] = NULL;
405: fem->g3Funcs[0] = g3_uu; /* < \nabla v, \nabla u > */
406: fem->f0BdFuncs[0] = f0_bd_zero;
407: fem->f1BdFuncs[0] = f1_bd_zero;
408: switch (user->dim) {
409: case 2:
410: user->exactFuncs[0] = quadratic_u_2d;
411: if (user->bcType == NEUMANN) {
412: fem->f0BdFuncs[0] = f0_bd_u;
413: }
414: break;
415: case 3:
416: user->exactFuncs[0] = quadratic_u_3d;
417: break;
418: default:
419: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
420: }
421: DMPlexSetFEMIntegration(dm, FEMIntegrateResidualBatch, FEMIntegrateBdResidualBatch, FEMIntegrateJacobianActionBatch, FEMIntegrateJacobianBatch);
422: return(0);
423: }
427: /*
428: FormJacobianAction - Form the global Jacobian action Y = JX from the global input X
430: Input Parameters:
431: + mat - The Jacobian shell matrix
432: - X - Global input vector
434: Output Parameter:
435: . Y - Local output vector
437: Note:
438: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
439: like a GPU, or vectorize on a multicore machine.
441: .seealso: FormJacobianActionLocal()
442: */
443: PetscErrorCode FormJacobianAction(Mat J, Vec X, Vec Y)
444: {
445: JacActionCtx *ctx;
446: DM dm;
447: Vec localX, localY;
448: PetscInt N, n;
452: #if 0
453: /* Needs petscimpl.h */
457: #endif
458: MatShellGetContext(J, &ctx);
459: dm = ctx->dm;
461: /* determine whether X = localX */
462: DMGetLocalVector(dm, &localX);
463: DMGetLocalVector(dm, &localY);
464: VecGetSize(X, &N);
465: VecGetSize(localX, &n);
467: if (n != N) { /* X != localX */
468: VecSet(localX, 0.0);
469: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
470: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
471: } else {
472: DMRestoreLocalVector(dm, &localX);
473: localX = X;
474: }
475: DMPlexComputeJacobianActionFEM(dm, J, localX, localY, ctx->user);
476: if (n != N) {
477: DMRestoreLocalVector(dm, &localX);
478: }
479: VecSet(Y, 0.0);
480: DMLocalToGlobalBegin(dm, localY, ADD_VALUES, Y);
481: DMLocalToGlobalEnd(dm, localY, ADD_VALUES, Y);
482: DMRestoreLocalVector(dm, &localY);
483: if (0) {
484: Vec r;
485: PetscReal norm;
487: VecDuplicate(X, &r);
488: MatMult(ctx->J, X, r);
489: VecAXPY(r, -1.0, Y);
490: VecNorm(r, NORM_2, &norm);
491: if (norm > 1.0e-8) {
492: PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Input:\n");
493: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
494: PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Result:\n");
495: VecView(Y, PETSC_VIEWER_STDOUT_WORLD);
496: PetscPrintf(PETSC_COMM_WORLD, "Difference:\n");
497: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
498: SETERRQ1(PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "The difference with assembled multiply is too large %g", norm);
499: }
500: VecDestroy(&r);
501: }
502: return(0);
503: }
507: int main(int argc, char **argv)
508: {
509: SNES snes; /* nonlinear solver */
510: Vec u,r; /* solution, residual vectors */
511: Mat A,J; /* Jacobian matrix */
512: MatNullSpace nullSpace; /* May be necessary for Neumann conditions */
513: AppCtx user; /* user-defined work context */
514: JacActionCtx userJ; /* context for Jacobian MF action */
515: PetscInt its; /* iterations for convergence */
516: PetscReal error = 0.0; /* L_2 error in the solution */
517: const PetscInt numComponents = NUM_BASIS_COMPONENTS_TOTAL;
520: PetscInitialize(&argc, &argv, NULL, help);
521: ProcessOptions(PETSC_COMM_WORLD, &user);
522: SNESCreate(PETSC_COMM_WORLD, &snes);
523: CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);
524: SNESSetDM(snes, user.dm);
526: SetupExactSolution(user.dm, &user);
527: SetupQuadrature(&user);
528: SetupSection(user.dm, &user);
530: DMCreateGlobalVector(user.dm, &u);
531: VecDuplicate(u, &r);
533: DMCreateMatrix(user.dm, MATAIJ, &J);
534: if (user.jacobianMF) {
535: PetscInt M, m, N, n;
537: MatGetSize(J, &M, &N);
538: MatGetLocalSize(J, &m, &n);
539: MatCreate(PETSC_COMM_WORLD, &A);
540: MatSetSizes(A, m, n, M, N);
541: MatSetType(A, MATSHELL);
542: MatSetUp(A);
543: MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
545: userJ.dm = user.dm;
546: userJ.J = J;
547: userJ.user = &user;
549: DMCreateLocalVector(user.dm, &userJ.u);
550: DMPlexProjectFunctionLocal(user.dm, numComponents, user.exactFuncs, INSERT_BC_VALUES, userJ.u);
551: MatShellSetContext(A, &userJ);
552: } else {
553: A = J;
554: }
555: if (user.bcType == NEUMANN) {
556: MatNullSpaceCreate(PetscObjectComm((PetscObject) user.dm), PETSC_TRUE, 0, NULL, &nullSpace);
557: MatSetNullSpace(J, nullSpace);
558: if (A != J) {
559: MatSetNullSpace(A, nullSpace);
560: }
561: }
563: DMSNESSetFunctionLocal(user.dm, (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);
564: DMSNESSetJacobianLocal(user.dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*))DMPlexComputeJacobianFEM,&user);
565: SNESSetJacobian(snes, A, J, NULL, NULL);
567: SNESSetFromOptions(snes);
569: DMPlexProjectFunction(user.dm, numComponents, user.exactFuncs, INSERT_ALL_VALUES, u);
570: if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
571: if (user.runType == RUN_FULL) {
572: PetscScalar (*initialGuess[numComponents])(const PetscReal x[]);
573: PetscInt c;
575: for (c = 0; c < numComponents; ++c) initialGuess[c] = zero;
576: DMPlexProjectFunction(user.dm, numComponents, initialGuess, INSERT_VALUES, u);
577: if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
578: if (user.debug) {
579: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
580: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
581: }
582: SNESSolve(snes, NULL, u);
583: SNESGetIterationNumber(snes, &its);
584: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
585: DMPlexComputeL2Diff(user.dm, user.fem.quad, user.exactFuncs, u, &error);
586: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g\n", error);
587: if (user.showSolution) {
588: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
589: VecChop(u, 3.0e-9);
590: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
591: }
592: } else {
593: PetscReal res = 0.0;
595: /* Check discretization error */
596: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
597: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
598: DMPlexComputeL2Diff(user.dm, user.fem.quad, user.exactFuncs, u, &error);
599: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
600: /* Check residual */
601: SNESComputeFunction(snes, u, r);
602: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
603: VecChop(r, 1.0e-10);
604: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
605: VecNorm(r, NORM_2, &res);
606: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
607: /* Check Jacobian */
608: {
609: Vec b;
610: MatStructure flag;
612: SNESComputeJacobian(snes, u, &A, &A, &flag);
613: VecDuplicate(u, &b);
614: VecSet(r, 0.0);
615: SNESComputeFunction(snes, r, b);
616: MatMult(A, u, r);
617: VecAXPY(r, 1.0, b);
618: VecDestroy(&b);
619: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
620: VecChop(r, 1.0e-10);
621: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
622: VecNorm(r, NORM_2, &res);
623: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
624: }
625: }
627: if (user.runType == RUN_FULL) {
628: PetscViewer viewer;
629: Vec uLocal;
631: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
632: PetscViewerSetType(viewer, PETSCVIEWERVTK);
633: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
634: PetscViewerFileSetName(viewer, "ex12_sol.vtk");
636: DMGetLocalVector(user.dm, &uLocal);
637: DMGlobalToLocalBegin(user.dm, u, INSERT_VALUES, uLocal);
638: DMGlobalToLocalEnd(user.dm, u, INSERT_VALUES, uLocal);
640: PetscObjectReference((PetscObject) user.dm); /* Needed because viewer destroys the DM */
641: PetscObjectReference((PetscObject) uLocal); /* Needed because viewer destroys the Vec */
642: PetscViewerVTKAddField(viewer, (PetscObject) user.dm, DMPlexVTKWriteAll, PETSC_VTK_POINT_FIELD, (PetscObject) uLocal);
643: DMRestoreLocalVector(user.dm, &uLocal);
644: PetscViewerDestroy(&viewer);
645: }
647: if (user.bcType == NEUMANN) {
648: MatNullSpaceDestroy(&nullSpace);
649: }
650: if (user.jacobianMF) {
651: VecDestroy(&userJ.u);
652: }
653: if (A != J) {
654: MatDestroy(&A);
655: }
656: MatDestroy(&J);
657: VecDestroy(&u);
658: VecDestroy(&r);
659: SNESDestroy(&snes);
660: DMDestroy(&user.dm);
661: PetscFinalize();
662: return 0;
663: }