Actual source code: ex12.c

petsc-master 2015-08-27
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_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:   /* Problem definition */
 29:   BCType        bcType;
 30:   CoeffType     variableCoefficient;
 31:   PetscErrorCode (**exactFuncs)(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 32:   /* Solver */
 33:   PC            pcmg;              /* This is needed for error monitoring */
 34: } AppCtx;

 36: PetscErrorCode zero(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 37: {
 38:   u[0] = 0.0;
 39:   return 0;
 40: }

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

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

 48:   so that

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

 52:   For Neumann conditions, we have

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

 59:   Which we can express as

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

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

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

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

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

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

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

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

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

131:   so that

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

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

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

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

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

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

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

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

195:   so that

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

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

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

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

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

246:     u = x^2 + y^2 + z^2
247:     f = 6

249:   so that

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

253:   For Neumann conditions, we have

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

262:   Which we can express as

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

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

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

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

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

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

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

333: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
334: {
335:   PetscInt       dim             = user->dim;
336:   const char    *filename        = user->filename;
337:   PetscBool      interpolate     = user->interpolate;
338:   PetscReal      refinementLimit = user->refinementLimit;
339:   size_t         len;

343:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
344:   PetscStrlen(filename, &len);
345:   if (!len) {
346:     DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
347:     PetscObjectSetName((PetscObject) *dm, "Mesh");
348:   } else {
349:     DMPlexCreateFromFile(comm, filename, interpolate, dm);
350:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
351:   }
352:   {
353:     DM refinedMesh     = NULL;
354:     DM distributedMesh = NULL;

356:     /* Refine mesh using a volume constraint */
357:     DMPlexSetRefinementLimit(*dm, refinementLimit);
358:     DMRefine(*dm, comm, &refinedMesh);
359:     if (refinedMesh) {
360:       const char *name;

362:       PetscObjectGetName((PetscObject) *dm,         &name);
363:       PetscObjectSetName((PetscObject) refinedMesh,  name);
364:       DMDestroy(dm);
365:       *dm  = refinedMesh;
366:     }
367:     /* Distribute mesh over processes */
368:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
369:     if (distributedMesh) {
370:       DMDestroy(dm);
371:       *dm  = distributedMesh;
372:     }
373:   }
374:   if (user->bcType == NEUMANN) {
375:     DMLabel label;

377:     DMPlexCreateLabel(*dm, "boundary");
378:     DMPlexGetLabel(*dm, "boundary", &label);
379:     DMPlexMarkBoundaryFaces(*dm, label);
380:   }
381:   DMSetFromOptions(*dm);
382:   /* Must have boundary marker for Dirichlet conditions */
383:   {
384:     DM cdm = *dm;

386:     while (cdm) {
387:       PetscBool hasBdLabel;

389:       DMPlexHasLabel(cdm, "marker", &hasBdLabel);
390:       if (!hasBdLabel) {
391:         DMLabel label;

393:         DMPlexCreateLabel(cdm, "marker");
394:         DMPlexGetLabel(cdm, "marker", &label);
395:         DMPlexMarkBoundaryFaces(cdm, label);
396:         DMPlexLabelComplete(cdm, label);
397:       }
398:       DMPlexGetCoarseDM(cdm, &cdm);
399:     }
400:   }

402:   DMViewFromOptions(*dm, NULL, "-dm_view");
403:   if (user->viewHierarchy) {
404:     DM       cdm = *dm;
405:     PetscInt i   = 0;
406:     char     buf[256];

408:     while (cdm) {DMPlexGetCoarseDM(cdm, &cdm); ++i;}
409:     cdm = *dm;
410:     while (cdm) {
411:       PetscViewer viewer;

413:       --i;
414:       PetscSNPrintf(buf, 256, "ex12-%d.h5", i);
415:       PetscViewerHDF5Open(comm, buf, FILE_MODE_WRITE, &viewer);
416:       DMView(cdm, viewer);
417:       PetscViewerDestroy(&viewer);
418:       DMPlexGetCoarseDM(cdm, &cdm);
419:     }
420:   }
421:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
422:   return(0);
423: }

427: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
428: {
429:   PetscDS        prob;

433:   DMGetDS(dm, &prob);
434:   switch (user->variableCoefficient) {
435:   case COEFF_NONE:
436:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
437:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
438:     break;
439:   case COEFF_ANALYTIC:
440:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
441:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
442:     break;
443:   case COEFF_FIELD:
444:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
445:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
446:     break;
447:   case COEFF_NONLINEAR:
448:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
449:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
450:     break;
451:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
452:   }
453:   switch (user->dim) {
454:   case 2:
455:     user->exactFuncs[0] = quadratic_u_2d;
456:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
457:     break;
458:   case 3:
459:     user->exactFuncs[0] = quadratic_u_3d;
460:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
461:     break;
462:   default:
463:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
464:   }
465:   return(0);
466: }

470: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
471: {
472:   PetscErrorCode (*matFuncs[1])(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {nu_2d};
473:   Vec            nu;

477:   DMCreateLocalVector(dmAux, &nu);
478:   DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);
479:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
480:   VecDestroy(&nu);
481:   return(0);
482: }

486: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
487: {
488:   DM             cdm   = dm;
489:   const PetscInt dim   = user->dim;
490:   const PetscInt id    = 1;
491:   PetscFE        feAux = NULL;
492:   PetscFE        feBd  = NULL;
493:   PetscFE        feCh  = NULL;
494:   PetscFE        fe;
495:   PetscDS        prob;

499:   /* Create finite element */
500:   PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, NULL, -1, &fe);
501:   PetscObjectSetName((PetscObject) fe, "potential");
502:   if (user->bcType == NEUMANN) {
503:     PetscFECreateDefault(dm, dim-1, 1, PETSC_TRUE, "bd_", -1, &feBd);
504:     PetscObjectSetName((PetscObject) feBd, "potential");
505:   }
506:   if (user->variableCoefficient == COEFF_FIELD) {
507:     PetscQuadrature q;

509:     PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "mat_", -1, &feAux);
510:     PetscFEGetQuadrature(fe, &q);
511:     PetscFESetQuadrature(feAux, q);
512:   }
513:   if (user->check) {PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "ch_", -1, &feCh);}
514:   /* Set discretization and boundary conditions for each mesh */
515:   while (cdm) {
516:     DMGetDS(cdm, &prob);
517:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
518:     PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
519:     if (feAux) {
520:       DM      dmAux;
521:       PetscDS probAux;

523:       DMClone(cdm, &dmAux);
524:       DMPlexCopyCoordinates(cdm, dmAux);
525:       DMGetDS(dmAux, &probAux);
526:       PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
527:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
528:       SetupMaterial(cdm, dmAux, user);
529:       DMDestroy(&dmAux);
530:     }
531:     if (feCh) {
532:       DM      dmCh;
533:       PetscDS probCh;

535:       DMClone(cdm, &dmCh);
536:       DMPlexCopyCoordinates(cdm, dmCh);
537:       DMGetDS(dmCh, &probCh);
538:       PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
539:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
540:       DMDestroy(&dmCh);
541:     }
542:     SetupProblem(cdm, user);
543:     DMPlexAddBoundary(cdm, user->bcType == DIRICHLET ? PETSC_TRUE : PETSC_FALSE, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, 0, NULL, (void (*)()) user->exactFuncs[0], 1, &id, user);
544:     DMPlexGetCoarseDM(cdm, &cdm);
545:   }
546:   PetscFEDestroy(&fe);
547:   PetscFEDestroy(&feBd);
548:   PetscFEDestroy(&feAux);
549:   PetscFEDestroy(&feCh);
550:   return(0);
551: }

