Actual source code: ex62.c

petsc-3.4.4 2014-03-13
  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 a Python script to generate a tabulation of the finite element basis
 18: functions at quadrature points, which we put in a C header file. The generic
 19: command would be:

 21:     bin/pythonscripts/PetscGenerateFEMQuadrature.py dim order dim 1 laplacian dim order 1 1 gradient src/snes/examples/tutorials/ex62.h

 23: We can currently generate an arbitrary order Lagrange element. The underlying
 24: FIAT code is capable of handling more exotic elements, but these have not been
 25: tested with this code.

 27: Field Data:

 29:   Sieve data is organized by point, and the closure operation just stacks up the
 30: data from each sieve point in the closure. Thus, for a P_2-P_1 Stokes element, we
 31: have

 33:   cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
 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} p_{v_0} u_{v_1} v_{v_1} p_{v_1} u_{v_2} v_{v_2} p_{v_2}]

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

 40:   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}]

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

 45: Next Steps:

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

 51: For tensor product meshes, see SNES ex67, ex72
 52: */

 54: #include <petscdmplex.h>
 55: #include <petscsnes.h>

 57: /*------------------------------------------------------------------------------
 58:   This code can be generated using 'bin/pythonscripts/PetscGenerateFEMQuadrature.py dim order dim 1 laplacian dim order 1 1 gradient src/snes/examples/tutorials/ex62.h'
 59:  -----------------------------------------------------------------------------*/
 60: #include "ex62.h"

 62: typedef enum {NEUMANN, DIRICHLET} BCType;
 63: typedef enum {RUN_FULL, RUN_TEST} RunType;

 65: typedef struct {
 66:   DM            dm;                /* REQUIRED in order to use SNES evaluation functions */
 67:   PetscFEM      fem;               /* REQUIRED to use DMPlexComputeResidualFEM() */
 68:   PetscInt      debug;             /* The debugging level */
 69:   PetscMPIInt   rank;              /* The process rank */
 70:   PetscMPIInt   numProcs;          /* The number of processes */
 71:   RunType       runType;           /* Whether to run tests, or solve the full problem */
 72:   PetscBool     jacobianMF;        /* Whether to calculate the Jacobian action on the fly */
 73:   PetscLogEvent createMeshEvent;
 74:   PetscBool     showInitial, showSolution;
 75:   /* Domain and mesh definition */
 76:   PetscInt      dim;               /* The topological mesh dimension */
 77:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 78:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 79:   char          partitioner[2048]; /* The graph partitioner */
 80:   /* GPU partitioning */
 81:   PetscInt      numBatches;        /* The number of cell batches per kernel */
 82:   PetscInt      numBlocks;         /* The number of concurrent blocks per kernel */
 83:   /* Element quadrature */
 84:   PetscInt        order[NUM_FIELDS];
 85:   PetscQuadrature q[NUM_FIELDS];
 86:   /* Problem definition */
 87:   void (*f0Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]); /* f0_u(x,y,z), and f0_p(x,y,z) */
 88:   void (*f1Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]); /* f1_u(x,y,z), and f1_p(x,y,z) */
 89:   void (*g0Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], 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) */
 90:   void (*g1Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], 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) */
 91:   void (*g2Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], 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) */
 92:   void (*g3Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], 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) */
 93:   void (*exactFuncs[NUM_BASIS_COMPONENTS_TOTAL])(const PetscReal x[], PetscScalar *u); /* The exact solution function u(x,y,z), v(x,y,z), and p(x,y,z) */
 94:   BCType bcType;
 95: } AppCtx;

 97: void zero(const PetscReal coords[], PetscScalar *u)
 98: {
 99:   *u = 0.0;
100: }

