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