Actual source code: ex12.c

petsc-3.7.7 2017-09-25
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: typedef enum {NEUMANN, DIRICHLET, NONE} BCType;
 13: typedef enum {RUN_FULL, RUN_EXACT, RUN_TEST, RUN_PERF} RunType;
 14: typedef enum {COEFF_NONE, COEFF_ANALYTIC, COEFF_FIELD, COEFF_NONLINEAR} CoeffType;

 16: typedef struct {
 17:   PetscInt       debug;             /* The debugging level */
 18:   RunType        runType;           /* Whether to run tests, or solve the full problem */
 19:   PetscBool      jacobianMF;        /* Whether to calculate the Jacobian action on the fly */
 20:   PetscLogEvent  createMeshEvent;
 21:   PetscBool      showInitial, showSolution, restart, check;
 22:   /* Domain and mesh definition */
 23:   PetscInt       dim;               /* The topological mesh dimension */
 24:   char           filename[2048];    /* The optional ExodusII file */
 25:   PetscBool      interpolate;       /* Generate intermediate mesh elements */
 26:   PetscReal      refinementLimit;   /* The largest allowable cell volume */
 27:   PetscBool      viewHierarchy;     /* Whether to view the hierarchy */
 28:   PetscBool      simplex;           /* Simplicial mesh */
 29:   /* Problem definition */
 30:   BCType         bcType;
 31:   CoeffType      variableCoefficient;
 32:   PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 33:   /* Solver */
 34:   PC            pcmg;              /* This is needed for error monitoring */
 35: } AppCtx;

 37: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 38: {
 39:   u[0] = 0.0;
 40:   return 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: static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 65: {
 66:   *u = x[0]*x[0] + x[1]*x[1];
 67:   return 0;
 68: }

 70: void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 71:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 72:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 73:           PetscReal t, const PetscReal x[], PetscScalar f0[])
 74: {
 75:   f0[0] = 4.0;
 76: }

 78: void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 79:              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 80:              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 81:              PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 82: {
 83:   PetscInt d;
 84:   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += -n[d]*2.0*x[d];
 85: }

 87: void f0_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 88:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 89:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 90:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 91: {
 92:   f0[0] = 0.0;
 93: }

 95: void f1_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 96:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 97:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 98:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
 99: {
100:   PetscInt comp;
101:   for (comp = 0; comp < dim; ++comp) f1[comp] = 0.0;
102: }

104: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
105: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
106:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
107:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
108:           PetscReal t, const PetscReal x[], PetscScalar f1[])
109: {
110:   PetscInt d;
111:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
112: }

114: /* < \nabla v, \nabla u + {\nabla u}^T >
115:    This just gives \nabla u, give the perdiagonal for the transpose */
116: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
117:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
118:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
119:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
120: {
121:   PetscInt d;
122:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
123: }

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

128:     u  = x^2 + y^2
129:     f  = 6 (x + y)
130:     nu = (x + y)

132:   so that

134:     -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
135: */
136: static PetscErrorCode nu_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
137: {
138:   *u = x[0] + x[1];
139:   return 0;
140: }

142: void f0_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
143:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
144:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
145:                    PetscReal t, const PetscReal x[], PetscScalar f0[])
146: {
147:   f0[0] = 6.0*(x[0] + x[1]);
148: }

150: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
151: void f1_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
152:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
153:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
154:                    PetscReal t, const PetscReal x[], PetscScalar f1[])
155: {
156:   PetscInt d;
157:   for (d = 0; d < dim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
158: }

160: void f1_field_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
161:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
162:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
163:                 PetscReal t, const PetscReal x[], PetscScalar f1[])
164: {
165:   PetscInt d;
166:   for (d = 0; d < dim; ++d) f1[d] = a[0]*u_x[d];
167: }