102: /*
103:   In 2D we use exact solution:

105:     u = x^2 + y^2
106:     v = 2 x^2 - 2xy
107:     p = x + y - 1
108:     f_x = f_y = 3

110:   so that

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

120: void quadratic_v_2d(const PetscReal x[], PetscScalar *v)
121: {
122:   *v = 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
123: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

217:   so that

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

227: void quadratic_v_3d(const PetscReal x[], PetscScalar *v)
228: {
229:   *v = x[1]*x[1] + x[2]*x[2];
230: }

232: void quadratic_w_3d(const PetscReal x[], PetscScalar *w)
233: {
234:   *w = x[0]*x[0] + x[1]*x[1] - 2.0*(x[0] + x[1])*x[2];
235: }

237: void linear_p_3d(const PetscReal x[], PetscScalar *p)
238: {
239:   *p = x[0] + x[1] + x[2] - 1.5;
240: }

244: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
245: {
246:   const char    *bcTypes[2]  = {"neumann", "dirichlet"};
247:   const char    *runTypes[2] = {"full", "test"};
248:   PetscInt       bc, run, n;

252:   options->debug           = 0;
253:   options->runType         = RUN_FULL;
254:   options->dim             = 2;
255:   options->interpolate     = PETSC_FALSE;
256:   options->refinementLimit = 0.0;
257:   options->bcType          = DIRICHLET;
258:   options->numBatches      = 1;
259:   options->numBlocks       = 1;
260:   options->jacobianMF      = PETSC_FALSE;
261:   options->showInitial     = PETSC_FALSE;
262:   options->showSolution    = PETSC_TRUE;
263:   options->order[0]        = 1;
264:   options->order[1]        = 1;

266:   options->fem.quad    = (PetscQuadrature*) &options->q;
267:   options->fem.f0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f0Funcs;
268:   options->fem.f1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f1Funcs;
269:   options->fem.g0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g0Funcs;
270:   options->fem.g1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g1Funcs;
271:   options->fem.g2Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g2Funcs;
272:   options->fem.g3Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g3Funcs;
273:   options->fem.bcFuncs = (void (**)(const PetscReal[], PetscScalar *)) &options->exactFuncs;

275:   MPI_Comm_size(comm, &options->numProcs);
276:   MPI_Comm_rank(comm, &options->rank);
277:   PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
278:   PetscOptionsInt("-debug", "The debugging level", "ex62.c", options->debug, &options->debug, NULL);
279:   run  = options->runType;
280:   PetscOptionsEList("-run_type", "The run type", "ex62.c", runTypes, 2, runTypes[options->runType], &run, NULL);

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

284:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex62.c", options->dim, &options->dim, NULL);
285:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex62.c", options->interpolate, &options->interpolate, NULL);
286:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex62.c", options->refinementLimit, &options->refinementLimit, NULL);
287:   PetscStrcpy(options->partitioner, "chaco");
288:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
289:   bc   = options->bcType;
290:   PetscOptionsEList("-bc_type","Type of boundary condition","ex62.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);

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

294:   PetscOptionsInt("-gpu_batches", "The number of cell batches per kernel", "ex62.c", options->numBatches, &options->numBatches, NULL);
295:   PetscOptionsInt("-gpu_blocks", "The number of concurrent blocks per kernel", "ex62.c", options->numBlocks, &options->numBlocks, NULL);
296:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex62.c", options->jacobianMF, &options->jacobianMF, NULL);
297:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex62.c", options->showInitial, &options->showInitial, NULL);
298:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex62.c", options->showSolution, &options->showSolution, NULL);
299:   n    = NUM_FIELDS;
300:   PetscOptionsIntArray("-order", "The FEM order", "ex62.c", options->order, &n, NULL);
301:   PetscOptionsEnd();

303:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
304:   return(0);
305: }

309: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
310: {
311:   Vec            lv;
312:   PetscInt       p;
313:   PetscMPIInt    rank, numProcs;

317:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
318:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
319:   DMGetLocalVector(dm, &lv);
320:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
321:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
322:   PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
323:   for (p = 0; p < numProcs; ++p) {
324:     if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
325:     PetscBarrier((PetscObject) dm);
326:   }
327:   DMRestoreLocalVector(dm, &lv);
328:   return(0);
329: }

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

342:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
343:   DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
344:   {
345:     DM refinedMesh     = NULL;
346:     DM distributedMesh = NULL;

348:     /* Refine mesh using a volume constraint */
349:     DMPlexSetRefinementLimit(*dm, refinementLimit);
350:     DMRefine(*dm, comm, &refinedMesh);
351:     if (refinedMesh) {
352:       DMDestroy(dm);
353:       *dm  = refinedMesh;
354:     }
355:     /* Distribute mesh over processes */
356:     DMPlexDistribute(*dm, partitioner, 0, &distributedMesh);
357:     if (distributedMesh) {
358:       DMDestroy(dm);
359:       *dm  = distributedMesh;
360:     }
361:   }
362:   DMSetFromOptions(*dm);
363:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
364:   user->dm = *dm;
365:   return(0);
366: }

