Actual source code: ex12.c

petsc-dev 2014-07-30
Report Typos and Errors
  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\
  4: This example supports discretized auxiliary fields (conductivity) as well as\n\
  5: multilevel nonlienar solvers.\n\n\n";

  7: #include <petscdmplex.h>
  8: #include <petscsnes.h>
  9: #include <petscds.h>
 10: #include <petscviewerhdf5.h>

 12: PetscInt spatialDim = 0;

 14: typedef enum {NEUMANN, DIRICHLET, NONE} BCType;
 15: typedef enum {RUN_FULL, RUN_TEST, RUN_PERF} RunType;
 16: typedef enum {COEFF_NONE, COEFF_ANALYTIC, COEFF_FIELD, COEFF_NONLINEAR} CoeffType;

 18: typedef struct {
 19:   PetscInt      debug;             /* The debugging level */
 20:   RunType       runType;           /* Whether to run tests, or solve the full problem */
 21:   PetscBool     jacobianMF;        /* Whether to calculate the Jacobian action on the fly */
 22:   PetscLogEvent createMeshEvent;
 23:   PetscBool     showInitial, showSolution, restart, check;
 24:   PetscViewer   checkpoint;
 25:   /* Domain and mesh definition */
 26:   PetscInt      dim;               /* The topological mesh dimension */
 27:   char          filename[2048];    /* The optional ExodusII file */
 28:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 29:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 30:   PetscBool     refinementUniform; /* Uniformly refine the mesh */
 31:   PetscInt      refinementRounds;  /* The number of uniform refinements */
 32:   char          partitioner[2048]; /* The graph partitioner */
 33:   /* Problem definition */
 34:   BCType        bcType;
 35:   CoeffType     variableCoefficient;
 36:   void       (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
 37: } AppCtx;

 39: void zero(const PetscReal coords[], PetscScalar *u, void *ctx)
 40: {
 41:   *u = 0.0;
 42: }

 44: /*
 45:   In 2D for Dirichlet conditions, we use exact solution:

 47:     u = x^2 + y^2
 48:     f = 4

 50:   so that

 52:     -\Delta u + f = -4 + 4 = 0

 54:   For Neumann conditions, we have

 56:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (bottom)
 57:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top)
 58:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
 59:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

 61:   Which we can express as

 63:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y} \cdot \hat n = 2 (x + y)
 64: */
 65: void quadratic_u_2d(const PetscReal x[], PetscScalar *u, void *ctx)
 66: {
 67:   *u = x[0]*x[0] + x[1]*x[1];
 68: }

 70: void f0_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
 71: {
 72:   f0[0] = 4.0;
 73: }

 75: void f0_bd_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 76: {
 77:   PetscInt  d;
 78:   for (d = 0, f0[0] = 0.0; d < spatialDim; ++d) f0[0] += -n[d]*2.0*x[d];
 79: }

 81: void f0_bd_zero(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 82: {
 83:   f0[0] = 0.0;
 84: }

 86: void f1_bd_zero(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[])
 87: {
 88:   PetscInt comp;
 89:   for (comp = 0; comp < spatialDim; ++comp) f1[comp] = 0.0;
 90: }

 92: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
 93: void f1_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
 94: {
 95:   PetscInt d;

 97:   for (d = 0; d < spatialDim; ++d) f1[d] = u_x[d];
 98: }

100: /* < \nabla v, \nabla u + {\nabla u}^T >
101:    This just gives \nabla u, give the perdiagonal for the transpose */
102: void g3_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
103: {
104:   PetscInt d;

106:   for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = 1.0;
107: }

109: /*
110:   In 2D for Dirichlet conditions with a variable coefficient, we use exact solution:

112:     u  = x^2 + y^2
113:     f  = 6 (x + y)
114:     nu = (x + y)

116:   so that

118:     -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
119: */
120: void nu_2d(const PetscReal x[], PetscScalar *u, void *ctx)
121: {
122:   *u = x[0] + x[1];
123: }

125: void f0_analytic_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
126: {
127:   f0[0] = 6.0*(x[0] + x[1]);
128: }

130: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
131: void f1_analytic_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
132: {
133:   PetscInt d;

135:   for (d = 0; d < spatialDim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
136: }
137: void f1_field_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
138: {
139:   PetscInt d;

141:   for (d = 0; d < spatialDim; ++d) f1[d] = a[0]*u_x[d];
142: }

144: /* < \nabla v, \nabla u + {\nabla u}^T >
145:    This just gives \nabla u, give the perdiagonal for the transpose */
146: void g3_analytic_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
147: {
148:   PetscInt d;

150:   for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = x[0] + x[1];
151: }
152: void g3_field_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
153: {
154:   PetscInt d;

156:   for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = a[0];
157: }

159: /*
160:   In 2D for Dirichlet conditions with a nonlinear coefficient (p-Laplacian with p = 4), we use exact solution:

162:     u  = x^2 + y^2
163:     f  = 16 (x^2 + y^2)
164:     nu = 1/2 |grad u|^2

166:   so that

168:     -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
169: */
170: void f0_analytic_nonlinear_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
171: {
172:   f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
173: }

175: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
176: void f1_analytic_nonlinear_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
177: {
178:   PetscScalar nu = 0.0;
179:   PetscInt    d;
180:   for (d = 0; d < spatialDim; ++d) nu += u_x[d]*u_x[d];
181:   for (d = 0; d < spatialDim; ++d) f1[d] = 0.5*nu*u_x[d];
182: }

184: /*
185:   grad (u + eps w) - grad u = eps grad w

187:   1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
188: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
189: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
190: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
191: */
192: void g3_analytic_nonlinear_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
193: {
194:   PetscScalar nu = 0.0;
195:   PetscInt    d, e;
196:   for (d = 0; d < spatialDim; ++d) nu += u_x[d]*u_x[d];
197:   for (d = 0; d < spatialDim; ++d) {
198:     g3[d*spatialDim+d] = 0.5*nu;
199:     for (e = 0; e < spatialDim; ++e) {
200:       g3[d*spatialDim+e] += u_x[d]*u_x[e];
201:     }
202:   }
203: }

205: /*
206:   In 3D for Dirichlet conditions we use exact solution:

208:     u = x^2 + y^2 + z^2
209:     f = 6

211:   so that

213:     -\Delta u + f = -6 + 6 = 0

215:   For Neumann conditions, we have

217:     -\nabla u \cdot -\hat z |_{z=0} =  (2z)|_{z=0} =  0 (bottom)
218:     -\nabla u \cdot  \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
219:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (front)
220:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
221:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
222:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

224:   Which we can express as

226:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
227: */
228: void quadratic_u_3d(const PetscReal x[], PetscScalar *u, void *ctx)
229: {
230:   *u = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
231: }

235: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
236: {
237:   const char    *bcTypes[3]  = {"neumann", "dirichlet", "none"};
238:   const char    *runTypes[3] = {"full", "test", "perf"};
239:   const char    *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
240:   PetscInt       bc, run, coeff;
241:   PetscBool      flg;

245:   options->debug               = 0;
246:   options->runType             = RUN_FULL;
247:   options->dim                 = 2;
248:   options->filename[0]         = '\0';
249:   options->interpolate         = PETSC_FALSE;
250:   options->refinementLimit     = 0.0;
251:   options->refinementUniform   = PETSC_FALSE;
252:   options->refinementRounds    = 1;
253:   options->bcType              = DIRICHLET;
254:   options->variableCoefficient = COEFF_NONE;
255:   options->jacobianMF          = PETSC_FALSE;
256:   options->showInitial         = PETSC_FALSE;
257:   options->showSolution        = PETSC_FALSE;
258:   options->restart             = PETSC_FALSE;
259:   options->check               = PETSC_FALSE;
260:   options->checkpoint          = NULL;

262:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
263:   PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, NULL);
264:   run  = options->runType;
265:   PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 3, runTypes[options->runType], &run, NULL);

267:   options->runType = (RunType) run;

269:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
270:   spatialDim = options->dim;
271:   PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
272: #if !defined(PETSC_HAVE_EXODUSII)
273:   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
274: #endif
275:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
276:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
277:   PetscOptionsBool("-refinement_uniform", "Uniformly refine the mesh", "ex52.c", options->refinementUniform, &options->refinementUniform, NULL);
278:   PetscOptionsInt("-refinement_rounds", "The number of uniform refinements", "ex52.c", options->refinementRounds, &options->refinementRounds, NULL);
279:   PetscStrcpy(options->partitioner, "chaco");
280:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
281:   bc   = options->bcType;
282:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
283:   options->bcType = (BCType) bc;
284:   coeff = options->variableCoefficient;
285:   PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
286:   options->variableCoefficient = (CoeffType) coeff;

288:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
289:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
290:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
291:   PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
292:   PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
293:   PetscOptionsEnd();

295:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);

297:   if (options->restart) {
298:     PetscViewerCreate(comm, &options->checkpoint);
299:     PetscViewerSetType(options->checkpoint, PETSCVIEWERHDF5);
300:     PetscViewerFileSetMode(options->checkpoint, FILE_MODE_READ);
301:     PetscViewerFileSetName(options->checkpoint, options->filename);
302:   }
303:   return(0);
304: }

