Actual source code: dmplexts.c
petsc-3.5.4 2015-05-23
1: #include <petscdmplex.h> /*I "petscdmplex.h" I*/
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
3: #include <petscfv.h>
5: PETSC_STATIC_INLINE void WaxpyD(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];}
6: PETSC_STATIC_INLINE PetscReal DotD(PetscInt dim, const PetscScalar *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d])*y[d]; return sum;}
7: PETSC_STATIC_INLINE PetscReal DotRealD(PetscInt dim, const PetscReal *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;}
8: PETSC_STATIC_INLINE PetscReal NormD(PetscInt dim, const PetscReal *x) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*x[d]; return PetscSqrtReal(sum);}
10: typedef struct {
11: PetscBool setupGeom; /* Flag for geometry setup */
12: PetscBool setupGrad; /* Flag for gradient calculation setup */
13: Vec facegeom; /* FaceGeom struct for each face */
14: Vec cellgeom; /* CellGeom struct for each cell */
15: DM dmGrad; /* Layout for the gradient data */
16: PetscReal minradius; /* Minimum distance from centroid to face */
17: void (*riemann)(const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[], void *);
18: void *rhsfunctionlocalctx;
19: } DMTS_Plex;
23: static PetscErrorCode DMTSDestroy_Plex(DMTS dmts)
24: {
25: DMTS_Plex *dmplexts = (DMTS_Plex *) dmts->data;
29: DMDestroy(&dmplexts->dmGrad);
30: VecDestroy(&dmplexts->facegeom);
31: VecDestroy(&dmplexts->cellgeom);
32: PetscFree(dmts->data);
33: return(0);
34: }
38: static PetscErrorCode DMTSDuplicate_Plex(DMTS olddmts, DMTS dmts)
39: {
43: PetscNewLog(dmts, (DMTS_Plex **) &dmts->data);
44: if (olddmts->data) {PetscMemcpy(dmts->data, olddmts->data, sizeof(DMTS_Plex));}
45: return(0);
46: }
50: static PetscErrorCode DMPlexTSGetContext(DM dm, DMTS dmts, DMTS_Plex **dmplexts)
51: {
55: *dmplexts = NULL;
56: if (!dmts->data) {
57: PetscNewLog(dm, (DMTS_Plex **) &dmts->data);
58: dmts->ops->destroy = DMTSDestroy_Plex;
59: dmts->ops->duplicate = DMTSDuplicate_Plex;
60: }
61: *dmplexts = (DMTS_Plex *) dmts->data;
62: return(0);
63: }
67: /*@
68: DMPlexTSGetGeometry - Return precomputed geometric data
70: Input Parameter:
71: . dm - The DM
73: Output Parameters:
74: + facegeom - The values precomputed from face geometry
75: . cellgeom - The values precomputed from cell geometry
76: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
78: Level: developer
80: .seealso: DMPlexTSSetRHSFunctionLocal()
81: @*/
82: PetscErrorCode DMPlexTSGetGeometry(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
83: {
84: DMTS dmts;
85: DMTS_Plex *dmplexts;
89: DMGetDMTS(dm, &dmts);
90: DMPlexTSGetContext(dm, dmts, &dmplexts);
91: if (facegeom) *facegeom = dmplexts->facegeom;
92: if (cellgeom) *cellgeom = dmplexts->cellgeom;
93: if (minRadius) *minRadius = dmplexts->minradius;
94: return(0);
95: }
99: static PetscErrorCode DMPlexTSSetupGeometry(DM dm, PetscFV fvm, DMTS_Plex *dmplexts)
100: {
101: DM dmFace, dmCell;
102: DMLabel ghostLabel;
103: PetscSection sectionFace, sectionCell;
104: PetscSection coordSection;
105: Vec coordinates;
106: PetscReal minradius;
107: PetscScalar *fgeom, *cgeom;
108: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
112: if (dmplexts->setupGeom) return(0);
113: DMPlexGetDimension(dm, &dim);
114: DMGetCoordinateSection(dm, &coordSection);
115: DMGetCoordinatesLocal(dm, &coordinates);
116: /* Make cell centroids and volumes */
117: DMClone(dm, &dmCell);
118: DMSetCoordinateSection(dmCell, coordSection);
119: DMSetCoordinatesLocal(dmCell, coordinates);
120: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);
121: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
122: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
123: PetscSectionSetChart(sectionCell, cStart, cEnd);
124: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));}
125: PetscSectionSetUp(sectionCell);
126: DMSetDefaultSection(dmCell, sectionCell);
127: PetscSectionDestroy(§ionCell);
128: DMCreateLocalVector(dmCell, &dmplexts->cellgeom);
129: VecGetArray(dmplexts->cellgeom, &cgeom);
130: for (c = cStart; c < cEndInterior; ++c) {
131: CellGeom *cg;
133: DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
134: PetscMemzero(cg, sizeof(*cg));
135: DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
136: }
137: /* Compute face normals and minimum cell radius */
138: DMClone(dm, &dmFace);
139: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);
140: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
141: PetscSectionSetChart(sectionFace, fStart, fEnd);
142: for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));}
143: PetscSectionSetUp(sectionFace);
144: DMSetDefaultSection(dmFace, sectionFace);
145: PetscSectionDestroy(§ionFace);
146: DMCreateLocalVector(dmFace, &dmplexts->facegeom);
147: VecGetArray(dmplexts->facegeom, &fgeom);
148: DMPlexGetLabel(dm, "ghost", &ghostLabel);
149: minradius = PETSC_MAX_REAL;
150: for (f = fStart; f < fEnd; ++f) {
151: FaceGeom *fg;
152: PetscReal area;
153: PetscInt ghost, d;
155: DMLabelGetValue(ghostLabel, f, &ghost);
156: if (ghost >= 0) continue;
157: DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
158: DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
159: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
160: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
161: {
162: CellGeom *cL, *cR;
163: const PetscInt *cells;
164: PetscReal *lcentroid, *rcentroid;
165: PetscReal v[3];
167: DMPlexGetSupport(dm, f, &cells);
168: DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
169: DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
170: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
171: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
172: WaxpyD(dim, -1, lcentroid, rcentroid, v);
173: if (DotRealD(dim, fg->normal, v) < 0) {
174: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
175: }
176: if (DotRealD(dim, fg->normal, v) <= 0) {
177: if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) v[0], (double) v[1]);
178: if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) fg->normal[2], (double) v[0], (double) v[1], (double) v[2]);
179: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
180: }
181: if (cells[0] < cEndInterior) {
182: WaxpyD(dim, -1, fg->centroid, cL->centroid, v);
183: minradius = PetscMin(minradius, NormD(dim, v));
184: }
185: if (cells[1] < cEndInterior) {
186: WaxpyD(dim, -1, fg->centroid, cR->centroid, v);
187: minradius = PetscMin(minradius, NormD(dim, v));
188: }
189: }
190: }
191: MPI_Allreduce(&minradius, &dmplexts->minradius, 1, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm));
192: /* Compute centroids of ghost cells */
193: for (c = cEndInterior; c < cEnd; ++c) {
194: FaceGeom *fg;
195: const PetscInt *cone, *support;
196: PetscInt coneSize, supportSize, s;
198: DMPlexGetConeSize(dmCell, c, &coneSize);
199: if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
200: DMPlexGetCone(dmCell, c, &cone);
201: DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
202: if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize);
203: DMPlexGetSupport(dmCell, cone[0], &support);
204: DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
205: for (s = 0; s < 2; ++s) {
206: /* Reflect ghost centroid across plane of face */
207: if (support[s] == c) {
208: const CellGeom *ci;
209: CellGeom *cg;
210: PetscReal c2f[3], a;
212: DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
213: WaxpyD(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
214: a = DotRealD(dim, c2f, fg->normal)/DotRealD(dim, fg->normal, fg->normal);
215: DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
216: WaxpyD(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
217: cg->volume = ci->volume;
218: }
219: }
220: }
221: VecRestoreArray(dmplexts->facegeom, &fgeom);
222: VecRestoreArray(dmplexts->cellgeom, &cgeom);
223: DMDestroy(&dmCell);
224: DMDestroy(&dmFace);
225: dmplexts->setupGeom = PETSC_TRUE;
226: return(0);
227: }
231: static PetscErrorCode BuildGradientReconstruction(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
232: {
233: DMLabel ghostLabel;
234: PetscScalar *dx, *grad, **gref;
235: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
239: DMPlexGetDimension(dm, &dim);
240: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
241: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
242: DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
243: PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
244: DMPlexGetLabel(dm, "ghost", &ghostLabel);
245: PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
246: for (c = cStart; c < cEndInterior; c++) {
247: const PetscInt *faces;
248: PetscInt numFaces, usedFaces, f, d;
249: const CellGeom *cg;
250: PetscBool boundary;
251: PetscInt ghost;
253: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
254: DMPlexGetConeSize(dm, c, &numFaces);
255: DMPlexGetCone(dm, c, &faces);
256: if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
257: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
258: const CellGeom *cg1;
259: FaceGeom *fg;
260: const PetscInt *fcells;
261: PetscInt ncell, side;
263: DMLabelGetValue(ghostLabel, faces[f], &ghost);
264: DMPlexIsBoundaryPoint(dm, faces[f], &boundary);
265: if ((ghost >= 0) || boundary) continue;
266: DMPlexGetSupport(dm, faces[f], &fcells);
267: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
268: ncell = fcells[!side]; /* the neighbor */
269: DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
270: DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
271: for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
272: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
273: }
274: if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
275: PetscFVComputeGradient(fvm, usedFaces, dx, grad);
276: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
277: DMLabelGetValue(ghostLabel, faces[f], &ghost);
278: DMPlexIsBoundaryPoint(dm, faces[f], &boundary);
279: if ((ghost >= 0) || boundary) continue;
280: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
281: ++usedFaces;
282: }
283: }
284: PetscFree3(dx, grad, gref);
285: return(0);
286: }
290: static PetscErrorCode DMPlexTSSetupGradient(DM dm, PetscFV fvm, DMTS_Plex *dmplexts)
291: {
292: DM dmFace, dmCell;
293: PetscScalar *fgeom, *cgeom;
294: PetscSection sectionGrad;
295: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
299: if (dmplexts->setupGrad) return(0);
300: DMPlexGetDimension(dm, &dim);
301: PetscFVGetNumComponents(fvm, &pdim);
302: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
303: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
304: /* Construct the interpolant corresponding to each face from the leat-square solution over the cell neighborhood */
305: VecGetDM(dmplexts->facegeom, &dmFace);
306: VecGetDM(dmplexts->cellgeom, &dmCell);
307: VecGetArray(dmplexts->facegeom, &fgeom);
308: VecGetArray(dmplexts->cellgeom, &cgeom);
309: BuildGradientReconstruction(dm, fvm, dmFace, fgeom, dmCell, cgeom);
310: VecRestoreArray(dmplexts->facegeom, &fgeom);
311: VecRestoreArray(dmplexts->cellgeom, &cgeom);
312: /* Create storage for gradients */
313: DMClone(dm, &dmplexts->dmGrad);
314: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);
315: PetscSectionSetChart(sectionGrad, cStart, cEnd);
316: for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
317: PetscSectionSetUp(sectionGrad);
318: DMSetDefaultSection(dmplexts->dmGrad, sectionGrad);
319: PetscSectionDestroy(§ionGrad);
320: dmplexts->setupGrad = PETSC_TRUE;
321: return(0);
322: }
326: static PetscErrorCode DMPlexInsertBoundaryValuesFVM_Static(DM dm, PetscFV fvm, PetscReal time, Vec locX, Vec Grad, DMTS_Plex *dmplexts)
327: {
328: Vec faceGeometry = dmplexts->facegeom;
329: Vec cellGeometry = dmplexts->cellgeom;
330: DM dmFace, dmCell;
331: const PetscScalar *facegeom, *cellgeom, *grad;
332: PetscScalar *x, *fx;
333: PetscInt numBd, b, dim, pdim, fStart, fEnd;
334: PetscErrorCode ierr;
337: DMPlexGetDimension(dm, &dim);
338: PetscFVGetNumComponents(fvm, &pdim);
339: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
340: DMPlexGetNumBoundary(dm, &numBd);
341: if (Grad) {
342: VecGetDM(cellGeometry, &dmCell);
343: VecGetArrayRead(cellGeometry, &cellgeom);
344: DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
345: VecGetArrayRead(Grad, &grad);
346: }
347: VecGetDM(faceGeometry, &dmFace);
348: VecGetArrayRead(faceGeometry, &facegeom);
349: VecGetArray(locX, &x);
350: for (b = 0; b < numBd; ++b) {
351: PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*);
352: DMLabel label;
353: const char *labelname;
354: const PetscInt *ids;
355: PetscInt numids, i;
356: void *ctx;
358: DMPlexGetBoundary(dm, b, NULL, NULL, &labelname, NULL, (void (**)()) &func, &numids, &ids, &ctx);
359: DMPlexGetLabel(dm, labelname, &label);
360: for (i = 0; i < numids; ++i) {
361: IS faceIS;
362: const PetscInt *faces;
363: PetscInt numFaces, f;
365: DMLabelGetStratumIS(label, ids[i], &faceIS);
366: if (!faceIS) continue; /* No points with that id on this process */
367: ISGetLocalSize(faceIS, &numFaces);
368: ISGetIndices(faceIS, &faces);
369: for (f = 0; f < numFaces; ++f) {
370: const PetscInt face = faces[f], *cells;
371: const FaceGeom *fg;
373: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
374: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
375: DMPlexGetSupport(dm, face, &cells);
376: if (Grad) {
377: const CellGeom *cg;
378: const PetscScalar *cx, *cgrad;
379: PetscScalar *xG;
380: PetscReal dx[3];
381: PetscInt d;
383: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
384: DMPlexPointLocalRead(dm, cells[0], x, &cx);
385: DMPlexPointLocalRead(dmplexts->dmGrad, cells[0], grad, &cgrad);
386: DMPlexPointLocalRef(dm, cells[1], x, &xG);
387: WaxpyD(dim, -1, cg->centroid, fg->centroid, dx);
388: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DotD(dim, &cgrad[d*dim], dx);
389: (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
390: } else {
391: const PetscScalar *xI;
392: PetscScalar *xG;
394: DMPlexPointLocalRead(dm, cells[0], x, &xI);
395: DMPlexPointLocalRef(dm, cells[1], x, &xG);
396: (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
397: }
398: }
399: ISRestoreIndices(faceIS, &faces);
400: ISDestroy(&faceIS);
401: }
402: }
403: VecRestoreArrayRead(faceGeometry, &facegeom);
404: VecRestoreArray(locX, &x);
405: if (Grad) {
406: DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
407: VecRestoreArrayRead(Grad, &grad);
408: }
409: return(0);
410: }
414: static PetscErrorCode TSComputeRHSFunction_DMPlex(TS ts, PetscReal time, Vec X, Vec F, void *ctx)
415: {
416: DM dm;
417: DMTS_Plex *dmplexts = (DMTS_Plex *) ctx;
418: void (*riemann)(const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[], void *) = dmplexts->riemann;
419: PetscFV fvm;
420: PetscLimiter lim;
421: Vec faceGeometry = dmplexts->facegeom;
422: Vec cellGeometry = dmplexts->cellgeom;
423: Vec Grad = NULL, locGrad, locX;
424: DM dmFace, dmCell;
425: DMLabel ghostLabel;
426: PetscCellGeometry fgeom, cgeom;
427: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
428: PetscScalar *grad, *f, *uL, *uR, *fluxL, *fluxR;
429: PetscReal *centroid, *normal, *vol, *cellPhi;
430: PetscBool computeGradients;
431: PetscInt Nf, dim, pdim, fStart, fEnd, numFaces = 0, face, iface, cell, cStart, cEnd, cEndInterior;
432: PetscErrorCode ierr;
438: TSGetDM(ts, &dm);
439: DMGetLocalVector(dm, &locX);
440: VecZeroEntries(locX);
441: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
442: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
443: VecZeroEntries(F);
444: DMPlexGetDimension(dm, &dim);
445: DMGetNumFields(dm, &Nf);
446: DMGetField(dm, 0, (PetscObject *) &fvm);
447: PetscFVGetLimiter(fvm, &lim);
448: PetscFVGetNumComponents(fvm, &pdim);
449: PetscFVGetComputeGradients(fvm, &computeGradients);
450: if (computeGradients) {
451: DMGetGlobalVector(dmplexts->dmGrad, &Grad);
452: VecZeroEntries(Grad);
453: VecGetArray(Grad, &grad);
454: }
455: DMPlexGetLabel(dm, "ghost", &ghostLabel);
456: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
457: VecGetDM(faceGeometry, &dmFace);
458: VecGetDM(cellGeometry, &dmCell);
459: VecGetArrayRead(faceGeometry, &facegeom);
460: VecGetArrayRead(cellGeometry, &cellgeom);
461: VecGetArrayRead(locX, &x);
462: /* Count faces and reconstruct gradients */
463: for (face = fStart; face < fEnd; ++face) {
464: const PetscInt *cells;
465: const FaceGeom *fg;
466: const PetscScalar *cx[2];
467: PetscScalar *cgrad[2];
468: PetscBool boundary;
469: PetscInt ghost, c, pd, d;
471: DMLabelGetValue(ghostLabel, face, &ghost);
472: if (ghost >= 0) continue;
473: ++numFaces;
474: if (!computeGradients) continue;
475: DMPlexIsBoundaryPoint(dm, face, &boundary);
476: if (boundary) continue;
477: DMPlexGetSupport(dm, face, &cells);
478: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
479: for (c = 0; c < 2; ++c) {
480: DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
481: DMPlexPointGlobalRef(dmplexts->dmGrad, cells[c], grad, &cgrad[c]);
482: }
483: for (pd = 0; pd < pdim; ++pd) {
484: PetscScalar delta = cx[1][pd] - cx[0][pd];
486: for (d = 0; d < dim; ++d) {
487: if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta;
488: if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta;
489: }
490: }
491: }
492: /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
493: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
494: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
495: DMGetWorkArray(dm, pdim, PETSC_REAL, &cellPhi);
496: for (cell = computeGradients && lim ? cStart : cEnd; cell < cEndInterior; ++cell) {
497: const PetscInt *faces;
498: const PetscScalar *cx;
499: const CellGeom *cg;
500: PetscScalar *cgrad;
501: PetscInt coneSize, f, pd, d;
503: DMPlexGetConeSize(dm, cell, &coneSize);
504: DMPlexGetCone(dm, cell, &faces);
505: DMPlexPointLocalRead(dm, cell, x, &cx);
506: DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
507: DMPlexPointGlobalRef(dmplexts->dmGrad, cell, grad, &cgrad);
508: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
509: /* Limiter will be minimum value over all neighbors */
510: for (d = 0; d < pdim; ++d) cellPhi[d] = PETSC_MAX_REAL;
511: for (f = 0; f < coneSize; ++f) {
512: const PetscScalar *ncx;
513: const CellGeom *ncg;
514: const PetscInt *fcells;
515: PetscInt face = faces[f], ncell, ghost;
516: PetscReal v[3];
517: PetscBool boundary;
519: DMLabelGetValue(ghostLabel, face, &ghost);
520: DMPlexIsBoundaryPoint(dm, face, &boundary);
521: if ((ghost >= 0) || boundary) continue;
522: DMPlexGetSupport(dm, face, &fcells);
523: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
524: DMPlexPointLocalRead(dm, ncell, x, &ncx);
525: DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
526: WaxpyD(dim, -1, cg->centroid, ncg->centroid, v);
527: for (d = 0; d < pdim; ++d) {
528: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
529: PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / DotD(dim, &cgrad[d*dim], v);
531: PetscLimiterLimit(lim, flim, &phi);
532: cellPhi[d] = PetscMin(cellPhi[d], phi);
533: }
534: }
535: /* Apply limiter to gradient */
536: for (pd = 0; pd < pdim; ++pd)
537: /* Scalar limiter applied to each component separately */
538: for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
539: }
540: DMRestoreWorkArray(dm, pdim, PETSC_REAL, &cellPhi);
541: DMPlexInsertBoundaryValuesFVM_Static(dm, fvm, time, locX, Grad, dmplexts);
542: if (computeGradients) {
543: VecRestoreArray(Grad, &grad);
544: DMGetLocalVector(dmplexts->dmGrad, &locGrad);
545: DMGlobalToLocalBegin(dmplexts->dmGrad, Grad, INSERT_VALUES, locGrad);
546: DMGlobalToLocalEnd(dmplexts->dmGrad, Grad, INSERT_VALUES, locGrad);
547: DMRestoreGlobalVector(dmplexts->dmGrad, &Grad);
548: VecGetArrayRead(locGrad, &lgrad);
549: }
550: PetscMalloc7(numFaces*dim,¢roid,numFaces*dim,&normal,numFaces*2,&vol,numFaces*pdim,&uL,numFaces*pdim,&uR,numFaces*pdim,&fluxL,numFaces*pdim,&fluxR);
551: /* Read out values */
552: for (face = fStart, iface = 0; face < fEnd; ++face) {
553: const PetscInt *cells;
554: const FaceGeom *fg;
555: const CellGeom *cgL, *cgR;
556: const PetscScalar *xL, *xR, *gL, *gR;
557: PetscInt ghost, d;
559: DMLabelGetValue(ghostLabel, face, &ghost);
560: if (ghost >= 0) continue;
561: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
562: DMPlexGetSupport(dm, face, &cells);
563: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
564: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
565: DMPlexPointLocalRead(dm, cells[0], x, &xL);
566: DMPlexPointLocalRead(dm, cells[1], x, &xR);
567: if (computeGradients) {
568: PetscReal dxL[3], dxR[3];
570: DMPlexPointLocalRead(dmplexts->dmGrad, cells[0], lgrad, &gL);
571: DMPlexPointLocalRead(dmplexts->dmGrad, cells[1], lgrad, &gR);
572: WaxpyD(dim, -1, cgL->centroid, fg->centroid, dxL);
573: WaxpyD(dim, -1, cgR->centroid, fg->centroid, dxR);
574: for (d = 0; d < pdim; ++d) {
575: uL[iface*pdim+d] = xL[d] + DotD(dim, &gL[d*dim], dxL);
576: uR[iface*pdim+d] = xR[d] + DotD(dim, &gR[d*dim], dxR);
577: }
578: } else {
579: for (d = 0; d < pdim; ++d) {
580: uL[iface*pdim+d] = xL[d];
581: uR[iface*pdim+d] = xR[d];
582: }
583: }
584: for (d = 0; d < dim; ++d) {
585: centroid[iface*dim+d] = fg->centroid[d];
586: normal[iface*dim+d] = fg->normal[d];
587: }
588: vol[iface*2+0] = cgL->volume;
589: vol[iface*2+1] = cgR->volume;
590: ++iface;
591: }
592: if (computeGradients) {
593: VecRestoreArrayRead(locGrad,&lgrad);
594: DMRestoreLocalVector(dmplexts->dmGrad, &locGrad);
595: }
596: VecRestoreArrayRead(locX, &x);
597: VecRestoreArrayRead(faceGeometry, &facegeom);
598: VecRestoreArrayRead(cellGeometry, &cellgeom);
599: fgeom.v0 = centroid;
600: fgeom.n = normal;
601: cgeom.vol = vol;
602: /* Riemann solve */
603: PetscFVIntegrateRHSFunction(fvm, numFaces, Nf, &fvm, 0, fgeom, cgeom, uL, uR, riemann, fluxL, fluxR, dmplexts->rhsfunctionlocalctx);
604: /* Insert fluxes */
605: VecGetArray(F, &f);
606: for (face = fStart, iface = 0; face < fEnd; ++face) {
607: const PetscInt *cells;
608: PetscScalar *fL, *fR;
609: PetscInt ghost, d;
611: DMLabelGetValue(ghostLabel, face, &ghost);
612: if (ghost >= 0) continue;
613: DMPlexGetSupport(dm, face, &cells);
614: DMPlexPointGlobalRef(dm, cells[0], f, &fL);
615: DMPlexPointGlobalRef(dm, cells[1], f, &fR);
616: for (d = 0; d < pdim; ++d) {
617: if (fL) fL[d] -= fluxL[iface*pdim+d];
618: if (fR) fR[d] += fluxR[iface*pdim+d];
619: }
620: ++iface;
621: }
622: VecRestoreArray(F, &f);
623: PetscFree7(centroid,normal,vol,uL,uR,fluxL,fluxR);
624: DMRestoreLocalVector(dm, &locX);
625: return(0);
626: }
630: /*@C
631: DMPlexTSSetRHSFunctionLocal - set a local residual evaluation function
633: Logically Collective
635: Input Arguments:
636: + dm - DM to associate callback with
637: . riemann - Riemann solver
638: - ctx - optional context for Riemann solve
640: Calling sequence for riemann:
642: $ riemann(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
644: + x - The coordinates at a point on the interface
645: . n - The normal vector to the interface
646: . uL - The state vector to the left of the interface
647: . uR - The state vector to the right of the interface
648: . flux - output array of flux through the interface
649: - ctx - optional user context
651: Level: beginner
653: .seealso: DMTSSetRHSFunctionLocal()
654: @*/
655: PetscErrorCode DMPlexTSSetRHSFunctionLocal(DM dm, void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx), void *ctx)
656: {
657: DMTS dmts;
658: DMTS_Plex *dmplexts;
659: PetscFV fvm;
660: PetscInt Nf;
661: PetscBool computeGradients;
666: DMGetDMTSWrite(dm, &dmts);
667: DMPlexTSGetContext(dm, dmts, &dmplexts);
668: dmplexts->riemann = riemann;
669: dmplexts->rhsfunctionlocalctx = ctx;
670: DMGetNumFields(dm, &Nf);
671: DMGetField(dm, 0, (PetscObject *) &fvm);
672: DMPlexTSSetupGeometry(dm, fvm, dmplexts);
673: PetscFVGetComputeGradients(fvm, &computeGradients);
674: if (computeGradients) {DMPlexTSSetupGradient(dm, fvm, dmplexts);}
675: DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMPlex, dmplexts);
676: return(0);
677: }