Actual source code: ex12.c

petsc-master 2016-08-23
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, quiet, nonzInit;
 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: static PetscErrorCode ecks(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 44: {
 45:   u[0] = x[0];
 46:   return 0;
 47: }

 49: /*
 50:   In 2D for Dirichlet conditions, we use exact solution:

 52:     u = x^2 + y^2
 53:     f = 4

 55:   so that

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

 59:   For Neumann conditions, we have

 61:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (bottom)
 62:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top)
 63:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
 64:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

 66:   Which we can express as

 68:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y} \cdot \hat n = 2 (x + y)
 69: */
 70: static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 71: {
 72:   *u = x[0]*x[0] + x[1]*x[1];
 73:   return 0;
 74: }

 76: void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 77:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 78:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 79:           PetscReal t, const PetscReal x[], PetscScalar f0[])
 80: {
 81:   f0[0] = 4.0;
 82: }

 84: void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 85:              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 86:              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 87:              PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 88: {
 89:   PetscInt d;
 90:   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += -n[d]*2.0*x[d];
 91: }

 93: void f0_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 94:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 95:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 96:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 97: {
 98:   f0[0] = 0.0;
 99: }

101: void f1_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
102:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
103:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
104:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
105: {
106:   PetscInt comp;
107:   for (comp = 0; comp < dim; ++comp) f1[comp] = 0.0;
108: }

110: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
111: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
112:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
113:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
114:           PetscReal t, const PetscReal x[], PetscScalar f1[])
115: {
116:   PetscInt d;
117:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
118: }

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

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

134:     u  = x^2 + y^2
135:     f  = 6 (x + y)
136:     nu = (x + y)

138:   so that

140:     -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
141: */
142: static PetscErrorCode nu_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
143: {
144:   *u = x[0] + x[1];
145:   return 0;
146: }

148: void f0_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
149:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
150:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
151:                    PetscReal t, const PetscReal x[], PetscScalar f0[])
152: {
153:   f0[0] = 6.0*(x[0] + x[1]);
154: }

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

166: void f1_field_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
167:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
168:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
169:                 PetscReal t, const PetscReal x[], PetscScalar f1[])
170: {
171:   PetscInt d;
172:   for (d = 0; d < dim; ++d) f1[d] = a[0]*u_x[d];
173: }

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

186: void g3_field_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
187:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
188:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
189:                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
190: {
191:   PetscInt d;
192:   for (d = 0; d < dim; ++d) g3[d*dim+d] = a[0];
193: }

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

198:     u  = x^2 + y^2
199:     f  = 16 (x^2 + y^2)
200:     nu = 1/2 |grad u|^2

202:   so that

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

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

226: /*
227:   grad (u + eps w) - grad u = eps grad w

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

250: /*
251:   In 3D for Dirichlet conditions we use exact solution:

253:     u = 2/3 (x^2 + y^2 + z^2)
254:     f = 4

256:   so that

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

260:   For Neumann conditions, we have

262:     -\nabla u \cdot -\hat z |_{z=0} =  (2z)|_{z=0} =  0 (bottom)
263:     -\nabla u \cdot  \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
264:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (front)
265:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
266:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
267:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

269:   Which we can express as

271:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
272: */
273: static PetscErrorCode quadratic_u_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
274: {
275:   *u = 2.0*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/3.0;
276:   return 0;
277: }