308: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
309: {
310:   PetscInt       dim               = user->dim;
311:   const char    *filename          = user->filename;
312:   PetscBool      interpolate       = user->interpolate;
313:   PetscReal      refinementLimit   = user->refinementLimit;
314:   PetscBool      refinementUniform = user->refinementUniform;
315:   PetscInt       refinementRounds  = user->refinementRounds;
316:   const char    *partitioner       = user->partitioner;
317:   size_t         len;

321:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
322:   PetscStrlen(filename, &len);
323:   if (!len) {
324:     DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
325:     PetscObjectSetName((PetscObject) *dm, "Mesh");
326:   } else if (user->checkpoint) {
327:     DMCreate(comm, dm);
328:     DMSetType(*dm, DMPLEX);
329:     DMLoad(*dm, user->checkpoint);
330:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
331:   } else {
332:     PetscMPIInt rank;

334:     MPI_Comm_rank(comm, &rank);
335:     DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
336:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
337:     /* Must have boundary marker for Dirichlet conditions */
338:   }
339:   {
340:     DM refinedMesh     = NULL;
341:     DM distributedMesh = NULL;

343:     /* Refine mesh using a volume constraint */
344:     DMPlexSetRefinementLimit(*dm, refinementLimit);
345:     DMRefine(*dm, comm, &refinedMesh);
346:     if (refinedMesh) {
347:       const char *name;

349:       PetscObjectGetName((PetscObject) *dm,         &name);
350:       PetscObjectSetName((PetscObject) refinedMesh,  name);
351:       DMDestroy(dm);
352:       *dm  = refinedMesh;
353:     }
354:     /* Distribute mesh over processes */
355:     DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
356:     if (distributedMesh) {
357:       DMDestroy(dm);
358:       *dm  = distributedMesh;
359:     }
360:     /* Use regular refinement in parallel */
361:     if (refinementUniform) {
362:       PetscInt r;

364:       DMPlexSetRefinementUniform(*dm, refinementUniform);
365:       for (r = 0; r < refinementRounds; ++r) {
366:         DMRefine(*dm, comm, &refinedMesh);
367:         if (refinedMesh) {
368:           DMDestroy(dm);
369:           *dm  = refinedMesh;
370:         }
371:       }
372:     }
373:   }
374:   DMSetFromOptions(*dm);
375:   if (user->bcType == NEUMANN) {
376:     DMLabel label;

378:     DMPlexCreateLabel(*dm, "boundary");
379:     DMPlexGetLabel(*dm, "boundary", &label);
380:     DMPlexMarkBoundaryFaces(*dm, label);
381:   }
382:   DMViewFromOptions(*dm, NULL, "-dm_view");
383:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
384:   return(0);
385: }

