Actual source code: ex62.c

petsc-3.5.1 2014-08-01
Report Typos and Errors
  1: static char help[] = "Stokes Problem in 2d and 3d with simplicial finite elements.\n\
  2: We solve the Stokes problem in a rectangular\n\
  3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";

  5: /*
  6: The isoviscous Stokes problem, which we discretize using the finite
  7: element method on an unstructured mesh. The weak form equations are

  9:   < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > + < v, f > = 0
 10:   < q, \nabla\cdot u >                                                    = 0

 12: We start with homogeneous Dirichlet conditions. We will expand this as the set
 13: of test problems is developed.

 15: Discretization:

 17: We use PetscFE to generate a tabulation of the finite element basis functions
 18: at quadrature points. We can currently generate an arbitrary order Lagrange
 19: element.

 21: Field Data:

 23:   DMPLEX data is organized by point, and the closure operation just stacks up the
 24: data from each sieve point in the closure. Thus, for a P_2-P_1 Stokes element, we
 25: have

 27:   cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
 28:   x     = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} p_{v_0} u_{v_1} v_{v_1} p_{v_1} u_{v_2} v_{v_2} p_{v_2}]

 30: The problem here is that we would like to loop over each field separately for
 31: integration. Therefore, the closure visitor in DMPlexVecGetClosure() reorders
 32: the data so that each field is contiguous

 34:   x'    = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} u_{v_1} v_{v_1} u_{v_2} v_{v_2} p_{v_0} p_{v_1} p_{v_2}]

 36: Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
 37: puts it into the Sieve ordering.

 39: Next Steps:

 41: - Refine and show convergence of correct order automatically (use femTest.py)
 42: - Fix InitialGuess for arbitrary disc (means making dual application work again)
 43: - Redo slides from GUCASTutorial for this new example

 45: For tensor product meshes, see SNES ex67, ex72
 46: */

 48: #include <petscdmplex.h>
 49: #include <petscsnes.h>
 50: #include <petscds.h>

 52: PetscInt spatialDim = 0;

 54: typedef enum {NEUMANN, DIRICHLET} BCType;
 55: typedef enum {RUN_FULL, RUN_TEST} RunType;

 57: typedef struct {
 58:   PetscInt      debug;             /* The debugging level */
 59:   RunType       runType;           /* Whether to run tests, or solve the full problem */
 60:   PetscLogEvent createMeshEvent;
 61:   PetscBool     showInitial, showSolution, showError;
 62:   /* Domain and mesh definition */
 63:   PetscInt      dim;               /* The topological mesh dimension */
 64:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 65:   PetscBool     simplex;           /* Use simplices or tensor product cells */
 66:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 67:   char          partitioner[2048]; /* The graph partitioner */
 68:   /* Problem definition */
 69:   BCType        bcType;
 70:   void       (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
 71: } AppCtx;

 73: void zero_scalar(const PetscReal coords[], PetscScalar *u, void *ctx)
 74: {
 75:   u[0] = 0.0;
 76: }
 77: void zero_vector(const PetscReal coords[], PetscScalar *u, void *ctx)
 78: {
 79:   PetscInt d;
 80:   for (d = 0; d < spatialDim; ++d) u[d] = 0.0;
 81: }

 83: /*
 84:   In 2D we use exact solution:

 86:     u = x^2 + y^2
 87:     v = 2 x^2 - 2xy
 88:     p = x + y - 1
 89:     f_x = f_y = 3

 91:   so that

 93:     -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
 94:     \nabla \cdot u           = 2x - 2x                    = 0
 95: */
 96: void quadratic_u_2d(const PetscReal x[], PetscScalar *u, void *ctx)
 97: {
 98:   u[0] = x[0]*x[0] + x[1]*x[1];
 99:   u[1] = 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
100: }

102: void linear_p_2d(const PetscReal x[], PetscScalar *p, void *ctx)
103: {
104:   *p = x[0] + x[1] - 1.0;
105: }
106: void constant_p(const PetscReal x[], PetscScalar *p, void *ctx)
107: {
108:   *p = 1.0;
109: }

111: void f0_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
112: {
113:   PetscInt c;
114:   for (c = 0; c < spatialDim; ++c) f0[c] = 3.0;
115: }

117: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z}
118:    u[Ncomp]          = {p} */
119: void f1_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
120: {
121:   const PetscInt dim   = spatialDim;
122:   const PetscInt Ncomp = spatialDim;
123:   PetscInt       comp, d;

125:   for (comp = 0; comp < Ncomp; ++comp) {
126:     for (d = 0; d < dim; ++d) {
127:       /* f1[comp*dim+d] = 0.5*(gradU[comp*dim+d] + gradU[d*dim+comp]); */
128:       f1[comp*dim+d] = u_x[comp*dim+d];
129:     }
130:     f1[comp*dim+comp] -= u[Ncomp];
131:   }
132: }