281: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
282: {
283:   const char    *bcTypes[3]  = {"neumann", "dirichlet", "none"};
284:   const char    *runTypes[4] = {"full", "exact", "test", "perf"};
285:   const char    *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
286:   PetscInt       bc, run, coeff;
287:   PetscBool      flg;

291:   options->debug               = 0;
292:   options->runType             = RUN_FULL;
293:   options->dim                 = 2;
294:   options->filename[0]         = '\0';
295:   options->interpolate         = PETSC_FALSE;
296:   options->refinementLimit     = 0.0;
297:   options->bcType              = DIRICHLET;
298:   options->variableCoefficient = COEFF_NONE;
299:   options->jacobianMF          = PETSC_FALSE;
300:   options->showInitial         = PETSC_FALSE;
301:   options->showSolution        = PETSC_FALSE;
302:   options->restart             = PETSC_FALSE;
303:   options->check               = PETSC_FALSE;
304:   options->viewHierarchy       = PETSC_FALSE;
305:   options->simplex             = PETSC_TRUE;
306:   options->quiet               = PETSC_FALSE;
307:   options->nonzInit            = PETSC_FALSE;

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

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

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

330:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
331:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
332:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
333:   PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
334:   PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
335:   PetscOptionsBool("-dm_view_hierarchy", "View the coarsened hierarchy", "ex12.c", options->viewHierarchy, &options->viewHierarchy, NULL);
336:   PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex12.c", options->simplex, &options->simplex, NULL);
337:   PetscOptionsBool("-quiet", "Don't print any vecs", "ex12.c", options->quiet, &options->quiet, NULL);
338:   PetscOptionsBool("-nonzero_initial_guess", "nonzero intial guess", "ex12.c", options->nonzInit, &options->nonzInit, NULL);
339:   PetscOptionsEnd();
340:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
341:   return(0);
342: }

346: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
347: {
348:   DMLabel        label;

352:   DMCreateLabel(dm, name);
353:   DMGetLabel(dm, name, &label);
354:   DMPlexMarkBoundaryFaces(dm, label);
355:   DMPlexLabelComplete(dm, label);
356:   return(0);
357: }

361: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
362: {
363:   PetscInt       dim             = user->dim;
364:   const char    *filename        = user->filename;
365:   PetscBool      interpolate     = user->interpolate;
366:   PetscReal      refinementLimit = user->refinementLimit;
367:   size_t         len;

371:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
372:   PetscStrlen(filename, &len);
373:   if (!len) {
374:     if (user->simplex) {
375:       DMPlexCreateBoxMesh(comm, dim, dim == 2 ? 2 : 1, interpolate, dm);
376:     }
377:     else {
378:       PetscInt cells[3] = {1, 1, 1}; /* coarse mesh is one cell; refine from there */

380:       DMPlexCreateHexBoxMesh(comm, dim, cells,  DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
381:     }
382:     PetscObjectSetName((PetscObject) *dm, "Mesh");
383:   } else {
384:     DMPlexCreateFromFile(comm, filename, interpolate, dm);
385:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
386:   }
387:   {
388:     DM refinedMesh     = NULL;
389:     DM distributedMesh = NULL;

391:     /* Refine mesh using a volume constraint */
392:     DMPlexSetRefinementLimit(*dm, refinementLimit);
393:     DMRefine(*dm, comm, &refinedMesh);
394:     if (refinedMesh) {
395:       const char *name;

397:       PetscObjectGetName((PetscObject) *dm,         &name);
398:       PetscObjectSetName((PetscObject) refinedMesh,  name);
399:       DMDestroy(dm);
400:       *dm  = refinedMesh;
401:     }
402:     /* Distribute mesh over processes */
403:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
404:     if (distributedMesh) {
405:       DMDestroy(dm);
406:       *dm  = distributedMesh;
407:     }
408:   }
409:   if (user->bcType == NEUMANN) {
410:     DMLabel   label;

412:     DMCreateLabel(*dm, "boundary");
413:     DMGetLabel(*dm, "boundary", &label);
414:     DMPlexMarkBoundaryFaces(*dm, label);
415:   } else if (user->bcType == DIRICHLET) {
416:     PetscBool hasLabel;

418:     DMHasLabel(*dm,"marker",&hasLabel);
419:     if (!hasLabel) {CreateBCLabel(*dm, "marker");}
420:   }
421:   {
422:     char      convType[256];
423:     PetscBool flg;

425:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
426:     PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);
427:     PetscOptionsEnd();
428:     if (flg) {
429:       DM dmConv;

431:       DMConvert(*dm,convType,&dmConv);
432:       if (dmConv) {
433:         DMDestroy(dm);
434:         *dm  = dmConv;
435:       }
436:     }
437:   }
438:   DMSetFromOptions(*dm);
439:   DMViewFromOptions(*dm, NULL, "-dm_view");
440:   if (user->viewHierarchy) {
441:     DM       cdm = *dm;
442:     PetscInt i   = 0;
443:     char     buf[256];

445:     while (cdm) {
446:       DMSetUp(cdm);
447:       DMGetCoarseDM(cdm, &cdm);
448:       ++i;
449:     }
450:     cdm = *dm;
451:     while (cdm) {
452:       PetscViewer       viewer;
453:       PetscBool   isHDF5, isVTK;

455:       --i;
456:       PetscViewerCreate(comm,&viewer);
457:       PetscViewerSetType(viewer,PETSCVIEWERHDF5);
458:       PetscViewerSetOptionsPrefix(viewer,"hierarchy_");
459:       PetscViewerSetFromOptions(viewer);
460:       PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
461:       PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
462:       if (isHDF5) {
463:         PetscSNPrintf(buf, 256, "ex12-%d.h5", i);
464:       }
465:       else if (isVTK) {
466:         PetscSNPrintf(buf, 256, "ex12-%d.vtu", i);
467:         PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
468:       }
469:       else {
470:         PetscSNPrintf(buf, 256, "ex12-%d", i);
471:       }
472:       PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
473:       PetscViewerFileSetName(viewer,buf);
474:       DMView(cdm, viewer);
475:       PetscViewerDestroy(&viewer);
476:       DMGetCoarseDM(cdm, &cdm);
477:     }
478:   }
479:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
480:   return(0);
481: }