389: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
390: {
391:   PetscDS        prob;

395:   DMGetDS(dm, &prob);
396:   switch (user->variableCoefficient) {
397:   case COEFF_NONE:
398:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
399:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
400:     break;
401:   case COEFF_ANALYTIC:
402:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
403:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
404:     break;
405:   case COEFF_FIELD:
406:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
407:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
408:     break;
409:   case COEFF_NONLINEAR:
410:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
411:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
412:     break;
413:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
414:   }
415:   switch (user->dim) {
416:   case 2:
417:     user->exactFuncs[0] = quadratic_u_2d;
418:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
419:     break;
420:   case 3:
421:     user->exactFuncs[0] = quadratic_u_3d;
422:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
423:     break;
424:   default:
425:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
426:   }
427:   return(0);
428: }

432: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
433: {
434:   void (*matFuncs[1])(const PetscReal x[], PetscScalar *u, void *ctx) = {nu_2d};
435:   Vec            nu;

439:   DMCreateLocalVector(dmAux, &nu);
440:   DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);
441:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
442:   VecDestroy(&nu);
443:   return(0);
444: }

448: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
449: {
450:   DM             cdm   = dm;
451:   const PetscInt dim   = user->dim;
452:   const PetscInt id    = 1;
453:   PetscFE        feAux = NULL;
454:   PetscFE        feBd  = NULL;
455:   PetscFE        feCh  = NULL;
456:   PetscFE        fe;
457:   PetscDS        prob;

461:   /* Create finite element */
462:   PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, NULL, -1, &fe);
463:   PetscObjectSetName((PetscObject) fe, "potential");
464:   if (user->bcType == NEUMANN) {
465:     PetscFECreateDefault(dm, dim-1, 1, PETSC_TRUE, "bd_", -1, &feBd);
466:     PetscObjectSetName((PetscObject) feBd, "potential");
467:   }
468:   if (user->variableCoefficient == COEFF_FIELD) {
469:     PetscQuadrature q;

471:     PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "mat_", -1, &feAux);
472:     PetscFEGetQuadrature(fe, &q);
473:     PetscFESetQuadrature(feAux, q);
474:   }
475:   if (user->check) {PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "ch_", -1, &feCh);}
476:   /* Set discretization and boundary conditions for each mesh */
477:   while (cdm) {
478:     DMGetDS(cdm, &prob);
479:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
480:     PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
481:     if (feAux) {
482:       DM      dmAux;
483:       PetscDS probAux;

485:       DMClone(cdm, &dmAux);
486:       DMPlexCopyCoordinates(cdm, dmAux);
487:       DMGetDS(dmAux, &probAux);
488:       PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
489:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
490:       SetupMaterial(cdm, dmAux, user);
491:       DMDestroy(&dmAux);
492:     }
493:     if (feCh) {
494:       DM      dmCh;
495:       PetscDS probCh;

497:       DMClone(cdm, &dmCh);
498:       DMPlexCopyCoordinates(cdm, dmCh);
499:       DMGetDS(dmCh, &probCh);
500:       PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
501:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
502:       DMDestroy(&dmCh);
503:     }
504:     SetupProblem(cdm, user);
505:     DMPlexAddBoundary(cdm, user->bcType == DIRICHLET, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1, &id, user);
506:     DMPlexGetCoarseDM(cdm, &cdm);
507:   }
508:   PetscFEDestroy(&fe);
509:   PetscFEDestroy(&feBd);
510:   PetscFEDestroy(&feAux);
511:   PetscFEDestroy(&feCh);
512:   return(0);
513: }