134: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z} */
135: void f0_p(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
136: {
137:   const PetscInt dim = spatialDim;
138:   PetscInt       d;

140:   f0[0] = 0.0;
141:   for (d = 0; d < dim; ++d) f0[0] += u_x[d*dim+d];
142: }

144: void f1_p(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
145: {
146:   PetscInt d;
147:   for (d = 0; d < spatialDim; ++d) f1[d] = 0.0;
148: }

150: /* < q, \nabla\cdot u >
151:    NcompI = 1, NcompJ = dim */
152: void g1_pu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[])
153: {
154:   const PetscInt dim = spatialDim;
155:   PetscInt       d;

157:   for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
158: }

160: /* -< \nabla\cdot v, p >
161:     NcompI = dim, NcompJ = 1 */
162: void g2_up(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[])
163: {
164:   const PetscInt dim = spatialDim;
165:   PetscInt       d;

167:   for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
168: }

170: /* < \nabla v, \nabla u + {\nabla u}^T >
171:    This just gives \nabla u, give the perdiagonal for the transpose */
172: void g3_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
173: {
174:   const PetscInt dim   = spatialDim;
175:   const PetscInt Ncomp = spatialDim;
176:   PetscInt       compI, d;

178:   for (compI = 0; compI < Ncomp; ++compI) {
179:     for (d = 0; d < dim; ++d) {
180:       g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
181:     }
182:   }
183: }

185: /*
186:   In 3D we use exact solution:

188:     u = x^2 + y^2
189:     v = y^2 + z^2
190:     w = x^2 + y^2 - 2(x+y)z
191:     p = x + y + z - 3/2
192:     f_x = f_y = f_z = 3

194:   so that

196:     -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
197:     \nabla \cdot u           = 2x + 2y - 2(x + y)                   = 0
198: */
199: void quadratic_u_3d(const PetscReal x[], PetscScalar *u, void *ctx)
200: {
201:   u[0] = x[0]*x[0] + x[1]*x[1];
202:   u[1] = x[1]*x[1] + x[2]*x[2];
203:   u[2] = x[0]*x[0] + x[1]*x[1] - 2.0*(x[0] + x[1])*x[2];
204: }

206: void linear_p_3d(const PetscReal x[], PetscScalar *p, void *ctx)
207: {
208:   *p = x[0] + x[1] + x[2] - 1.5;
209: }

213: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
214: {
215:   const char    *bcTypes[2]  = {"neumann", "dirichlet"};
216:   const char    *runTypes[2] = {"full", "test"};
217:   PetscInt       bc, run;

221:   options->debug           = 0;
222:   options->runType         = RUN_FULL;
223:   options->dim             = 2;
224:   options->interpolate     = PETSC_FALSE;
225:   options->simplex         = PETSC_TRUE;
226:   options->refinementLimit = 0.0;
227:   options->bcType          = DIRICHLET;
228:   options->showInitial     = PETSC_FALSE;
229:   options->showSolution    = PETSC_TRUE;
230:   options->showError       = PETSC_FALSE;

232:   PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
233:   PetscOptionsInt("-debug", "The debugging level", "ex62.c", options->debug, &options->debug, NULL);
234:   run  = options->runType;
235:   PetscOptionsEList("-run_type", "The run type", "ex62.c", runTypes, 2, runTypes[options->runType], &run, NULL);

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

239:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex62.c", options->dim, &options->dim, NULL);
240:   spatialDim = options->dim;
241:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex62.c", options->interpolate, &options->interpolate, NULL);
242:   PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex62.c", options->simplex, &options->simplex, NULL);
243:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex62.c", options->refinementLimit, &options->refinementLimit, NULL);
244:   PetscStrcpy(options->partitioner, "chaco");
245:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
246:   bc   = options->bcType;
247:   PetscOptionsEList("-bc_type","Type of boundary condition","ex62.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);

249:   options->bcType = (BCType) bc;

251:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex62.c", options->showInitial, &options->showInitial, NULL);
252:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex62.c", options->showSolution, &options->showSolution, NULL);
253:   PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL);
254:   PetscOptionsEnd();

256:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
257:   return(0);
258: }

262: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
263: {
264:   Vec            lv;
265:   PetscInt       p;
266:   PetscMPIInt    rank, numProcs;

270:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
271:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
272:   DMGetLocalVector(dm, &lv);
273:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
274:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
275:   PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
276:   for (p = 0; p < numProcs; ++p) {
277:     if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
278:     PetscBarrier((PetscObject) dm);
279:   }
280:   DMRestoreLocalVector(dm, &lv);
281:   return(0);
282: }

286: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
287: {
288:   DMLabel        label;
289:   PetscInt       dim             = user->dim;
290:   PetscBool      interpolate     = user->interpolate;
291:   PetscReal      refinementLimit = user->refinementLimit;
292:   const char     *partitioner    = user->partitioner;
293:   const PetscInt cells[3]        = {3, 3, 3};

297:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
298:   if (user->simplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
299:   else               {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);}
300:   DMPlexGetLabel(*dm, "marker", &label);
301:   if (label) {DMPlexLabelComplete(*dm, label);}
302:   {
303:     DM refinedMesh     = NULL;
304:     DM distributedMesh = NULL;

306:     /* Refine mesh using a volume constraint */
307:     DMPlexSetRefinementLimit(*dm, refinementLimit);
308:     if (user->simplex) {DMRefine(*dm, comm, &refinedMesh);}
309:     if (refinedMesh) {
310:       DMDestroy(dm);
311:       *dm  = refinedMesh;
312:     }
313:     /* Distribute mesh over processes */
314:     DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
315:     if (distributedMesh) {
316:       DMDestroy(dm);
317:       *dm  = distributedMesh;
318:     }
319:   }
320:   DMSetFromOptions(*dm);
321:   DMViewFromOptions(*dm, NULL, "-dm_view");
322:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
323:   return(0);
324: }

328: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
329: {
330:   PetscDS        prob;

334:   DMGetDS(dm, &prob);
335:   PetscDSSetResidual(prob, 0, f0_u, f1_u);
336:   PetscDSSetResidual(prob, 1, f0_p, f1_p);
337:   PetscDSSetJacobian(prob, 0, 0, NULL, NULL,  NULL,  g3_uu);
338:   PetscDSSetJacobian(prob, 0, 1, NULL, NULL,  g2_up, NULL);
339:   PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu, NULL,  NULL);
340:   switch (user->dim) {
341:   case 2:
342:     user->exactFuncs[0] = quadratic_u_2d;
343:     user->exactFuncs[1] = linear_p_2d;
344:     break;
345:   case 3:
346:     user->exactFuncs[0] = quadratic_u_3d;
347:     user->exactFuncs[1] = linear_p_3d;
348:     break;
349:   default:
350:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
351:   }
352:   return(0);
353: }

