Actual source code: ex3.c
petsc-3.7.7 2017-09-25
1: static char help[] = "Check that a DM can accurately represent and interpolate functions of a given polynomial order\n\n";
3: #include <petscdmplex.h>
4: #include <petscdm.h>
5: #include <petscdmda.h>
6: #include <petscfe.h>
7: #include <petscds.h>
8: #include <petscksp.h>
9: #include <petscsnes.h>
11: typedef struct {
12: PetscInt debug; /* The debugging level */
13: /* Domain and mesh definition */
14: PetscInt dim; /* The topological mesh dimension */
15: PetscBool simplex; /* Flag for simplex or tensor product mesh */
16: PetscBool useDA; /* Flag DMDA tensor product mesh */
17: PetscBool interpolate; /* Generate intermediate mesh elements */
18: PetscReal refinementLimit; /* The largest allowable cell volume */
19: /* Element definition */
20: PetscInt qorder; /* Order of the quadrature */
21: PetscInt numComponents; /* Number of field components */
22: PetscFE fe; /* The finite element */
23: /* Testing space */
24: PetscInt porder; /* Order of polynomials to test */
25: PetscBool convergence; /* Test for order of convergence */
26: PetscBool convRefine; /* Test for convergence using refinement, otherwise use coarsening */
27: PetscBool constraints; /* Test local constraints */
28: PetscBool tree; /* Test tree routines */
29: PetscBool testFEjacobian; /* Test finite element Jacobian assembly */
30: PetscBool testFVgrad; /* Test finite difference gradient routine */
31: PetscBool testInjector; /* Test finite element injection routines */
32: PetscInt treeCell; /* Cell to refine in tree test */
33: PetscReal constants[3]; /* Constant values for each dimension */
34: } AppCtx;
36: /* u = 1 */
37: PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
38: {
39: AppCtx *user = (AppCtx *) ctx;
40: PetscInt d;
41: for (d = 0; d < user->dim; ++d) u[d] = user->constants[d];
42: return 0;
43: }
44: PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
45: {
46: AppCtx *user = (AppCtx *) ctx;
47: PetscInt d;
48: for (d = 0; d < user->dim; ++d) u[d] = 0.0;
49: return 0;
50: }
52: /* u = x */
53: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
54: {
55: AppCtx *user = (AppCtx *) ctx;
56: PetscInt d;
57: for (d = 0; d < user->dim; ++d) u[d] = coords[d];
58: return 0;
59: }
60: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
61: {
62: AppCtx *user = (AppCtx *) ctx;
63: PetscInt d, e;
64: for (d = 0; d < user->dim; ++d) {
65: u[d] = 0.0;
66: for (e = 0; e < user->dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
67: }
68: return 0;
69: }
71: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
72: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
73: {
74: AppCtx *user = (AppCtx *) ctx;
75: if (user->dim > 2) {u[0] = coords[0]*coords[1]; u[1] = coords[1]*coords[2]; u[2] = coords[2]*coords[0];}
76: else if (user->dim > 1) {u[0] = coords[0]*coords[0]; u[1] = coords[0]*coords[1];}
77: else if (user->dim > 0) {u[0] = coords[0]*coords[0];}
78: return 0;
79: }
80: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
81: {
82: AppCtx *user = (AppCtx *) ctx;
83: if (user->dim > 2) {u[0] = coords[1]*n[0] + coords[0]*n[1]; u[1] = coords[2]*n[1] + coords[1]*n[2]; u[2] = coords[2]*n[0] + coords[0]*n[2];}
84: else if (user->dim > 1) {u[0] = 2.0*coords[0]*n[0]; u[1] = coords[1]*n[0] + coords[0]*n[1];}
85: else if (user->dim > 0) {u[0] = 2.0*coords[0]*n[0];}
86: return 0;
87: }
89: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
90: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
91: {
92: AppCtx *user = (AppCtx *) ctx;
93: if (user->dim > 2) {u[0] = coords[0]*coords[0]*coords[1]; u[1] = coords[1]*coords[1]*coords[2]; u[2] = coords[2]*coords[2]*coords[0];}
94: else if (user->dim > 1) {u[0] = coords[0]*coords[0]*coords[0]; u[1] = coords[0]*coords[0]*coords[1];}
95: else if (user->dim > 0) {u[0] = coords[0]*coords[0]*coords[0];}
96: return 0;
97: }
98: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
99: {
100: AppCtx *user = (AppCtx *) ctx;
101: if (user->dim > 2) {u[0] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1]; u[1] = 2.0*coords[1]*coords[2]*n[1] + coords[1]*coords[1]*n[2]; u[2] = 2.0*coords[2]*coords[0]*n[2] + coords[2]*coords[2]*n[0];}
102: else if (user->dim > 1) {u[0] = 3.0*coords[0]*coords[0]*n[0]; u[1] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1];}
103: else if (user->dim > 0) {u[0] = 3.0*coords[0]*coords[0]*n[0];}
104: return 0;
105: }
107: /* u = tanh(x) */
108: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
109: {
110: AppCtx *user = (AppCtx *) ctx;
111: PetscInt d;
112: for (d = 0; d < user->dim; ++d) u[d] = tanh(coords[d] - 0.5);
113: return 0;
114: }
115: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
116: {
117: AppCtx *user = (AppCtx *) ctx;
118: PetscInt d;
119: for (d = 0; d < user->dim; ++d) u[d] = 1.0/PetscSqr(cosh(coords[d] - 0.5)) * n[d];
120: return 0;
121: }
125: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
126: {
130: options->debug = 0;
131: options->dim = 2;
132: options->simplex = PETSC_TRUE;
133: options->useDA = PETSC_TRUE;
134: options->interpolate = PETSC_TRUE;
135: options->refinementLimit = 0.0;
136: options->qorder = 0;
137: options->numComponents = 1;
138: options->porder = 0;
139: options->convergence = PETSC_FALSE;
140: options->convRefine = PETSC_TRUE;
141: options->constraints = PETSC_FALSE;
142: options->tree = PETSC_FALSE;
143: options->treeCell = 0;
144: options->testFEjacobian = PETSC_FALSE;
145: options->testFVgrad = PETSC_FALSE;
146: options->testInjector = PETSC_FALSE;
148: PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
149: PetscOptionsInt("-debug", "The debugging level", "ex3.c", options->debug, &options->debug, NULL);
150: PetscOptionsInt("-dim", "The topological mesh dimension", "ex3.c", options->dim, &options->dim, NULL);
151: PetscOptionsBool("-simplex", "Flag for simplices or hexhedra", "ex3.c", options->simplex, &options->simplex, NULL);
152: PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL);
153: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex3.c", options->interpolate, &options->interpolate, NULL);
154: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex3.c", options->refinementLimit, &options->refinementLimit, NULL);
155: PetscOptionsInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL);
156: PetscOptionsInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL);
157: PetscOptionsInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL);
158: PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL);
159: PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL);
160: PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL);
161: PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL);
162: PetscOptionsInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL);
163: PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL);
164: PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL);
165: PetscOptionsBool("-test_injector","Test finite element injection", "ex3.c", options->testInjector, &options->testInjector,NULL);
166: PetscOptionsEnd();
168: return(0);
169: }
173: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
174: {
175: PetscInt dim = user->dim;
176: PetscBool interpolate = user->interpolate;
177: PetscReal refinementLimit = user->refinementLimit;
178: PetscBool isPlex;
182: if (user->simplex) {
183: DM refinedMesh = NULL;
185: DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
186: /* Refine mesh using a volume constraint */
187: DMPlexSetRefinementLimit(*dm, refinementLimit);
188: DMRefine(*dm, comm, &refinedMesh);
189: if (refinedMesh) {
190: DMDestroy(dm);
191: *dm = refinedMesh;
192: }
193: DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
194: } else {
195: if (user->constraints || user->tree || !user->useDA) {
196: PetscInt cells[3] = {2, 2, 2};
198: PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&cells[0],NULL);
199: PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&cells[1],NULL);
200: PetscOptionsGetInt(NULL,NULL,"-da_grid_z",&cells[2],NULL);
201: DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
202: } else {
203: switch (user->dim) {
204: case 2:
205: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);
206: DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
207: break;
208: default:
209: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %d", dim);
210: }
211: PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
212: }
213: }
214: PetscObjectTypeCompare((PetscObject)*dm,DMPLEX,&isPlex);
215: if (isPlex) {
216: DM distributedMesh = NULL;
217: if (user->tree) {
218: DM refTree;
219: DM ncdm = NULL;
221: DMPlexCreateDefaultReferenceTree(comm,user->dim,user->simplex,&refTree);
222: DMPlexSetReferenceTree(*dm,refTree);
223: DMDestroy(&refTree);
224: DMPlexTreeRefineCell(*dm,user->treeCell,&ncdm);
225: if (ncdm) {
226: DMDestroy(dm);
227: *dm = ncdm;
228: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
229: }
230: }
231: else {
232: DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
233: }
234: /* Distribute mesh over processes */
235: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
236: if (distributedMesh) {
237: DMDestroy(dm);
238: *dm = distributedMesh;
239: }
240: if (user->simplex) {
241: PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");
242: }
243: else {
244: PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
245: }
246: }
247: DMSetFromOptions(*dm);
248: DMViewFromOptions(*dm,NULL,"-dm_view");
249: return(0);
250: }
254: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux,
255: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
256: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
257: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
258: {
259: PetscInt d, e;
260: for (d = 0, e = 0; d < dim; d++, e+=dim+1) {
261: g0[e] = 1.;
262: }
263: }
265: /* < \nabla v, 1/2(\nabla u + {\nabla u}^T) > */
268: static void symmetric_gradient_inner_product(PetscInt dim, PetscInt Nf, PetscInt NfAux,
269: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
270: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
271: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar C[])
272: {
273: PetscInt compI, compJ, d, e;
275: for (compI = 0; compI < dim; ++compI) {
276: for (compJ = 0; compJ < dim; ++compJ) {
277: for (d = 0; d < dim; ++d) {
278: for (e = 0; e < dim; e++) {
279: if (d == e && d == compI && d == compJ) {
280: C[((compI*dim+compJ)*dim+d)*dim+e] = 1.0;
281: }
282: else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
283: C[((compI*dim+compJ)*dim+d)*dim+e] = 0.5;
284: }
285: else {
286: C[((compI*dim+compJ)*dim+d)*dim+e] = 0.0;
287: }
288: }
289: }
290: }
291: }
292: }
296: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
297: {
298: PetscBool isPlex;
302: if (!user->simplex && user->constraints) {
303: /* test local constraints */
304: DM coordDM;
305: PetscInt fStart, fEnd, f, vStart, vEnd, v;
306: PetscInt edgesx = 2, vertsx;
307: PetscInt edgesy = 2, vertsy;
308: PetscMPIInt size;
309: PetscInt numConst;
310: PetscSection aSec;
311: PetscInt *anchors;
312: PetscInt offset;
313: IS aIS;
314: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
316: MPI_Comm_size(comm,&size);
317: if (size > 1) SETERRQ(comm,PETSC_ERR_SUP,"Local constraint test can only be performed in serial");
319: /* we are going to test constraints by using them to enforce periodicity
320: * in one direction, and comparing to the existing method of enforcing
321: * periodicity */
323: /* first create the coordinate section so that it does not clone the
324: * constraints */
325: DMGetCoordinateDM(dm,&coordDM);
327: /* create the constrained-to-anchor section */
328: DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
329: DMPlexGetDepthStratum(dm,1,&fStart,&fEnd);
330: PetscSectionCreate(PETSC_COMM_SELF,&aSec);
331: PetscSectionSetChart(aSec,PetscMin(fStart,vStart),PetscMax(fEnd,vEnd));
333: /* define the constraints */
334: PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&edgesx,NULL);
335: PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&edgesy,NULL);
336: vertsx = edgesx + 1;
337: vertsy = edgesy + 1;
338: numConst = vertsy + edgesy;
339: PetscMalloc1(numConst,&anchors);
340: offset = 0;
341: for (v = vStart + edgesx; v < vEnd; v+= vertsx) {
342: PetscSectionSetDof(aSec,v,1);
343: anchors[offset++] = v - edgesx;
344: }
345: for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
346: PetscSectionSetDof(aSec,f,1);
347: anchors[offset++] = f - edgesx * edgesy;
348: }
349: PetscSectionSetUp(aSec);
350: ISCreateGeneral(PETSC_COMM_SELF,numConst,anchors,PETSC_OWN_POINTER,&aIS);
352: DMPlexSetAnchors(dm,aSec,aIS);
353: PetscSectionDestroy(&aSec);
354: ISDestroy(&aIS);
355: }
356: DMSetNumFields(dm, 1);
357: DMSetField(dm, 0, (PetscObject) user->fe);
358: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
359: if (!isPlex) {
360: PetscSection section;
361: const PetscInt *numDof;
362: PetscInt numComp;
364: PetscFEGetNumComponents(user->fe, &numComp);
365: PetscFEGetNumDof(user->fe, &numDof);
366: DMDACreateSection(dm, &numComp, numDof, NULL, §ion);
367: DMSetDefaultSection(dm, section);
368: PetscSectionDestroy(§ion);
369: }
370: if (!user->simplex && user->constraints) {
371: /* test getting local constraint matrix that matches section */
372: PetscSection aSec;
373: IS aIS;
375: DMPlexGetAnchors(dm,&aSec,&aIS);
376: if (aSec) {
377: PetscDS ds;
378: PetscSection cSec, section;
379: PetscInt cStart, cEnd, c, numComp;
380: Mat cMat, mass;
381: Vec local;
382: const PetscInt *anchors;
384: DMGetDefaultSection(dm,§ion);
385: /* this creates the matrix and preallocates the matrix structure: we
386: * just have to fill in the values */
387: DMGetDefaultConstraints(dm,&cSec,&cMat);
388: PetscSectionGetChart(cSec,&cStart,&cEnd);
389: ISGetIndices(aIS,&anchors);
390: PetscFEGetNumComponents(user->fe, &numComp);
391: for (c = cStart; c < cEnd; c++) {
392: PetscInt cDof;
394: /* is this point constrained? (does it have an anchor?) */
395: PetscSectionGetDof(aSec,c,&cDof);
396: if (cDof) {
397: PetscInt cOff, a, aDof, aOff, j;
398: if (cDof != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Found %d anchor points: should be just one",cDof);
400: /* find the anchor point */
401: PetscSectionGetOffset(aSec,c,&cOff);
402: a = anchors[cOff];
404: /* find the constrained dofs (row in constraint matrix) */
405: PetscSectionGetDof(cSec,c,&cDof);
406: PetscSectionGetOffset(cSec,c,&cOff);
408: /* find the anchor dofs (column in constraint matrix) */
409: PetscSectionGetDof(section,a,&aDof);
410: PetscSectionGetOffset(section,a,&aOff);
412: if (cDof != aDof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point and anchor have different number of dofs: %d, %d\n",cDof,aDof);
413: if (cDof % numComp) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point dofs not divisible by field components: %d, %d\n",cDof,numComp);
415: /* put in a simple equality constraint */
416: for (j = 0; j < cDof; j++) {
417: MatSetValue(cMat,cOff+j,aOff+j,1.,INSERT_VALUES);
418: }
419: }
420: }
421: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
422: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
423: ISRestoreIndices(aIS,&anchors);
425: /* Now that we have constructed the constraint matrix, any FE matrix
426: * that we construct will apply the constraints during construction */
428: DMCreateMatrix(dm,&mass);
429: /* get a dummy local variable to serve as the solution */
430: DMGetLocalVector(dm,&local);
431: DMGetDS(dm,&ds);
432: /* set the jacobian to be the mass matrix */
433: PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL);
434: /* build the mass matrix */
435: DMPlexSNESComputeJacobianFEM(dm,local,mass,mass,NULL);
436: MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
437: MatDestroy(&mass);
438: DMRestoreLocalVector(dm,&local);
439: #if 0
440: {
441: /* compare this to periodicity with DMDA: this is broken right now
442: * because DMCreateMatrix() doesn't respect the default section that I
443: * set */
444: DM dmda;
445: PetscSection section;
446: const PetscInt *numDof;
447: PetscInt numComp;
449: /* periodic x */
450: DMDACreate2d(PetscObjectComm((PetscObject)dm), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dmda);
451: DMDASetVertexCoordinates(dmda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
454: PetscFEGetNumComponents(user->fe, &numComp);
455: PetscFEGetNumDof(user->fe, &numDof);
456: DMDACreateSection(dmda, &numComp, numDof, NULL, §ion);
457: DMSetDefaultSection(dmda, section);
458: PetscSectionDestroy(§ion);
459: DMCreateMatrix(dmda,&mass);
460: /* there isn't a DMDA equivalent of DMPlexSNESComputeJacobianFEM()
461: * right now, but we can at least verify the nonzero structure */
462: MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
463: MatDestroy(&mass);
464: DMDestroy(&dmda);
465: }
466: #endif
467: }
468: }
469: return(0);
470: }
474: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
475: {
476: PetscBool isPlex;
480: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
481: if (isPlex) {
482: Vec local;
483: Mat E;
484: MatNullSpace sp;
485: PetscBool isNullSpace;
486: PetscDS ds;
488: if (user->numComponents != user->dim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "The number of components %d must be equal to the dimension %d for this test", user->numComponents, user->dim);
489: DMGetDS(dm,&ds);
490: PetscDSSetJacobian(ds,0,0,NULL,NULL,NULL,symmetric_gradient_inner_product);
491: DMCreateMatrix(dm,&E);
492: DMGetLocalVector(dm,&local);
493: DMPlexSNESComputeJacobianFEM(dm,local,E,E,NULL);
494: DMPlexCreateRigidBody(dm,&sp);
495: MatNullSpaceTest(sp,E,&isNullSpace);
496: if (isNullSpace) {
497: PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: PASS\n");
498: }
499: else {
500: PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: FAIL\n");
501: }
502: MatNullSpaceDestroy(&sp);
503: MatDestroy(&E);
504: DMRestoreLocalVector(dm,&local);
505: }
506: return(0);
507: }
511: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
512: {
513: DM refTree;
514: PetscMPIInt rank;
518: DMPlexGetReferenceTree(dm,&refTree);
519: if (refTree) {
520: Mat inj;
522: DMPlexComputeInjectorReferenceTree(refTree,&inj);
523: PetscObjectSetName((PetscObject)inj,"Reference Tree Injector");
524: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
525: if (!rank) {
526: MatView(inj,PETSC_VIEWER_STDOUT_SELF);
527: }
528: MatDestroy(&inj);
529: }
530: return(0);
531: }
535: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
536: {
537: MPI_Comm comm;
538: DM dmRedist, dmfv, dmgrad, dmCell, refTree;
539: PetscFV fv;
540: PetscInt nvecs, v, cStart, cEnd, cEndInterior;
541: PetscMPIInt size;
542: Vec cellgeom, grad, locGrad;
543: const PetscScalar *cgeom;
544: PetscReal allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;
548: comm = PetscObjectComm((PetscObject)dm);
549: /* duplicate DM, give dup. a FV discretization */
550: DMPlexSetAdjacencyUseCone(dm,PETSC_TRUE);
551: DMPlexSetAdjacencyUseClosure(dm,PETSC_FALSE);
552: MPI_Comm_size(comm,&size);
553: dmRedist = NULL;
554: if (size > 1) {
555: DMPlexDistributeOverlap(dm,1,NULL,&dmRedist);
556: }
557: if (!dmRedist) {
558: dmRedist = dm;
559: PetscObjectReference((PetscObject)dmRedist);
560: }
561: PetscFVCreate(comm,&fv);
562: PetscFVSetType(fv,PETSCFVLEASTSQUARES);
563: PetscFVSetNumComponents(fv,user->numComponents);
564: PetscFVSetSpatialDimension(fv,user->dim);
565: PetscFVSetFromOptions(fv);
566: PetscFVSetUp(fv);
567: DMPlexConstructGhostCells(dmRedist,NULL,NULL,&dmfv);
568: DMDestroy(&dmRedist);
569: DMSetNumFields(dmfv,1);
570: DMSetField(dmfv, 0, (PetscObject) fv);
571: DMPlexGetReferenceTree(dm,&refTree);
572: if (refTree) {
573: PetscDS ds;
574: DMGetDS(dmfv,&ds);
575: DMSetDS(refTree,ds);
576: }
577: DMPlexSNESGetGradientDM(dmfv, fv, &dmgrad);
578: DMPlexGetHeightStratum(dmfv,0,&cStart,&cEnd);
579: nvecs = user->dim * (user->dim+1) / 2;
580: DMPlexSNESGetGeometryFVM(dmfv,NULL,&cellgeom,NULL);
581: VecGetDM(cellgeom,&dmCell);
582: VecGetArrayRead(cellgeom,&cgeom);
583: DMGetGlobalVector(dmgrad,&grad);
584: DMGetLocalVector(dmgrad,&locGrad);
585: DMPlexGetHybridBounds(dmgrad,&cEndInterior,NULL,NULL,NULL);
586: cEndInterior = (cEndInterior < 0) ? cEnd: cEndInterior;
587: for (v = 0; v < nvecs; v++) {
588: Vec locX;
589: PetscInt c;
590: PetscScalar trueGrad[3][3] = {{0.}};
591: const PetscScalar *gradArray;
592: PetscReal maxDiff, maxDiffGlob;
594: DMGetLocalVector(dmfv,&locX);
595: /* get the local projection of the rigid body mode */
596: for (c = cStart; c < cEnd; c++) {
597: PetscFVCellGeom *cg;
598: PetscScalar cx[3] = {0.,0.,0.};
600: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
601: if (v < user->dim) {
602: cx[v] = 1.;
603: }
604: else {
605: PetscInt w = v - user->dim;
607: cx[(w + 1) % user->dim] = cg->centroid[(w + 2) % user->dim];
608: cx[(w + 2) % user->dim] = -cg->centroid[(w + 1) % user->dim];
609: }
610: DMPlexVecSetClosure(dmfv,NULL,locX,c,cx,INSERT_ALL_VALUES);
611: }
612: /* TODO: this isn't in any header */
613: DMPlexReconstructGradientsFVM(dmfv,locX,grad);
614: DMGlobalToLocalBegin(dmgrad,grad,INSERT_VALUES,locGrad);
615: DMGlobalToLocalEnd(dmgrad,grad,INSERT_VALUES,locGrad);
616: VecGetArrayRead(locGrad,&gradArray);
617: /* compare computed gradient to exact gradient */
618: if (v >= user->dim) {
619: PetscInt w = v - user->dim;
621: trueGrad[(w + 1) % user->dim][(w + 2) % user->dim] = 1.;
622: trueGrad[(w + 2) % user->dim][(w + 1) % user->dim] = -1.;
623: }
624: maxDiff = 0.;
625: for (c = cStart; c < cEndInterior; c++) {
626: PetscScalar *compGrad;
627: PetscInt i, j, k;
628: PetscReal FrobDiff = 0.;
630: DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad);
632: for (i = 0, k = 0; i < user->dim; i++) {
633: for (j = 0; j < user->dim; j++, k++) {
634: PetscScalar diff = compGrad[k] - trueGrad[i][j];
635: FrobDiff += PetscRealPart(diff * PetscConj(diff));
636: }
637: }
638: FrobDiff = PetscSqrtReal(FrobDiff);
639: maxDiff = PetscMax(maxDiff,FrobDiff);
640: }
641: MPI_Allreduce(&maxDiff,&maxDiffGlob,1,MPIU_REAL,MPIU_MAX,comm);
642: allVecMaxDiff = PetscMax(allVecMaxDiff,maxDiffGlob);
643: VecRestoreArrayRead(locGrad,&gradArray);
644: DMRestoreLocalVector(dmfv,&locX);
645: }
646: if (allVecMaxDiff < fvTol) {
647: PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: PASS\n");
648: }
649: else {
650: PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n",fvTol,allVecMaxDiff);
651: }
652: DMRestoreLocalVector(dmgrad,&locGrad);
653: DMRestoreGlobalVector(dmgrad,&grad);
654: VecRestoreArrayRead(cellgeom,&cgeom);
655: DMDestroy(&dmfv);
656: PetscFVDestroy(&fv);
657: return(0);
658: }
662: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *),
663: PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
664: void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
665: {
666: Vec u;
667: PetscReal n[3] = {1.0, 1.0, 1.0};
671: DMGetGlobalVector(dm, &u);
672: /* Project function into FE function space */
673: DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
674: /* Compare approximation to exact in L_2 */
675: DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error);
676: DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer);
677: DMRestoreGlobalVector(dm, &u);
678: return(0);
679: }
683: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
684: {
685: PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
686: PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
687: void *exactCtxs[3];
688: MPI_Comm comm;
689: PetscReal error, errorDer, tol = PETSC_SMALL;
690: PetscErrorCode ierr;
693: exactCtxs[0] = user;
694: exactCtxs[1] = user;
695: exactCtxs[2] = user;
696: user->constants[0] = 1.0;
697: user->constants[1] = 2.0;
698: user->constants[2] = 3.0;
699: PetscObjectGetComm((PetscObject)dm, &comm);
700: /* Setup functions to approximate */
701: switch (order) {
702: case 0:
703: exactFuncs[0] = constant;
704: exactFuncDers[0] = constantDer;
705: break;
706: case 1:
707: exactFuncs[0] = linear;
708: exactFuncDers[0] = linearDer;
709: break;
710: case 2:
711: exactFuncs[0] = quadratic;
712: exactFuncDers[0] = quadraticDer;
713: break;
714: case 3:
715: exactFuncs[0] = cubic;
716: exactFuncDers[0] = cubicDer;
717: break;
718: default:
719: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %d order %d", user->dim, order);
720: }
721: ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
722: /* Report result */
723: if (error > tol) {PetscPrintf(comm, "Function tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol,(double) error);}
724: else {PetscPrintf(comm, "Function tests pass for order %D at tolerance %g\n", order, (double)tol);}
725: if (errorDer > tol) {PetscPrintf(comm, "Function tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
726: else {PetscPrintf(comm, "Function tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
727: return(0);
728: }
732: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
733: {
734: PetscErrorCode (*exactFuncs[1]) (PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
735: PetscErrorCode (*exactFuncDers[1]) (PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
736: PetscReal n[3] = {1.0, 1.0, 1.0};
737: void *exactCtxs[3];
738: DM rdm, idm, fdm;
739: Mat Interp;
740: Vec iu, fu, scaling;
741: MPI_Comm comm;
742: PetscInt dim = user->dim;
743: PetscReal error, errorDer, tol = 1.0e-10;
744: PetscBool isPlex, isDA;
745: PetscErrorCode ierr;
748: exactCtxs[0] = user;
749: exactCtxs[1] = user;
750: exactCtxs[2] = user;
751: user->constants[0] = 1.0;
752: user->constants[1] = 2.0;
753: user->constants[2] = 3.0;
754: PetscObjectGetComm((PetscObject)dm,&comm);
755: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
756: PetscObjectTypeCompare((PetscObject) dm, DMDA, &isDA);
757: DMRefine(dm, comm, &rdm);
758: DMSetCoarseDM(rdm, dm);
759: DMPlexSetRegularRefinement(rdm, user->convRefine);
760: if (user->tree && isPlex) {
761: DM refTree;
762: DMPlexGetReferenceTree(dm,&refTree);
763: DMPlexSetReferenceTree(rdm,refTree);
764: }
765: if (!user->simplex && !user->constraints) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
766: SetupSection(rdm, user);
767: /* Setup functions to approximate */
768: switch (order) {
769: case 0:
770: exactFuncs[0] = constant;
771: exactFuncDers[0] = constantDer;
772: break;
773: case 1:
774: exactFuncs[0] = linear;
775: exactFuncDers[0] = linearDer;
776: break;
777: case 2:
778: exactFuncs[0] = quadratic;
779: exactFuncDers[0] = quadraticDer;
780: break;
781: default:
782: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %D order %D", dim, order);
783: }
784: idm = checkRestrict ? rdm : dm;
785: fdm = checkRestrict ? dm : rdm;
786: DMGetGlobalVector(idm, &iu);
787: DMGetGlobalVector(fdm, &fu);
788: DMSetApplicationContext(dm, user);
789: DMSetApplicationContext(rdm, user);
790: DMCreateInterpolation(dm, rdm, &Interp, &scaling);
791: /* Project function into initial FE function space */
792: DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);
793: /* Interpolate function into final FE function space */
794: if (checkRestrict) {MatRestrict(Interp, iu, fu);VecPointwiseMult(fu, scaling, fu);}
795: else {MatInterpolate(Interp, iu, fu);}
796: /* Compare approximation to exact in L_2 */
797: DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error);
798: DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer);
799: /* Report result */
800: if (error > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol, (double)error);}
801: else {PetscPrintf(comm, "Interpolation tests pass for order %D at tolerance %g\n", order, (double)tol);}
802: if (errorDer > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
803: else {PetscPrintf(comm, "Interpolation tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
804: DMRestoreGlobalVector(idm, &iu);
805: DMRestoreGlobalVector(fdm, &fu);
806: MatDestroy(&Interp);
807: VecDestroy(&scaling);
808: DMDestroy(&rdm);
809: return(0);
810: }
814: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
815: {
816: DM odm = dm, rdm = NULL, cdm = NULL;
817: PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
818: PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
819: void *exactCtxs[3];
820: PetscInt r, c, cStart, cEnd;
821: PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
822: double p;
823: PetscErrorCode ierr;
826: if (!user->convergence) return(0);
827: exactCtxs[0] = user;
828: exactCtxs[1] = user;
829: exactCtxs[2] = user;
830: PetscObjectReference((PetscObject) odm);
831: if (!user->convRefine) {
832: for (r = 0; r < Nr; ++r) {
833: DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
834: DMDestroy(&odm);
835: odm = rdm;
836: }
837: SetupSection(odm, user);
838: }
839: ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user);
840: if (user->convRefine) {
841: for (r = 0; r < Nr; ++r) {
842: DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
843: if (!user->simplex) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
844: SetupSection(rdm, user);
845: ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
846: p = PetscLog2Real(errorOld/error);
847: PetscPrintf(PetscObjectComm((PetscObject) dm), "Function convergence rate at refinement %D: %.2g\n", r, (double)p);
848: p = PetscLog2Real(errorDerOld/errorDer);
849: PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at refinement %D: %.2g\n", r, (double)p);
850: DMDestroy(&odm);
851: odm = rdm;
852: errorOld = error;
853: errorDerOld = errorDer;
854: }
855: } else {
856: /* ComputeLongestEdge(dm, &lenOld); */
857: DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd);
858: lenOld = cEnd - cStart;
859: for (c = 0; c < Nr; ++c) {
860: DMCoarsen(odm, PetscObjectComm((PetscObject) dm), &cdm);
861: if (!user->simplex) {DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
862: SetupSection(cdm, user);
863: ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
864: /* ComputeLongestEdge(cdm, &len); */
865: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
866: len = cEnd - cStart;
867: rel = error/errorOld;
868: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
869: PetscPrintf(PetscObjectComm((PetscObject) dm), "Function convergence rate at coarsening %D: %.2g\n", c, (double)p);
870: rel = errorDer/errorDerOld;
871: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
872: PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at coarsening %D: %.2g\n", c, (double)p);
873: DMDestroy(&odm);
874: odm = cdm;
875: errorOld = error;
876: errorDerOld = errorDer;
877: lenOld = len;
878: }
879: }
880: DMDestroy(&odm);
881: return(0);
882: }
886: int main(int argc, char **argv)
887: {
888: DM dm;
889: AppCtx user; /* user-defined work context */
892: PetscInitialize(&argc, &argv, NULL, help);
893: ProcessOptions(PETSC_COMM_WORLD, &user);
894: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
895: PetscFECreateDefault(dm, user.dim, user.numComponents, user.simplex, NULL, user.qorder, &user.fe);
896: SetupSection(dm, &user);
897: if (user.testFEjacobian) {
898: TestFEJacobian(dm, &user);
899: }
900: if (user.testFVgrad) {
901: TestFVGrad(dm, &user);
902: }
903: if (user.testInjector) {
904: TestInjector(dm, &user);
905: }
906: CheckFunctions(dm, user.porder, &user);
907: if (user.dim == 2 && user.simplex == PETSC_TRUE && user.tree == PETSC_FALSE) {
908: CheckInterpolation(dm, PETSC_FALSE, user.porder, &user);
909: CheckInterpolation(dm, PETSC_TRUE, user.porder, &user);
910: }
911: CheckConvergence(dm, 3, &user);
912: PetscFEDestroy(&user.fe);
913: DMDestroy(&dm);
914: PetscFinalize();
915: return 0;
916: }