517: int main(int argc, char **argv)
518: {
519:   DM             dm;          /* Problem specification */
520:   SNES           snes;        /* nonlinear solver */
521:   Vec            u,r;         /* solution, residual vectors */
522:   Mat            A,J;         /* Jacobian matrix */
523:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
524:   AppCtx         user;        /* user-defined work context */
525:   JacActionCtx   userJ;       /* context for Jacobian MF action */
526:   PetscInt       its;         /* iterations for convergence */
527:   PetscReal      error = 0.0; /* L_2 error in the solution */

530:   PetscInitialize(&argc, &argv, NULL, help);
531:   ProcessOptions(PETSC_COMM_WORLD, &user);
532:   SNESCreate(PETSC_COMM_WORLD, &snes);
533:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
534:   SNESSetDM(snes, dm);
535:   DMSetApplicationContext(dm, &user);

537:   PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
538:   SetupDiscretization(dm, &user);

540:   DMCreateGlobalVector(dm, &u);
541:   PetscObjectSetName((PetscObject) u, "potential");
542:   VecDuplicate(u, &r);

544:   DMSetMatType(dm,MATAIJ);
545:   DMCreateMatrix(dm, &J);
546:   if (user.jacobianMF) {
547:     PetscInt M, m, N, n;

549:     MatGetSize(J, &M, &N);
550:     MatGetLocalSize(J, &m, &n);
551:     MatCreate(PETSC_COMM_WORLD, &A);
552:     MatSetSizes(A, m, n, M, N);
553:     MatSetType(A, MATSHELL);
554:     MatSetUp(A);
555: #if 0
556:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
557: #endif

559:     userJ.dm   = dm;
560:     userJ.J    = J;
561:     userJ.user = &user;

563:     DMCreateLocalVector(dm, &userJ.u);
564:     DMPlexProjectFunctionLocal(dm, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
565:     MatShellSetContext(A, &userJ);
566:   } else {
567:     A = J;
568:   }
569:   if (user.bcType == NEUMANN) {
570:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
571:     MatSetNullSpace(J, nullSpace);
572:     if (A != J) {
573:       MatSetNullSpace(A, nullSpace);
574:     }
575:   }

577:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);
578:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);
579:   SNESSetJacobian(snes, A, J, NULL, NULL);