357: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
358: {
359:   DM             cdm   = dm;
360:   const PetscInt dim   = user->dim;
361:   const PetscInt id    = 1;
362:   PetscFE        fe[2];
363:   PetscSpace     P;
364:   PetscDS        prob;
365:   PetscInt       order;

369:   /* Create finite element */
370:   PetscFECreateDefault(dm, dim, dim, user->simplex, "vel_", -1, &fe[0]);
371:   PetscObjectSetName((PetscObject) fe[0], "velocity");
372:   PetscFEGetBasisSpace(fe[0], &P);
373:   PetscSpaceGetOrder(P, &order);
374:   PetscFECreateDefault(dm, dim, 1, user->simplex, "pres_", order, &fe[1]);
375:   PetscObjectSetName((PetscObject) fe[1], "pressure");
376:   /* Set discretization and boundary conditions for each mesh */
377:   while (cdm) {
378:     DMGetDS(cdm, &prob);
379:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);
380:     PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);
381:     SetupProblem(cdm, user);
382:     DMPlexAddBoundary(cdm, user->bcType == DIRICHLET, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1, &id, user);
383:     DMPlexGetCoarseDM(cdm, &cdm);
384:   }
385:   PetscFEDestroy(&fe[0]);
386:   PetscFEDestroy(&fe[1]);
387:   return(0);
388: }

392: PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, Vec *v, MatNullSpace *nullSpace)
393: {
394:   Vec            vec;
395:   void         (*funcs[2])(const PetscReal x[], PetscScalar *u, void* ctx) = {zero_vector, constant_p};

399:   DMGetGlobalVector(dm, &vec);
400:   DMPlexProjectFunction(dm, funcs, NULL, INSERT_ALL_VALUES, vec);
401:   VecNormalize(vec, NULL);
402:   if (user->debug) {
403:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
404:     VecView(vec, PETSC_VIEWER_STDOUT_WORLD);
405:   }
406:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
407:   if (v) {
408:     DMCreateGlobalVector(dm, v);
409:     VecCopy(vec, *v);
410:   }
411:   DMRestoreGlobalVector(dm, &vec);
412:   /* New style for field null spaces */
413:   {
414:     PetscObject  pressure;
415:     MatNullSpace nullSpacePres;

417:     DMGetField(dm, 1, &pressure);
418:     MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);
419:     PetscObjectCompose(pressure, "nullspace", (PetscObject) nullSpacePres);
420:     MatNullSpaceDestroy(&nullSpacePres);
421:   }
422:   return(0);
423: }