485: static PetscErrorCode SetupProblem(PetscDS prob, AppCtx *user)
486: {
487:   const PetscInt id = 1;

491:   switch (user->variableCoefficient) {
492:   case COEFF_NONE:
493:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
494:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
495:     break;
496:   case COEFF_ANALYTIC:
497:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
498:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
499:     break;
500:   case COEFF_FIELD:
501:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
502:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
503:     break;
504:   case COEFF_NONLINEAR:
505:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
506:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
507:     break;
508:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
509:   }
510:   switch (user->dim) {
511:   case 2:
512:     user->exactFuncs[0] = quadratic_u_2d;
513:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
514:     break;
515:   case 3:
516:     user->exactFuncs[0] = quadratic_u_3d;
517:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
518:     break;
519:   default:
520:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
521:   }
522:   PetscDSAddBoundary(prob, user->bcType == DIRICHLET ? PETSC_TRUE : PETSC_FALSE, "wall", user->bcType == DIRICHLET ? "marker" : "boundary", 0, 0, NULL, (void (*)()) user->exactFuncs[0], 1, &id, user);
523:   return(0);
524: }

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

535:   DMCreateLocalVector(dmAux, &nu);
536:   DMProjectFunctionLocal(dmAux, 0.0, matFuncs, NULL, INSERT_ALL_VALUES, nu);
537:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
538:   VecDestroy(&nu);
539:   return(0);
540: }

