Actual source code: ex12.c

petsc-master 2014-11-28
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 nonlinear 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:   /* Problem definition */
 31:   BCType        bcType;
 32:   CoeffType     variableCoefficient;
 33:   void       (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
 34: } AppCtx;

 36: void zero(const PetscReal coords[], PetscScalar *u, void *ctx)
 37: {
 38:   *u = 0.0;
 39: }

 41: /*
 42:   In 2D for Dirichlet conditions, we use exact solution:

 44:     u = x^2 + y^2
 45:     f = 4

 47:   so that

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

 51:   For Neumann conditions, we have

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

 58:   Which we can express as

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

 67: 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[])
 68: {
 69:   f0[0] = 4.0;
 70: }

 72: 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[])
 73: {
 74:   PetscInt  d;
 75:   for (d = 0, f0[0] = 0.0; d < spatialDim; ++d) f0[0] += -n[d]*2.0*x[d];
 76: }

 78: 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[])
 79: {
 80:   f0[0] = 0.0;
 81: }

 83: 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[])
 84: {
 85:   PetscInt comp;
 86:   for (comp = 0; comp < spatialDim; ++comp) f1[comp] = 0.0;
 87: }

 89: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
 90: 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[])
 91: {
 92:   PetscInt d;

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

 97: /* < \nabla v, \nabla u + {\nabla u}^T >
 98:    This just gives \nabla u, give the perdiagonal for the transpose */
 99: 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[])
100: {
101:   PetscInt d;

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

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

109:     u  = x^2 + y^2
110:     f  = 6 (x + y)
111:     nu = (x + y)

113:   so that

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

122: 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[])
123: {
124:   f0[0] = 6.0*(x[0] + x[1]);
125: }

127: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
128: 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[])
129: {
130:   PetscInt d;

132:   for (d = 0; d < spatialDim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
133: }
134: 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[])
135: {
136:   PetscInt d;

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

141: /* < \nabla v, \nabla u + {\nabla u}^T >
142:    This just gives \nabla u, give the perdiagonal for the transpose */
143: 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[])
144: {
145:   PetscInt d;

147:   for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = x[0] + x[1];
148: }
149: 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[])
150: {
151:   PetscInt d;

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

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

159:     u  = x^2 + y^2
160:     f  = 16 (x^2 + y^2)
161:     nu = 1/2 |grad u|^2

163:   so that

165:     -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
166: */
167: 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[])
168: {
169:   f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
170: }

172: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
173: 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[])
174: {
175:   PetscScalar nu = 0.0;
176:   PetscInt    d;
177:   for (d = 0; d < spatialDim; ++d) nu += u_x[d]*u_x[d];
178:   for (d = 0; d < spatialDim; ++d) f1[d] = 0.5*nu*u_x[d];
179: }

181: /*
182:   grad (u + eps w) - grad u = eps grad w

184:   1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
185: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
186: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
187: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
188: */
189: 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[])
190: {
191:   PetscScalar nu = 0.0;
192:   PetscInt    d, e;
193:   for (d = 0; d < spatialDim; ++d) nu += u_x[d]*u_x[d];
194:   for (d = 0; d < spatialDim; ++d) {
195:     g3[d*spatialDim+d] = 0.5*nu;
196:     for (e = 0; e < spatialDim; ++e) {
197:       g3[d*spatialDim+e] += u_x[d]*u_x[e];
198:     }
199:   }
200: }