370: PetscErrorCode SetupQuadrature(AppCtx *user)
371: {
372:   PetscReal     *x, *w;
373:   const PetscInt dim = user->dim;
374:   PetscInt       order, numPoints, p, d;

378:   /* Velocity discretization */
379:   order     = PetscMax(user->order[0], user->order[1]);
380:   numPoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
381:   if (numPoints != NUM_QUADRATURE_POINTS_0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of quadrature points: %d != %d", numPoints, NUM_QUADRATURE_POINTS_0);
382:   PetscDTGaussJacobiQuadrature(dim, order, -1.0, 1.0, &x, &w);
383:   for (p = 0; p < numPoints; ++p) {
384:     for (d = 0; d < dim; ++d) {
385:       if (fabs(x[p*dim+d] - points_0[p*dim+d]) > 1.0e-10) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid point %d, component %d: %g != %g", p, d, x[p*dim+d], points_0[p*dim+d]);
386:     }
387:     if (fabs(w[p] - weights_0[p]) > 1.0e-10) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid weight %d: %g != %g", p, w[p], weights_0[p]);
388:   }
389:   user->fem.quad[0].numQuadPoints = numPoints;
390:   user->fem.quad[0].quadPoints    = x;
391:   user->fem.quad[0].quadWeights   = w;
392:   user->fem.quad[0].numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
393:   user->fem.quad[0].numComponents = NUM_BASIS_COMPONENTS_0;
394:   user->fem.quad[0].basis         = Basis_0;
395:   user->fem.quad[0].basisDer      = BasisDerivatives_0;
396:   /* Pressure discretization */
397:   order     = PetscMax(user->order[0], user->order[1]);
398:   numPoints = dim > 1 ? dim > 2 ? order*PetscSqr(order) : PetscSqr(order) : order;
399:   if (numPoints != NUM_QUADRATURE_POINTS_1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of quadrature points: %d != %d", numPoints, NUM_QUADRATURE_POINTS_1);
400:   PetscDTGaussJacobiQuadrature(dim, order, -1.0, 1.0, &x, &w);
401:   for (p = 0; p < numPoints; ++p) {
402:     for (d = 0; d < dim; ++d) {
403:       if (fabs(x[p*dim+d] - points_1[p*dim+d]) > 1.0e-10) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid point %d, component %d: %g != %g", p, d, x[p*dim+d], points_1[p*dim+d]);
404:     }
405:     if (fabs(w[p] - weights_1[p]) > 1.0e-10) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid weight %d: %g != %g", p, w[p], weights_1[p]);
406:   }
407:   user->fem.quad[1].numQuadPoints = numPoints;
408:   user->fem.quad[1].quadPoints    = x;
409:   user->fem.quad[1].quadWeights   = w;
410:   user->fem.quad[1].numBasisFuncs = NUM_BASIS_FUNCTIONS_1;
411:   user->fem.quad[1].numComponents = NUM_BASIS_COMPONENTS_1;
412:   user->fem.quad[1].basis         = Basis_1;
413:   user->fem.quad[1].basisDer      = BasisDerivatives_1;
414:   return(0);
415: }

419: PetscErrorCode DestroyQuadrature(AppCtx *user)
420: {

424:   PetscFree(user->fem.quad[0].quadPoints);
425:   PetscFree(user->fem.quad[0].quadWeights);
426:   PetscFree(user->fem.quad[1].quadPoints);
427:   PetscFree(user->fem.quad[1].quadWeights);
428:   return(0);
429: }

433: /*
434:   There is a problem here with uninterpolated meshes. The index in numDof[] is not dimension in this case,
435:   but sieve depth.
436: */
437: PetscErrorCode SetupSection(DM dm, AppCtx *user)
438: {
439:   PetscSection   section;
440:   const PetscInt numFields           = NUM_FIELDS;
441:   PetscInt       dim                 = user->dim;
442:   PetscInt       numBC               = 0;
443:   PetscInt       numComp[NUM_FIELDS] = {NUM_BASIS_COMPONENTS_0, NUM_BASIS_COMPONENTS_1};
444:   PetscInt       bcFields[1]         = {0};
445:   IS             bcPoints[1]         = {NULL};
446:   PetscInt       numDof[NUM_FIELDS*(SPATIAL_DIM_0+1)];
447:   PetscInt       f, d;

451:   if (dim != SPATIAL_DIM_0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_0);
452:   if (dim != SPATIAL_DIM_1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_1);
453:   for (d = 0; d <= dim; ++d) {
454:     numDof[0*(dim+1)+d] = numDof_0[d];
455:     numDof[1*(dim+1)+d] = numDof_1[d];
456:   }
457:   for (f = 0; f < numFields; ++f) {
458:     for (d = 1; d < dim; ++d) {
459:       if ((numDof[f*(dim+1)+d] > 0) && !user->interpolate) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
460:     }
461:   }
462:   if (user->bcType == DIRICHLET) {
463:     numBC = 1;
464:     DMPlexGetStratumIS(dm, "marker", 1, &bcPoints[0]);
465:   }
466:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, &section);
467:   PetscSectionSetFieldName(section, 0, "velocity");
468:   PetscSectionSetFieldName(section, 1, "pressure");
469:   DMSetDefaultSection(dm, section);
470:   if (user->bcType == DIRICHLET) {
471:     ISDestroy(&bcPoints[0]);
472:   }
473:   return(0);
474: }

478: PetscErrorCode SetupExactSolution(DM dm, AppCtx *user)
479: {
480:   PetscFEM       *fem = &user->fem;

484:   fem->f0Funcs[0] = f0_u;
485:   fem->f0Funcs[1] = f0_p;
486:   fem->f1Funcs[0] = f1_u;
487:   fem->f1Funcs[1] = f1_p;
488:   fem->g0Funcs[0] = NULL;
489:   fem->g0Funcs[1] = NULL;
490:   fem->g0Funcs[2] = NULL;
491:   fem->g0Funcs[3] = NULL;
492:   fem->g1Funcs[0] = NULL;
493:   fem->g1Funcs[1] = NULL;
494:   fem->g1Funcs[2] = g1_pu;      /* < q, \nabla\cdot v > */
495:   fem->g1Funcs[3] = NULL;
496:   fem->g2Funcs[0] = NULL;
497:   fem->g2Funcs[1] = g2_up;      /* < \nabla\cdot v, p > */
498:   fem->g2Funcs[2] = NULL;
499:   fem->g2Funcs[3] = NULL;
500:   fem->g3Funcs[0] = g3_uu;      /* < \nabla v, \nabla u + {\nabla u}^T > */
501:   fem->g3Funcs[1] = NULL;
502:   fem->g3Funcs[2] = NULL;
503:   fem->g3Funcs[3] = NULL;
504:   switch (user->dim) {
505:   case 2:
506:     user->exactFuncs[0] = quadratic_u_2d;
507:     user->exactFuncs[1] = quadratic_v_2d;
508:     user->exactFuncs[2] = linear_p_2d;
509:     break;
510:   case 3:
511:     user->exactFuncs[0] = quadratic_u_3d;
512:     user->exactFuncs[1] = quadratic_v_3d;
513:     user->exactFuncs[2] = quadratic_w_3d;
514:     user->exactFuncs[3] = linear_p_3d;
515:     break;
516:   default:
517:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
518:   }
519:   DMPlexSetFEMIntegration(dm, FEMIntegrateResidualBatch, NULL, FEMIntegrateJacobianActionBatch, FEMIntegrateJacobianBatch);
520:   return(0);
521: }

525: PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, MatNullSpace *nullSpace)
526: {
527:   Vec            vec, localVec;

531:   DMGetGlobalVector(dm, &vec);
532:   DMGetLocalVector(dm, &localVec);
533:   VecSet(vec,  0.0);
534:   /* Put a constant in for all pressures
535:      Could change this to project the constant function onto the pressure space (when that is finished) */
536:   {
537:     PetscSection section;
538:     PetscInt     pStart, pEnd, p;
539:     PetscScalar  *a;

541:     DMGetDefaultSection(dm, &section);
542:     PetscSectionGetChart(section, &pStart, &pEnd);
543:     VecGetArray(localVec, &a);
544:     for (p = pStart; p < pEnd; ++p) {
545:       PetscInt fDim, off, d;

547:       PetscSectionGetFieldDof(section, p, 1, &fDim);
548:       PetscSectionGetFieldOffset(section, p, 1, &off);
549:       for (d = 0; d < fDim; ++d) a[off+d] = 1.0;
550:     }
551:     VecRestoreArray(localVec, &a);
552:   }
553:   DMLocalToGlobalBegin(dm, localVec, INSERT_VALUES, vec);
554:   DMLocalToGlobalEnd(dm, localVec, INSERT_VALUES, vec);
555:   DMRestoreLocalVector(dm, &localVec);
556:   VecNormalize(vec, NULL);
557:   if (user->debug) {
558:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
559:     VecView(vec, PETSC_VIEWER_STDOUT_WORLD);
560:   }
561:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
562:   DMRestoreGlobalVector(dm, &vec);
563:   /* New style for field null spaces */
564:   {
565:     PetscObject  pressure;
566:     MatNullSpace nullSpacePres;

568:     DMGetField(dm, 1, &pressure);
569:     MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);
570:     PetscObjectCompose(pressure, "nullspace", (PetscObject) nullSpacePres);
571:     MatNullSpaceDestroy(&nullSpacePres);
572:   }
573:   return(0);
574: }

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