553:  #include petsc/private/petscimpl.h

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

560:   Collective on KSP

562:   Input Parameters:
563: + ksp   - the KSP
564: . its   - iteration number
565: . rnorm - 2-norm, preconditioned residual value (may be estimated).
566: - ctx   - monitor context

568:   Level: intermediate

570: .keywords: KSP, default, monitor, residual
571: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorDefault()
572: @*/
573: PetscErrorCode KSPMonitorError(KSP ksp, PetscInt its, PetscReal rnorm, void *ctx)
574: {
575:   AppCtx        *user = (AppCtx *) ctx;
576:   DM             dm, edm;
577:   PetscDS        prob;
578:   PetscFE        fe;
579:   Vec            du, r;
580:   PetscInt       level = 0;
581:   PetscBool      hasLevel;
582:   PetscViewer    viewer;
583:   char           buf[256];

587:   KSPGetDM(ksp, &dm);
588:   /* Create FE space for the error */
589:   PetscFECreateDefault(dm, user->dim, 1, PETSC_TRUE, "error_", -1, &fe);
590:   PetscSNPrintf(buf, 256, "%D", its);
591:   PetscObjectSetName((PetscObject) fe, buf);
592:   /* Create DM for error */
593:   DMClone(dm, &edm);
594:   DMPlexCopyCoordinates(dm, edm);
595:   DMGetDS(edm, &prob);
596:   PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
597:   /* Calculate solution */
598:   {
599:     PC       pc = user->pcmg; /* The MG PC */
600:     DM       fdm,  cdm;
601:     KSP      fksp, cksp;
602:     Vec      fu,   cu;
603:     PetscInt levels, l;

605:     KSPBuildSolution(ksp, NULL, &du);
606:     PetscObjectComposedDataGetInt((PetscObject) ksp, PetscMGLevelId, level, hasLevel);
607:     PCMGGetLevels(pc, &levels);
608:     PCMGGetSmoother(pc, levels-1, &fksp);
609:     KSPBuildSolution(fksp, NULL, &fu);
610:     for (l = levels-1; l > level; --l) {
611:       Mat R;
612:       Vec s;

614:       PCMGGetSmoother(pc, l-1, &cksp);
615:       KSPGetDM(cksp, &cdm);
616:       DMGetGlobalVector(cdm, &cu);
617:       PCMGGetRestriction(pc, l, &R);
618:       PCMGGetRScale(pc, l, &s);
619:       MatRestrict(R, fu, cu);
620:       VecPointwiseMult(cu, cu, s);
621:       if (l < levels-1) {DMRestoreGlobalVector(fdm, &fu);}
622:       fdm  = cdm;
623:       fu   = cu;
624:     }
625:     if (levels-1 > level) {
626:       VecAXPY(du, 1.0, cu);
627:       DMRestoreGlobalVector(cdm, &cu);
628:     }
629:   }
630:   /* Calculate error */
631:   DMGetGlobalVector(edm, &r);
632:   PetscObjectSetName((PetscObject) r, "solution error");
633:   DMPlexComputeL2DiffVec(edm, user->exactFuncs, NULL, du, r);
634:   /* View error */
635:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
636:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
637:   VecView(r, viewer);
638:   /* Cleanup */
639:   PetscViewerDestroy(&viewer);
640:   DMRestoreGlobalVector(edm, &r);
641:   DMDestroy(&edm);
642:   PetscFEDestroy(&fe);
643:   return(0);
644: }

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