202: /*
203:   In 3D for Dirichlet conditions we use exact solution:

205:     u = x^2 + y^2 + z^2
206:     f = 6

208:   so that

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

212:   For Neumann conditions, we have

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

221:   Which we can express as

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

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

242:   options->debug               = 0;
243:   options->runType             = RUN_FULL;
244:   options->dim                 = 2;
245:   options->filename[0]         = '\0';
246:   options->interpolate         = PETSC_FALSE;
247:   options->refinementLimit     = 0.0;
248:   options->bcType              = DIRICHLET;
249:   options->variableCoefficient = COEFF_NONE;
250:   options->jacobianMF          = PETSC_FALSE;
251:   options->showInitial         = PETSC_FALSE;
252:   options->showSolution        = PETSC_FALSE;
253:   options->restart             = PETSC_FALSE;
254:   options->check               = PETSC_FALSE;
255:   options->checkpoint          = NULL;

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

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

264:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
265:   spatialDim = options->dim;
266:   PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
267: #if !defined(PETSC_HAVE_EXODUSII)
268:   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
269: #endif
270:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
271:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
272:   bc   = options->bcType;
273:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
274:   options->bcType = (BCType) bc;
275:   coeff = options->variableCoefficient;
276:   PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
277:   options->variableCoefficient = (CoeffType) coeff;

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

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

288:   if (options->restart) {
289:     PetscViewerCreate(comm, &options->checkpoint);
290:     PetscViewerSetType(options->checkpoint, PETSCVIEWERHDF5);
291:     PetscViewerFileSetMode(options->checkpoint, FILE_MODE_READ);
292:     PetscViewerFileSetName(options->checkpoint, options->filename);
293:   }
294:   return(0);
295: }

299: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
300: {
301:   PetscInt       dim             = user->dim;
302:   const char    *filename        = user->filename;
303:   PetscBool      interpolate     = user->interpolate;
304:   PetscReal      refinementLimit = user->refinementLimit;
305:   size_t         len;

309:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
310:   PetscStrlen(filename, &len);
311:   if (!len) {
312:     DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
313:     PetscObjectSetName((PetscObject) *dm, "Mesh");
314:   } else if (user->checkpoint) {
315:     DMCreate(comm, dm);
316:     DMSetType(*dm, DMPLEX);
317:     DMLoad(*dm, user->checkpoint);
318:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
319:   } else {
320:     PetscMPIInt rank;

322:     MPI_Comm_rank(comm, &rank);
323:     DMPlexCreateFromFile(comm, filename, interpolate, dm);
324:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
325:     /* Must have boundary marker for Dirichlet conditions */
326:   }
327:   {
328:     DM refinedMesh     = NULL;
329:     DM distributedMesh = NULL;

331:     /* Refine mesh using a volume constraint */
332:     DMPlexSetRefinementLimit(*dm, refinementLimit);
333:     DMRefine(*dm, comm, &refinedMesh);
334:     if (refinedMesh) {
335:       const char *name;

337:       PetscObjectGetName((PetscObject) *dm,         &name);
338:       PetscObjectSetName((PetscObject) refinedMesh,  name);
339:       DMDestroy(dm);
340:       *dm  = refinedMesh;
341:     }
342:     /* Distribute mesh over processes */
343:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
344:     if (distributedMesh) {
345:       DMDestroy(dm);
346:       *dm  = distributedMesh;
347:     } else {
348:       DMSetFromOptions(*dm);
349:     }
350:   }
351:   if (user->bcType == NEUMANN) {
352:     DMLabel label;

354:     DMPlexCreateLabel(*dm, "boundary");
355:     DMPlexGetLabel(*dm, "boundary", &label);
356:     DMPlexMarkBoundaryFaces(*dm, label);
357:   }
358:   DMViewFromOptions(*dm, NULL, "-dm_view");
359:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
360:   return(0);
361: }

365: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
366: {
367:   PetscDS        prob;

371:   DMGetDS(dm, &prob);
372:   switch (user->variableCoefficient) {
373:   case COEFF_NONE:
374:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
375:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
376:     break;
377:   case COEFF_ANALYTIC:
378:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
379:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
380:     break;
381:   case COEFF_FIELD:
382:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
383:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
384:     break;
385:   case COEFF_NONLINEAR:
386:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
387:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
388:     break;
389:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
390:   }
391:   switch (user->dim) {
392:   case 2:
393:     user->exactFuncs[0] = quadratic_u_2d;
394:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
395:     break;
396:   case 3:
397:     user->exactFuncs[0] = quadratic_u_3d;
398:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
399:     break;
400:   default:
401:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
402:   }
403:   return(0);
404: }

