Actual source code: ex12.c

petsc-master 2017-02-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, 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:   PetscBool      fieldBC;
 34:   void           (**exactFields)(PetscInt, PetscInt, PetscInt,
 35:                                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
 36:                                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
 37:                                  PetscReal, const PetscReal[], PetscScalar[]);
 38:   /* Solver */
 39:   PC            pcmg;              /* This is needed for error monitoring */
 40: } AppCtx;

 42: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 43: {
 44:   u[0] = 0.0;
 45:   return 0;
 46: }

 48: static PetscErrorCode ecks(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 49: {
 50:   u[0] = x[0];
 51:   return 0;
 52: }

 54: /*
 55:   In 2D for Dirichlet conditions, we use exact solution:

 57:     u = x^2 + y^2
 58:     f = 4

 60:   so that

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

 64:   For Neumann conditions, we have

 66:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (bottom)
 67:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top)
 68:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
 69:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

 71:   Which we can express as

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

 81: static void quadratic_u_field_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 82:                                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 83:                                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 84:                                  PetscReal t, const PetscReal x[], PetscScalar uexact[])
 85: {
 86:   uexact[0] = a[0];
 87: }

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

 97: void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 98:              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 99:              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
100:              PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
101: {
102:   PetscInt d;
103:   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += -n[d]*2.0*x[d];
104: }

106: void f0_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
107:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
108:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
109:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
110: {
111:   f0[0] = 0.0;
112: }

114: void f1_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
115:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
116:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
117:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
118: {
119:   PetscInt comp;
120:   for (comp = 0; comp < dim; ++comp) f1[comp] = 0.0;
121: }

123: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
124: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
125:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
126:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
127:           PetscReal t, const PetscReal x[], PetscScalar f1[])
128: {
129:   PetscInt d;
130:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
131: }

133: /* < \nabla v, \nabla u + {\nabla u}^T >
134:    This just gives \nabla u, give the perdiagonal for the transpose */
135: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
136:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
137:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
138:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
139: {
140:   PetscInt d;
141:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
142: }

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

147:     u  = x^2 + y^2
148:     f  = 6 (x + y)
149:     nu = (x + y)

151:   so that

153:     -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
154: */
155: static PetscErrorCode nu_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
156: {
157:   *u = x[0] + x[1];
158:   return 0;
159: }

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

169: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
170: void f1_analytic_u(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, const PetscReal x[], PetscScalar f1[])
174: {
175:   PetscInt d;
176:   for (d = 0; d < dim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
177: }

179: void f1_field_u(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, const PetscReal x[], PetscScalar f1[])
183: {
184:   PetscInt d;
185:   for (d = 0; d < dim; ++d) f1[d] = a[0]*u_x[d];
186: }

188: /* < \nabla v, \nabla u + {\nabla u}^T >
189:    This just gives \nabla u, give the perdiagonal for the transpose */
190: void g3_analytic_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
191:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
192:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
193:                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
194: {
195:   PetscInt d;
196:   for (d = 0; d < dim; ++d) g3[d*dim+d] = x[0] + x[1];
197: }

199: void g3_field_uu(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, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
203: {
204:   PetscInt d;
205:   for (d = 0; d < dim; ++d) g3[d*dim+d] = a[0];
206: }

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

211:     u  = x^2 + y^2
212:     f  = 16 (x^2 + y^2)
213:     nu = 1/2 |grad u|^2

215:   so that

217:     -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
218: */
219: void f0_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
220:                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
221:                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
222:                              PetscReal t, const PetscReal x[], PetscScalar f0[])
223: {
224:   f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
225: }

227: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
228: void f1_analytic_nonlinear_u(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, const PetscReal x[], PetscScalar f1[])
232: {
233:   PetscScalar nu = 0.0;
234:   PetscInt    d;
235:   for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
236:   for (d = 0; d < dim; ++d) f1[d] = 0.5*nu*u_x[d];
237: }

239: /*
240:   grad (u + eps w) - grad u = eps grad w

242:   1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
243: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
244: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
245: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
246: */
247: void g3_analytic_nonlinear_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
248:                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
249:                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
250:                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
251: {
252:   PetscScalar nu = 0.0;
253:   PetscInt    d, e;
254:   for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
255:   for (d = 0; d < dim; ++d) {
256:     g3[d*dim+d] = 0.5*nu;
257:     for (e = 0; e < dim; ++e) {
258:       g3[d*dim+e] += u_x[d]*u_x[e];
259:     }
260:   }
261: }

263: /*
264:   In 3D for Dirichlet conditions we use exact solution:

266:     u = 2/3 (x^2 + y^2 + z^2)
267:     f = 4

269:   so that

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

273:   For Neumann conditions, we have

275:     -\nabla u \cdot -\hat z |_{z=0} =  (2z)|_{z=0} =  0 (bottom)
276:     -\nabla u \cdot  \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
277:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (front)
278:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
279:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
280:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

282:   Which we can express as

284:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
285: */
286: static PetscErrorCode quadratic_u_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
287: {
288:   *u = 2.0*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/3.0;
289:   return 0;
290: }

292: static void quadratic_u_field_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
293:                                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
294:                                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
295:                                  PetscReal t, const PetscReal x[], PetscScalar uexact[])
296: {
297:   uexact[0] = a[0];
298: }

300: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
301: {
302:   const char    *bcTypes[3]  = {"neumann", "dirichlet", "none"};
303:   const char    *runTypes[4] = {"full", "exact", "test", "perf"};
304:   const char    *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
305:   PetscInt       bc, run, coeff;
306:   PetscBool      flg;

310:   options->debug               = 0;
311:   options->runType             = RUN_FULL;
312:   options->dim                 = 2;
313:   options->filename[0]         = '\0';
314:   options->interpolate         = PETSC_FALSE;
315:   options->refinementLimit     = 0.0;
316:   options->bcType              = DIRICHLET;
317:   options->variableCoefficient = COEFF_NONE;
318:   options->fieldBC             = PETSC_FALSE;
319:   options->jacobianMF          = PETSC_FALSE;
320:   options->showInitial         = PETSC_FALSE;
321:   options->showSolution        = PETSC_FALSE;
322:   options->restart             = PETSC_FALSE;
323:   options->check               = PETSC_FALSE;
324:   options->viewHierarchy       = PETSC_FALSE;
325:   options->simplex             = PETSC_TRUE;
326:   options->quiet               = PETSC_FALSE;
327:   options->nonzInit            = PETSC_FALSE;

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

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

336:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
337:   PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
338: #if !defined(PETSC_HAVE_EXODUSII)
339:   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
340: #endif
341:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
342:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
343:   bc   = options->bcType;
344:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
345:   options->bcType = (BCType) bc;
346:   coeff = options->variableCoefficient;
347:   PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
348:   options->variableCoefficient = (CoeffType) coeff;

350:   PetscOptionsBool("-field_bc", "Use a field representation for the BC", "ex12.c", options->fieldBC, &options->fieldBC, NULL);
351:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
352:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
353:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
354:   PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
355:   PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
356:   PetscOptionsBool("-dm_view_hierarchy", "View the coarsened hierarchy", "ex12.c", options->viewHierarchy, &options->viewHierarchy, NULL);
357:   PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex12.c", options->simplex, &options->simplex, NULL);
358:   PetscOptionsBool("-quiet", "Don't print any vecs", "ex12.c", options->quiet, &options->quiet, NULL);
359:   PetscOptionsBool("-nonzero_initial_guess", "nonzero intial guess", "ex12.c", options->nonzInit, &options->nonzInit, NULL);
360:   PetscOptionsEnd();
361:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
362:   return(0);
363: }

365: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
366: {
367:   DMLabel        label;

371:   DMCreateLabel(dm, name);
372:   DMGetLabel(dm, name, &label);
373:   DMPlexMarkBoundaryFaces(dm, label);
374:   DMPlexLabelComplete(dm, label);
375:   return(0);
376: }

378: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
379: {
380:   PetscInt       dim             = user->dim;
381:   const char    *filename        = user->filename;
382:   PetscBool      interpolate     = user->interpolate;
383:   PetscReal      refinementLimit = user->refinementLimit;
384:   size_t         len;

388:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
389:   PetscStrlen(filename, &len);
390:   if (!len) {
391:     if (user->simplex) {
392:       DMPlexCreateBoxMesh(comm, dim, dim == 2 ? 2 : 1, interpolate, dm);
393:     }
394:     else {
395:       PetscInt cells[3] = {1, 1, 1}; /* coarse mesh is one cell; refine from there */

397:       DMPlexCreateHexBoxMesh(comm, dim, cells,  DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
398:     }
399:     PetscObjectSetName((PetscObject) *dm, "Mesh");
400:   } else {
401:     DMPlexCreateFromFile(comm, filename, interpolate, dm);
402:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
403:   }
404:   {
405:     DM refinedMesh     = NULL;
406:     DM distributedMesh = NULL;

408:     /* Refine mesh using a volume constraint */
409:     DMPlexSetRefinementLimit(*dm, refinementLimit);
410:     DMRefine(*dm, comm, &refinedMesh);
411:     if (refinedMesh) {
412:       const char *name;

414:       PetscObjectGetName((PetscObject) *dm,         &name);
415:       PetscObjectSetName((PetscObject) refinedMesh,  name);
416:       DMDestroy(dm);
417:       *dm  = refinedMesh;
418:     }
419:     /* Distribute mesh over processes */
420:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
421:     if (distributedMesh) {
422:       DMDestroy(dm);
423:       *dm  = distributedMesh;
424:     }
425:   }
426:   if (user->bcType == NEUMANN) {
427:     DMLabel   label;

429:     DMCreateLabel(*dm, "boundary");
430:     DMGetLabel(*dm, "boundary", &label);
431:     DMPlexMarkBoundaryFaces(*dm, label);
432:   } else if (user->bcType == DIRICHLET) {
433:     PetscBool hasLabel;

435:     DMHasLabel(*dm,"marker",&hasLabel);
436:     if (!hasLabel) {CreateBCLabel(*dm, "marker");}
437:   }
438:   {
439:     char      convType[256];
440:     PetscBool flg;

442:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
443:     PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);
444:     PetscOptionsEnd();
445:     if (flg) {
446:       DM dmConv;

448:       DMConvert(*dm,convType,&dmConv);
449:       if (dmConv) {
450:         DMDestroy(dm);
451:         *dm  = dmConv;
452:       }
453:     }
454:   }
455:   DMSetFromOptions(*dm);
456:   DMViewFromOptions(*dm, NULL, "-dm_view");
457:   if (user->viewHierarchy) {
458:     DM       cdm = *dm;
459:     PetscInt i   = 0;
460:     char     buf[256];

462:     while (cdm) {
463:       DMSetUp(cdm);
464:       DMGetCoarseDM(cdm, &cdm);
465:       ++i;
466:     }
467:     cdm = *dm;
468:     while (cdm) {
469:       PetscViewer       viewer;
470:       PetscBool   isHDF5, isVTK;

472:       --i;
473:       PetscViewerCreate(comm,&viewer);
474:       PetscViewerSetType(viewer,PETSCVIEWERHDF5);
475:       PetscViewerSetOptionsPrefix(viewer,"hierarchy_");
476:       PetscViewerSetFromOptions(viewer);
477:       PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
478:       PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
479:       if (isHDF5) {
480:         PetscSNPrintf(buf, 256, "ex12-%d.h5", i);
481:       }
482:       else if (isVTK) {
483:         PetscSNPrintf(buf, 256, "ex12-%d.vtu", i);
484:         PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
485:       }
486:       else {
487:         PetscSNPrintf(buf, 256, "ex12-%d", i);
488:       }
489:       PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
490:       PetscViewerFileSetName(viewer,buf);
491:       DMView(cdm, viewer);
492:       PetscViewerDestroy(&viewer);
493:       DMGetCoarseDM(cdm, &cdm);
494:     }
495:   }
496:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
497:   return(0);
498: }

500: static PetscErrorCode SetupProblem(PetscDS prob, AppCtx *user)
501: {
502:   const PetscInt id = 1;

506:   switch (user->variableCoefficient) {
507:   case COEFF_NONE:
508:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
509:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
510:     break;
511:   case COEFF_ANALYTIC:
512:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
513:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
514:     break;
515:   case COEFF_FIELD:
516:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
517:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
518:     break;
519:   case COEFF_NONLINEAR:
520:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
521:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
522:     break;
523:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
524:   }
525:   switch (user->dim) {
526:   case 2:
527:     user->exactFuncs[0]  = quadratic_u_2d;
528:     user->exactFields[0] = quadratic_u_field_2d;
529:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
530:     break;
531:   case 3:
532:     user->exactFuncs[0]  = quadratic_u_3d;
533:     user->exactFields[0] = quadratic_u_field_3d;
534:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
535:     break;
536:   default:
537:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
538:   }
539:   PetscDSAddBoundary(prob, user->bcType == DIRICHLET ? (user->fieldBC ? DM_BC_ESSENTIAL_FIELD : DM_BC_ESSENTIAL) : DM_BC_NATURAL,
540:                             "wall", user->bcType == DIRICHLET ? "marker" : "boundary", 0, 0, NULL,
541:                             user->fieldBC ? (void (*)()) user->exactFields[0] : (void (*)()) user->exactFuncs[0], 1, &id, user);
542:   return(0);
543: }

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

552:   DMCreateLocalVector(dmAux, &nu);
553:   DMProjectFunctionLocal(dmAux, 0.0, matFuncs, NULL, INSERT_ALL_VALUES, nu);
554:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
555:   VecDestroy(&nu);
556:   return(0);
557: }