581:   Input Parameters:
582: + mat - The Jacobian shell matrix
583: - X  - Global input vector

585:   Output Parameter:
586: . Y  - Local output vector

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

592: .seealso: FormJacobianActionLocal()
593: */
594: PetscErrorCode FormJacobianAction(Mat J, Vec X,  Vec Y)
595: {
596:   JacActionCtx   *ctx;
597:   DM             dm;
598:   Vec            localX, localY;
599:   PetscInt       N, n;

603: #if 0
604:   /* Needs petscimpl.h */
608: #endif
609:   MatShellGetContext(J, &ctx);
610:   dm   = ctx->dm;

612:   /* determine whether X = localX */
613:   DMGetLocalVector(dm, &localX);
614:   DMGetLocalVector(dm, &localY);
615:   VecGetSize(X, &N);
616:   VecGetSize(localX, &n);

618:   if (n != N) { /* X != localX */
619:     VecSet(localX, 0.0);
620:     DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
621:     DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
622:   } else {
623:     DMRestoreLocalVector(dm, &localX);
624:     localX = X;
625:   }
626:   DMPlexComputeJacobianActionFEM(dm, J, localX, localY, ctx->user);
627:   if (n != N) {
628:     DMRestoreLocalVector(dm, &localX);
629:   }
630:   VecSet(Y, 0.0);
631:   DMLocalToGlobalBegin(dm, localY, ADD_VALUES, Y);
632:   DMLocalToGlobalEnd(dm, localY, ADD_VALUES, Y);
633:   DMRestoreLocalVector(dm, &localY);
634:   if (0) {
635:     Vec       r;
636:     PetscReal norm;

638:     VecDuplicate(X, &r);
639:     MatMult(ctx->J, X, r);
640:     VecAXPY(r, -1.0, Y);
641:     VecNorm(r, NORM_2, &norm);
642:     if (norm > 1.0e-8) {
643:       PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Input:\n");
644:       VecView(X, PETSC_VIEWER_STDOUT_WORLD);
645:       PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Result:\n");
646:       VecView(Y, PETSC_VIEWER_STDOUT_WORLD);
647:       PetscPrintf(PETSC_COMM_WORLD, "Difference:\n");
648:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
649:       SETERRQ1(PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "The difference with assembled multiply is too large %g", norm);
650:     }
651:     VecDestroy(&r);
652:   }
653:   return(0);
654: }