408: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
409: {
410:   void (*matFuncs[1])(const PetscReal x[], PetscScalar *u, void *ctx) = {nu_2d};
411:   Vec            nu;

415:   DMCreateLocalVector(dmAux, &nu);
416:   DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);
417:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
418:   VecDestroy(&nu);
419:   return(0);
420: }

424: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
425: {
426:   DM             cdm   = dm;
427:   const PetscInt dim   = user->dim;
428:   const PetscInt id    = 1;
429:   PetscFE        feAux = NULL;
430:   PetscFE        feBd  = NULL;
431:   PetscFE        feCh  = NULL;
432:   PetscFE        fe;
433:   PetscDS        prob;

437:   /* Create finite element */
438:   PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, NULL, -1, &fe);
439:   PetscObjectSetName((PetscObject) fe, "potential");
440:   if (user->bcType == NEUMANN) {
441:     PetscFECreateDefault(dm, dim-1, 1, PETSC_TRUE, "bd_", -1, &feBd);
442:     PetscObjectSetName((PetscObject) feBd, "potential");
443:   }
444:   if (user->variableCoefficient == COEFF_FIELD) {
445:     PetscQuadrature q;

447:     PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "mat_", -1, &feAux);
448:     PetscFEGetQuadrature(fe, &q);
449:     PetscFESetQuadrature(feAux, q);
450:   }
451:   if (user->check) {PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "ch_", -1, &feCh);}
452:   /* Set discretization and boundary conditions for each mesh */
453:   while (cdm) {
454:     DMGetDS(cdm, &prob);
455:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
456:     PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
457:     if (feAux) {
458:       DM      dmAux;
459:       PetscDS probAux;

461:       DMClone(cdm, &dmAux);
462:       DMPlexCopyCoordinates(cdm, dmAux);
463:       DMGetDS(dmAux, &probAux);
464:       PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
465:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
466:       SetupMaterial(cdm, dmAux, user);
467:       DMDestroy(&dmAux);
468:     }
469:     if (feCh) {
470:       DM      dmCh;
471:       PetscDS probCh;

473:       DMClone(cdm, &dmCh);
474:       DMPlexCopyCoordinates(cdm, dmCh);
475:       DMGetDS(dmCh, &probCh);
476:       PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
477:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
478:       DMDestroy(&dmCh);
479:     }
480:     SetupProblem(cdm, user);
481:     DMPlexAddBoundary(cdm, user->bcType == DIRICHLET, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1, &id, user);
482:     DMPlexGetCoarseDM(cdm, &cdm);
483:   }
484:   PetscFEDestroy(&fe);
485:   PetscFEDestroy(&feBd);
486:   PetscFEDestroy(&feAux);
487:   PetscFEDestroy(&feCh);
488:   return(0);
489: }