559: static PetscErrorCode SetupBC(DM dm, DM dmAux, AppCtx *user)
560: {
561:   PetscErrorCode (*bcFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
562:   Vec            uexact;
563:   PetscInt       dim;

567:   DMGetDimension(dm, &dim);
568:   if (dim == 2) bcFuncs[0] = quadratic_u_2d;
569:   else          bcFuncs[0] = quadratic_u_3d;
570:   DMCreateLocalVector(dmAux, &uexact);
571:   DMProjectFunctionLocal(dmAux, 0.0, bcFuncs, NULL, INSERT_ALL_VALUES, uexact);
572:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) uexact);
573:   VecDestroy(&uexact);
574:   return(0);
575: }

577: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
578: {
579:   DM             cdm   = dm;
580:   const PetscInt dim   = user->dim;
581:   PetscFE        feAux = NULL;
582:   PetscFE        feCh  = NULL;
583:   PetscFE        fe;
584:   PetscDS        prob, probAux = NULL, probCh = NULL;
585:   PetscBool      simplex = user->simplex;

589:   /* Create finite element */
590:   PetscFECreateDefault(dm, dim, 1, simplex, NULL, -1, &fe);
591:   PetscObjectSetName((PetscObject) fe, "potential");
592:   if (user->variableCoefficient == COEFF_FIELD) {
593:     PetscQuadrature q;

595:     PetscFECreateDefault(dm, dim, 1, simplex, "mat_", -1, &feAux);
596:     PetscFEGetQuadrature(fe, &q);
597:     PetscFESetQuadrature(feAux, q);
598:     PetscDSCreate(PetscObjectComm((PetscObject)dm),&probAux);
599:     PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
600:   } else if (user->fieldBC) {
601:     PetscQuadrature q;

603:     PetscFECreateDefault(dm, dim, 1, simplex, "bc_", -1, &feAux);
604:     PetscFEGetQuadrature(fe, &q);
605:     PetscFESetQuadrature(feAux, q);
606:     PetscDSCreate(PetscObjectComm((PetscObject)dm),&probAux);
607:     PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
608:   }
609:   if (user->check) {
610:     PetscFECreateDefault(dm, dim, 1, simplex, "ch_", -1, &feCh);
611:     PetscDSCreate(PetscObjectComm((PetscObject)dm),&probCh);
612:     PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
613:   }
614:   /* Set discretization and boundary conditions for each mesh */
615:   DMGetDS(dm, &prob);
616:   PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
617:   SetupProblem(prob, user);
618:   while (cdm) {
619:     DM coordDM;

621:     DMSetDS(cdm,prob);
622:     DMGetCoordinateDM(cdm,&coordDM);
623:     if (feAux) {
624:       DM      dmAux;

626:       DMClone(cdm, &dmAux);
627:       DMSetCoordinateDM(dmAux, coordDM);
628:       DMSetDS(dmAux, probAux);
629:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
630:       if (user->fieldBC) {SetupBC(cdm, dmAux, user);}
631:       else               {SetupMaterial(cdm, dmAux, user);}
632:       DMDestroy(&dmAux);
633:     }
634:     if (feCh) {
635:       DM      dmCh;

637:       DMClone(cdm, &dmCh);
638:       DMSetCoordinateDM(dmCh, coordDM);
639:       DMSetDS(dmCh, probCh);
640:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
641:       DMDestroy(&dmCh);
642:     }
643:     if (user->bcType == DIRICHLET) {
644:       PetscBool hasLabel;

646:       DMHasLabel(cdm, "marker", &hasLabel);
647:       if (!hasLabel) {CreateBCLabel(cdm, "marker");}
648:     }
649:     DMGetCoarseDM(cdm, &cdm);
650:   }
651:   PetscFEDestroy(&fe);
652:   PetscFEDestroy(&feAux);
653:   PetscFEDestroy(&feCh);
654:   PetscDSDestroy(&probAux);
655:   PetscDSDestroy(&probCh);
656:   return(0);
657: }

