Actual source code: ex31.c

petsc-3.4.4 2014-03-13
  1: static char help[] = "Stokes Problem with Temperature 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: TODO for Mantle Convection:
  7:   - Variable viscosity
  8:   - Free-slip boundary condition on upper surface
  9:   - Parse Citcom input

 11: The isoviscous Stokes problem, which we discretize using the finite
 12: element method on an unstructured mesh. The weak form equations are

 14:   < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > + < v, f > = 0
 15:   < q, \nabla\cdot v >                                                    = 0
 16:   < \nabla t, \nabla T>                                                   = q_T

 18: Boundary Conditions:

 20:   -bc_type dirichlet   Dirichlet conditions on the entire boundary, coming from the exact solution functions

 22:   -dm_plex_separate_marker
 23:   -bc_type freeslip    Dirichlet conditions on the normal velocity at each boundary

 25: Discretization:

 27: We use a Python script to generate a tabulation of the finite element basis
 28: functions at quadrature points, which we put in a C header file. The generic
 29: command would be:

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

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

 37: Field Data:

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

 43:   cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
 44:   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}]

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

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

 52: Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
 53: puts it into the Sieve ordering.
 54: */

 56: #include <petscdmplex.h>
 57: #include <petscsnes.h>

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

 64: typedef enum {DIRICHLET, FREE_SLIP} BCType;
 65: typedef enum {RUN_FULL, RUN_TEST} RunType;
 66: typedef enum {FORCING_CONSTANT, FORCING_LINEAR, FORCING_CUBIC} ForcingType;

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

100: void zero(const PetscReal coords[], PetscScalar *u)
101: {
102:   *u = 0.0;
103: }

105: /*
106:   In 2D, for constant forcing,

108:     f_x = f_y = 3

110:   we use the exact solution,

112:     u = x^2 + y^2
113:     v = 2 x^2 - 2xy
114:     p = x + y - 1
115:     T = x + y

117:   so that

119:     -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
120:     \nabla \cdot u           = 2x - 2x                    = 0
121:     -\Delta T + q_T          = 0
122: */
123: void quadratic_u_2d(const PetscReal x[], PetscScalar *u)
124: {
125:   *u = x[0]*x[0] + x[1]*x[1];
126: };

128: void quadratic_v_2d(const PetscReal x[], PetscScalar *v)
129: {
130:   *v = 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
131: };

133: void linear_p_2d(const PetscReal x[], PetscScalar *p)
134: {
135:   *p = x[0] + x[1] - 1.0;
136: };

138: void linear_T_2d(const PetscReal x[], PetscScalar *T)
139: {
140:   *T = x[0] + x[1];
141: };

143: /*
144:   In 2D, for linear forcing,

146:     f_x =  3 - 8y
147:     f_y = -5 + 8x

149:   we use the exact solution,

151:     u =  2 x (x-1) (1 - 2 y)
152:     v = -2 y (y-1) (1 - 2 x)
153:     p = x + y - 1
154:     T = x + y

156:   so that

158:     -\Delta u + \nabla p + f = <-4+8y, 4-8x> + <1, 1> + <3-8y, 8x-5> = 0
159:     \nabla \cdot u           = (4x-2) (1-2y) - (4y-2) (1-2x)         = 0
160:     -\Delta T + q_T          = 0
161: */
162: void cubic_u_2d(const PetscReal x[], PetscScalar *u)
163: {
164:   *u = 2.0*x[0]*(x[0]-1.0)*(1.0 - 2.0*x[1]);
165: };

167: void cubic_v_2d(const PetscReal x[], PetscScalar *v)
168: {
169:   *v = -2.0*x[1]*(x[1]-1.0)*(1.0 - 2.0*x[0]);
170: };

172: /*
173:   Let \sigma = (\nabla u + \nabla u^T) = < \sigma_{ij} >, where \sigma_{ij} = \sigma_{ji}
174:   Then at the top and bottom (t = <1,0>),
175:     <\sigma_{00}, \sigma_{01}> = 0 so \sigma_{00} = A(x,y) y(1-y) \sigma_{01} = B(x,y) y(1-y)
176:   Using the left and right (t = <0,1>),
177:     <\sigma_{10}, \sigma_{11}> = 0 so \sigma_{10} = C(x,y) x(1-x) \sigma_{11} = D(x,y) x(1-x)
178:   Which means
179:     \sigma_{00} = A(x,y) y(1-y)        = 2 u_x
180:     \sigma_{01} = E(x,y) x(1-x) y(1-y) = u_y + v_x
181:     \sigma_{11} = D(x,y) x(1-x)        = 2 v_y
182:   Also we have
183:     u(x=0,1) = 0 ==> u = A'(x,y) x(1-x)
184:     v(y=0,1) = 0 ==> v = D'(x,y) y(1-y)
185:   Thus we need
186:     \int x - x^2 = x^2/2 - x^3/3 + C ==> 3 x^2 - 2 x^3 + 1 = 0 at x=0,1
187:   so that
188:     u =  (3 x^2 - 2 x^3 + 1) y(1-y)
189:     v = -(3 y^2 - 2 y^3 + 1) x(1-x)
190:     u_x =  6 x(1-x) y(1-y)
191:     v_y = -6 x(1-x) y(1-y)
192:     u_xx =  6 (1-2x) y(1-y)
193:     v_yy = -6 (1-2y) x(1-x)

195:   In 2D, for cubic forcing,

197:     f_x = -1 + 6 (1-2x) y(1-y)
198:     f_y = -1 - 6 (1-2y) x(1-x)

200:   we use the exact solution,

202:     u =  (3 x^2 - 2 x^3 + 1) y(1-y)
203:     v = -(3 y^2 - 2 y^3 + 1) x(1-x)
204:     p = x + y - 1
205:     T = x + y

207:   so that

209:     -\Delta u + \nabla p + f = <-6 (1-2x) y(1-y), 6 (1-2y) x(1-x)> + <1, 1> + <-1 + 6 (1-2x) y(1-y), -1 - 6 (1-2y) x(1-x)> = 0
210:     \nabla \cdot u           = 6 x(1-x) y(1-y) -6 (1-2y) x(1-x) = 0
211:     -\Delta T + q_T          = 0
212: */
213: void quintic_u_2d(const PetscReal x[], PetscScalar *u)
214: {
215:   *u = (3.0*x[0]*x[0] - 2.0*x[0]*x[0]*x[0] + 1.0)*x[1]*(1.0-x[1]);
216: };

