Actual source code: plexfvm.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petsc/private/petscfvimpl.h>
7: static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter lim, PetscInt dim, PetscInt dof, PetscInt cell, PetscInt field, PetscInt face, PetscInt fStart, PetscInt fEnd,
8: PetscReal *cellPhi, const PetscScalar *x, const PetscScalar *cellgeom, const PetscFVCellGeom *cg, const PetscScalar *cx, const PetscScalar *cgrad)
9: {
10: const PetscInt *children;
11: PetscInt numChildren;
12: PetscErrorCode ierr;
15: DMPlexGetTreeChildren(dm,face,&numChildren,&children);
16: if (numChildren) {
17: PetscInt c;
19: for (c = 0; c < numChildren; c++) {
20: PetscInt childFace = children[c];
22: if (childFace >= fStart && childFace < fEnd) {
23: DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,dof,cell,field,childFace,fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
24: }
25: }
26: } else {
27: PetscScalar *ncx;
28: PetscFVCellGeom *ncg;
29: const PetscInt *fcells;
30: PetscInt ncell, d;
31: PetscReal v[3];
33: DMPlexGetSupport(dm, face, &fcells);
34: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
35: if (field >= 0) {
36: DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx);
37: } else {
38: DMPlexPointLocalRead(dm, ncell, x, &ncx);
39: }
40: DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
41: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
42: for (d = 0; d < dof; ++d) {
43: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
44: PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v);
45: PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom;
47: PetscLimiterLimit(lim, flim, &phi);
48: cellPhi[d] = PetscMin(cellPhi[d], phi);
49: }
50: }
51: return(0);
52: }
54: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscFV fvm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
55: {
56: DM dmFace, dmCell, dmGrad;
57: DMLabel ghostLabel;
58: PetscDS prob;
59: PetscLimiter lim;
60: const PetscScalar *facegeom, *cellgeom, *x;
61: PetscScalar *gr;
62: PetscReal *cellPhi;
63: PetscInt dim, face, cell, field, dof, cStart, cEnd, cEndInterior, nFields;
64: PetscErrorCode ierr;
67: DMGetDimension(dm, &dim);
68: DMGetDS(dm, &prob);
69: PetscDSGetNumFields(prob, &nFields);
70: PetscDSGetFieldIndex(prob, (PetscObject) fvm, &field);
71: PetscDSGetFieldSize(prob, field, &dof);
72: DMGetLabel(dm, "ghost", &ghostLabel);
73: PetscFVGetLimiter(fvm, &lim);
74: VecGetDM(faceGeometry, &dmFace);
75: VecGetArrayRead(faceGeometry, &facegeom);
76: VecGetDM(cellGeometry, &dmCell);
77: VecGetArrayRead(cellGeometry, &cellgeom);
78: VecGetArrayRead(locX, &x);
79: VecGetDM(grad, &dmGrad);
80: VecZeroEntries(grad);
81: VecGetArray(grad, &gr);
82: /* Reconstruct gradients */
83: for (face = fStart; face < fEnd; ++face) {
84: const PetscInt *cells;
85: PetscFVFaceGeom *fg;
86: PetscScalar *cx[2];
87: PetscScalar *cgrad[2];
88: PetscBool boundary;
89: PetscInt ghost, c, pd, d, numChildren, numCells;
91: DMLabelGetValue(ghostLabel, face, &ghost);
92: DMIsBoundaryPoint(dm, face, &boundary);
93: DMPlexGetTreeChildren(dm, face, &numChildren, NULL);
94: if (ghost >= 0 || boundary || numChildren) continue;
95: DMPlexGetSupportSize(dm, face, &numCells);
96: if (numCells != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "facet %d has %d support points: expected 2",face,numCells);
97: DMPlexGetSupport(dm, face, &cells);
98: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
99: for (c = 0; c < 2; ++c) {
100: if (nFields > 1) {
101: DMPlexPointLocalFieldRead(dm, cells[c], field, x, &cx[c]);
102: } else {
103: DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
104: }
105: DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
106: }
107: for (pd = 0; pd < dof; ++pd) {
108: PetscScalar delta = cx[1][pd] - cx[0][pd];
110: for (d = 0; d < dim; ++d) {
111: if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta;
112: if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta;
113: }
114: }
115: }
116: /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
117: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
118: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
119: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
120: DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi);
121: for (cell = (dmGrad && lim) ? cStart : cEnd; cell < cEndInterior; ++cell) {
122: const PetscInt *faces;
123: PetscScalar *cx;
124: PetscFVCellGeom *cg;
125: PetscScalar *cgrad;
126: PetscInt coneSize, f, pd, d;
128: DMPlexGetConeSize(dm, cell, &coneSize);
129: DMPlexGetCone(dm, cell, &faces);
130: if (nFields > 1) {
131: DMPlexPointLocalFieldRead(dm, cell, field, x, &cx);
132: }
133: else {
134: DMPlexPointLocalRead(dm, cell, x, &cx);
135: }
136: DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
137: DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
138: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
139: /* Limiter will be minimum value over all neighbors */
140: for (d = 0; d < dof; ++d) cellPhi[d] = PETSC_MAX_REAL;
141: for (f = 0; f < coneSize; ++f) {
142: DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,dof,cell,nFields > 1 ? field : -1,faces[f],fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
143: }
144: /* Apply limiter to gradient */
145: for (pd = 0; pd < dof; ++pd)
146: /* Scalar limiter applied to each component separately */
147: for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
148: }
149: DMRestoreWorkArray(dm, dof, MPIU_REAL, &cellPhi);
150: VecRestoreArrayRead(faceGeometry, &facegeom);
151: VecRestoreArrayRead(cellGeometry, &cellgeom);
152: VecRestoreArrayRead(locX, &x);
153: VecRestoreArray(grad, &gr);
154: return(0);
155: }
157: /*@
158: DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.
160: Input Parameters:
161: + dm - the mesh
162: - locX - the local representation of the vector
164: Output Parameter:
165: . grad - the global representation of the gradient
167: Level: developer
169: .seealso: DMPlexSNESGetGradientDM()
170: @*/
171: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
172: {
173: PetscDS prob;
174: PetscInt Nf, f, fStart, fEnd;
175: PetscBool useFVM = PETSC_FALSE;
176: PetscFV fvm = NULL;
177: Vec faceGeometryFVM, cellGeometryFVM;
178: PetscFVCellGeom *cgeomFVM = NULL;
179: PetscFVFaceGeom *fgeomFVM = NULL;
180: DM dmGrad = NULL;
181: PetscErrorCode ierr;
184: DMGetDS(dm, &prob);
185: PetscDSGetNumFields(prob, &Nf);
186: for (f = 0; f < Nf; ++f) {
187: PetscObject obj;
188: PetscClassId id;
190: PetscDSGetDiscretization(prob, f, &obj);
191: PetscObjectGetClassId(obj, &id);
192: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
193: }
194: if (!useFVM) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm does not have a finite volume discretization");
195: DMPlexGetDataFVM(dm, fvm, &cellGeometryFVM, &faceGeometryFVM, &dmGrad);
196: if (!dmGrad) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm's finite volume discretization does not reconstruct gradients");
197: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
198: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
199: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
200: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
201: return(0);
202: }