Actual source code: plexfvm.c

  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, PetscReal *cellPhi, const PetscScalar *x, const PetscScalar *cellgeom, const PetscFVCellGeom *cg, const PetscScalar *cx, const PetscScalar *cgrad)
  8: {
  9:   const PetscInt *children;
 10:   PetscInt        numChildren;

 12:   DMPlexGetTreeChildren(dm, face, &numChildren, &children);
 13:   if (numChildren) {
 14:     PetscInt c;

 16:     for (c = 0; c < numChildren; c++) {
 17:       PetscInt childFace = children[c];

 19:       if (childFace >= fStart && childFace < fEnd) DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, field, childFace, fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad);
 20:     }
 21:   } else {
 22:     PetscScalar     *ncx;
 23:     PetscFVCellGeom *ncg;
 24:     const PetscInt  *fcells;
 25:     PetscInt         ncell, d;
 26:     PetscReal        v[3];

 28:     DMPlexGetSupport(dm, face, &fcells);
 29:     ncell = cell == fcells[0] ? fcells[1] : fcells[0];
 30:     if (field >= 0) {
 31:       DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx);
 32:     } else {
 33:       DMPlexPointLocalRead(dm, ncell, x, &ncx);
 34:     }
 35:     DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
 36:     DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
 37:     for (d = 0; d < dof; ++d) {
 38:       /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
 39:       PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v);
 40:       PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom;

 42:       PetscLimiterLimit(lim, flim, &phi);
 43:       cellPhi[d] = PetscMin(cellPhi[d], phi);
 44:     }
 45:   }
 46:   return 0;
 47: }

 49: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscFV fvm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
 50: {
 51:   DM                 dmFace, dmCell, dmGrad;
 52:   DMLabel            ghostLabel;
 53:   PetscDS            prob;
 54:   PetscLimiter       lim;
 55:   const PetscScalar *facegeom, *cellgeom, *x;
 56:   PetscScalar       *gr;
 57:   PetscReal         *cellPhi;
 58:   PetscInt           dim, face, cell, field, dof, cStart, cEnd, nFields;

 60:   DMGetDimension(dm, &dim);
 61:   DMGetDS(dm, &prob);
 62:   PetscDSGetNumFields(prob, &nFields);
 63:   PetscDSGetFieldIndex(prob, (PetscObject)fvm, &field);
 64:   PetscDSGetFieldSize(prob, field, &dof);
 65:   DMGetLabel(dm, "ghost", &ghostLabel);
 66:   PetscFVGetLimiter(fvm, &lim);
 67:   VecGetDM(faceGeometry, &dmFace);
 68:   VecGetArrayRead(faceGeometry, &facegeom);
 69:   VecGetDM(cellGeometry, &dmCell);
 70:   VecGetArrayRead(cellGeometry, &cellgeom);
 71:   VecGetArrayRead(locX, &x);
 72:   VecGetDM(grad, &dmGrad);
 73:   VecZeroEntries(grad);
 74:   VecGetArray(grad, &gr);
 75:   /* Reconstruct gradients */
 76:   for (face = fStart; face < fEnd; ++face) {
 77:     const PetscInt  *cells;
 78:     PetscFVFaceGeom *fg;
 79:     PetscScalar     *cx[2];
 80:     PetscScalar     *cgrad[2];
 81:     PetscBool        boundary;
 82:     PetscInt         ghost, c, pd, d, numChildren, numCells;

 84:     DMLabelGetValue(ghostLabel, face, &ghost);
 85:     DMIsBoundaryPoint(dm, face, &boundary);
 86:     DMPlexGetTreeChildren(dm, face, &numChildren, NULL);
 87:     if (ghost >= 0 || boundary || numChildren) continue;
 88:     DMPlexGetSupportSize(dm, face, &numCells);
 90:     DMPlexGetSupport(dm, face, &cells);
 91:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
 92:     for (c = 0; c < 2; ++c) {
 93:       if (nFields > 1) {
 94:         DMPlexPointLocalFieldRead(dm, cells[c], field, x, &cx[c]);
 95:       } else {
 96:         DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
 97:       }
 98:       DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
 99:     }
100:     for (pd = 0; pd < dof; ++pd) {
101:       PetscScalar delta = cx[1][pd] - cx[0][pd];

103:       for (d = 0; d < dim; ++d) {
104:         if (cgrad[0]) cgrad[0][pd * dim + d] += fg->grad[0][d] * delta;
105:         if (cgrad[1]) cgrad[1][pd * dim + d] -= fg->grad[1][d] * delta;
106:       }
107:     }
108:   }
109:   /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
110:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
111:   DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi);
112:   for (cell = (dmGrad && lim) ? cStart : cEnd; cell < cEnd; ++cell) {
113:     const PetscInt  *faces;
114:     PetscScalar     *cx;
115:     PetscFVCellGeom *cg;
116:     PetscScalar     *cgrad;
117:     PetscInt         coneSize, f, pd, d;

119:     DMPlexGetConeSize(dm, cell, &coneSize);
120:     DMPlexGetCone(dm, cell, &faces);
121:     if (nFields > 1) {
122:       DMPlexPointLocalFieldRead(dm, cell, field, x, &cx);
123:     } else {
124:       DMPlexPointLocalRead(dm, cell, x, &cx);
125:     }
126:     DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
127:     DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
128:     if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
129:     /* Limiter will be minimum value over all neighbors */
130:     for (d = 0; d < dof; ++d) cellPhi[d] = PETSC_MAX_REAL;
131:     for (f = 0; f < coneSize; ++f) DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, nFields > 1 ? field : -1, faces[f], fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad);
132:     /* Apply limiter to gradient */
133:     for (pd = 0; pd < dof; ++pd) /* Scalar limiter applied to each component separately */
134:       for (d = 0; d < dim; ++d) cgrad[pd * dim + d] *= cellPhi[pd];
135:   }
136:   DMRestoreWorkArray(dm, dof, MPIU_REAL, &cellPhi);
137:   VecRestoreArrayRead(faceGeometry, &facegeom);
138:   VecRestoreArrayRead(cellGeometry, &cellgeom);
139:   VecRestoreArrayRead(locX, &x);
140:   VecRestoreArray(grad, &gr);
141:   return 0;
142: }

