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