427: int main(int argc, char **argv)
428: {
429:   SNES           snes;                 /* nonlinear solver */
430:   DM             dm;                   /* problem definition */
431:   Vec            u,r;                  /* solution, residual vectors */
432:   Mat            A,J;                  /* Jacobian matrix */
433:   MatNullSpace   nullSpace;            /* May be necessary for pressure */
434:   AppCtx         user;                 /* user-defined work context */
435:   PetscInt       its;                  /* iterations for convergence */
436:   PetscReal      error         = 0.0;  /* L_2 error in the solution */
437:   PetscReal      ferrors[2];

440:   PetscInitialize(&argc, &argv, NULL, help);
441:   ProcessOptions(PETSC_COMM_WORLD, &user);
442:   SNESCreate(PETSC_COMM_WORLD, &snes);
443:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
444:   SNESSetDM(snes, dm);
445:   DMSetApplicationContext(dm, &user);

447:   PetscMalloc(2 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
448:   SetupDiscretization(dm, &user);
449:   DMPlexCreateClosureIndex(dm, NULL);

451:   DMCreateGlobalVector(dm, &u);
452:   VecDuplicate(u, &r);

454:   DMSetMatType(dm,MATAIJ);
455:   DMCreateMatrix(dm, &J);
456:   A = J;
457:   CreatePressureNullSpace(dm, &user, NULL, &nullSpace);
458:   MatSetNullSpace(J, nullSpace);
459:   if (A != J) {
460:     MatSetNullSpace(A, nullSpace);
461:   }

463:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexSNESComputeResidualFEM,&user);
464:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*))DMPlexSNESComputeJacobianFEM,&user);
465:   SNESSetJacobian(snes, A, J, NULL, NULL);

467:   SNESSetFromOptions(snes);

469:   DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
470:   if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}
471:   if (user.runType == RUN_FULL) {
472:     void (*initialGuess[2])(const PetscReal x[], PetscScalar *u, void* ctx) = {zero_vector, zero_scalar};

474:     DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
475:     if (user.debug) {
476:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
477:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
478:     }
479:     SNESSolve(snes, NULL, u);
480:     SNESGetIterationNumber(snes, &its);
481:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
482:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
483:     DMPlexComputeL2FieldDiff(dm, user.exactFuncs, NULL, u, ferrors);
484:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g [%.3g, %.3g]\n", error, ferrors[0], ferrors[1]);
485:     if (user.showError) {
486:       Vec r;
487:       DMGetGlobalVector(dm, &r);
488:       DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, r);
489:       VecAXPY(r, -1.0, u);
490:       PetscPrintf(PETSC_COMM_WORLD, "Solution Error\n");
491:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
492:       DMRestoreGlobalVector(dm, &r);
493:     }
494:     if (user.showSolution) {
495:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
496:       VecChop(u, 3.0e-9);
497:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
498:     }
499:   } else {
500:     PetscReal res = 0.0;

502:     /* Check discretization error */
503:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
504:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
505:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
506:     if (error >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
507:     else                  {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
508:     /* Check residual */
509:     SNESComputeFunction(snes, u, r);
510:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
511:     VecChop(r, 1.0e-10);
512:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
513:     VecNorm(r, NORM_2, &res);
514:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
515:     /* Check Jacobian */
516:     {
517:       Vec          b;
518:       PetscBool    isNull;

520:       SNESComputeJacobian(snes, u, A, A);
521:       MatNullSpaceTest(nullSpace, J, &isNull);
522:       //if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
523:       VecDuplicate(u, &b);
524:       VecSet(r, 0.0);
525:       SNESComputeFunction(snes, r, b);
526:       MatMult(A, u, r);
527:       VecAXPY(r, 1.0, b);
528:       VecDestroy(&b);
529:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
530:       VecChop(r, 1.0e-10);
531:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
532:       VecNorm(r, NORM_2, &res);
533:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
534:     }
535:   }
536:   VecViewFromOptions(u, NULL, "-sol_vec_view");

538:   MatNullSpaceDestroy(&nullSpace);
539:   if (A != J) {MatDestroy(&A);}
540:   MatDestroy(&J);
541:   VecDestroy(&u);
542:   VecDestroy(&r);
543:   SNESDestroy(&snes);
544:   DMDestroy(&dm);
545:   PetscFree(user.exactFuncs);
546:   PetscFinalize();
547:   return 0;
548: }