Actual source code: ex62.c
petsc-3.4.5 2014-06-29
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, §ion);
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, §ion);
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: }