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