169: /* < \nabla v, \nabla u + {\nabla u}^T >
170:    This just gives \nabla u, give the perdiagonal for the transpose */
171: void g3_analytic_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
172:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
173:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
174:                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
175: {
176:   PetscInt d;
177:   for (d = 0; d < dim; ++d) g3[d*dim+d] = x[0] + x[1];
178: }

180: void g3_field_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
181:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
182:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
183:                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
184: {
185:   PetscInt d;
186:   for (d = 0; d < dim; ++d) g3[d*dim+d] = a[0];
187: }

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

192:     u  = x^2 + y^2
193:     f  = 16 (x^2 + y^2)
194:     nu = 1/2 |grad u|^2

196:   so that

198:     -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
199: */
200: void f0_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
201:                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
202:                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
203:                              PetscReal t, const PetscReal x[], PetscScalar f0[])
204: {
205:   f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
206: }

208: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
209: void f1_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
210:                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
211:                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
212:                              PetscReal t, const PetscReal x[], PetscScalar f1[])
213: {
214:   PetscScalar nu = 0.0;
215:   PetscInt    d;
216:   for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
217:   for (d = 0; d < dim; ++d) f1[d] = 0.5*nu*u_x[d];
218: }

220: /*
221:   grad (u + eps w) - grad u = eps grad w

223:   1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
224: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
225: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
226: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
227: */
228: void g3_analytic_nonlinear_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
229:                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
230:                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
231:                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
232: {
233:   PetscScalar nu = 0.0;
234:   PetscInt    d, e;
235:   for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
236:   for (d = 0; d < dim; ++d) {
237:     g3[d*dim+d] = 0.5*nu;
238:     for (e = 0; e < dim; ++e) {
239:       g3[d*dim+e] += u_x[d]*u_x[e];
240:     }
241:   }
242: }

244: /*
245:   In 3D for Dirichlet conditions we use exact solution:

247:     u = 2/3 (x^2 + y^2 + z^2)
248:     f = 4

250:   so that

252:     -\Delta u + f = -2/3 * 6 + 4 = 0

254:   For Neumann conditions, we have

256:     -\nabla u \cdot -\hat z |_{z=0} =  (2z)|_{z=0} =  0 (bottom)
257:     -\nabla u \cdot  \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
258:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (front)
259:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
260:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
261:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

263:   Which we can express as

265:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
266: */
267: static PetscErrorCode quadratic_u_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
268: {
269:   *u = 2.0*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/3.0;
270:   return 0;
271: }

275: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
276: {
277:   const char    *bcTypes[3]  = {"neumann", "dirichlet", "none"};
278:   const char    *runTypes[4] = {"full", "exact", "test", "perf"};
279:   const char    *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
280:   PetscInt       bc, run, coeff;
281:   PetscBool      flg;

285:   options->debug               = 0;
286:   options->runType             = RUN_FULL;
287:   options->dim                 = 2;
288:   options->filename[0]         = '\0';
289:   options->interpolate         = PETSC_FALSE;
290:   options->refinementLimit     = 0.0;
291:   options->bcType              = DIRICHLET;
292:   options->variableCoefficient = COEFF_NONE;
293:   options->jacobianMF          = PETSC_FALSE;
294:   options->showInitial         = PETSC_FALSE;
295:   options->showSolution        = PETSC_FALSE;
296:   options->restart             = PETSC_FALSE;
297:   options->check               = PETSC_FALSE;
298:   options->viewHierarchy       = PETSC_FALSE;
299:   options->simplex             = PETSC_TRUE;

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

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

308:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
309:   PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
310: #if !defined(PETSC_HAVE_EXODUSII)
311:   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
312: #endif
313:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
314:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
315:   bc   = options->bcType;
316:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
317:   options->bcType = (BCType) bc;
318:   coeff = options->variableCoefficient;
319:   PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
320:   options->variableCoefficient = (CoeffType) coeff;

322:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
323:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
324:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
325:   PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
326:   PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
327:   PetscOptionsBool("-dm_view_hierarchy", "View the coarsened hierarchy", "ex12.c", options->viewHierarchy, &options->viewHierarchy, NULL);
328:   PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex12.c", options->simplex, &options->simplex, NULL);
329:   PetscOptionsEnd();
330:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
331:   return(0);
332: }