581:   SNESSetFromOptions(snes);

583:   DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
584:   if (user.checkpoint) {
585: #if defined(PETSC_HAVE_HDF5)
586:     PetscViewerHDF5PushGroup(user.checkpoint, "/fields");
587:     VecLoad(u, user.checkpoint);
588:     PetscViewerHDF5PopGroup(user.checkpoint);
589: #endif
590:   }
591:   PetscViewerDestroy(&user.checkpoint);
592:   if (user.showInitial) {
593:     Vec lv;
594:     DMGetLocalVector(dm, &lv);
595:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
596:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
597:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
598:     DMRestoreLocalVector(dm, &lv);
599:   }
600:   if (user.runType == RUN_FULL) {
601:     void (*initialGuess[1])(const PetscReal x[], PetscScalar *, void *ctx) = {zero};

603:     DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
604:     if (user.debug) {
605:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
606:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
607:     }
608:     SNESSolve(snes, NULL, u);
609:     SNESGetIterationNumber(snes, &its);
610:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
611:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
612:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
613:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
614:     if (user.showSolution) {
615:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
616:       VecChop(u, 3.0e-9);
617:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
618:     }
619:   } else if (user.runType == RUN_PERF) {
620:     PetscReal res = 0.0;

622:     SNESComputeFunction(snes, u, r);
623:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
624:     VecChop(r, 1.0e-10);
625:     VecNorm(r, NORM_2, &res);
626:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
627:   } else {
628:     PetscReal res = 0.0;

630:     /* Check discretization error */
631:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
632:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
633:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
634:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
635:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
636:     /* Check residual */
637:     SNESComputeFunction(snes, u, r);
638:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
639:     VecChop(r, 1.0e-10);
640:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
641:     VecNorm(r, NORM_2, &res);
642:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
643:     /* Check Jacobian */
644:     {
645:       Vec          b;

647:       SNESComputeJacobian(snes, u, A, A);
648:       VecDuplicate(u, &b);
649:       VecSet(r, 0.0);
650:       SNESComputeFunction(snes, r, b);
651:       MatMult(A, u, r);
652:       VecAXPY(r, 1.0, b);
653:       VecDestroy(&b);
654:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
655:       VecChop(r, 1.0e-10);
656:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
657:       VecNorm(r, NORM_2, &res);
658:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
659:     }
660:   }
661:   VecViewFromOptions(u, NULL, "-vec_view");

663:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
664:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
665:   if (A != J) {MatDestroy(&A);}
666:   MatDestroy(&J);
667:   VecDestroy(&u);
668:   VecDestroy(&r);
669:   SNESDestroy(&snes);
670:   DMDestroy(&dm);
671:   PetscFree(user.exactFuncs);
672:   PetscFinalize();
673:   return 0;
674: }