144: /*@
145:   DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.

147:   Input Parameters:
148: + dm - the mesh
149: - locX - the local representation of the vector

151:   Output Parameter:
152: . grad - the global representation of the gradient

154:   Level: developer

156: .seealso: [](chapter_unstructured), `DM`, `Vec`, `DMPlexGetGradientDM()`
157: @*/
158: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
159: {
160:   PetscDS          prob;
161:   PetscInt         Nf, f, fStart, fEnd;
162:   PetscBool        useFVM = PETSC_FALSE;
163:   PetscFV          fvm    = NULL;
164:   Vec              faceGeometryFVM, cellGeometryFVM;
165:   PetscFVCellGeom *cgeomFVM = NULL;
166:   PetscFVFaceGeom *fgeomFVM = NULL;
167:   DM               dmGrad   = NULL;

169:   DMGetDS(dm, &prob);
170:   PetscDSGetNumFields(prob, &Nf);
171:   for (f = 0; f < Nf; ++f) {
172:     PetscObject  obj;
173:     PetscClassId id;

175:     PetscDSGetDiscretization(prob, f, &obj);
176:     PetscObjectGetClassId(obj, &id);
177:     if (id == PETSCFV_CLASSID) {
178:       useFVM = PETSC_TRUE;
179:       fvm    = (PetscFV)obj;
180:     }
181:   }
183:   DMPlexGetDataFVM(dm, fvm, &cellGeometryFVM, &faceGeometryFVM, &dmGrad);
185:   VecGetArrayRead(faceGeometryFVM, (const PetscScalar **)&fgeomFVM);
186:   VecGetArrayRead(cellGeometryFVM, (const PetscScalar **)&cgeomFVM);
187:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
188:   DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
189:   return 0;
190: }