336: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
337: {
338:   DMLabel        label;

342:   DMCreateLabel(dm, name);
343:   DMGetLabel(dm, name, &label);
344:   DMPlexMarkBoundaryFaces(dm, label);
345:   DMPlexLabelComplete(dm, label);
346:   return(0);
347: }

351: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
352: {
353:   PetscInt       dim             = user->dim;
354:   const char    *filename        = user->filename;
355:   PetscBool      interpolate     = user->interpolate;
356:   PetscReal      refinementLimit = user->refinementLimit;
357:   size_t         len;

361:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
362:   PetscStrlen(filename, &len);
363:   if (!len) {
364:     if (user->simplex) {
365:       DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
366:     }
367:     else {
368:       PetscInt cells[3] = {1, 1, 1}; /* coarse mesh is one cell; refine from there */

370:       DMPlexCreateHexBoxMesh(comm, dim, cells,  DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
371:     }
372:     PetscObjectSetName((PetscObject) *dm, "Mesh");
373:   } else {
374:     DMPlexCreateFromFile(comm, filename, interpolate, dm);
375:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
376:   }
377:   {
378:     DM refinedMesh     = NULL;
379:     DM distributedMesh = NULL;

381:     /* Refine mesh using a volume constraint */
382:     DMPlexSetRefinementLimit(*dm, refinementLimit);
383:     DMRefine(*dm, comm, &refinedMesh);
384:     if (refinedMesh) {
385:       const char *name;

387:       PetscObjectGetName((PetscObject) *dm,         &name);
388:       PetscObjectSetName((PetscObject) refinedMesh,  name);
389:       DMDestroy(dm);
390:       *dm  = refinedMesh;
391:     }
392:     /* Distribute mesh over processes */
393:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
394:     if (distributedMesh) {
395:       DMDestroy(dm);
396:       *dm  = distributedMesh;
397:     }
398:   }
399:   if (user->bcType == NEUMANN) {
400:     DMLabel   label;

402:     DMCreateLabel(*dm, "boundary");
403:     DMGetLabel(*dm, "boundary", &label);
404:     DMPlexMarkBoundaryFaces(*dm, label);
405:   } else if (user->bcType == DIRICHLET) {
406:     PetscBool hasLabel;

408:     DMHasLabel(*dm,"marker",&hasLabel);
409:     if (!hasLabel) {CreateBCLabel(*dm, "marker");}
410:   }
411:   {
412:     char      convType[256];
413:     PetscBool flg;

415:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
416:     PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);
417:     PetscOptionsEnd();
418:     if (flg) {
419:       DM dmConv;

421:       DMConvert(*dm,convType,&dmConv);
422:       if (dmConv) {
423:         DMDestroy(dm);
424:         *dm  = dmConv;
425:       }
426:     }
427:   }
428:   DMSetFromOptions(*dm);
429:   DMViewFromOptions(*dm, NULL, "-dm_view");
430:   if (user->viewHierarchy) {
431:     DM       cdm = *dm;
432:     PetscInt i   = 0;
433:     char     buf[256];

435:     while (cdm) {
436:       DMSetUp(cdm);
437:       DMGetCoarseDM(cdm, &cdm);
438:       ++i;
439:     }
440:     cdm = *dm;
441:     while (cdm) {
442:       PetscViewer       viewer;
443:       PetscBool   isHDF5, isVTK;

445:       --i;
446:       PetscViewerCreate(comm,&viewer);
447:       PetscViewerSetType(viewer,PETSCVIEWERHDF5);
448:       PetscViewerSetOptionsPrefix(viewer,"hierarchy_");
449:       PetscViewerSetFromOptions(viewer);
450:       PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
451:       PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
452:       if (isHDF5) {
453:         PetscSNPrintf(buf, 256, "ex12-%d.h5", i);
454:       }
455:       else if (isVTK) {
456:         PetscSNPrintf(buf, 256, "ex12-%d.vtu", i);
457:         PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
458:       }
459:       else {
460:         PetscSNPrintf(buf, 256, "ex12-%d", i);
461:       }
462:       PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
463:       PetscViewerFileSetName(viewer,buf);
464:       DMView(cdm, viewer);
465:       PetscViewerDestroy(&viewer);
466:       DMGetCoarseDM(cdm, &cdm);
467:     }
468:   }
469:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
470:   return(0);
471: }