544: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
545: {
546:   DM             cdm   = dm;
547:   const PetscInt dim   = user->dim;
548:   PetscFE        feAux = NULL;
549:   PetscFE        feBd  = NULL;
550:   PetscFE        feCh  = NULL;
551:   PetscFE        fe;
552:   PetscDS        prob, probAux = NULL, probCh = NULL;
553:   PetscBool      simplex = user->simplex;

557:   /* Create finite element */
558:   PetscFECreateDefault(dm, dim, 1, simplex, NULL, -1, &fe);
559:   PetscObjectSetName((PetscObject) fe, "potential");
560:   if (user->bcType == NEUMANN) {
561:     PetscFECreateDefault(dm, dim-1, 1, simplex, "bd_", -1, &feBd);
562:     PetscObjectSetName((PetscObject) feBd, "potential");
563:   }
564:   if (user->variableCoefficient == COEFF_FIELD) {
565:     PetscQuadrature q;

567:     PetscFECreateDefault(dm, dim, 1, simplex, "mat_", -1, &feAux);
568:     PetscFEGetQuadrature(fe, &q);
569:     PetscFESetQuadrature(feAux, q);
570:     PetscDSCreate(PetscObjectComm((PetscObject)dm),&probAux);
571:     PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
572:   }
573:   if (user->check) {
574:     PetscFECreateDefault(dm, dim, 1, simplex, "ch_", -1, &feCh);
575:     PetscDSCreate(PetscObjectComm((PetscObject)dm),&probCh);
576:     PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
577:   }
578:   /* Set discretization and boundary conditions for each mesh */
579:   DMGetDS(dm, &prob);
580:   PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
581:   PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
582:   SetupProblem(prob, user);
583:   while (cdm) {
584:     DM coordDM;

586:     DMSetDS(cdm,prob);
587:     DMGetCoordinateDM(cdm,&coordDM);
588:     if (feAux) {
589:       DM      dmAux;

591:       DMClone(cdm, &dmAux);
592:       DMSetCoordinateDM(dmAux, coordDM);
593:       DMSetDS(dmAux, probAux);
594:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
595:       SetupMaterial(cdm, dmAux, user);
596:       DMDestroy(&dmAux);
597:     }
598:     if (feCh) {
599:       DM      dmCh;

601:       DMClone(cdm, &dmCh);
602:       DMSetCoordinateDM(dmCh, coordDM);
603:       DMSetDS(dmCh, probCh);
604:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
605:       DMDestroy(&dmCh);
606:     }
607:     if (user->bcType == DIRICHLET) {
608:       PetscBool hasLabel;

610:       DMHasLabel(cdm, "marker", &hasLabel);
611:       if (!hasLabel) {CreateBCLabel(cdm, "marker");}
612:     }
613:     DMGetCoarseDM(cdm, &cdm);
614:   }
615:   PetscFEDestroy(&fe);
616:   PetscFEDestroy(&feBd);
617:   PetscFEDestroy(&feAux);
618:   PetscFEDestroy(&feCh);
619:   PetscDSDestroy(&probAux);
620:   PetscDSDestroy(&probCh);
621:   return(0);
622: }

624:  #include petsc/private/petscimpl.h

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

631:   Collective on KSP

633:   Input Parameters:
634: + ksp   - the KSP
635: . its   - iteration number
636: . rnorm - 2-norm, preconditioned residual value (may be estimated).
637: - ctx   - monitor context

639:   Level: intermediate

641: .keywords: KSP, default, monitor, residual
642: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorDefault()
643: @*/
644: static PetscErrorCode KSPMonitorError(KSP ksp, PetscInt its, PetscReal rnorm, void *ctx)
645: {
646:   AppCtx        *user = (AppCtx *) ctx;
647:   DM             dm;
648:   Vec            du, r;
649:   PetscInt       level = 0;
650:   PetscBool      hasLevel;
651:   PetscViewer    viewer;
652:   char           buf[256];

656:   KSPGetDM(ksp, &dm);
657:   /* Calculate solution */
658:   {
659:     PC        pc = user->pcmg; /* The MG PC */
660:     DM        fdm,  cdm;
661:     KSP       fksp, cksp;
662:     Vec       fu,   cu;
663:     PetscInt  levels, l;

665:     KSPBuildSolution(ksp, NULL, &du);
666:     PetscObjectComposedDataGetInt((PetscObject) ksp, PetscMGLevelId, level, hasLevel);
667:     PCMGGetLevels(pc, &levels);
668:     PCMGGetSmoother(pc, levels-1, &fksp);
669:     KSPBuildSolution(fksp, NULL, &fu);
670:     for (l = levels-1; l > level; --l) {
671:       Mat R;
672:       Vec s;

674:       PCMGGetSmoother(pc, l-1, &cksp);
675:       KSPGetDM(cksp, &cdm);
676:       DMGetGlobalVector(cdm, &cu);
677:       PCMGGetRestriction(pc, l, &R);
678:       PCMGGetRScale(pc, l, &s);
679:       MatRestrict(R, fu, cu);
680:       VecPointwiseMult(cu, cu, s);
681:       if (l < levels-1) {DMRestoreGlobalVector(fdm, &fu);}
682:       fdm  = cdm;
683:       fu   = cu;
684:     }
685:     if (levels-1 > level) {
686:       VecAXPY(du, 1.0, cu);
687:       DMRestoreGlobalVector(cdm, &cu);
688:     }
689:   }
690:   /* Calculate error */
691:   DMGetGlobalVector(dm, &r);
692:   DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
693:   VecAXPY(r,-1.0,du);
694:   PetscObjectSetName((PetscObject) r, "solution error");
695:   /* View error */
696:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
697:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
698:   VecView(r, viewer);
699:   /* Cleanup */
700:   PetscViewerDestroy(&viewer);
701:   DMRestoreGlobalVector(dm, &r);
702:   return(0);
703: }

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