659: #include "petsc/private/petscimpl.h"

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

664:   Collective on KSP

666:   Input Parameters:
667: + ksp   - the KSP
668: . its   - iteration number
669: . rnorm - 2-norm, preconditioned residual value (may be estimated).
670: - ctx   - monitor context

672:   Level: intermediate

674: .keywords: KSP, default, monitor, residual
675: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorDefault()
676: @*/
677: static PetscErrorCode KSPMonitorError(KSP ksp, PetscInt its, PetscReal rnorm, void *ctx)
678: {
679:   AppCtx        *user = (AppCtx *) ctx;
680:   DM             dm;
681:   Vec            du, r;
682:   PetscInt       level = 0;
683:   PetscBool      hasLevel;
684:   PetscViewer    viewer;
685:   char           buf[256];

689:   KSPGetDM(ksp, &dm);
690:   /* Calculate solution */
691:   {
692:     PC        pc = user->pcmg; /* The MG PC */
693:     DM        fdm,  cdm;
694:     KSP       fksp, cksp;
695:     Vec       fu,   cu;
696:     PetscInt  levels, l;

698:     KSPBuildSolution(ksp, NULL, &du);
699:     PetscObjectComposedDataGetInt((PetscObject) ksp, PetscMGLevelId, level, hasLevel);
700:     PCMGGetLevels(pc, &levels);
701:     PCMGGetSmoother(pc, levels-1, &fksp);
702:     KSPBuildSolution(fksp, NULL, &fu);
703:     for (l = levels-1; l > level; --l) {
704:       Mat R;
705:       Vec s;

707:       PCMGGetSmoother(pc, l-1, &cksp);
708:       KSPGetDM(cksp, &cdm);
709:       DMGetGlobalVector(cdm, &cu);
710:       PCMGGetRestriction(pc, l, &R);
711:       PCMGGetRScale(pc, l, &s);
712:       MatRestrict(R, fu, cu);
713:       VecPointwiseMult(cu, cu, s);
714:       if (l < levels-1) {DMRestoreGlobalVector(fdm, &fu);}
715:       fdm  = cdm;
716:       fu   = cu;
717:     }
718:     if (levels-1 > level) {
719:       VecAXPY(du, 1.0, cu);
720:       DMRestoreGlobalVector(cdm, &cu);
721:     }
722:   }
723:   /* Calculate error */
724:   DMGetGlobalVector(dm, &r);
725:   DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
726:   VecAXPY(r,-1.0,du);
727:   PetscObjectSetName((PetscObject) r, "solution error");
728:   /* View error */
729:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
730:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
731:   VecView(r, viewer);
732:   /* Cleanup */
733:   PetscViewerDestroy(&viewer);
734:   DMRestoreGlobalVector(dm, &r);
735:   return(0);
736: }

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