475: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
476: {
477:   PetscDS        prob;

481:   DMGetDS(dm, &prob);
482:   switch (user->variableCoefficient) {
483:   case COEFF_NONE:
484:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
485:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
486:     break;
487:   case COEFF_ANALYTIC:
488:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
489:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
490:     break;
491:   case COEFF_FIELD:
492:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
493:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
494:     break;
495:   case COEFF_NONLINEAR:
496:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
497:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
498:     break;
499:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
500:   }
501:   switch (user->dim) {
502:   case 2:
503:     user->exactFuncs[0] = quadratic_u_2d;
504:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
505:     break;
506:   case 3:
507:     user->exactFuncs[0] = quadratic_u_3d;
508:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
509:     break;
510:   default:
511:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
512:   }
513:   return(0);
514: }

518: static PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
519: {
520:   PetscErrorCode (*matFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {nu_2d};
521:   Vec            nu;

525:   DMCreateLocalVector(dmAux, &nu);
526:   DMProjectFunctionLocal(dmAux, 0.0, matFuncs, NULL, INSERT_ALL_VALUES, nu);
527:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
528:   VecDestroy(&nu);
529:   return(0);
530: }

534: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
535: {
536:   DM             cdm   = dm;
537:   const PetscInt dim   = user->dim;
538:   const PetscInt id    = 1;
539:   PetscFE        feAux = NULL;
540:   PetscFE        feBd  = NULL;
541:   PetscFE        feCh  = NULL;
542:   PetscFE        fe;
543:   PetscDS        prob;
544:   PetscBool      simplex = user->simplex;

548:   /* Create finite element */
549:   PetscFECreateDefault(dm, dim, 1, simplex, NULL, -1, &fe);
550:   PetscObjectSetName((PetscObject) fe, "potential");
551:   if (user->bcType == NEUMANN) {
552:     PetscFECreateDefault(dm, dim-1, 1, simplex, "bd_", -1, &feBd);
553:     PetscObjectSetName((PetscObject) feBd, "potential");
554:   }
555:   if (user->variableCoefficient == COEFF_FIELD) {
556:     PetscQuadrature q;

558:     PetscFECreateDefault(dm, dim, 1, simplex, "mat_", -1, &feAux);
559:     PetscFEGetQuadrature(fe, &q);
560:     PetscFESetQuadrature(feAux, q);
561:   }
562:   if (user->check) {PetscFECreateDefault(dm, dim, 1, simplex, "ch_", -1, &feCh);}
563:   /* Set discretization and boundary conditions for each mesh */
564:   while (cdm) {
565:     DM coordDM;

567:     DMGetDS(cdm, &prob);
568:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
569:     PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
570:     DMGetCoordinateDM(cdm,&coordDM);
571:     if (feAux) {
572:       DM      dmAux;
573:       PetscDS probAux;

575:       DMClone(cdm, &dmAux);
576:       DMSetCoordinateDM(dmAux, coordDM);
577:       DMGetDS(dmAux, &probAux);
578:       PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
579:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
580:       SetupMaterial(cdm, dmAux, user);
581:       DMDestroy(&dmAux);
582:     }
583:     if (feCh) {
584:       DM      dmCh;
585:       PetscDS probCh;

587:       DMClone(cdm, &dmCh);
588:       DMSetCoordinateDM(dmCh, coordDM);
589:       DMGetDS(dmCh, &probCh);
590:       PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
591:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
592:       DMDestroy(&dmCh);
593:     }
594:     if (user->bcType == DIRICHLET) {
595:       PetscBool hasLabel;

597:       DMHasLabel(cdm, "marker", &hasLabel);
598:       if (!hasLabel) {CreateBCLabel(cdm, "marker");}
599:     }
600:     SetupProblem(cdm, user);
601:     DMAddBoundary(cdm, user->bcType == DIRICHLET ? PETSC_TRUE : PETSC_FALSE, "wall", user->bcType == DIRICHLET ? "marker" : "boundary", 0, 0, NULL, (void (*)()) user->exactFuncs[0], 1, &id, user);
602:     DMGetCoarseDM(cdm, &cdm);
603:   }
604:   PetscFEDestroy(&fe);
605:   PetscFEDestroy(&feBd);
606:   PetscFEDestroy(&feAux);
607:   PetscFEDestroy(&feCh);
608:   return(0);
609: }

611:  #include petsc/private/petscimpl.h

615: /*@C
616:   KSPMonitorError - Outputs the error at each iteration of an iterative solver.

618:   Collective on KSP

620:   Input Parameters:
621: + ksp   - the KSP
622: . its   - iteration number
623: . rnorm - 2-norm, preconditioned residual value (may be estimated).
624: - ctx   - monitor context

626:   Level: intermediate

628: .keywords: KSP, default, monitor, residual
629: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorDefault()
630: @*/
631: static PetscErrorCode KSPMonitorError(KSP ksp, PetscInt its, PetscReal rnorm, void *ctx)
632: {
633:   AppCtx        *user = (AppCtx *) ctx;
634:   DM             dm;
635:   Vec            du, r;
636:   PetscInt       level = 0;
637:   PetscBool      hasLevel;
638:   PetscViewer    viewer;
639:   char           buf[256];

643:   KSPGetDM(ksp, &dm);
644:   /* Calculate solution */
645:   {
646:     PC        pc = user->pcmg; /* The MG PC */
647:     DM        fdm,  cdm;
648:     KSP       fksp, cksp;
649:     Vec       fu,   cu;
650:     PetscInt  levels, l;

652:     KSPBuildSolution(ksp, NULL, &du);
653:     PetscObjectComposedDataGetInt((PetscObject) ksp, PetscMGLevelId, level, hasLevel);
654:     PCMGGetLevels(pc, &levels);
655:     PCMGGetSmoother(pc, levels-1, &fksp);
656:     KSPBuildSolution(fksp, NULL, &fu);
657:     for (l = levels-1; l > level; --l) {
658:       Mat R;
659:       Vec s;

661:       PCMGGetSmoother(pc, l-1, &cksp);
662:       KSPGetDM(cksp, &cdm);
663:       DMGetGlobalVector(cdm, &cu);
664:       PCMGGetRestriction(pc, l, &R);
665:       PCMGGetRScale(pc, l, &s);
666:       MatRestrict(R, fu, cu);
667:       VecPointwiseMult(cu, cu, s);
668:       if (l < levels-1) {DMRestoreGlobalVector(fdm, &fu);}
669:       fdm  = cdm;
670:       fu   = cu;
671:     }
672:     if (levels-1 > level) {
673:       VecAXPY(du, 1.0, cu);
674:       DMRestoreGlobalVector(cdm, &cu);
675:     }
676:   }
677:   /* Calculate error */
678:   DMGetGlobalVector(dm, &r);
679:   DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
680:   VecAXPY(r,-1.0,du);
681:   PetscObjectSetName((PetscObject) r, "solution error");
682:   /* View error */
683:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
684:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
685:   VecView(r, viewer);
686:   /* Cleanup */
687:   PetscViewerDestroy(&viewer);
688:   DMRestoreGlobalVector(dm, &r);
689:   return(0);
690: }

694: /*@C
695:   SNESMonitorError - Outputs the error at each iteration of an iterative solver.

697:   Collective on SNES

699:   Input Parameters:
700: + snes  - the SNES
701: . its   - iteration number
702: . rnorm - 2-norm of residual
703: - ctx   - user context

705:   Level: intermediate

707: .keywords: SNES, nonlinear, default, monitor, norm
708: .seealso: SNESMonitorDefault(), SNESMonitorSet(), SNESMonitorSolution()
709: @*/
710: static PetscErrorCode SNESMonitorError(SNES snes, PetscInt its, PetscReal rnorm, void *ctx)
711: {
712:   AppCtx        *user = (AppCtx *) ctx;
713:   DM             dm;
714:   Vec            u, r;
715:   PetscInt       level;
716:   PetscBool      hasLevel;
717:   PetscViewer    viewer;
718:   char           buf[256];

722:   SNESGetDM(snes, &dm);
723:   /* Calculate error */
724:   SNESGetSolution(snes, &u);
725:   DMGetGlobalVector(dm, &r);
726:   PetscObjectSetName((PetscObject) r, "solution error");
727:   DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
728:   VecAXPY(r, -1.0, u);
729:   /* View error */
730:   PetscObjectComposedDataGetInt((PetscObject) snes, PetscMGLevelId, level, hasLevel);
731:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
732:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
733:   VecView(r, viewer);
734:   /* Cleanup */
735:   PetscViewerDestroy(&viewer);
736:   DMRestoreGlobalVector(dm, &r);
737:   return(0);
738: }

742: int main(int argc, char **argv)
743: {
744:   DM             dm;          /* Problem specification */
745:   SNES           snes;        /* nonlinear solver */
746:   Vec            u,r;         /* solution, residual vectors */
747:   Mat            A,J;         /* Jacobian matrix */
748:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
749:   AppCtx         user;        /* user-defined work context */
750:   JacActionCtx   userJ;       /* context for Jacobian MF action */
751:   PetscInt       its;         /* iterations for convergence */
752:   PetscReal      error = 0.0; /* L_2 error in the solution */
753:   PetscBool      isFAS;

756:   PetscInitialize(&argc, &argv, NULL, help);
757:   ProcessOptions(PETSC_COMM_WORLD, &user);
758:   SNESCreate(PETSC_COMM_WORLD, &snes);
759:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
760:   SNESSetDM(snes, dm);
761:   DMSetApplicationContext(dm, &user);

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

766:   DMCreateGlobalVector(dm, &u);
767:   PetscObjectSetName((PetscObject) u, "potential");
768:   VecDuplicate(u, &r);

770:   DMSetMatType(dm,MATAIJ);
771:   DMCreateMatrix(dm, &J);
772:   if (user.jacobianMF) {
773:     PetscInt M, m, N, n;

775:     MatGetSize(J, &M, &N);
776:     MatGetLocalSize(J, &m, &n);
777:     MatCreate(PETSC_COMM_WORLD, &A);
778:     MatSetSizes(A, m, n, M, N);
779:     MatSetType(A, MATSHELL);
780:     MatSetUp(A);
781: #if 0
782:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
783: #endif

785:     userJ.dm   = dm;
786:     userJ.J    = J;
787:     userJ.user = &user;

789:     DMCreateLocalVector(dm, &userJ.u);
790:     DMProjectFunctionLocal(dm, 0.0, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
791:     MatShellSetContext(A, &userJ);
792:   } else {
793:     A = J;
794:   }
795:   if (user.bcType == NEUMANN) {
796:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
797:     MatSetNullSpace(A, nullSpace);
798:   }

800:   DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
801:   SNESSetJacobian(snes, A, J, NULL, NULL);

803:   SNESSetFromOptions(snes);

805:   DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
806:   if (user.restart) {
807: #if defined(PETSC_HAVE_HDF5)
808:     PetscViewer viewer;

810:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
811:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
812:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
813:     PetscViewerFileSetName(viewer, user.filename);
814:     PetscViewerHDF5PushGroup(viewer, "/fields");
815:     VecLoad(u, viewer);
816:     PetscViewerHDF5PopGroup(viewer);
817:     PetscViewerDestroy(&viewer);
818: #endif
819:   }
820:   if (user.showInitial) {
821:     Vec lv;
822:     DMGetLocalVector(dm, &lv);
823:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
824:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
825:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
826:     DMRestoreLocalVector(dm, &lv);
827:   }
828:   if (user.viewHierarchy) {
829:     SNES      lsnes;
830:     KSP       ksp;
831:     PC        pc;
832:     PetscInt  numLevels, l;
833:     PetscBool isMG;

835:     PetscObjectTypeCompare((PetscObject) snes, SNESFAS, &isFAS);
836:     if (isFAS) {
837:       SNESFASGetLevels(snes, &numLevels);
838:       for (l = 0; l < numLevels; ++l) {
839:         SNESFASGetCycleSNES(snes, l, &lsnes);
840:         SNESMonitorSet(lsnes, SNESMonitorError, &user, NULL);
841:       }
842:     } else {
843:       SNESGetKSP(snes, &ksp);
844:       KSPGetPC(ksp, &pc);
845:       PetscObjectTypeCompare((PetscObject) pc, PCMG, &isMG);
846:       if (isMG) {
847:         user.pcmg = pc;
848:         PCMGGetLevels(pc, &numLevels);
849:         for (l = 0; l < numLevels; ++l) {
850:           PCMGGetSmootherDown(pc, l, &ksp);
851:           KSPMonitorSet(ksp, KSPMonitorError, &user, NULL);
852:         }
853:       }
854:     }
855:   }
856:   if (user.runType == RUN_FULL || user.runType == RUN_EXACT) {
857:     PetscErrorCode (*initialGuess[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {zero};

859:     if (user.runType == RUN_FULL) {
860:       DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
861:     }
862:     if (user.debug) {
863:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
864:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
865:     }
866:     SNESSolve(snes, NULL, u);
867:     SNESGetIterationNumber(snes, &its);
868:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
869:     DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
870:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
871:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
872:     if (user.showSolution) {
873:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
874:       VecChop(u, 3.0e-9);
875:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
876:     }
877:   } else if (user.runType == RUN_PERF) {
878:     PetscReal res = 0.0;

880:     SNESComputeFunction(snes, u, r);
881:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
882:     VecChop(r, 1.0e-10);
883:     VecNorm(r, NORM_2, &res);
884:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
885:   } else {
886:     PetscReal res = 0.0;

888:     /* Check discretization error */
889:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
890:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
891:     DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
892:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
893:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
894:     /* Check residual */
895:     SNESComputeFunction(snes, u, r);
896:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
897:     VecChop(r, 1.0e-10);
898:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
899:     VecNorm(r, NORM_2, &res);
900:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
901:     /* Check Jacobian */
902:     {
903:       Vec          b;

905:       SNESComputeJacobian(snes, u, A, A);
906:       VecDuplicate(u, &b);
907:       VecSet(r, 0.0);
908:       SNESComputeFunction(snes, r, b);
909:       MatMult(A, u, r);
910:       VecAXPY(r, 1.0, b);
911:       VecDestroy(&b);
912:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
913:       VecChop(r, 1.0e-10);
914:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
915:       VecNorm(r, NORM_2, &res);
916:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
917:     }
918:   }
919:   VecViewFromOptions(u, NULL, "-vec_view");

921:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
922:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
923:   if (A != J) {MatDestroy(&A);}
924:   MatDestroy(&J);
925:   VecDestroy(&u);
926:   VecDestroy(&r);
927:   SNESDestroy(&snes);
928:   DMDestroy(&dm);
929:   PetscFree(user.exactFuncs);
930:   PetscFinalize();
931:   return 0;
932: }