218: void quintic_v_2d(const PetscReal x[], PetscScalar *v)
219: {
220:   *v = -(3.0*x[1]*x[1] - 2.0*x[1]*x[1]*x[1] + 1.0)*x[0]*(1.0-x[0]);
221: };

223: void f0_u_constant(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
224: {
225:   const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
226:   PetscInt       comp;

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

231: void f0_u_linear_2d(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
232: {
233:   f0[0] =  3.0 - 8.0*x[1];
234:   f0[1] = -5.0 + 8.0*x[0];
235: }

237: void f0_u_cubic_2d(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
238: {
239:   f0[0] = -1.0 + 6.0*(1.0 - 2.0*x[0])*x[1]*(1.0 - x[1]);
240:   f0[1] = -1.0 - 6.0*(1.0 - 2.0*x[1])*x[0]*(1.0 - x[0]);
241: }

243: /* 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}
244:    u[Ncomp]          = {p} */
245: void f1_u(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
246: {
247:   const PetscInt dim   = SPATIAL_DIM_0;
248:   const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
249:   PetscInt       comp, d;

251:   for (comp = 0; comp < Ncomp; ++comp) {
252:     for (d = 0; d < dim; ++d) {
253:       /* f1[comp*dim+d] = 0.5*(gradU[comp*dim+d] + gradU[d*dim+comp]); */
254:       f1[comp*dim+d] = gradU[comp*dim+d];
255:     }
256:     f1[comp*dim+comp] -= u[Ncomp];
257:   }
258: }

260: /* 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} */
261: void f0_p(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
262: {
263:   const PetscInt dim = SPATIAL_DIM_0;
264:   PetscInt       d;

266:   f0[0] = 0.0;
267:   for (d = 0; d < dim; ++d) f0[0] += gradU[d*dim+d];
268: }

270: void f1_p(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
271: {
272:   const PetscInt dim = SPATIAL_DIM_0;
273:   PetscInt       d;

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

278: void f0_T(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
279: {
280:   f0[0] = 0.0;
281: }

283: void f1_T(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
284: {
285:   const PetscInt dim = SPATIAL_DIM_2;
286:   const PetscInt off = SPATIAL_DIM_0*NUM_BASIS_COMPONENTS_0+SPATIAL_DIM_1*NUM_BASIS_COMPONENTS_1;
287:   PetscInt       d;

289:   for (d = 0; d < dim; ++d) f1[d] = gradU[off+d];
290: }

292: /* < v_t, I t > */
293: void g0_TT(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[])
294: {
295:   g0[0] = 1.0;
296: }

298: /* < q, \nabla\cdot v >
299:    NcompI = 1, NcompJ = dim */
300: void g1_pu(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[])
301: {
302:   const PetscInt dim = SPATIAL_DIM_0;
303:   PetscInt       d;

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

308: /* -< \nabla\cdot v, p >
309:     NcompI = dim, NcompJ = 1 */
310: void g2_up(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[])
311: {
312:   const PetscInt dim = SPATIAL_DIM_0;
313:   PetscInt       d;

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

318: /* < \nabla v, \nabla u + {\nabla u}^T >
319:    This just gives \nabla u, give the perdiagonal for the transpose */
320: void g3_uu(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[])
321: {
322:   const PetscInt dim   = SPATIAL_DIM_0;
323:   const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
324:   PetscInt       compI, d;

326:   for (compI = 0; compI < Ncomp; ++compI) {
327:     for (d = 0; d < dim; ++d) {
328:       g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
329:     }
330:   }
331: }

333: /* < \nabla t, \nabla T + {\nabla u}^T >
334:    This just gives \nabla T, give the perdiagonal for the transpose */
335: void g3_TT(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[])
336: {
337:   const PetscInt dim   = SPATIAL_DIM_2;
338:   const PetscInt Ncomp = NUM_BASIS_COMPONENTS_2;
339:   PetscInt       compI, d;

341:   for (compI = 0; compI < Ncomp; ++compI) {
342:     for (d = 0; d < dim; ++d) {
343:       g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
344:     }
345:   }
346: }

348: /*
349:   In 3D we use exact solution:

351:     u = x^2 + y^2
352:     v = y^2 + z^2
353:     w = x^2 + y^2 - 2(x+y)z
354:     p = x + y + z - 3/2
355:     f_x = f_y = f_z = 3

357:   so that

359:     -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
360:     \nabla \cdot u           = 2x + 2y - 2(x + y)                   = 0
361: */
362: void quadratic_u_3d(const PetscReal x[], PetscScalar *u)
363: {
364:   *u = x[0]*x[0] + x[1]*x[1];
365: };

367: void quadratic_v_3d(const PetscReal x[], PetscScalar *v)
368: {
369:   *v = x[1]*x[1] + x[2]*x[2];
370: };

372: void quadratic_w_3d(const PetscReal x[], PetscScalar *w)
373: {
374:   *w = x[0]*x[0] + x[1]*x[1] - 2.0*(x[0] + x[1])*x[2];
375: };

377: void linear_p_3d(const PetscReal x[], PetscScalar *p)
378: {
379:   *p = x[0] + x[1] + x[2] - 1.5;
380: };

384: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
385: {
386:   const char     *bcTypes[2]      = {"dirichlet", "freeslip"};
387:   const char     *forcingTypes[3] = {"constant", "linear", "cubic"};
388:   const char     *runTypes[2]     = {"full", "test"};
389:   PetscInt       bc, forcing, run;

393:   options->debug           = 0;
394:   options->runType         = RUN_FULL;
395:   options->dim             = 2;
396:   options->interpolate     = PETSC_FALSE;
397:   options->refinementLimit = 0.0;
398:   options->bcType          = DIRICHLET;
399:   options->forcingType     = FORCING_CONSTANT;
400:   options->numBatches      = 1;
401:   options->numBlocks       = 1;
402:   options->jacobianMF      = PETSC_FALSE;
403:   options->showInitial     = PETSC_FALSE;
404:   options->showSolution    = PETSC_TRUE;

406:   options->fem.quad    = (PetscQuadrature*) &options->q;
407:   options->fem.f0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f0Funcs;
408:   options->fem.f1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f1Funcs;
409:   options->fem.g0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g0Funcs;
410:   options->fem.g1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g1Funcs;
411:   options->fem.g2Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g2Funcs;
412:   options->fem.g3Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g3Funcs;

414:   MPI_Comm_size(comm, &options->numProcs);
415:   MPI_Comm_rank(comm, &options->rank);
416:   PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
417:   PetscOptionsInt("-debug", "The debugging level", "ex31.c", options->debug, &options->debug, NULL);
418:   run  = options->runType;
419:   PetscOptionsEList("-run_type", "The run type", "ex31.c", runTypes, 2, runTypes[options->runType], &run, NULL);

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

423:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex31.c", options->dim, &options->dim, NULL);
424:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex31.c", options->interpolate, &options->interpolate, NULL);
425:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex31.c", options->refinementLimit, &options->refinementLimit, NULL);
426:   PetscStrcpy(options->partitioner, "chaco");
427:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
428:   bc   = options->bcType;
429:   PetscOptionsEList("-bc_type","Type of boundary condition","ex31.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);

431:   options->bcType = (BCType) bc;
432:   forcing         = options->forcingType;

434:   PetscOptionsEList("-forcing_type","Type of forcing function","ex31.c",forcingTypes,3,forcingTypes[options->forcingType],&forcing,NULL);

436:   options->forcingType = (ForcingType) forcing;

438:   PetscOptionsInt("-gpu_batches", "The number of cell batches per kernel", "ex31.c", options->numBatches, &options->numBatches, NULL);
439:   PetscOptionsInt("-gpu_blocks", "The number of concurrent blocks per kernel", "ex31.c", options->numBlocks, &options->numBlocks, NULL);
440:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex31.c", options->jacobianMF, &options->jacobianMF, NULL);
441:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex31.c", options->showInitial, &options->showInitial, NULL);
442:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex31.c", options->showSolution, &options->showSolution, NULL);
443:   PetscOptionsEnd();

445:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
446:   return(0);
447: };

451: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
452: {
453:   Vec            lv;
454:   PetscInt       p;
455:   PetscMPIInt    rank, numProcs;

459:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
460:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
461:   DMGetLocalVector(dm, &lv);
462:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
463:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
464:   PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
465:   for (p = 0; p < numProcs; ++p) {
466:     if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
467:     PetscBarrier((PetscObject) dm);
468:   }
469:   DMRestoreLocalVector(dm, &lv);
470:   return(0);
471: }

475: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
476: {
477:   PetscInt       dim             = user->dim;
478:   PetscBool      interpolate     = user->interpolate;
479:   PetscReal      refinementLimit = user->refinementLimit;
480:   const char     *partitioner    = user->partitioner;

484:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
485:   DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
486:   {
487:     DM refinedMesh     = NULL;
488:     DM distributedMesh = NULL;

490:     /* Refine mesh using a volume constraint */
491:     DMPlexSetRefinementLimit(*dm, refinementLimit);
492:     DMRefine(*dm, comm, &refinedMesh);
493:     if (refinedMesh) {
494:       DMDestroy(dm);
495:       *dm  = refinedMesh;
496:     }
497:     /* Distribute mesh over processes */
498:     DMPlexDistribute(*dm, partitioner, 0, &distributedMesh);
499:     if (distributedMesh) {
500:       DMDestroy(dm);
501:       *dm  = distributedMesh;
502:     }
503:   }
504:   DMSetFromOptions(*dm);
505:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
506:   user->dm = *dm;
507:   return(0);
508: }

512: PetscErrorCode PointOnBoundary_2D(const PetscScalar coords[], PetscBool onBd[])
513: {
514:   const PetscInt  corner = 0, bottom = 1, right = 2, top = 3, left = 4;
515:   const PetscReal eps    = 1.0e-10;

518:   onBd[bottom] = PetscAbsScalar(coords[1])       < eps ? PETSC_TRUE : PETSC_FALSE;
519:   onBd[right]  = PetscAbsScalar(coords[0] - 1.0) < eps ? PETSC_TRUE : PETSC_FALSE;
520:   onBd[top]    = PetscAbsScalar(coords[1] - 1.0) < eps ? PETSC_TRUE : PETSC_FALSE;
521:   onBd[left]   = PetscAbsScalar(coords[0])       < eps ? PETSC_TRUE : PETSC_FALSE;
522:   onBd[corner] = onBd[bottom] + onBd[right] + onBd[top] + onBd[left] > 1 ? PETSC_TRUE : PETSC_FALSE;
523:   return(0);
524: }

528: PetscErrorCode CreateBoundaryPointIS_Square(DM dm, PetscInt *numBoundaries, PetscInt **numBoundaryConstraints, IS **boundaryPoints, IS **constraintIndices)
529: {
530:   MPI_Comm       comm;
531:   PetscSection   coordSection;
532:   Vec            coordinates;
533:   PetscScalar    *coords;
534:   PetscInt       vStart, vEnd;
535:   IS             bcPoints;
536:   const PetscInt *points;
537:   const PetscInt corner               = 0, bottom = 1, right = 2, top = 3, left = 4;
538:   PetscInt       numBoundaryPoints[5] = {0, 0, 0, 0, 0}, bd, numPoints, p;
539:   PetscInt       *bdPoints[5], *idx;

543:   PetscObjectGetComm((PetscObject)dm,&comm);
544:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
545:   /* boundary 0: corners
546:      boundary 1: bottom
547:      boundary 2: right
548:      boundary 3: top
549:      boundary 4: left
550:   */
551:   *numBoundaries = 5;

553:   PetscMalloc(*numBoundaries * sizeof(PetscInt), numBoundaryConstraints);
554:   PetscMalloc(*numBoundaries * sizeof(IS), boundaryPoints);
555:   PetscMalloc(*numBoundaries * sizeof(IS), constraintIndices);

557:   /* Set number of constraints for each boundary */
558:   (*numBoundaryConstraints)[corner] = 2;
559:   (*numBoundaryConstraints)[bottom] = 1;
560:   (*numBoundaryConstraints)[right]  = 1;
561:   (*numBoundaryConstraints)[top]    = 1;
562:   (*numBoundaryConstraints)[left]   = 1;
563:   /* Set local constraint indices for each boundary */
564:   PetscMalloc((*numBoundaryConstraints)[corner] * sizeof(PetscInt), &idx);
565:   idx[0] = 0; idx[1] = 1;
566:   ISCreateGeneral(comm, (*numBoundaryConstraints)[corner], idx, PETSC_OWN_POINTER, &(*constraintIndices)[corner]);
567:   PetscMalloc((*numBoundaryConstraints)[bottom] * sizeof(PetscInt), &idx);
568:   idx[0] = 1;
569:   ISCreateGeneral(comm, (*numBoundaryConstraints)[bottom], idx, PETSC_OWN_POINTER, &(*constraintIndices)[bottom]);
570:   PetscMalloc((*numBoundaryConstraints)[right] * sizeof(PetscInt), &idx);
571:   idx[0] = 0;
572:   ISCreateGeneral(comm, (*numBoundaryConstraints)[right], idx, PETSC_OWN_POINTER, &(*constraintIndices)[right]);
573:   PetscMalloc((*numBoundaryConstraints)[top] * sizeof(PetscInt), &idx);
574:   idx[0] = 1;
575:   ISCreateGeneral(comm, (*numBoundaryConstraints)[top], idx, PETSC_OWN_POINTER, &(*constraintIndices)[top]);
576:   PetscMalloc((*numBoundaryConstraints)[left] * sizeof(PetscInt), &idx);
577:   idx[0] = 0;
578:   ISCreateGeneral(comm, (*numBoundaryConstraints)[left], idx, PETSC_OWN_POINTER, &(*constraintIndices)[left]);

580:   /* Count points on each boundary */
581:   DMPlexGetCoordinateSection(dm, &coordSection);
582:   DMGetCoordinatesLocal(dm, &coordinates);
583:   VecGetArray(coordinates, &coords);
584:   DMPlexGetStratumIS(dm, "marker", 1, &bcPoints);
585:   ISGetLocalSize(bcPoints, &numPoints);
586:   ISGetIndices(bcPoints, &points);
587:   for (p = 0; p < numPoints; ++p) {
588:     PetscBool onBd[5];
589:     PetscInt  off, bd;

591:     if ((points[p] >= vStart) && (points[p] < vEnd)) {
592:       PetscSectionGetOffset(coordSection, points[p], &off);
593:       PointOnBoundary_2D(&coords[off], onBd);
594:     } else {
595:       PetscInt *closure = NULL;
596:       PetscInt closureSize, q, r;

598:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
599:       /* Compress out non-vertices */
600:       for (q = 0, r = 0; q < closureSize*2; q += 2) {
601:         if ((closure[q] >= vStart) && (closure[q] < vEnd)) {
602:           closure[r] = closure[q];
603:           ++r;
604:         }
605:       }
606:       closureSize = r;
607:       for (q = 0; q < closureSize; ++q) {
608:         PetscSectionGetOffset(coordSection, closure[q], &off);
609:         PointOnBoundary_2D(&coords[off], onBd);
610:         if (!onBd[corner]) break;
611:       }
612:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
613:       if (q == closureSize) SETERRQ1(comm, PETSC_ERR_PLIB, "Cannot handle face %d which has every vertex on a corner", points[p]);
614:     }

616:     for (bd = 0; bd < 5; ++bd) {
617:       if (onBd[bd]) {
618:         ++numBoundaryPoints[bd];
619:         break;
620:       }
621:     }
622:   }
623:   /* Set points on each boundary */
624:   for (bd = 0; bd < 5; ++bd) {
625:     PetscMalloc(numBoundaryPoints[bd] * sizeof(PetscInt), &bdPoints[bd]);
626:     numBoundaryPoints[bd] = 0;
627:   }
628:   for (p = 0; p < numPoints; ++p) {
629:     PetscBool onBd[5];
630:     PetscInt  off, bd;

632:     if ((points[p] >= vStart) && (points[p] < vEnd)) {
633:       PetscSectionGetOffset(coordSection, points[p], &off);
634:       PointOnBoundary_2D(&coords[off], onBd);
635:     } else {
636:       PetscInt *closure = NULL;
637:       PetscInt closureSize, q, r;

639:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
640:       /* Compress out non-vertices */
641:       for (q = 0, r = 0; q < closureSize*2; q += 2) {
642:         if ((closure[q] >= vStart) && (closure[q] < vEnd)) {
643:           closure[r] = closure[q];
644:           ++r;
645:         }
646:       }
647:       closureSize = r;
648:       for (q = 0; q < closureSize; ++q) {
649:         PetscSectionGetOffset(coordSection, closure[q], &off);
650:         PointOnBoundary_2D(&coords[off], onBd);
651:         if (!onBd[corner]) break;
652:       }
653:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
654:       if (q == closureSize) SETERRQ1(comm, PETSC_ERR_PLIB, "Cannot handle face %d which has every vertex on a corner", points[p]);
655:     }

657:     for (bd = 0; bd < 5; ++bd) {
658:       if (onBd[bd]) {
659:         bdPoints[bd][numBoundaryPoints[bd]++] = points[p];
660:         break;
661:       }
662:     }
663:   }
664:   VecRestoreArray(coordinates, &coords);
665:   ISRestoreIndices(bcPoints, &points);
666:   ISDestroy(&bcPoints);
667:   for (bd = 0; bd < 5; ++bd) {
668:     ISCreateGeneral(comm, numBoundaryPoints[bd], bdPoints[bd], PETSC_OWN_POINTER, &(*boundaryPoints)[bd]);
669:   }
670:   return(0);
671: }

675: PetscErrorCode CreateBoundaryPointIS_Cube(DM dm, PetscInt *numBoundaries, PetscInt **numBoundaryConstraints, IS **boundaryPoints, IS **constraintIndices)
676: {
678:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Just lazy");
679:   return(0);
680: }

684: /* This will only work for the square/cube, but I think the interface is robust */
685: PetscErrorCode CreateBoundaryPointIS(DM dm, PetscInt *numBoundaries, PetscInt **numBoundaryConstraints, IS **boundaryPoints, IS **constraintIndices)
686: {
687:   PetscInt       dim;

691:   DMPlexGetDimension(dm, &dim);
692:   switch (dim) {
693:   case 2:
694:     CreateBoundaryPointIS_Square(dm, numBoundaries, numBoundaryConstraints, boundaryPoints, constraintIndices);
695:     break;
696:   case 3:
697:     CreateBoundaryPointIS_Cube(dm, numBoundaries, numBoundaryConstraints, boundaryPoints, constraintIndices);
698:     break;
699:   default:
700:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "No boundary creatin routine for dimension %d", dim);
701:   }
702:   return(0);
703: }

707: PetscErrorCode SetupQuadrature(AppCtx *user)
708: {
710:   user->fem.quad[0].numQuadPoints = NUM_QUADRATURE_POINTS_0;
711:   user->fem.quad[0].quadPoints    = points_0;
712:   user->fem.quad[0].quadWeights   = weights_0;
713:   user->fem.quad[0].numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
714:   user->fem.quad[0].numComponents = NUM_BASIS_COMPONENTS_0;
715:   user->fem.quad[0].basis         = Basis_0;
716:   user->fem.quad[0].basisDer      = BasisDerivatives_0;
717:   user->fem.quad[1].numQuadPoints = NUM_QUADRATURE_POINTS_1;
718:   user->fem.quad[1].quadPoints    = points_1;
719:   user->fem.quad[1].quadWeights   = weights_1;
720:   user->fem.quad[1].numBasisFuncs = NUM_BASIS_FUNCTIONS_1;
721:   user->fem.quad[1].numComponents = NUM_BASIS_COMPONENTS_1;
722:   user->fem.quad[1].basis         = Basis_1;
723:   user->fem.quad[1].basisDer      = BasisDerivatives_1;
724:   user->fem.quad[2].numQuadPoints = NUM_QUADRATURE_POINTS_2;
725:   user->fem.quad[2].quadPoints    = points_2;
726:   user->fem.quad[2].quadWeights   = weights_2;
727:   user->fem.quad[2].numBasisFuncs = NUM_BASIS_FUNCTIONS_2;
728:   user->fem.quad[2].numComponents = NUM_BASIS_COMPONENTS_2;
729:   user->fem.quad[2].basis         = Basis_2;
730:   user->fem.quad[2].basisDer      = BasisDerivatives_2;
731:   return(0);
732: }

736: /*
737:   There is a problem here with uninterpolated meshes. The index in numDof[] is not dimension in this case,
738:   but sieve depth.
739: */
740: PetscErrorCode SetupSection(DM dm, AppCtx *user)
741: {
742:   PetscSection   section;
743:   const PetscInt numFields           = NUM_FIELDS;
744:   PetscInt       dim                 = user->dim;
745:   PetscInt       numBC               = 0;
746:   PetscInt       numComp[NUM_FIELDS] = {NUM_BASIS_COMPONENTS_0, NUM_BASIS_COMPONENTS_1, NUM_BASIS_COMPONENTS_2};
747:   PetscInt       bcFields[2]         = {0, 2};
748:   IS             bcPoints[2]         = {NULL, NULL};
749:   PetscInt       numDof[NUM_FIELDS*(SPATIAL_DIM_0+1)];
750:   PetscInt       f, d;
751:   PetscBool      view;

755:   if (dim != SPATIAL_DIM_0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_0);
756:   if (dim != SPATIAL_DIM_1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_1);
757:   for (d = 0; d <= dim; ++d) {
758:     numDof[0*(dim+1)+d] = numDof_0[d];
759:     numDof[1*(dim+1)+d] = numDof_1[d];
760:     numDof[2*(dim+1)+d] = numDof_2[d];
761:   }
762:   for (f = 0; f < numFields; ++f) {
763:     for (d = 1; d < dim; ++d) {
764:       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.");
765:     }
766:   }
767:   if (user->bcType == FREE_SLIP) {
768:     PetscInt numBoundaries, b;
769:     PetscInt *numBoundaryConstraints;
770:     IS       *boundaryPoints, *constraintIndices;

772:     DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, &section);
773:     /* Velocity conditions */
774:     CreateBoundaryPointIS(dm, &numBoundaries, &numBoundaryConstraints, &boundaryPoints, &constraintIndices);
775:     for (b = 0; b < numBoundaries; ++b) {
776:       DMPlexCreateSectionBCDof(dm, 1, &bcFields[0], &boundaryPoints[b], numBoundaryConstraints[b], section);
777:     }
778:     /* Temperature conditions */
779:     DMPlexGetStratumIS(dm, "marker", 1, &bcPoints[0]);
780:     DMPlexCreateSectionBCDof(dm, 1, &bcFields[1], &bcPoints[0], PETSC_DETERMINE, section);
781:     PetscSectionSetUp(section);
782:     for (b = 0; b < numBoundaries; ++b) {
783:       DMPlexCreateSectionBCIndicesField(dm, bcFields[0], boundaryPoints[b], constraintIndices[b], section);
784:     }
785:     DMPlexCreateSectionBCIndicesField(dm, bcFields[1], bcPoints[0], NULL, section);
786:     DMPlexCreateSectionBCIndices(dm, section);
787:   } else {
788:     if (user->bcType == DIRICHLET) {
789:       numBC       = 2;
790:       DMPlexGetStratumIS(dm, "marker", 1, &bcPoints[0]);
791:       bcPoints[1] = bcPoints[0];
792:       PetscObjectReference((PetscObject) bcPoints[1]);
793:     }
794:     DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, &section);
795:   }
796:   PetscSectionSetFieldName(section, 0, "velocity");
797:   PetscSectionSetFieldName(section, 1, "pressure");
798:   PetscSectionSetFieldName(section, 2, "temperature");
799:   DMSetDefaultSection(dm, section);
800:   ISDestroy(&bcPoints[0]);
801:   ISDestroy(&bcPoints[1]);
802:   PetscOptionsHasName(((PetscObject) dm)->prefix, "-section_view", &view);
803:   if ((user->bcType == FREE_SLIP) && view) {
804:     PetscSection s, gs;

806:     DMGetDefaultSection(dm, &s);
807:     PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD);
808:     DMGetDefaultGlobalSection(dm, &gs);
809:     PetscSectionView(gs, PETSC_VIEWER_STDOUT_WORLD);
810:   }
811:   return(0);
812: }

816: PetscErrorCode SetupExactSolution(DM dm, AppCtx *user)
817: {
818:   PetscFEM       *fem = &user->fem;

822:   switch (user->forcingType) {
823:   case FORCING_CONSTANT:
824:     if (user->bcType == FREE_SLIP) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Constant forcing is incompatible with freeslip boundary conditions");
825:     fem->f0Funcs[0] = f0_u_constant;
826:     break;
827:   case FORCING_LINEAR:
828:     switch (user->bcType) {
829:     case DIRICHLET:
830:     case FREE_SLIP:
831:       switch (user->dim) {
832:       case 2:
833:         fem->f0Funcs[0] = f0_u_linear_2d;
834:         break;
835:       default:
836:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
837:       }
838:       break;
839:     default:
840:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
841:     }
842:     break;
843:   case FORCING_CUBIC:
844:     switch (user->bcType) {
845:     case DIRICHLET:
846:     case FREE_SLIP:
847:       switch (user->dim) {
848:       case 2:
849:         fem->f0Funcs[0] = f0_u_cubic_2d;
850:         break;
851:       default:
852:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
853:       }
854:       break;
855:     default:
856:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
857:     }
858:     break;
859:   }
860:   fem->f0Funcs[1] = f0_p;
861:   fem->f0Funcs[2] = f0_T;
862:   fem->f1Funcs[0] = f1_u;
863:   fem->f1Funcs[1] = f1_p;
864:   fem->f1Funcs[2] = f1_T;
865:   fem->g0Funcs[0] = NULL;
866:   fem->g0Funcs[1] = NULL;
867:   fem->g0Funcs[2] = NULL;
868:   fem->g0Funcs[3] = NULL;
869:   fem->g0Funcs[4] = NULL;
870:   fem->g0Funcs[5] = NULL;
871:   fem->g0Funcs[6] = NULL;
872:   fem->g0Funcs[7] = NULL;
873:   fem->g0Funcs[8] = NULL;
874:   fem->g1Funcs[0] = NULL;
875:   fem->g1Funcs[1] = NULL;
876:   fem->g1Funcs[2] = NULL;
877:   fem->g1Funcs[3] = g1_pu;      /* < q, \nabla\cdot v > */
878:   fem->g1Funcs[4] = NULL;
879:   fem->g1Funcs[5] = NULL;
880:   fem->g1Funcs[6] = NULL;
881:   fem->g1Funcs[7] = NULL;
882:   fem->g1Funcs[8] = NULL;
883:   fem->g2Funcs[0] = NULL;
884:   fem->g2Funcs[1] = g2_up;      /* < \nabla\cdot v, p > */
885:   fem->g2Funcs[2] = NULL;
886:   fem->g2Funcs[3] = NULL;
887:   fem->g2Funcs[4] = NULL;
888:   fem->g2Funcs[5] = NULL;
889:   fem->g2Funcs[6] = NULL;
890:   fem->g2Funcs[7] = NULL;
891:   fem->g2Funcs[8] = NULL;
892:   fem->g3Funcs[0] = g3_uu;      /* < \nabla v, \nabla u + {\nabla u}^T > */
893:   fem->g3Funcs[1] = NULL;
894:   fem->g3Funcs[2] = NULL;
895:   fem->g3Funcs[3] = NULL;
896:   fem->g3Funcs[4] = NULL;
897:   fem->g3Funcs[5] = NULL;
898:   fem->g3Funcs[6] = NULL;
899:   fem->g3Funcs[7] = NULL;
900:   fem->g3Funcs[8] = g3_TT;      /* < \nabla t, \nabla T + {\nabla T}^T > */
901:   switch (user->forcingType) {
902:   case FORCING_CONSTANT:
903:     switch (user->bcType) {
904:     case DIRICHLET:
905:       switch (user->dim) {
906:       case 2:
907:         user->exactFuncs[0] = quadratic_u_2d;
908:         user->exactFuncs[1] = quadratic_v_2d;
909:         user->exactFuncs[2] = linear_p_2d;
910:         user->exactFuncs[3] = linear_T_2d;
911:         break;
912:       case 3:
913:         user->exactFuncs[0] = quadratic_u_3d;
914:         user->exactFuncs[1] = quadratic_v_3d;
915:         user->exactFuncs[2] = quadratic_w_3d;
916:         user->exactFuncs[3] = linear_p_3d;
917:         user->exactFuncs[4] = linear_T_2d;
918:         break;
919:       default:
920:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
921:       }
922:       break;
923:     default:
924:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
925:     }
926:     break;
927:   case FORCING_LINEAR:
928:     switch (user->bcType) {
929:     case DIRICHLET:
930:       switch (user->dim) {
931:       case 2:
932:         user->exactFuncs[0] = cubic_u_2d;
933:         user->exactFuncs[1] = cubic_v_2d;
934:         user->exactFuncs[2] = linear_p_2d;
935:         user->exactFuncs[3] = linear_T_2d;
936:         break;
937:       default:
938:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
939:       }
940:       break;
941:     default:
942:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
943:     }
944:     break;
945:   case FORCING_CUBIC:
946:     switch (user->bcType) {
947:     case DIRICHLET:
948:     case FREE_SLIP:
949:       switch (user->dim) {
950:       case 2:
951:         user->exactFuncs[0] = quintic_u_2d;
952:         user->exactFuncs[1] = quintic_v_2d;
953:         user->exactFuncs[2] = linear_p_2d;
954:         user->exactFuncs[3] = linear_T_2d;
955:         break;
956:       default:
957:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
958:       }
959:       break;
960:     default:
961:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
962:     }
963:     break;
964:   default:
965:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid forcing type %d", user->forcingType);
966:   }
967:   DMPlexSetFEMIntegration(dm, FEMIntegrateResidualBatch, FEMIntegrateJacobianActionBatch, FEMIntegrateJacobianBatch);
968:   return(0);
969: }

973: /*
974:   . field - The field whose diagonal block (of the Jacobian) has this null space
975: */
976: PetscErrorCode CreateNullSpaces(DM dm, PetscInt field, MatNullSpace *nullSpace)
977: {
978:   AppCtx         *user;
979:   Vec            nullVec, localNullVec;
980:   PetscSection   section;
981:   PetscScalar    *a;
982:   PetscInt       pressure = field;
983:   PetscInt       pStart, pEnd, p;

987:   DMGetApplicationContext(dm, (void**) &user);
988:   DMGetGlobalVector(dm, &nullVec);
989:   DMGetLocalVector(dm, &localNullVec);
990:   VecSet(nullVec, 0.0);
991:   /* Put a constant in for all pressures */
992:   DMGetDefaultSection(dm, &section);
993:   PetscSectionGetChart(section, &pStart, &pEnd);
994:   VecGetArray(localNullVec, &a);
995:   for (p = pStart; p < pEnd; ++p) {
996:     PetscInt fDim, off, d;

998:     PetscSectionGetFieldDof(section, p, pressure, &fDim);
999:     PetscSectionGetFieldOffset(section, p, pressure, &off);
1000:     for (d = 0; d < fDim; ++d) a[off+d] = 1.0;
1001:   }
1002:   VecRestoreArray(localNullVec, &a);
1003:   DMLocalToGlobalBegin(dm, localNullVec, INSERT_VALUES, nullVec);
1004:   DMLocalToGlobalEnd(dm, localNullVec, INSERT_VALUES, nullVec);
1005:   DMRestoreLocalVector(dm, &localNullVec);
1006:   VecNormalize(nullVec, NULL);
1007:   if (user->debug) {
1008:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
1009:     VecView(nullVec, PETSC_VIEWER_STDOUT_WORLD);
1010:   }
1011:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &nullVec, nullSpace);
1012:   DMRestoreGlobalVector(dm, &nullVec);
1013:   return(0);
1014: }

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

1021:   Input Parameters:
1022: + mat - The Jacobian shell matrix
1023: - X  - Global input vector

1025:   Output Parameter:
1026: . Y  - Local output vector

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

1032: .seealso: FormJacobianActionLocal()
1033: */
1034: PetscErrorCode FormJacobianAction(Mat J, Vec X,  Vec Y)
1035: {
1036:   JacActionCtx   *ctx;
1037:   DM             dm;
1038:   Vec            dummy, localX, localY;
1039:   PetscInt       N, n;

1046:   MatShellGetContext(J, &ctx);
1047:   dm   = ctx->dm;

1049:   /* determine whether X = localX */
1050:   DMGetLocalVector(dm, &dummy);
1051:   DMGetLocalVector(dm, &localX);
1052:   DMGetLocalVector(dm, &localY);
1053:   /* TODO: THIS dummy restore is necessary here so that the first available local vector has boundary conditions in it
1054:    I think the right thing to do is have the user put BC into a local vector and give it to us
1055:   */
1056:   DMRestoreLocalVector(dm, &dummy);
1057:   VecGetSize(X, &N);
1058:   VecGetSize(localX, &n);

1060:   if (n != N) { /* X != localX */
1061:     VecSet(localX, 0.0);
1062:     DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1063:     DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1064:   } else {
1065:     DMRestoreLocalVector(dm, &localX);
1066:     localX = X;
1067:   }
1068:   DMPlexComputeJacobianActionFEM(dm, J, localX, localY, ctx->user);
1069:   if (n != N) {
1070:     DMRestoreLocalVector(dm, &localX);
1071:   }
1072:   VecSet(Y, 0.0);
1073:   DMLocalToGlobalBegin(dm, localY, ADD_VALUES, Y);
1074:   DMLocalToGlobalEnd(dm, localY, ADD_VALUES, Y);
1075:   DMRestoreLocalVector(dm, &localY);
1076:   if (0) {
1077:     Vec       r;
1078:     PetscReal norm;

1080:     VecDuplicate(X, &r);
1081:     MatMult(ctx->J, X, r);
1082:     VecAXPY(r, -1.0, Y);
1083:     VecNorm(r, NORM_2, &norm);
1084:     if (norm > 1.0e-8) {
1085:       PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Input:\n");
1086:       VecView(X, PETSC_VIEWER_STDOUT_WORLD);
1087:       PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Result:\n");
1088:       VecView(Y, PETSC_VIEWER_STDOUT_WORLD);
1089:       PetscPrintf(PETSC_COMM_WORLD, "Difference:\n");
1090:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
1091:       SETERRQ1(PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "The difference with assembled multiply is too large %g", norm);
1092:     }
1093:     VecDestroy(&r);
1094:   }
1095:   return(0);
1096: }

1100: int main(int argc, char **argv)
1101: {
1102:   MPI_Comm       comm;
1103:   SNES           snes;                 /* nonlinear solver */
1104:   Vec            u,r;                  /* solution, residual vectors */
1105:   Mat            A,J;                  /* Jacobian matrix */
1106:   MatNullSpace   nullSpace = 0;            /* May be necessary for pressure */
1107:   AppCtx         user;                 /* user-defined work context */
1108:   JacActionCtx   userJ;                /* context for Jacobian MF action */
1109:   PetscInt       its;                  /* iterations for convergence */
1110:   PetscReal      error         = 0.0;  /* L_2 error in the solution */
1111:   const PetscInt numComponents = NUM_BASIS_COMPONENTS_TOTAL;

1114:   PetscInitialize(&argc, &argv, NULL, help);
1115:   comm = PETSC_COMM_WORLD;
1116:   ProcessOptions(comm, &user);
1117:   SNESCreate(comm, &snes);
1118:   CreateMesh(comm, &user, &user.dm);
1119:   SNESSetDM(snes, user.dm);
1120:   DMSetApplicationContext(user.dm, &user);

1122:   SetupExactSolution(user.dm, &user);
1123:   SetupQuadrature(&user);
1124:   SetupSection(user.dm, &user);

1126:   DMCreateGlobalVector(user.dm, &u);
1127:   PetscObjectSetName((PetscObject) u, "solution");
1128:   VecDuplicate(u, &r);

1130:   DMCreateMatrix(user.dm, MATAIJ, &J);
1131:   if (user.jacobianMF) {
1132:     PetscInt M, m, N, n;

1134:     MatGetSize(J, &M, &N);
1135:     MatGetLocalSize(J, &m, &n);
1136:     MatCreate(comm, &A);
1137:     MatSetSizes(A, m, n, M, N);
1138:     MatSetType(A, MATSHELL);
1139:     MatSetUp(A);
1140:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) FormJacobianAction);

1142:     userJ.dm   = user.dm;
1143:     userJ.J    = J;
1144:     userJ.user = &user;

1146:     DMCreateLocalVector(user.dm, &userJ.u);
1147:     MatShellSetContext(A, &userJ);
1148:   } else {
1149:     A = J;
1150:   }
1151:   DMSetNullSpaceConstructor(user.dm, 1, CreateNullSpaces);

1153:   DMSNESSetFunctionLocal(user.dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);
1154:   DMSNESSetJacobianLocal(user.dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*))DMPlexComputeJacobianFEM,&user);

1156:   SNESSetFromOptions(snes);

1158:   {
1159:     KSP               ksp; PC pc; Vec crd_vec;
1160:     const PetscScalar *v;
1161:     PetscInt          i,k,j,mlocal;
1162:     PetscReal         *coords;

1164:     SNESGetKSP(snes, &ksp);
1165:     KSPGetPC(ksp, &pc);
1166:     DMGetCoordinatesLocal(user.dm, &crd_vec);
1167:     VecGetLocalSize(crd_vec,&mlocal);
1168:     PetscMalloc(SPATIAL_DIM_0*mlocal*sizeof(*coords),&coords);
1169:     VecGetArrayRead(crd_vec,&v);
1170:     for (k=j=0; j<mlocal; j++) {
1171:       for (i=0; i<SPATIAL_DIM_0; i++,k++) {
1172:         coords[k] = PetscRealPart(v[k]);
1173:       }
1174:     }
1175:     VecRestoreArrayRead(crd_vec,&v);
1176:     PCSetCoordinates(pc, SPATIAL_DIM_0, mlocal, coords);
1177:     PetscFree(coords);
1178:   }

1180:   DMPlexProjectFunction(user.dm, numComponents, user.exactFuncs, INSERT_ALL_VALUES, u);
1181:   if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
1182:   if (user.runType == RUN_FULL) {
1183:     PetscScalar (*initialGuess[numComponents])(const PetscReal x[]);
1184:     PetscInt c;

1186:     for (c = 0; c < numComponents; ++c) initialGuess[c] = zero;
1187:     DMPlexProjectFunction(user.dm, numComponents, initialGuess, INSERT_VALUES, u);
1188:     if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
1189:     if (user.debug) {
1190:       PetscPrintf(comm, "Initial guess\n");
1191:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
1192:     }
1193:     SNESSolve(snes, NULL, u);
1194:     SNESGetIterationNumber(snes, &its);
1195:     PetscPrintf(comm, "Number of SNES iterations = %D\n", its);
1196:     DMPlexComputeL2Diff(user.dm, user.q, user.exactFuncs, u, &error);
1197:     PetscPrintf(comm, "L_2 Error: %.3g\n", error);
1198:     if (user.showSolution) {
1199:       PetscPrintf(comm, "Solution\n");
1200:       VecChop(u, 3.0e-9);
1201:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
1202:     }
1203:   } else {
1204:     PetscReal res = 0.0;

1206:     /* Check discretization error */
1207:     PetscPrintf(comm, "Initial guess\n");
1208:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
1209:     DMPlexComputeL2Diff(user.dm, user.q, user.exactFuncs, u, &error);
1210:     PetscPrintf(comm, "L_2 Error: %g\n", error);
1211:     /* Check residual */
1212:     SNESComputeFunction(snes, u, r);
1213:     PetscPrintf(comm, "Initial Residual\n");
1214:     VecChop(r, 1.0e-10);
1215:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
1216:     VecNorm(r, NORM_2, &res);
1217:     PetscPrintf(comm, "L_2 Residual: %g\n", res);
1218:     /* Check Jacobian */
1219:     {
1220:       Vec          b;
1221:       MatStructure flag;
1222:       MatNullSpace nullSpace2;
1223:       PetscBool    isNull;

1225:       CreateNullSpaces(user.dm, 1, &nullSpace2);
1226:       MatNullSpaceTest(nullSpace2, J, &isNull);
1227:       if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1228:       MatNullSpaceDestroy(&nullSpace2);

1230:       SNESComputeJacobian(snes, u, &A, &A, &flag);
1231:       VecDuplicate(u, &b);
1232:       VecSet(r, 0.0);
1233:       SNESComputeFunction(snes, r, b);
1234:       MatMult(A, u, r);
1235:       VecAXPY(r, 1.0, b);
1236:       VecDestroy(&b);
1237:       PetscPrintf(comm, "Au - b = Au + F(0)\n");
1238:       VecChop(r, 1.0e-10);
1239:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
1240:       VecNorm(r, NORM_2, &res);
1241:       PetscPrintf(comm, "Linear L_2 Residual: %g\n", res);
1242:     }
1243:   }

1245:   if (user.runType == RUN_FULL) {
1246:     PetscContainer c;
1247:     PetscSection   section;
1248:     Vec            sol;
1249:     PetscViewer    viewer;
1250:     const char     *name;

1252:     PetscViewerCreate(comm, &viewer);
1253:     PetscViewerSetType(viewer, PETSCVIEWERVTK);
1254:     PetscViewerFileSetName(viewer, "ex31_sol.vtk");
1255:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
1256:     DMGetLocalVector(user.dm, &sol);
1257:     PetscObjectGetName((PetscObject) u, &name);
1258:     PetscObjectSetName((PetscObject) sol, name);
1259:     DMGlobalToLocalBegin(user.dm, u, INSERT_VALUES, sol);
1260:     DMGlobalToLocalEnd(user.dm, u, INSERT_VALUES, sol);
1261:     DMGetDefaultSection(user.dm, &section);
1262:     PetscObjectReference((PetscObject) user.dm); /* Needed because viewer destroys the DM */
1263:     PetscViewerVTKAddField(viewer, (PetscObject) user.dm, DMPlexVTKWriteAll, PETSC_VTK_POINT_FIELD, (PetscObject) sol);
1264:     PetscObjectReference((PetscObject) sol); /* Needed because viewer destroys the Vec */
1265:     PetscContainerCreate(comm, &c);
1266:     PetscContainerSetPointer(c, section);
1267:     PetscObjectCompose((PetscObject) sol, "section", (PetscObject) c);
1268:     PetscContainerDestroy(&c);
1269:     DMRestoreLocalVector(user.dm, &sol);
1270:     PetscViewerDestroy(&viewer);
1271:   }

1273:   MatNullSpaceDestroy(&nullSpace);
1274:   if (user.jacobianMF) {
1275:     VecDestroy(&userJ.u);
1276:   }
1277:   if (A != J) {
1278:     MatDestroy(&A);
1279:   }
1280:   MatDestroy(&J);
1281:   VecDestroy(&u);
1282:   VecDestroy(&r);
1283:   SNESDestroy(&snes);
1284:   DMDestroy(&user.dm);
1285:   PetscFinalize();
1286:   return 0;
1287: }