493: int main(int argc, char **argv)
494: {
495:   DM             dm;          /* Problem specification */
496:   SNES           snes;        /* nonlinear solver */
497:   Vec            u,r;         /* solution, residual vectors */
498:   Mat            A,J;         /* Jacobian matrix */
499:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
500:   AppCtx         user;        /* user-defined work context */
501:   JacActionCtx   userJ;       /* context for Jacobian MF action */
502:   PetscInt       its;         /* iterations for convergence */
503:   PetscReal      error = 0.0; /* L_2 error in the solution */

506:   PetscInitialize(&argc, &argv, NULL, help);
507:   ProcessOptions(PETSC_COMM_WORLD, &user);
508:   SNESCreate(PETSC_COMM_WORLD, &snes);
509:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
510:   SNESSetDM(snes, dm);
511:   DMSetApplicationContext(dm, &user);

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

516:   DMCreateGlobalVector(dm, &u);
517:   PetscObjectSetName((PetscObject) u, "potential");
518:   VecDuplicate(u, &r);

520:   DMSetMatType(dm,MATAIJ);
521:   DMCreateMatrix(dm, &J);
522:   if (user.jacobianMF) {
523:     PetscInt M, m, N, n;

525:     MatGetSize(J, &M, &N);
526:     MatGetLocalSize(J, &m, &n);
527:     MatCreate(PETSC_COMM_WORLD, &A);
528:     MatSetSizes(A, m, n, M, N);
529:     MatSetType(A, MATSHELL);
530:     MatSetUp(A);
531: #if 0
532:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
533: #endif

535:     userJ.dm   = dm;
536:     userJ.J    = J;
537:     userJ.user = &user;

539:     DMCreateLocalVector(dm, &userJ.u);
540:     DMPlexProjectFunctionLocal(dm, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
541:     MatShellSetContext(A, &userJ);
542:   } else {
543:     A = J;
544:   }
545:   if (user.bcType == NEUMANN) {
546:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
547:     MatSetNullSpace(J, nullSpace);
548:     if (A != J) {
549:       MatSetNullSpace(A, nullSpace);
550:     }
551:   }

553:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);
554:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);
555:   SNESSetJacobian(snes, A, J, NULL, NULL);

557:   SNESSetFromOptions(snes);

559:   DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
560:   if (user.checkpoint) {
561: #if defined(PETSC_HAVE_HDF5)
562:     PetscViewerHDF5PushGroup(user.checkpoint, "/fields");
563:     VecLoad(u, user.checkpoint);
564:     PetscViewerHDF5PopGroup(user.checkpoint);
565: #endif
566:   }
567:   PetscViewerDestroy(&user.checkpoint);
568:   if (user.showInitial) {
569:     Vec lv;
570:     DMGetLocalVector(dm, &lv);
571:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
572:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
573:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
574:     DMRestoreLocalVector(dm, &lv);
575:   }
576:   if (user.runType == RUN_FULL) {
577:     void (*initialGuess[1])(const PetscReal x[], PetscScalar *, void *ctx) = {zero};

579:     DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
580:     if (user.debug) {
581:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
582:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
583:     }
584:     SNESSolve(snes, NULL, u);
585:     SNESGetIterationNumber(snes, &its);
586:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
587:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
588:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
589:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
590:     if (user.showSolution) {
591:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
592:       VecChop(u, 3.0e-9);
593:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
594:     }
595:   } else if (user.runType == RUN_PERF) {
596:     PetscReal res = 0.0;

598:     SNESComputeFunction(snes, u, r);
599:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
600:     VecChop(r, 1.0e-10);
601:     VecNorm(r, NORM_2, &res);
602:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
603:   } else {
604:     PetscReal res = 0.0;

606:     /* Check discretization error */
607:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
608:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
609:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
610:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
611:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
612:     /* Check residual */
613:     SNESComputeFunction(snes, u, r);
614:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
615:     VecChop(r, 1.0e-10);
616:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
617:     VecNorm(r, NORM_2, &res);
618:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
619:     /* Check Jacobian */
620:     {
621:       Vec          b;

623:       SNESComputeJacobian(snes, u, A, A);
624:       VecDuplicate(u, &b);
625:       VecSet(r, 0.0);
626:       SNESComputeFunction(snes, r, b);
627:       MatMult(A, u, r);
628:       VecAXPY(r, 1.0, b);
629:       VecDestroy(&b);
630:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
631:       VecChop(r, 1.0e-10);
632:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
633:       VecNorm(r, NORM_2, &res);
634:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
635:     }
636:   }
637:   VecViewFromOptions(u, NULL, "-vec_view");

639:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
640:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
641:   if (A != J) {MatDestroy(&A);}
642:   MatDestroy(&J);
643:   VecDestroy(&u);
644:   VecDestroy(&r);
645:   SNESDestroy(&snes);
646:   DMDestroy(&dm);
647:   PetscFree(user.exactFuncs);
648:   PetscFinalize();
649:   return 0;
650: }