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,
      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, 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:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
    118:   DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi);
    119:   for (cell = (dmGrad && lim) ? cStart : cEnd; cell < cEnd; ++cell) {
    120:     const PetscInt        *faces;
    121:     PetscScalar           *cx;
    122:     PetscFVCellGeom       *cg;
    123:     PetscScalar           *cgrad;
    124:     PetscInt               coneSize, f, pd, d;

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

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

    158:   Input Parameters:
    159: + dm - the mesh
    160: - locX - the local representation of the vector

    162:   Output Parameter:
    163: . grad - the global representation of the gradient

    165:   Level: developer

    167: .seealso: DMPlexGetGradientDM()
    168: @*/
    169: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
    170: {
    171:   PetscDS          prob;
    172:   PetscInt         Nf, f, fStart, fEnd;
    173:   PetscBool        useFVM = PETSC_FALSE;
    174:   PetscFV          fvm = NULL;
    175:   Vec              faceGeometryFVM, cellGeometryFVM;
    176:   PetscFVCellGeom  *cgeomFVM   = NULL;
    177:   PetscFVFaceGeom  *fgeomFVM   = NULL;
    178:   DM               dmGrad = NULL;
    179:   PetscErrorCode   ierr;

    182:   DMGetDS(dm, &prob);
    183:   PetscDSGetNumFields(prob, &Nf);
    184:   for (f = 0; f < Nf; ++f) {
    185:     PetscObject  obj;
    186:     PetscClassId id;

    188:     PetscDSGetDiscretization(prob, f, &obj);
    189:     PetscObjectGetClassId(obj, &id);
    190:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
    191:   }
    192:   if (!useFVM) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm does not have a finite volume discretization");
    193:   DMPlexGetDataFVM(dm, fvm, &cellGeometryFVM, &faceGeometryFVM, &dmGrad);
    194:   if (!dmGrad) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm's finite volume discretization does not reconstruct gradients");
    195:   VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
    196:   VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
    197:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
    198:   DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
    199:   return(0);
    200: }