651:   Collective on SNES

653:   Input Parameters:
654: + snes  - the SNES
655: . its   - iteration number
656: . rnorm - 2-norm of residual
657: - ctx   - user context

659:   Level: intermediate

661: .keywords: SNES, nonlinear, default, monitor, norm
662: .seealso: SNESMonitorDefault(), SNESMonitorSet(), SNESMonitorSolution()
663: @*/
664: PetscErrorCode SNESMonitorError(SNES snes, PetscInt its, PetscReal rnorm, void *ctx)
665: {
666:   AppCtx        *user = (AppCtx *) ctx;
667:   DM             dm, edm;
668:   PetscDS        prob;
669:   PetscFE        fe;
670:   Vec            u, r;
671:   PetscInt       level;
672:   PetscBool      hasLevel;
673:   PetscViewer    viewer;
674:   char           buf[256];

678:   SNESGetDM(snes, &dm);
679:   /* Create FE space for the error */
680:   PetscFECreateDefault(dm, user->dim, 1, PETSC_TRUE, "error_", -1, &fe);
681:   PetscSNPrintf(buf, 256, "%D", its);
682:   PetscObjectSetName((PetscObject) fe, buf);
683:   /* Create DM for error */
684:   DMClone(dm, &edm);
685:   DMPlexCopyCoordinates(dm, edm);
686:   DMGetDS(edm, &prob);
687:   PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
688:   /* Calculate error */
689:   SNESGetSolution(snes, &u);
690:   DMGetGlobalVector(edm, &r);
691:   PetscObjectSetName((PetscObject) r, "solution error");
692:   DMPlexComputeL2DiffVec(edm, user->exactFuncs, NULL, u, r);
693:   /* View error */
694:   PetscObjectComposedDataGetInt((PetscObject) snes, PetscMGLevelId, level, hasLevel);
695:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
696:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
697:   VecView(r, viewer);
698:   /* Cleanup */
699:   PetscViewerDestroy(&viewer);
700:   DMRestoreGlobalVector(edm, &r);
701:   DMDestroy(&edm);
702:   PetscFEDestroy(&fe);
703:   return(0);
704: }