658: int main(int argc, char **argv)
659: {
660:   SNES           snes;                 /* nonlinear solver */
661:   Vec            u,r;                  /* solution, residual vectors */
662:   Mat            A,J;                  /* Jacobian matrix */
663:   MatNullSpace   nullSpace;            /* May be necessary for pressure */
664:   AppCtx         user;                 /* user-defined work context */
665:   JacActionCtx   userJ;                /* context for Jacobian MF action */
666:   PetscInt       its;                  /* iterations for convergence */
667:   PetscReal      error         = 0.0;  /* L_2 error in the solution */
668:   const PetscInt numComponents = NUM_BASIS_COMPONENTS_TOTAL;

671:   PetscInitialize(&argc, &argv, NULL, help);
672:   ProcessOptions(PETSC_COMM_WORLD, &user);
673:   SNESCreate(PETSC_COMM_WORLD, &snes);
674:   CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);
675:   SNESSetDM(snes, user.dm);

677:   SetupExactSolution(user.dm, &user);
678:   SetupQuadrature(&user);
679:   SetupSection(user.dm, &user);

681:   DMCreateGlobalVector(user.dm, &u);
682:   VecDuplicate(u, &r);

684:   DMCreateMatrix(user.dm, MATAIJ, &J);
685:   if (user.jacobianMF) {
686:     PetscInt M, m, N, n;

688:     MatGetSize(J, &M, &N);
689:     MatGetLocalSize(J, &m, &n);
690:     MatCreate(PETSC_COMM_WORLD, &A);
691:     MatSetSizes(A, m, n, M, N);
692:     MatSetType(A, MATSHELL);
693:     MatSetUp(A);
694:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);

