Actual source code: ex62.c

petsc-dev 2014-04-14
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 v >                                                    = 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 <petscdt.h>
 50: #include <petscfe.h>
 51: #include <petscsnes.h>

 53: #define NUM_FIELDS 2
 54: PetscInt spatialDim = 0;

 56: typedef enum {NEUMANN, DIRICHLET} BCType;
 57: typedef enum {RUN_FULL, RUN_TEST} RunType;

 59: typedef struct {
 60:   PetscFEM      fem;               /* REQUIRED to use DMPlexComputeResidualFEM() */
 61:   PetscInt      debug;             /* The debugging level */
 62:   PetscMPIInt   rank;              /* The process rank */
 63:   PetscMPIInt   numProcs;          /* The number of processes */
 64:   RunType       runType;           /* Whether to run tests, or solve the full problem */
 65:   PetscBool     jacobianMF;        /* Whether to calculate the Jacobian action on the fly */
 66:   PetscLogEvent createMeshEvent;
 67:   PetscBool     showInitial, showSolution;
 68:   /* Domain and mesh definition */
 69:   PetscInt      dim;               /* The topological mesh dimension */
 70:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 71:   PetscBool     simplex;           /* Use simplices or tensor product cells */
 72:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 73:   char          partitioner[2048]; /* The graph partitioner */
 74:   /* GPU partitioning */
 75:   PetscInt      numBatches;        /* The number of cell batches per kernel */
 76:   PetscInt      numBlocks;         /* The number of concurrent blocks per kernel */
 77:   /* Element quadrature */
 78:   PetscFE       fe[NUM_FIELDS];    /* Element definitions for each field */
 79:   PetscQuadrature q[NUM_FIELDS];
 80:   /* Problem definition */
 81:   void (*f0Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[]); /* f0_u(x,y,z), and f0_p(x,y,z) */
 82:   void (*f1Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[]); /* f1_u(x,y,z), and f1_p(x,y,z) */
 83:   void (*g0Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g0[]); /* g0_uu(x,y,z), g0_up(x,y,z), g0_pu(x,y,z), and g0_pp(x,y,z) */
 84:   void (*g1Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g1[]); /* g1_uu(x,y,z), g1_up(x,y,z), g1_pu(x,y,z), and g1_pp(x,y,z) */
 85:   void (*g2Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g2[]); /* g2_uu(x,y,z), g2_up(x,y,z), g2_pu(x,y,z), and g2_pp(x,y,z) */
 86:   void (*g3Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]); /* g3_uu(x,y,z), g3_up(x,y,z), g3_pu(x,y,z), and g3_pp(x,y,z) */
 87:   void (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx); /* The exact solution function u(x,y,z), v(x,y,z), and p(x,y,z) */
 88:   void (*initialGuess[NUM_FIELDS])(const PetscReal x[], PetscScalar *u, void* ctx);
 89:   BCType bcType;
 90: } AppCtx;

 92: void zero_1d(const PetscReal coords[], PetscScalar *u, void *ctx)
 93: {
 94:   u[0] = 0.0;
 95: }
 96: void zero_2d(const PetscReal coords[], PetscScalar *u, void *ctx)
 97: {
 98:   u[0] = 0.0; u[1] = 0.0;
 99: }
100: void zero_3d(const PetscReal coords[], PetscScalar *u, void *ctx)
101: {
102:   u[0] = 0.0; u[1] = 0.0; u[2] = 0.0;
103: }

105: /*
106:   In 2D we use exact solution:

108:     u = x^2 + y^2
109:     v = 2 x^2 - 2xy
110:     p = x + y - 1
111:     f_x = f_y = 3

113:   so that

115:     -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
116:     \nabla \cdot u           = 2x - 2x                    = 0
117: */
118: void quadratic_u_2d(const PetscReal x[], PetscScalar *u, void *ctx)
119: {
120:   u[0] = x[0]*x[0] + x[1]*x[1];
121:   u[1] = 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
122: }