708: int main(int argc, char **argv)
709: {
710:   DM             dm;          /* Problem specification */
711:   SNES           snes;        /* nonlinear solver */
712:   Vec            u,r;         /* solution, residual vectors */
713:   Mat            A,J;         /* Jacobian matrix */
714:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
715:   AppCtx         user;        /* user-defined work context */
716:   JacActionCtx   userJ;       /* context for Jacobian MF action */
717:   PetscInt       its;         /* iterations for convergence */
718:   PetscReal      error = 0.0; /* L_2 error in the solution */
719:   PetscBool      isFAS;

722:   PetscInitialize(&argc, &argv, NULL, help);
723:   ProcessOptions(PETSC_COMM_WORLD, &user);
724:   SNESCreate(PETSC_COMM_WORLD, &snes);
725:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
726:   SNESSetDM(snes, dm);
727:   DMSetApplicationContext(dm, &user);

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

732:   DMCreateGlobalVector(dm, &u);
733:   PetscObjectSetName((PetscObject) u, "potential");
734:   VecDuplicate(u, &r);

736:   DMSetMatType(dm,MATAIJ);
737:   DMCreateMatrix(dm, &J);
738:   if (user.jacobianMF) {
739:     PetscInt M, m, N, n;

741:     MatGetSize(J, &M, &N);
742:     MatGetLocalSize(J, &m, &n);
743:     MatCreate(PETSC_COMM_WORLD, &A);
744:     MatSetSizes(A, m, n, M, N);
745:     MatSetType(A, MATSHELL);
746:     MatSetUp(A);
747: #if 0
748:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
749: #endif

751:     userJ.dm   = dm;
752:     userJ.J    = J;
753:     userJ.user = &user;

755:     DMCreateLocalVector(dm, &userJ.u);
756:     DMPlexProjectFunctionLocal(dm, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
757:     MatShellSetContext(A, &userJ);
758:   } else {
759:     A = J;
760:   }
761:   if (user.bcType == NEUMANN) {
762:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
763:     MatSetNullSpace(A, nullSpace);
764:   }

766:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);
767:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);
768:   SNESSetJacobian(snes, A, J, NULL, NULL);

770:   SNESSetFromOptions(snes);