696:     userJ.dm   = user.dm;
697:     userJ.J    = J;
698:     userJ.user = &user;

700:     DMCreateLocalVector(user.dm, &userJ.u);
701:     DMPlexProjectFunctionLocal(user.dm, numComponents, user.exactFuncs, INSERT_BC_VALUES, userJ.u);
702:     MatShellSetContext(A, &userJ);
703:   } else {
704:     A = J;
705:   }
706:   CreatePressureNullSpace(user.dm, &user, &nullSpace);
707:   MatSetNullSpace(J, nullSpace);
708:   if (A != J) {
709:     MatSetNullSpace(A, nullSpace);
710:   }

712:   DMSNESSetFunctionLocal(user.dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);
713:   DMSNESSetJacobianLocal(user.dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*))DMPlexComputeJacobianFEM,&user);
714:   SNESSetJacobian(snes, A, J, NULL, NULL);

716:   SNESSetFromOptions(snes);

718:   DMPlexProjectFunction(user.dm, numComponents, user.exactFuncs, INSERT_ALL_VALUES, u);
719:   if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
720:   if (user.runType == RUN_FULL) {
721:     PetscScalar (*initialGuess[numComponents])(const PetscReal x[]);
722:     PetscInt c;

724:     for (c = 0; c < numComponents; ++c) initialGuess[c] = zero;
725:     DMPlexProjectFunction(user.dm, numComponents, initialGuess, INSERT_VALUES, u);
726:     if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
727:     if (user.debug) {
728:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
729:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
730:     }
731:     SNESSolve(snes, NULL, u);
732:     SNESGetIterationNumber(snes, &its);
733:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
734:     DMPlexComputeL2Diff(user.dm, user.fem.quad, user.exactFuncs, u, &error);
735:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g\n", error);
736:     if (user.showSolution) {
737:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
738:       VecChop(u, 3.0e-9);
739:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
740:     }
741:   } else {
742:     PetscReal res = 0.0;

744:     /* Check discretization error */
745:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
746:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
747:     DMPlexComputeL2Diff(user.dm, user.fem.quad, user.exactFuncs, u, &error);
748:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
749:     /* Check residual */
750:     SNESComputeFunction(snes, u, r);
751:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
752:     VecChop(r, 1.0e-10);
753:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
754:     VecNorm(r, NORM_2, &res);
755:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
756:     /* Check Jacobian */
757:     {
758:       Vec          b;
759:       MatStructure flag;
760:       PetscBool    isNull;

762:       SNESComputeJacobian(snes, u, &A, &A, &flag);
763:       MatNullSpaceTest(nullSpace, J, &isNull);
764:       if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
765:       VecDuplicate(u, &b);
766:       VecSet(r, 0.0);
767:       SNESComputeFunction(snes, r, b);
768:       MatMult(A, u, r);
769:       VecAXPY(r, 1.0, b);
770:       VecDestroy(&b);
771:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
772:       VecChop(r, 1.0e-10);
773:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
774:       VecNorm(r, NORM_2, &res);
775:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
776:     }
777:   }