124: void linear_p_2d(const PetscReal x[], PetscScalar *p, void *ctx)
125: {
126:   *p = x[0] + x[1] - 1.0;
127: }

129: void f0_u(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[])
130: {
131:   const PetscInt Ncomp = spatialDim;
132:   PetscInt       comp;

134:   for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 3.0;
135: }

137: /* 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}
138:    u[Ncomp]          = {p} */
139: void f1_u(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[])
140: {
141:   const PetscInt dim   = spatialDim;
142:   const PetscInt Ncomp = spatialDim;
143:   PetscInt       comp, d;

145:   for (comp = 0; comp < Ncomp; ++comp) {
146:     for (d = 0; d < dim; ++d) {
147:       /* f1[comp*dim+d] = 0.5*(gradU[comp*dim+d] + gradU[d*dim+comp]); */
148:       f1[comp*dim+d] = gradU[comp*dim+d];
149:     }
150:     f1[comp*dim+comp] -= u[Ncomp];
151:   }
152: }

154: /* 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} */
155: void f0_p(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f0[])
156: {
157:   const PetscInt dim = spatialDim;
158:   PetscInt       d;

160:   f0[0] = 0.0;
161:   for (d = 0; d < dim; ++d) f0[0] += gradU[d*dim+d];
162: }

164: void f1_p(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar f1[])
165: {
166:   const PetscInt dim = spatialDim;
167:   PetscInt       d;

169:   for (d = 0; d < dim; ++d) f1[d] = 0.0;
170: }