772:   DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
773:   if (user.restart) {
774: #if defined(PETSC_HAVE_HDF5)
775:     PetscViewer viewer;

777:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
778:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
779:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
780:     PetscViewerFileSetName(viewer, user.filename);
781:     PetscViewerHDF5PushGroup(viewer, "/fields");
782:     VecLoad(u, viewer);
783:     PetscViewerHDF5PopGroup(viewer);
784:     PetscViewerDestroy(&viewer);
785: #endif
786:   }
787:   if (user.showInitial) {
788:     Vec lv;
789:     DMGetLocalVector(dm, &lv);
790:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
791:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
792:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
793:     DMRestoreLocalVector(dm, &lv);
794:   }
795:   if (user.viewHierarchy) {
796:     SNES      lsnes;
797:     KSP       ksp;
798:     PC        pc;
799:     PetscInt  numLevels, l;
800:     PetscBool isMG;

802:     PetscObjectTypeCompare((PetscObject) snes, SNESFAS, &isFAS);
803:     if (isFAS) {
804:       SNESFASGetLevels(snes, &numLevels);
805:       for (l = 0; l < numLevels; ++l) {
806:         SNESFASGetCycleSNES(snes, l, &lsnes);
807:         SNESMonitorSet(lsnes, SNESMonitorError, &user, NULL);
808:       }
809:     } else {
810:       SNESGetKSP(snes, &ksp);
811:       KSPGetPC(ksp, &pc);
812:       PetscObjectTypeCompare((PetscObject) pc, PCMG, &isMG);
813:       if (isMG) {
814:         user.pcmg = pc;
815:         PCMGGetLevels(pc, &numLevels);
816:         for (l = 0; l < numLevels; ++l) {
817:           PCMGGetSmootherDown(pc, l, &ksp);
818:           KSPMonitorSet(ksp, KSPMonitorError, &user, NULL);
819:         }
820:       }
821:     }
822:   }
823:   if (user.runType == RUN_FULL) {
824:     PetscErrorCode (*initialGuess[1])(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {zero};

826:     DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
827:     if (user.debug) {
828:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
829:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
830:     }
831:     SNESSolve(snes, NULL, u);
832:     SNESGetIterationNumber(snes, &its);
833:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
834:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
835:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
836:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
837:     if (user.showSolution) {
838:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
839:       VecChop(u, 3.0e-9);
840:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
841:     }
842:   } else if (user.runType == RUN_PERF) {
843:     PetscReal res = 0.0;

845:     SNESComputeFunction(snes, u, r);
846:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
847:     VecChop(r, 1.0e-10);
848:     VecNorm(r, NORM_2, &res);
849:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
850:   } else {
851:     PetscReal res = 0.0;

853:     /* Check discretization error */
854:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
855:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
856:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
857:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
858:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
859:     /* Check residual */
860:     SNESComputeFunction(snes, u, r);
861:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
862:     VecChop(r, 1.0e-10);
863:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
864:     VecNorm(r, NORM_2, &res);
865:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
866:     /* Check Jacobian */
867:     {
868:       Vec          b;

870:       SNESComputeJacobian(snes, u, A, A);
871:       VecDuplicate(u, &b);
872:       VecSet(r, 0.0);
873:       SNESComputeFunction(snes, r, b);
874:       MatMult(A, u, r);
875:       VecAXPY(r, 1.0, b);
876:       VecDestroy(&b);
877:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
878:       VecChop(r, 1.0e-10);
879:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
880:       VecNorm(r, NORM_2, &res);
881:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
882:     }
883:   }
884:   VecViewFromOptions(u, NULL, "-vec_view");

886:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
887:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
888:   if (A != J) {MatDestroy(&A);}
889:   MatDestroy(&J);
890:   VecDestroy(&u);
891:   VecDestroy(&r);
892:   SNESDestroy(&snes);
893:   DMDestroy(&dm);
894:   PetscFree(user.exactFuncs);
895:   PetscFinalize();
896:   return 0;
897: }