Actual source code: ex12.c

petsc-3.4.4 2014-03-13
  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, &section);
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: }