710:   Collective on SNES

712:   Input Parameters:
713: + snes  - the SNES
714: . its   - iteration number
715: . rnorm - 2-norm of residual
716: - ctx   - user context

718:   Level: intermediate

720: .keywords: SNES, nonlinear, default, monitor, norm
721: .seealso: SNESMonitorDefault(), SNESMonitorSet(), SNESMonitorSolution()
722: @*/
723: static PetscErrorCode SNESMonitorError(SNES snes, PetscInt its, PetscReal rnorm, void *ctx)
724: {
725:   AppCtx        *user = (AppCtx *) ctx;
726:   DM             dm;
727:   Vec            u, r;
728:   PetscInt       level;
729:   PetscBool      hasLevel;
730:   PetscViewer    viewer;
731:   char           buf[256];

735:   SNESGetDM(snes, &dm);
736:   /* Calculate error */
737:   SNESGetSolution(snes, &u);
738:   DMGetGlobalVector(dm, &r);
739:   PetscObjectSetName((PetscObject) r, "solution error");
740:   DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
741:   VecAXPY(r, -1.0, u);
742:   /* View error */
743:   PetscObjectComposedDataGetInt((PetscObject) snes, PetscMGLevelId, level, hasLevel);
744:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
745:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
746:   VecView(r, viewer);
747:   /* Cleanup */
748:   PetscViewerDestroy(&viewer);
749:   DMRestoreGlobalVector(dm, &r);
750:   return(0);
751: }