779:   if (user.runType == RUN_FULL) {
780:     PetscViewer viewer;
781:     Vec         uLocal;

783:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
784:     PetscViewerSetType(viewer, PETSCVIEWERVTK);
785:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
786:     PetscViewerFileSetName(viewer, "ex62_sol.vtk");

788:     DMGetLocalVector(user.dm, &uLocal);
789:     DMGlobalToLocalBegin(user.dm, u, INSERT_VALUES, uLocal);
790:     DMGlobalToLocalEnd(user.dm, u, INSERT_VALUES, uLocal);

792:     PetscObjectReference((PetscObject) user.dm); /* Needed because viewer destroys the DM */
793:     PetscObjectReference((PetscObject) uLocal); /* Needed because viewer destroys the Vec */
794:     PetscViewerVTKAddField(viewer, (PetscObject) user.dm, DMPlexVTKWriteAll, PETSC_VTK_POINT_FIELD, (PetscObject) uLocal);
795:     DMRestoreLocalVector(user.dm, &uLocal);
796:     PetscViewerDestroy(&viewer);
797:   }

799:   DestroyQuadrature(&user);
800:   MatNullSpaceDestroy(&nullSpace);
801:   if (user.jacobianMF) {
802:     VecDestroy(&userJ.u);
803:   }
804:   if (A != J) {
805:     MatDestroy(&A);
806:   }
807:   MatDestroy(&J);
808:   VecDestroy(&u);
809:   VecDestroy(&r);
810:   SNESDestroy(&snes);
811:   DMDestroy(&user.dm);
812:   PetscFinalize();
813:   return 0;
814: }