741:   Collective on SNES

743:   Input Parameters:
744: + snes  - the SNES
745: . its   - iteration number
746: . rnorm - 2-norm of residual
747: - ctx   - user context

749:   Level: intermediate

751: .keywords: SNES, nonlinear, default, monitor, norm
752: .seealso: SNESMonitorDefault(), SNESMonitorSet(), SNESMonitorSolution()
753: @*/
754: static PetscErrorCode SNESMonitorError(SNES snes, PetscInt its, PetscReal rnorm, void *ctx)
755: {
756:   AppCtx        *user = (AppCtx *) ctx;
757:   DM             dm;
758:   Vec            u, r;
759:   PetscInt       level;
760:   PetscBool      hasLevel;
761:   PetscViewer    viewer;
762:   char           buf[256];

766:   SNESGetDM(snes, &dm);
767:   /* Calculate error */
768:   SNESGetSolution(snes, &u);
769:   DMGetGlobalVector(dm, &r);
770:   PetscObjectSetName((PetscObject) r, "solution error");
771:   DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
772:   VecAXPY(r, -1.0, u);
773:   /* View error */
774:   PetscObjectComposedDataGetInt((PetscObject) snes, PetscMGLevelId, level, hasLevel);
775:   PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
776:   PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
777:   VecView(r, viewer);
778:   /* Cleanup */
779:   PetscViewerDestroy(&viewer);
780:   DMRestoreGlobalVector(dm, &r);
781:   return(0);
782: }