755: int main(int argc, char **argv)
756: {
757:   DM             dm;          /* Problem specification */
758:   SNES           snes;        /* nonlinear solver */
759:   Vec            u,r;         /* solution, residual vectors */
760:   Mat            A,J;         /* Jacobian matrix */
761:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
762:   AppCtx         user;        /* user-defined work context */
763:   JacActionCtx   userJ;       /* context for Jacobian MF action */
764:   PetscInt       its;         /* iterations for convergence */
765:   PetscReal      error = 0.0; /* L_2 error in the solution */
766:   PetscBool      isFAS;

769:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
770:   ProcessOptions(PETSC_COMM_WORLD, &user);
771:   SNESCreate(PETSC_COMM_WORLD, &snes);
772:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
773:   SNESSetDM(snes, dm);
774:   DMSetApplicationContext(dm, &user);

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

779:   DMCreateGlobalVector(dm, &u);
780:   PetscObjectSetName((PetscObject) u, "potential");
781:   VecDuplicate(u, &r);

783:   DMCreateMatrix(dm, &J);
784:   if (user.jacobianMF) {
785:     PetscInt M, m, N, n;

787:     MatGetSize(J, &M, &N);
788:     MatGetLocalSize(J, &m, &n);
789:     MatCreate(PETSC_COMM_WORLD, &A);
790:     MatSetSizes(A, m, n, M, N);
791:     MatSetType(A, MATSHELL);
792:     MatSetUp(A);
793: #if 0
794:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
795: #endif

797:     userJ.dm   = dm;
798:     userJ.J    = J;
799:     userJ.user = &user;

801:     DMCreateLocalVector(dm, &userJ.u);
802:     DMProjectFunctionLocal(dm, 0.0, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
803:     MatShellSetContext(A, &userJ);
804:   } else {
805:     A = J;
806:   }
807:   if (user.bcType == NEUMANN) {
808:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
809:     MatSetNullSpace(A, nullSpace);
810:   }

812:   DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
813:   SNESSetJacobian(snes, A, J, NULL, NULL);

815:   SNESSetFromOptions(snes);

817:   DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
818:   if (user.restart) {
819: #if defined(PETSC_HAVE_HDF5)
820:     PetscViewer viewer;

822:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
823:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
824:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
825:     PetscViewerFileSetName(viewer, user.filename);
826:     PetscViewerHDF5PushGroup(viewer, "/fields");
827:     VecLoad(u, viewer);
828:     PetscViewerHDF5PopGroup(viewer);
829:     PetscViewerDestroy(&viewer);
830: #endif
831:   }
832:   if (user.showInitial) {
833:     Vec lv;
834:     DMGetLocalVector(dm, &lv);
835:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
836:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
837:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
838:     DMRestoreLocalVector(dm, &lv);
839:   }
840:   if (user.viewHierarchy) {
841:     SNES      lsnes;
842:     KSP       ksp;
843:     PC        pc;
844:     PetscInt  numLevels, l;
845:     PetscBool isMG;

847:     PetscObjectTypeCompare((PetscObject) snes, SNESFAS, &isFAS);
848:     if (isFAS) {
849:       SNESFASGetLevels(snes, &numLevels);
850:       for (l = 0; l < numLevels; ++l) {
851:         SNESFASGetCycleSNES(snes, l, &lsnes);
852:         SNESMonitorSet(lsnes, SNESMonitorError, &user, NULL);
853:       }
854:     } else {
855:       SNESGetKSP(snes, &ksp);
856:       KSPGetPC(ksp, &pc);
857:       PetscObjectTypeCompare((PetscObject) pc, PCMG, &isMG);
858:       if (isMG) {
859:         user.pcmg = pc;
860:         PCMGGetLevels(pc, &numLevels);
861:         for (l = 0; l < numLevels; ++l) {
862:           PCMGGetSmootherDown(pc, l, &ksp);
863:           KSPMonitorSet(ksp, KSPMonitorError, &user, NULL);
864:         }
865:       }
866:     }
867:   }
868:   if (user.runType == RUN_FULL || user.runType == RUN_EXACT) {
869:     PetscErrorCode (*initialGuess[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {zero};

871:     if (user.nonzInit) initialGuess[0] = ecks;
872:     if (user.runType == RUN_FULL) {
873:       DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
874:     }
875:     if (user.debug) {
876:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
877:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
878:     }
879:     SNESSolve(snes, NULL, u);
880:     SNESGetIterationNumber(snes, &its);
881:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
882:     DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
883:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
884:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
885:     if (user.showSolution) {
886:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
887:       VecChop(u, 3.0e-9);
888:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
889:     }
890:   } else if (user.runType == RUN_PERF) {
891:     PetscReal res = 0.0;

893:     SNESComputeFunction(snes, u, r);
894:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
895:     VecChop(r, 1.0e-10);
896:     VecNorm(r, NORM_2, &res);
897:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
898:   } else {
899:     PetscReal res = 0.0;

901:     /* Check discretization error */
902:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
903:     if (!user.quiet) {VecView(u, PETSC_VIEWER_STDOUT_WORLD);}
904:     DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
905:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
906:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
907:     /* Check residual */
908:     SNESComputeFunction(snes, u, r);
909:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
910:     VecChop(r, 1.0e-10);
911:     if (!user.quiet) {VecView(r, PETSC_VIEWER_STDOUT_WORLD);}
912:     VecNorm(r, NORM_2, &res);
913:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
914:     /* Check Jacobian */
915:     {
916:       Vec          b;

918:       SNESComputeJacobian(snes, u, A, A);
919:       VecDuplicate(u, &b);
920:       VecSet(r, 0.0);
921:       SNESComputeFunction(snes, r, b);
922:       MatMult(A, u, r);
923:       VecAXPY(r, 1.0, b);
924:       VecDestroy(&b);
925:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
926:       VecChop(r, 1.0e-10);
927:       if (!user.quiet) {VecView(r, PETSC_VIEWER_STDOUT_WORLD);}
928:       VecNorm(r, NORM_2, &res);
929:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
930:     }
931:   }
932:   VecViewFromOptions(u, NULL, "-vec_view");

934:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
935:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
936:   if (A != J) {MatDestroy(&A);}
937:   MatDestroy(&J);
938:   VecDestroy(&u);
939:   VecDestroy(&r);
940:   SNESDestroy(&snes);
941:   DMDestroy(&dm);
942:   PetscFree(user.exactFuncs);
943:   PetscFinalize();
944:   return ierr;
945: }