Actual source code: ex12.c

petsc-dev 2014-08-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 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:   PetscBool     simplex;           /* Use simplices or tensor product cells */
 30:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 31:   char          partitioner[2048]; /* The graph partitioner */
 32:   /* Problem definition */
 33:   BCType        bcType;
 34:   CoeffType     variableCoefficient;
 35:   void       (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
 36: } AppCtx;

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

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

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

 49:   so that

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

 53:   For Neumann conditions, we have

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

 60:   Which we can express as

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

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

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

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

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

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

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

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

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

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

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

115:   so that

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

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

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

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

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

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

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

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

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

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

165:   so that

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

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

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

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

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

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

210:   so that

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

214:   For Neumann conditions, we have

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

223:   Which we can express as

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

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

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

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

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

267:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
268:   spatialDim = options->dim;
269:   PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
270: #if !defined(PETSC_HAVE_EXODUSII)
271:   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
272: #endif
273:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
274:   PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex12.c", options->simplex, &options->simplex, NULL);
275:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
276:   if (!options->simplex && options->refinementLimit > 0.0) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot refine tensor cells based on a volume criterion");
277:   PetscStrcpy(options->partitioner, "chaco");
278:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
279:   bc   = options->bcType;
280:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
281:   options->bcType = (BCType) bc;
282:   coeff = options->variableCoefficient;
283:   PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
284:   options->variableCoefficient = (CoeffType) coeff;

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

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

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

306: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
307: {
308:   DM             refinedMesh     = NULL;
309:   DM             distributedMesh = NULL;
310:   PetscInt       dim             = user->dim;
311:   const char    *filename        = user->filename;
312:   PetscBool      interpolate     = user->interpolate;
313:   PetscReal      refinementLimit = user->refinementLimit;
314:   const char    *partitioner     = user->partitioner;
315:   const PetscInt cells[3]        = {4, 4, 4};
316:   size_t         len;

320:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
321:   PetscStrlen(filename, &len);
322:   if (!len) {
323:     if (user->simplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
324:     else               {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 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:   /* Refine mesh using a volume constraint */
340:   if (refinementLimit > 0.0) {
341:     DMPlexSetRefinementLimit(*dm, refinementLimit);
342:     DMRefine(*dm, comm, &refinedMesh);
343:   }
344:   if (refinedMesh) {
345:     DMDestroy(dm);
346:     *dm  = refinedMesh;
347:   }
348:   /* Distribute mesh over processes */
349:   DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
350:   if (distributedMesh) {
351:     DMDestroy(dm);
352:     *dm  = distributedMesh;
353:   }
354:   DMSetFromOptions(*dm);
355:   if (user->bcType == NEUMANN) {
356:     DMLabel label;

358:     DMPlexCreateLabel(*dm, "boundary");
359:     DMPlexGetLabel(*dm, "boundary", &label);
360:     DMPlexMarkBoundaryFaces(*dm, label);
361:   }
362:   DMViewFromOptions(*dm, NULL, "-dm_view");
363:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
364:   return(0);
365: }

369: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
370: {
371:   PetscDS        prob;

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

412: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
413: {
414:   void (*matFuncs[1])(const PetscReal x[], PetscScalar *u, void *ctx) = {nu_2d};
415:   Vec            nu;

419:   DMCreateLocalVector(dmAux, &nu);
420:   DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);
421:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
422:   VecDestroy(&nu);
423:   return(0);
424: }

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

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

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

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

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

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

510:   PetscInitialize(&argc, &argv, NULL, help);
511:   ProcessOptions(PETSC_COMM_WORLD, &user);
512:   SNESCreate(PETSC_COMM_WORLD, &snes);
513:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
514:   SNESSetDM(snes, dm);
515:   DMSetApplicationContext(dm, &user);

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

520:   DMCreateGlobalVector(dm, &u);
521:   PetscObjectSetName((PetscObject) u, "potential");
522:   VecDuplicate(u, &r);

524:   DMSetMatType(dm,MATAIJ);
525:   DMCreateMatrix(dm, &J);
526:   if (user.jacobianMF) {
527:     PetscInt M, m, N, n;

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

539:     userJ.dm   = dm;
540:     userJ.J    = J;
541:     userJ.user = &user;

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

557:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);
558:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);
559:   SNESSetJacobian(snes, A, J, NULL, NULL);

561:   SNESSetFromOptions(snes);

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

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

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

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

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

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