784: int main(int argc, char **argv)
785: {
786:   DM             dm;          /* Problem specification */
787:   SNES           snes;        /* nonlinear solver */
788:   Vec            u,r;         /* solution, residual vectors */
789:   Mat            A,J;         /* Jacobian matrix */
790:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
791:   AppCtx         user;        /* user-defined work context */
792:   JacActionCtx   userJ;       /* context for Jacobian MF action */
793:   PetscInt       its;         /* iterations for convergence */
794:   PetscReal      error = 0.0; /* L_2 error in the solution */
795:   PetscBool      isFAS;

798:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
799:   ProcessOptions(PETSC_COMM_WORLD, &user);
800:   SNESCreate(PETSC_COMM_WORLD, &snes);
801:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
802:   SNESSetDM(snes, dm);
803:   DMSetApplicationContext(dm, &user);

805:   PetscMalloc2(1, &user.exactFuncs, 1, &user.exactFields);
806:   SetupDiscretization(dm, &user);

808:   DMCreateGlobalVector(dm, &u);
809:   PetscObjectSetName((PetscObject) u, "potential");
810:   VecDuplicate(u, &r);

812:   DMCreateMatrix(dm, &J);
813:   if (user.jacobianMF) {
814:     PetscInt M, m, N, n;

816:     MatGetSize(J, &M, &N);
817:     MatGetLocalSize(J, &m, &n);
818:     MatCreate(PETSC_COMM_WORLD, &A);
819:     MatSetSizes(A, m, n, M, N);
820:     MatSetType(A, MATSHELL);
821:     MatSetUp(A);
822: #if 0
823:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
824: #endif

826:     userJ.dm   = dm;
827:     userJ.J    = J;
828:     userJ.user = &user;

830:     DMCreateLocalVector(dm, &userJ.u);
831:     if (user.fieldBC) {DMProjectFieldLocal(dm, 0.0, userJ.u, user.exactFields, INSERT_BC_VALUES, userJ.u);}
832:     else              {DMProjectFunctionLocal(dm, 0.0, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);}
833:     MatShellSetContext(A, &userJ);
834:   } else {
835:     A = J;
836:   }
837:   if (user.bcType == NEUMANN) {
838:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
839:     MatSetNullSpace(A, nullSpace);
840:   }

842:   DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
843:   SNESSetJacobian(snes, A, J, NULL, NULL);

845:   SNESSetFromOptions(snes);

847:   if (user.fieldBC) {DMProjectField(dm, 0.0, u, user.exactFields, INSERT_ALL_VALUES, u);}
848:   else              {DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);}
849:   if (user.restart) {
850: #if defined(PETSC_HAVE_HDF5)
851:     PetscViewer viewer;

853:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
854:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
855:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
856:     PetscViewerFileSetName(viewer, user.filename);
857:     PetscViewerHDF5PushGroup(viewer, "/fields");
858:     VecLoad(u, viewer);
859:     PetscViewerHDF5PopGroup(viewer);
860:     PetscViewerDestroy(&viewer);
861: #endif
862:   }
863:   if (user.showInitial) {
864:     Vec lv;
865:     DMGetLocalVector(dm, &lv);
866:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
867:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
868:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
869:     DMRestoreLocalVector(dm, &lv);
870:   }
871:   if (user.viewHierarchy) {
872:     SNES      lsnes;
873:     KSP       ksp;
874:     PC        pc;
875:     PetscInt  numLevels, l;
876:     PetscBool isMG;

878:     PetscObjectTypeCompare((PetscObject) snes, SNESFAS, &isFAS);
879:     if (isFAS) {
880:       SNESFASGetLevels(snes, &numLevels);
881:       for (l = 0; l < numLevels; ++l) {
882:         SNESFASGetCycleSNES(snes, l, &lsnes);
883:         SNESMonitorSet(lsnes, SNESMonitorError, &user, NULL);
884:       }
885:     } else {
886:       SNESGetKSP(snes, &ksp);
887:       KSPGetPC(ksp, &pc);
888:       PetscObjectTypeCompare((PetscObject) pc, PCMG, &isMG);
889:       if (isMG) {
890:         user.pcmg = pc;
891:         PCMGGetLevels(pc, &numLevels);
892:         for (l = 0; l < numLevels; ++l) {
893:           PCMGGetSmootherDown(pc, l, &ksp);
894:           KSPMonitorSet(ksp, KSPMonitorError, &user, NULL);
895:         }
896:       }
897:     }
898:   }
899:   if (user.runType == RUN_FULL || user.runType == RUN_EXACT) {
900:     PetscErrorCode (*initialGuess[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {zero};

902:     if (user.nonzInit) initialGuess[0] = ecks;
903:     if (user.runType == RUN_FULL) {
904:       DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
905:     }
906:     if (user.debug) {
907:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
908:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
909:     }
910:     SNESSolve(snes, NULL, u);
911:     SNESGetIterationNumber(snes, &its);
912:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
913:     DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
914:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
915:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
916:     if (user.showSolution) {
917:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
918:       VecChop(u, 3.0e-9);
919:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
920:     }
921:   } else if (user.runType == RUN_PERF) {
922:     PetscReal res = 0.0;

924:     SNESComputeFunction(snes, u, r);
925:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
926:     VecChop(r, 1.0e-10);
927:     VecNorm(r, NORM_2, &res);
928:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
929:   } else {
930:     PetscReal res = 0.0;

932:     /* Check discretization error */
933:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
934:     if (!user.quiet) {VecView(u, PETSC_VIEWER_STDOUT_WORLD);}
935:     DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
936:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
937:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
938:     /* Check residual */
939:     SNESComputeFunction(snes, u, r);
940:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
941:     VecChop(r, 1.0e-10);
942:     if (!user.quiet) {VecView(r, PETSC_VIEWER_STDOUT_WORLD);}
943:     VecNorm(r, NORM_2, &res);
944:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
945:     /* Check Jacobian */
946:     {
947:       Vec          b;

949:       SNESComputeJacobian(snes, u, A, A);
950:       VecDuplicate(u, &b);
951:       VecSet(r, 0.0);
952:       SNESComputeFunction(snes, r, b);
953:       MatMult(A, u, r);
954:       VecAXPY(r, 1.0, b);
955:       VecDestroy(&b);
956:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
957:       VecChop(r, 1.0e-10);
958:       if (!user.quiet) {VecView(r, PETSC_VIEWER_STDOUT_WORLD);}
959:       VecNorm(r, NORM_2, &res);
960:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
961:     }
962:   }
963:   VecViewFromOptions(u, NULL, "-vec_view");

965:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
966:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
967:   if (A != J) {MatDestroy(&A);}
968:   MatDestroy(&J);
969:   VecDestroy(&u);
970:   VecDestroy(&r);
971:   SNESDestroy(&snes);
972:   DMDestroy(&dm);
973:   PetscFree2(user.exactFuncs, user.exactFields);
974:   PetscFinalize();
975:   return ierr;
976: }