172: /* < q, \nabla\cdot v >
173:    NcompI = 1, NcompJ = dim */
174: void g1_pu(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g1[])
175: {
176:   const PetscInt dim = spatialDim;
177:   PetscInt       d;

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

182: /* -< \nabla\cdot v, p >
183:     NcompI = dim, NcompJ = 1 */
184: void g2_up(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g2[])
185: {
186:   const PetscInt dim = spatialDim;
187:   PetscInt       d;

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

192: /* < \nabla v, \nabla u + {\nabla u}^T >
193:    This just gives \nabla u, give the perdiagonal for the transpose */
194: void g3_uu(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
195: {
196:   const PetscInt dim   = spatialDim;
197:   const PetscInt Ncomp = spatialDim;
198:   PetscInt       compI, d;

200:   for (compI = 0; compI < Ncomp; ++compI) {
201:     for (d = 0; d < dim; ++d) {
202:       g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
203:     }
204:   }
205: }

207: /*
208:   In 3D we use exact solution:

210:     u = x^2 + y^2
211:     v = y^2 + z^2
212:     w = x^2 + y^2 - 2(x+y)z
213:     p = x + y + z - 3/2
214:     f_x = f_y = f_z = 3

216:   so that

218:     -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
219:     \nabla \cdot u           = 2x + 2y - 2(x + y)                   = 0
220: */
221: void quadratic_u_3d(const PetscReal x[], PetscScalar *u, void *ctx)
222: {
223:   u[0] = x[0]*x[0] + x[1]*x[1];
224:   u[1] = x[1]*x[1] + x[2]*x[2];
225:   u[2] = x[0]*x[0] + x[1]*x[1] - 2.0*(x[0] + x[1])*x[2];
226: }

228: void linear_p_3d(const PetscReal x[], PetscScalar *p, void *ctx)
229: {
230:   *p = x[0] + x[1] + x[2] - 1.5;
231: }

235: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
236: {
237:   const char    *bcTypes[2]  = {"neumann", "dirichlet"};
238:   const char    *runTypes[2] = {"full", "test"};
239:   PetscInt       bc, run;

243:   options->debug           = 0;
244:   options->runType         = RUN_FULL;
245:   options->dim             = 2;
246:   options->interpolate     = PETSC_FALSE;
247:   options->simplex         = PETSC_TRUE;
248:   options->refinementLimit = 0.0;
249:   options->bcType          = DIRICHLET;
250:   options->numBatches      = 1;
251:   options->numBlocks       = 1;
252:   options->jacobianMF      = PETSC_FALSE;
253:   options->showInitial     = PETSC_FALSE;
254:   options->showSolution    = PETSC_TRUE;

256:   options->fem.f0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f0Funcs;
257:   options->fem.f1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f1Funcs;
258:   options->fem.g0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g0Funcs;
259:   options->fem.g1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g1Funcs;
260:   options->fem.g2Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g2Funcs;
261:   options->fem.g3Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g3Funcs;

263:   MPI_Comm_size(comm, &options->numProcs);
264:   MPI_Comm_rank(comm, &options->rank);
265:   PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
266:   PetscOptionsInt("-debug", "The debugging level", "ex62.c", options->debug, &options->debug, NULL);
267:   run  = options->runType;
268:   PetscOptionsEList("-run_type", "The run type", "ex62.c", runTypes, 2, runTypes[options->runType], &run, NULL);

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

272:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex62.c", options->dim, &options->dim, NULL);
273:   spatialDim = options->dim;
274:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex62.c", options->interpolate, &options->interpolate, NULL);
275:   PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex62.c", options->simplex, &options->simplex, NULL);
276:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex62.c", options->refinementLimit, &options->refinementLimit, NULL);
277:   PetscStrcpy(options->partitioner, "chaco");
278:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
279:   bc   = options->bcType;
280:   PetscOptionsEList("-bc_type","Type of boundary condition","ex62.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);

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

284:   PetscOptionsInt("-gpu_batches", "The number of cell batches per kernel", "ex62.c", options->numBatches, &options->numBatches, NULL);
285:   PetscOptionsInt("-gpu_blocks", "The number of concurrent blocks per kernel", "ex62.c", options->numBlocks, &options->numBlocks, NULL);
286:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex62.c", options->jacobianMF, &options->jacobianMF, NULL);
287:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex62.c", options->showInitial, &options->showInitial, NULL);
288:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex62.c", options->showSolution, &options->showSolution, NULL);
289:   PetscOptionsEnd();

291:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
292:   return(0);
293: }

297: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
298: {
299:   Vec            lv;
300:   PetscInt       p;
301:   PetscMPIInt    rank, numProcs;

305:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
306:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
307:   DMGetLocalVector(dm, &lv);
308:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
309:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
310:   PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
311:   for (p = 0; p < numProcs; ++p) {
312:     if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
313:     PetscBarrier((PetscObject) dm);
314:   }
315:   DMRestoreLocalVector(dm, &lv);
316:   return(0);
317: }

321: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
322: {
323:   DMLabel        label;
324:   PetscInt       dim             = user->dim;
325:   PetscBool      interpolate     = user->interpolate;
326:   PetscReal      refinementLimit = user->refinementLimit;
327:   const char     *partitioner    = user->partitioner;
328:   const PetscInt cells[3]        = {3, 3, 3};

332:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
333:   if (user->simplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
334:   else               {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);}
335:   DMPlexGetLabel(*dm, "marker", &label);
336:   if (label) {DMPlexLabelComplete(*dm, label);}
337:   {
338:     DM refinedMesh     = NULL;
339:     DM distributedMesh = NULL;

341:     /* Refine mesh using a volume constraint */
342:     DMPlexSetRefinementLimit(*dm, refinementLimit);
343:     DMRefine(*dm, comm, &refinedMesh);
344:     if (refinedMesh) {
345:       DMDestroy(dm);
346:       *dm  = refinedMesh;
347:     }
348:     /* Distribute mesh over processes */
349:     DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
350:     if (distributedMesh) {
351:       DMDestroy(dm);
352:       *dm  = distributedMesh;
353:     }
354:   }
355:   DMSetFromOptions(*dm);
356:   DMViewFromOptions(*dm, NULL, "-dm_view");
357:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
358:   return(0);
359: }

363: PetscErrorCode SetupElement(DM dm, AppCtx *user)
364: {
365:   const PetscInt  dim = user->dim, numFields = 2;
366:   const char     *prefix[2] = {"vel_", "pres_"};
367:   PetscInt        qorder    = 0, f;
368:   PetscErrorCode  ierr;

371:   for (f = 0; f < numFields; ++f) {
372:     PetscFE         fem;
373:     DM              K;
374:     PetscSpace      P;
375:     PetscDualSpace  Q;
376:     PetscInt        order;

378:     /* Create space */
379:     PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
380:     PetscObjectSetOptionsPrefix((PetscObject) P, prefix[f]);
381:     PetscSpaceSetFromOptions(P);
382:     PetscSpacePolynomialSetNumVariables(P, dim);
383:     PetscSpaceSetUp(P);
384:     PetscSpaceGetOrder(P, &order);
385:     qorder = PetscMax(qorder, order);
386:     /* Create dual space */
387:     PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
388:     PetscObjectSetOptionsPrefix((PetscObject) Q, prefix[f]);
389:     PetscDualSpaceCreateReferenceCell(Q, dim, PETSC_TRUE, &K);
390:     PetscDualSpaceSetDM(Q, K);
391:     DMDestroy(&K);
392:     PetscDualSpaceSetOrder(Q, order);
393:     PetscDualSpaceSetFromOptions(Q);
394:     PetscDualSpaceSetUp(Q);
395:     /* Create element */
396:     PetscFECreate(PetscObjectComm((PetscObject) dm), &fem);
397:     PetscObjectSetOptionsPrefix((PetscObject) fem, prefix[f]);
398:     PetscFESetFromOptions(fem);
399:     PetscFESetBasisSpace(fem, P);
400:     PetscFESetDualSpace(fem, Q);
401:     PetscFESetNumComponents(fem, f ? 1 : dim);
402:     PetscFESetUp(fem);

404:     PetscSpaceDestroy(&P);
405:     PetscDualSpaceDestroy(&Q);
406:     user->fe[f] = fem;
407:   }
408:   for (f = 0; f < numFields; ++f) {
409:     PetscQuadrature q;

411:     /* Create quadrature */
412:     PetscDTGaussJacobiQuadrature(dim, qorder, -1.0, 1.0, &q);
413:     PetscFESetQuadrature(user->fe[f], q);
414:     PetscQuadratureDestroy(&q);
415:   }
416:   user->fem.fe    = user->fe;
417:   user->fem.feBd  = NULL;
418:   user->fem.feAux = NULL;
419:   return(0);
420: }

424: PetscErrorCode DestroyElement(AppCtx *user)
425: {
426:   PetscInt       numFields = 2, f;

430:   for (f = 0; f < numFields; ++f) {
431:     PetscFEDestroy(&user->fe[f]);
432:   }
433:   return(0);
434: }

438: PetscErrorCode SetupSection(DM dm, AppCtx *user)
439: {
440:   const PetscInt id = 1;

444:   PetscObjectSetName((PetscObject) user->fe[0], "velocity");
445:   PetscObjectSetName((PetscObject) user->fe[1], "pressure");
446:   DMSetNumFields(dm, 2);
447:   DMSetField(dm, 0, user->fe[0]);
448:   DMSetField(dm, 1, user->fe[1]);
449:   DMPlexAddBoundary(dm, user->bcType == DIRICHLET, user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1, &id, user);
450:   return(0);
451: }

455: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
456: {
457:   PetscFEM *fem = &user->fem;

460:   fem->f0Funcs[0] = f0_u;
461:   fem->f0Funcs[1] = f0_p;
462:   fem->f1Funcs[0] = f1_u;
463:   fem->f1Funcs[1] = f1_p;
464:   fem->g0Funcs[0] = NULL;
465:   fem->g0Funcs[1] = NULL;
466:   fem->g0Funcs[2] = NULL;
467:   fem->g0Funcs[3] = NULL;
468:   fem->g1Funcs[0] = NULL;
469:   fem->g1Funcs[1] = NULL;
470:   fem->g1Funcs[2] = g1_pu;      /* < q, \nabla\cdot v > */
471:   fem->g1Funcs[3] = NULL;
472:   fem->g2Funcs[0] = NULL;
473:   fem->g2Funcs[1] = g2_up;      /* < \nabla\cdot v, p > */
474:   fem->g2Funcs[2] = NULL;
475:   fem->g2Funcs[3] = NULL;
476:   fem->g3Funcs[0] = g3_uu;      /* < \nabla v, \nabla u + {\nabla u}^T > */
477:   fem->g3Funcs[1] = NULL;
478:   fem->g3Funcs[2] = NULL;
479:   fem->g3Funcs[3] = NULL;
480:   switch (user->dim) {
481:   case 2:
482:     user->exactFuncs[0] = quadratic_u_2d;
483:     user->exactFuncs[1] = linear_p_2d;
484:     user->initialGuess[0] = zero_2d;
485:     user->initialGuess[1] = zero_1d;
486:     break;
487:   case 3:
488:     user->exactFuncs[0] = quadratic_u_3d;
489:     user->exactFuncs[1] = linear_p_3d;
490:     user->initialGuess[0] = zero_3d;
491:     user->initialGuess[1] = zero_1d;
492:     break;
493:   default:
494:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
495:   }
496:   return(0);
497: }

501: PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, MatNullSpace *nullSpace)
502: {
503:   Vec            vec, localVec;

507:   DMGetGlobalVector(dm, &vec);
508:   DMGetLocalVector(dm, &localVec);
509:   VecSet(vec,  0.0);
510:   /* Put a constant in for all pressures
511:      Could change this to project the constant function onto the pressure space (when that is finished) */
512:   {
513:     PetscSection section;
514:     PetscInt     pStart, pEnd, p;
515:     PetscScalar  *a;

517:     DMGetDefaultSection(dm, &section);
518:     PetscSectionGetChart(section, &pStart, &pEnd);
519:     VecGetArray(localVec, &a);
520:     for (p = pStart; p < pEnd; ++p) {
521:       PetscInt fDim, off, d;

523:       PetscSectionGetFieldDof(section, p, 1, &fDim);
524:       PetscSectionGetFieldOffset(section, p, 1, &off);
525:       for (d = 0; d < fDim; ++d) a[off+d] = 1.0;
526:     }
527:     VecRestoreArray(localVec, &a);
528:   }
529:   DMLocalToGlobalBegin(dm, localVec, INSERT_VALUES, vec);
530:   DMLocalToGlobalEnd(dm, localVec, INSERT_VALUES, vec);
531:   DMRestoreLocalVector(dm, &localVec);
532:   VecNormalize(vec, NULL);
533:   if (user->debug) {
534:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
535:     VecView(vec, PETSC_VIEWER_STDOUT_WORLD);
536:   }
537:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
538:   DMRestoreGlobalVector(dm, &vec);
539:   /* New style for field null spaces */
540:   {
541:     PetscObject  pressure;
542:     MatNullSpace nullSpacePres;

544:     DMGetField(dm, 1, &pressure);
545:     MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);
546:     PetscObjectCompose(pressure, "nullspace", (PetscObject) nullSpacePres);
547:     MatNullSpaceDestroy(&nullSpacePres);
548:   }
549:   return(0);
550: }

554: /*
555:   FormJacobianAction - Form the global Jacobian action Y = JX from the global input X

557:   Input Parameters:
558: + mat - The Jacobian shell matrix
559: - X  - Global input vector

561:   Output Parameter:
562: . Y  - Local output vector

564:   Note:
565:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
566:   like a GPU, or vectorize on a multicore machine.

568: .seealso: FormJacobianActionLocal()
569: */
570: PetscErrorCode FormJacobianAction(Mat J, Vec X,  Vec Y)
571: {
572:   JacActionCtx   *ctx;
573:   DM             dm;
574:   Vec            localX, localY;
575:   PetscInt       N, n;

579: #if 0
580:   /* Needs petscimpl.h */
584: #endif
585:   MatShellGetContext(J, &ctx);
586:   dm   = ctx->dm;

588:   /* determine whether X = localX */
589:   DMGetLocalVector(dm, &localX);
590:   DMGetLocalVector(dm, &localY);
591:   VecGetSize(X, &N);
592:   VecGetSize(localX, &n);

594:   if (n != N) { /* X != localX */
595:     VecSet(localX, 0.0);
596:     DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
597:     DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
598:   } else {
599:     DMRestoreLocalVector(dm, &localX);
600:     localX = X;
601:   }
602:   DMPlexComputeJacobianActionFEM(dm, J, localX, localY, ctx->user);
603:   if (n != N) {
604:     DMRestoreLocalVector(dm, &localX);
605:   }
606:   VecSet(Y, 0.0);
607:   DMLocalToGlobalBegin(dm, localY, ADD_VALUES, Y);
608:   DMLocalToGlobalEnd(dm, localY, ADD_VALUES, Y);
609:   DMRestoreLocalVector(dm, &localY);
610:   if (0) {
611:     Vec       r;
612:     PetscReal norm;

614:     VecDuplicate(X, &r);
615:     MatMult(ctx->J, X, r);
616:     VecAXPY(r, -1.0, Y);
617:     VecNorm(r, NORM_2, &norm);
618:     if (norm > 1.0e-8) {
619:       PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Input:\n");
620:       VecView(X, PETSC_VIEWER_STDOUT_WORLD);
621:       PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Result:\n");
622:       VecView(Y, PETSC_VIEWER_STDOUT_WORLD);
623:       PetscPrintf(PETSC_COMM_WORLD, "Difference:\n");
624:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
625:       SETERRQ1(PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "The difference with assembled multiply is too large %g", norm);
626:     }
627:     VecDestroy(&r);
628:   }
629:   return(0);
630: }

634: int main(int argc, char **argv)
635: {
636:   SNES           snes;                 /* nonlinear solver */
637:   DM             dm;                   /* problem definition */
638:   Vec            u,r;                  /* solution, residual vectors */
639:   Mat            A,J;                  /* Jacobian matrix */
640:   MatNullSpace   nullSpace;            /* May be necessary for pressure */
641:   AppCtx         user;                 /* user-defined work context */
642:   JacActionCtx   userJ;                /* context for Jacobian MF action */
643:   PetscInt       its;                  /* iterations for convergence */
644:   PetscReal      error         = 0.0;  /* L_2 error in the solution */
645:   PetscInt       numComponents = 0, order, f;
646:   PetscSpace     P;

649:   PetscInitialize(&argc, &argv, NULL, help);
650:   ProcessOptions(PETSC_COMM_WORLD, &user);
651:   SNESCreate(PETSC_COMM_WORLD, &snes);
652:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
653:   SNESSetDM(snes, dm);

655:   user.fem.fe    = user.fe;
656:   user.fem.feBd  = NULL;
657:   user.fem.feAux = NULL;
658:   PetscFECreateDefault(dm, user.dim, user.dim, user.simplex, "vel_", -1, &user.fe[0]);
659:   PetscFEGetBasisSpace(user.fe[0], &P);
660:   PetscSpaceGetOrder(P, &order);
661:   PetscFECreateDefault(dm, user.dim, 1, user.simplex, "pres_", order, &user.fe[1]);
662:   //SetupElement(dm, &user);
663:   for (f = 0; f < NUM_FIELDS; ++f) {
664:     PetscInt numComp;
665:     PetscFEGetNumComponents(user.fe[f], &numComp);
666:     numComponents += numComp;
667:   }
668:   PetscMalloc(NUM_FIELDS * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
669:   SetupProblem(dm, &user);
670:   SetupSection(dm, &user);
671:   DMPlexCreateClosureIndex(dm, NULL);

673:   DMCreateGlobalVector(dm, &u);
674:   VecDuplicate(u, &r);

676:   DMSetMatType(dm,MATAIJ);
677:   DMCreateMatrix(dm, &J);
678:   if (user.jacobianMF) {
679:     PetscInt M, m, N, n;

681:     MatGetSize(J, &M, &N);
682:     MatGetLocalSize(J, &m, &n);
683:     MatCreate(PETSC_COMM_WORLD, &A);
684:     MatSetSizes(A, m, n, M, N);
685:     MatSetType(A, MATSHELL);
686:     MatSetUp(A);
687:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);

689:     userJ.dm   = dm;
690:     userJ.J    = J;
691:     userJ.user = &user;

693:     DMCreateLocalVector(dm, &userJ.u);
694:     DMPlexProjectFunctionLocal(dm, user.fe, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
695:     MatShellSetContext(A, &userJ);
696:   } else {
697:     A = J;
698:   }
699:   CreatePressureNullSpace(dm, &user, &nullSpace);
700:   MatSetNullSpace(J, nullSpace);
701:   if (A != J) {
702:     MatSetNullSpace(A, nullSpace);
703:   }

705:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);
706:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*))DMPlexComputeJacobianFEM,&user);
707:   SNESSetJacobian(snes, A, J, NULL, NULL);

709:   SNESSetFromOptions(snes);

711:   DMPlexProjectFunction(dm, user.fe, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
712:   if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}
713:   if (user.runType == RUN_FULL) {
714:     DMPlexProjectFunction(dm, user.fe, user.initialGuess, NULL, INSERT_VALUES, u);
715:     if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}
716:     if (user.debug) {
717:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
718:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
719:     }
720:     SNESSolve(snes, NULL, u);
721:     SNESGetIterationNumber(snes, &its);
722:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
723:     DMPlexComputeL2Diff(dm, user.fe, user.exactFuncs, NULL, u, &error);
724:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g\n", error);
725:     if (user.showSolution) {
726:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
727:       VecChop(u, 3.0e-9);
728:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
729:     }
730:   } else {
731:     PetscReal res = 0.0;

733:     /* Check discretization error */
734:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
735:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
736:     DMPlexComputeL2Diff(dm, user.fe, user.exactFuncs, NULL, u, &error);
737:     if (error >= 1.0e-11) {
738:       PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
739:     } else {
740:       PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n", error);
741:     }
742:     /* Check residual */
743:     SNESComputeFunction(snes, u, r);
744:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
745:     VecChop(r, 1.0e-10);
746:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
747:     VecNorm(r, NORM_2, &res);
748:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
749:     /* Check Jacobian */
750:     {
751:       Vec          b;
752:       PetscBool    isNull;

754:       SNESComputeJacobian(snes, u, A, A);
755:       MatNullSpaceTest(nullSpace, J, &isNull);
756:       if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
757:       VecDuplicate(u, &b);
758:       VecSet(r, 0.0);
759:       SNESComputeFunction(snes, r, b);
760:       MatMult(A, u, r);
761:       VecAXPY(r, 1.0, b);
762:       VecDestroy(&b);
763:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
764:       VecChop(r, 1.0e-10);
765:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
766:       VecNorm(r, NORM_2, &res);
767:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
768:     }
769:   }
770:   VecViewFromOptions(u, NULL, "-sol_vec_view");

772:   PetscFree(user.exactFuncs);
773:   DestroyElement(&user);
774:   MatNullSpaceDestroy(&nullSpace);
775:   if (user.jacobianMF) {
776:     VecDestroy(&userJ.u);
777:   }
778:   if (A != J) {
779:     MatDestroy(&A);
780:   }
781:   MatDestroy(&J);
782:   VecDestroy(&u);
783:   VecDestroy(&r);
784:   SNESDestroy(&snes);
785:   DMDestroy(&dm);
786:   PetscFinalize();
787:   return 0;
788: }