1 #include <petscdmplex.h> /*I "petscdmplex.h" I*/ 2 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscds.h> 4 #include <petscfv.h> 5 6 PETSC_STATIC_INLINE void WaxpyD(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];} 7 PETSC_STATIC_INLINE PetscReal DotD(PetscInt dim, const PetscScalar *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d])*y[d]; return sum;} 8 PETSC_STATIC_INLINE PetscReal DotRealD(PetscInt dim, const PetscReal *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;} 9 PETSC_STATIC_INLINE PetscReal NormD(PetscInt dim, const PetscReal *x) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*x[d]; return PetscSqrtReal(sum);} 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "DMPlexTSGetGeometry" 13 /*@ 14 DMPlexTSGetGeometry - Return precomputed geometric data 15 16 Input Parameter: 17 . dm - The DM 18 19 Output Parameters: 20 + facegeom - The values precomputed from face geometry 21 . cellgeom - The values precomputed from cell geometry 22 - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell 23 24 Level: developer 25 26 .seealso: DMPlexTSSetRHSFunctionLocal() 27 @*/ 28 PetscErrorCode DMPlexTSGetGeometry(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius) 29 { 30 DMTS dmts; 31 PetscErrorCode ierr; 32 33 PetscFunctionBegin; 34 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 35 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 36 if (facegeom) {PetscValidPointer(facegeom, 2); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom", (PetscObject *) facegeom);CHKERRQ(ierr);} 37 if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_cellgeom", (PetscObject *) cellgeom);CHKERRQ(ierr);} 38 if (minRadius) {ierr = DMTSGetMinRadius(dm, minRadius);CHKERRQ(ierr);} 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "DMPlexTSGetGradientDM" 44 /*@ 45 DMPlexTSGetGradientDM - Return gradient data layout 46 47 Input Parameter: 48 . dm - The DM 49 50 Output Parameter: 51 . dmGrad - The layout for gradient values 52 53 Level: developer 54 55 .seealso: DMPlexTSGetGeometry(), DMPlexTSSetRHSFunctionLocal() 56 @*/ 57 PetscErrorCode DMPlexTSGetGradientDM(DM dm, DM *dmGrad) 58 { 59 DMTS dmts; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 64 PetscValidPointer(dmGrad,2); 65 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 66 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad", (PetscObject *) dmGrad);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "DMPlexTSSetupGeometry" 72 static PetscErrorCode DMPlexTSSetupGeometry(DM dm, PetscFV fvm, DMTS dmts) 73 { 74 DM dmFace, dmCell; 75 DMLabel ghostLabel; 76 PetscSection sectionFace, sectionCell; 77 PetscSection coordSection; 78 Vec coordinates, cellgeom, facegeom; 79 PetscScalar *fgeom, *cgeom; 80 PetscReal minradius, gminradius; 81 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 86 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 87 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 88 /* Make cell centroids and volumes */ 89 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 90 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 91 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 92 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 93 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 94 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 95 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 96 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 97 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 98 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 99 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 100 ierr = DMCreateLocalVector(dmCell, &cellgeom);CHKERRQ(ierr); 101 ierr = VecGetArray(cellgeom, &cgeom);CHKERRQ(ierr); 102 for (c = cStart; c < cEndInterior; ++c) { 103 CellGeom *cg; 104 105 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 106 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 107 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 108 } 109 /* Compute face normals and minimum cell radius */ 110 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 111 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 112 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 113 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 114 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 115 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 116 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 117 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 118 ierr = DMCreateLocalVector(dmFace, &facegeom);CHKERRQ(ierr); 119 ierr = VecGetArray(facegeom, &fgeom);CHKERRQ(ierr); 120 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 121 minradius = PETSC_MAX_REAL; 122 for (f = fStart; f < fEnd; ++f) { 123 FaceGeom *fg; 124 PetscReal area; 125 PetscInt ghost, d; 126 127 ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr); 128 if (ghost >= 0) continue; 129 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 130 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 131 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 132 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 133 { 134 CellGeom *cL, *cR; 135 const PetscInt *cells; 136 PetscReal *lcentroid, *rcentroid; 137 PetscReal v[3]; 138 139 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 140 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 141 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 142 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 143 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 144 WaxpyD(dim, -1, lcentroid, rcentroid, v); 145 if (DotRealD(dim, fg->normal, v) < 0) { 146 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 147 } 148 if (DotRealD(dim, fg->normal, v) <= 0) { 149 if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) v[0], (double) v[1]); 150 if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) fg->normal[2], (double) v[0], (double) v[1], (double) v[2]); 151 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 152 } 153 if (cells[0] < cEndInterior) { 154 WaxpyD(dim, -1, fg->centroid, cL->centroid, v); 155 minradius = PetscMin(minradius, NormD(dim, v)); 156 } 157 if (cells[1] < cEndInterior) { 158 WaxpyD(dim, -1, fg->centroid, cR->centroid, v); 159 minradius = PetscMin(minradius, NormD(dim, v)); 160 } 161 } 162 } 163 ierr = MPI_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 164 ierr = DMTSSetMinRadius(dm, gminradius);CHKERRQ(ierr); 165 /* Compute centroids of ghost cells */ 166 for (c = cEndInterior; c < cEnd; ++c) { 167 FaceGeom *fg; 168 const PetscInt *cone, *support; 169 PetscInt coneSize, supportSize, s; 170 171 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 172 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 173 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 174 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 175 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize); 176 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 177 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 178 for (s = 0; s < 2; ++s) { 179 /* Reflect ghost centroid across plane of face */ 180 if (support[s] == c) { 181 const CellGeom *ci; 182 CellGeom *cg; 183 PetscReal c2f[3], a; 184 185 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 186 WaxpyD(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 187 a = DotRealD(dim, c2f, fg->normal)/DotRealD(dim, fg->normal, fg->normal); 188 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 189 WaxpyD(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 190 cg->volume = ci->volume; 191 } 192 } 193 } 194 ierr = VecRestoreArray(facegeom, &fgeom);CHKERRQ(ierr); 195 ierr = VecRestoreArray(cellgeom, &cgeom);CHKERRQ(ierr); 196 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom", (PetscObject) facegeom);CHKERRQ(ierr); 197 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom", (PetscObject) cellgeom);CHKERRQ(ierr); 198 ierr = VecDestroy(&facegeom);CHKERRQ(ierr); 199 ierr = VecDestroy(&cellgeom);CHKERRQ(ierr); 200 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 201 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 205 #undef __FUNCT__ 206 #define __FUNCT__ "BuildGradientReconstruction" 207 static PetscErrorCode BuildGradientReconstruction(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 208 { 209 DMLabel ghostLabel; 210 PetscScalar *dx, *grad, **gref; 211 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 216 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 217 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 218 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 219 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 220 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 221 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 222 for (c = cStart; c < cEndInterior; c++) { 223 const PetscInt *faces; 224 PetscInt numFaces, usedFaces, f, d; 225 const CellGeom *cg; 226 PetscBool boundary; 227 PetscInt ghost; 228 229 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 230 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 231 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 232 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 233 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 234 const CellGeom *cg1; 235 FaceGeom *fg; 236 const PetscInt *fcells; 237 PetscInt ncell, side; 238 239 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 240 ierr = DMPlexIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 241 if ((ghost >= 0) || boundary) continue; 242 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 243 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 244 ncell = fcells[!side]; /* the neighbor */ 245 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 246 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 247 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 248 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 249 } 250 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 251 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 252 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 253 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 254 ierr = DMPlexIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 255 if ((ghost >= 0) || boundary) continue; 256 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 257 ++usedFaces; 258 } 259 } 260 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "DMPlexTSSetupGradient" 266 static PetscErrorCode DMPlexTSSetupGradient(DM dm, PetscFV fvm, DMTS dmts) 267 { 268 DM dmFace, dmCell, dmGrad; 269 Vec facegeom, cellgeom; 270 PetscScalar *fgeom, *cgeom; 271 PetscSection sectionGrad; 272 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 277 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 278 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 279 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 280 /* Construct the interpolant corresponding to each face from the leat-square solution over the cell neighborhood */ 281 ierr = DMPlexTSGetGeometry(dm, &facegeom, &cellgeom, NULL);CHKERRQ(ierr); 282 ierr = VecGetDM(facegeom, &dmFace);CHKERRQ(ierr); 283 ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr); 284 ierr = VecGetArray(facegeom, &fgeom);CHKERRQ(ierr); 285 ierr = VecGetArray(cellgeom, &cgeom);CHKERRQ(ierr); 286 ierr = BuildGradientReconstruction(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 287 ierr = VecRestoreArray(facegeom, &fgeom);CHKERRQ(ierr); 288 ierr = VecRestoreArray(cellgeom, &cgeom);CHKERRQ(ierr); 289 /* Create storage for gradients */ 290 ierr = DMClone(dm, &dmGrad);CHKERRQ(ierr); 291 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 292 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 293 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 294 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 295 ierr = DMSetDefaultSection(dmGrad, sectionGrad);CHKERRQ(ierr); 296 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 297 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad", (PetscObject) dmGrad);CHKERRQ(ierr); 298 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 299 PetscFunctionReturn(0); 300 } 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "DMPlexInsertBoundaryValuesFVM_Static" 304 static PetscErrorCode DMPlexInsertBoundaryValuesFVM_Static(DM dm, PetscFV fvm, PetscReal time, Vec locX, Vec Grad) 305 { 306 Vec faceGeometry, cellGeometry; 307 DM dmFace, dmCell, dmGrad; 308 const PetscScalar *facegeom, *cellgeom, *grad; 309 PetscScalar *x, *fx; 310 PetscInt numBd, b, dim, pdim, fStart, fEnd; 311 PetscErrorCode ierr; 312 313 PetscFunctionBegin; 314 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 315 ierr = DMPlexTSGetGeometry(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 316 ierr = DMPlexTSGetGradientDM(dm, &dmGrad);CHKERRQ(ierr); 317 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 318 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 319 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 320 if (Grad) { 321 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 322 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 323 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 324 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 325 } 326 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 327 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 328 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 329 for (b = 0; b < numBd; ++b) { 330 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*); 331 DMLabel label; 332 const char *labelname; 333 const PetscInt *ids; 334 PetscInt numids, i; 335 void *ctx; 336 337 ierr = DMPlexGetBoundary(dm, b, NULL, NULL, &labelname, NULL, (void (**)()) &func, &numids, &ids, &ctx);CHKERRQ(ierr); 338 ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr); 339 for (i = 0; i < numids; ++i) { 340 IS faceIS; 341 const PetscInt *faces; 342 PetscInt numFaces, f; 343 344 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 345 if (!faceIS) continue; /* No points with that id on this process */ 346 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 347 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 348 for (f = 0; f < numFaces; ++f) { 349 const PetscInt face = faces[f], *cells; 350 const FaceGeom *fg; 351 352 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 353 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 354 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 355 if (Grad) { 356 const CellGeom *cg; 357 const PetscScalar *cx, *cgrad; 358 PetscScalar *xG; 359 PetscReal dx[3]; 360 PetscInt d; 361 362 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 363 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 364 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 365 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 366 WaxpyD(dim, -1, cg->centroid, fg->centroid, dx); 367 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DotD(dim, &cgrad[d*dim], dx); 368 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 369 } else { 370 const PetscScalar *xI; 371 PetscScalar *xG; 372 373 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 374 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 375 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 376 } 377 } 378 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 379 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 380 } 381 } 382 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 383 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 384 if (Grad) { 385 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 386 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 387 } 388 PetscFunctionReturn(0); 389 } 390 391 #undef __FUNCT__ 392 #define __FUNCT__ "DMPlexTSComputeRHSFunctionFVM" 393 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *ctx) 394 { 395 PetscDS prob; 396 void (*riemann)(const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[], void *); 397 void *rctx; 398 PetscFV fvm; 399 PetscLimiter lim; 400 Vec faceGeometry, cellGeometry; 401 Vec Grad = NULL, locGrad; 402 DM dmFace, dmCell, dmGrad; 403 DMLabel ghostLabel; 404 PetscCellGeometry fgeom, cgeom; 405 const PetscScalar *facegeom, *cellgeom, *x, *lgrad; 406 PetscScalar *grad, *f, *uL, *uR, *fluxL, *fluxR; 407 PetscReal *centroid, *normal, *vol, *cellPhi; 408 PetscBool computeGradients; 409 PetscInt Nf, dim, pdim, fStart, fEnd, numFaces = 0, face, iface, cell, cStart, cEnd, cEndInterior; 410 PetscErrorCode ierr; 411 412 PetscFunctionBegin; 413 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 414 PetscValidHeaderSpecific(locX,VEC_CLASSID,3); 415 PetscValidHeaderSpecific(F,VEC_CLASSID,5); 416 ierr = DMPlexTSGetGeometry(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 417 ierr = DMPlexTSGetGradientDM(dm, &dmGrad);CHKERRQ(ierr); 418 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 419 ierr = PetscDSGetRiemannSolver(prob, 0, &riemann);CHKERRQ(ierr); 420 ierr = PetscDSGetContext(prob, 0, &rctx);CHKERRQ(ierr); 421 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 422 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 423 ierr = DMGetField(dm, 0, (PetscObject *) &fvm);CHKERRQ(ierr); 424 ierr = PetscFVGetLimiter(fvm, &lim);CHKERRQ(ierr); 425 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 426 ierr = PetscFVGetComputeGradients(fvm, &computeGradients);CHKERRQ(ierr); 427 if (computeGradients) { 428 ierr = DMGetGlobalVector(dmGrad, &Grad);CHKERRQ(ierr); 429 ierr = VecZeroEntries(Grad);CHKERRQ(ierr); 430 ierr = VecGetArray(Grad, &grad);CHKERRQ(ierr); 431 } 432 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 433 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 434 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 435 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 436 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 437 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 438 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 439 /* Count faces and reconstruct gradients */ 440 for (face = fStart; face < fEnd; ++face) { 441 const PetscInt *cells; 442 const FaceGeom *fg; 443 const PetscScalar *cx[2]; 444 PetscScalar *cgrad[2]; 445 PetscBool boundary; 446 PetscInt ghost, c, pd, d; 447 448 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 449 if (ghost >= 0) continue; 450 ++numFaces; 451 if (!computeGradients) continue; 452 ierr = DMPlexIsBoundaryPoint(dm, face, &boundary);CHKERRQ(ierr); 453 if (boundary) continue; 454 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 455 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 456 for (c = 0; c < 2; ++c) { 457 ierr = DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);CHKERRQ(ierr); 458 ierr = DMPlexPointGlobalRef(dmGrad, cells[c], grad, &cgrad[c]);CHKERRQ(ierr); 459 } 460 for (pd = 0; pd < pdim; ++pd) { 461 PetscScalar delta = cx[1][pd] - cx[0][pd]; 462 463 for (d = 0; d < dim; ++d) { 464 if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta; 465 if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta; 466 } 467 } 468 } 469 /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */ 470 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 471 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 472 ierr = DMGetWorkArray(dm, pdim, PETSC_REAL, &cellPhi);CHKERRQ(ierr); 473 for (cell = computeGradients && lim ? cStart : cEnd; cell < cEndInterior; ++cell) { 474 const PetscInt *faces; 475 const PetscScalar *cx; 476 const CellGeom *cg; 477 PetscScalar *cgrad; 478 PetscInt coneSize, f, pd, d; 479 480 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 481 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 482 ierr = DMPlexPointLocalRead(dm, cell, x, &cx);CHKERRQ(ierr); 483 ierr = DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);CHKERRQ(ierr); 484 ierr = DMPlexPointGlobalRef(dmGrad, cell, grad, &cgrad);CHKERRQ(ierr); 485 if (!cgrad) continue; /* Unowned overlap cell, we do not compute */ 486 /* Limiter will be minimum value over all neighbors */ 487 for (d = 0; d < pdim; ++d) cellPhi[d] = PETSC_MAX_REAL; 488 for (f = 0; f < coneSize; ++f) { 489 const PetscScalar *ncx; 490 const CellGeom *ncg; 491 const PetscInt *fcells; 492 PetscInt face = faces[f], ncell, ghost; 493 PetscReal v[3]; 494 PetscBool boundary; 495 496 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 497 ierr = DMPlexIsBoundaryPoint(dm, face, &boundary);CHKERRQ(ierr); 498 if ((ghost >= 0) || boundary) continue; 499 ierr = DMPlexGetSupport(dm, face, &fcells);CHKERRQ(ierr); 500 ncell = cell == fcells[0] ? fcells[1] : fcells[0]; 501 ierr = DMPlexPointLocalRead(dm, ncell, x, &ncx);CHKERRQ(ierr); 502 ierr = DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);CHKERRQ(ierr); 503 WaxpyD(dim, -1, cg->centroid, ncg->centroid, v); 504 for (d = 0; d < pdim; ++d) { 505 /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */ 506 PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / DotD(dim, &cgrad[d*dim], v); 507 508 ierr = PetscLimiterLimit(lim, flim, &phi);CHKERRQ(ierr); 509 cellPhi[d] = PetscMin(cellPhi[d], phi); 510 } 511 } 512 /* Apply limiter to gradient */ 513 for (pd = 0; pd < pdim; ++pd) 514 /* Scalar limiter applied to each component separately */ 515 for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd]; 516 } 517 ierr = DMRestoreWorkArray(dm, pdim, PETSC_REAL, &cellPhi);CHKERRQ(ierr); 518 ierr = DMPlexInsertBoundaryValuesFVM_Static(dm, fvm, time, locX, Grad);CHKERRQ(ierr); 519 if (computeGradients) { 520 ierr = VecRestoreArray(Grad, &grad);CHKERRQ(ierr); 521 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 522 ierr = DMGlobalToLocalBegin(dmGrad, Grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 523 ierr = DMGlobalToLocalEnd(dmGrad, Grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 524 ierr = DMRestoreGlobalVector(dmGrad, &Grad);CHKERRQ(ierr); 525 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 526 } 527 ierr = PetscMalloc7(numFaces*dim,¢roid,numFaces*dim,&normal,numFaces*2,&vol,numFaces*pdim,&uL,numFaces*pdim,&uR,numFaces*pdim,&fluxL,numFaces*pdim,&fluxR);CHKERRQ(ierr); 528 /* Read out values */ 529 for (face = fStart, iface = 0; face < fEnd; ++face) { 530 const PetscInt *cells; 531 const FaceGeom *fg; 532 const CellGeom *cgL, *cgR; 533 const PetscScalar *xL, *xR, *gL, *gR; 534 PetscInt ghost, d; 535 536 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 537 if (ghost >= 0) continue; 538 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 539 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 540 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 541 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 542 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xL);CHKERRQ(ierr); 543 ierr = DMPlexPointLocalRead(dm, cells[1], x, &xR);CHKERRQ(ierr); 544 if (computeGradients) { 545 PetscReal dxL[3], dxR[3]; 546 547 ierr = DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);CHKERRQ(ierr); 548 ierr = DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);CHKERRQ(ierr); 549 WaxpyD(dim, -1, cgL->centroid, fg->centroid, dxL); 550 WaxpyD(dim, -1, cgR->centroid, fg->centroid, dxR); 551 for (d = 0; d < pdim; ++d) { 552 uL[iface*pdim+d] = xL[d] + DotD(dim, &gL[d*dim], dxL); 553 uR[iface*pdim+d] = xR[d] + DotD(dim, &gR[d*dim], dxR); 554 } 555 } else { 556 for (d = 0; d < pdim; ++d) { 557 uL[iface*pdim+d] = xL[d]; 558 uR[iface*pdim+d] = xR[d]; 559 } 560 } 561 for (d = 0; d < dim; ++d) { 562 centroid[iface*dim+d] = fg->centroid[d]; 563 normal[iface*dim+d] = fg->normal[d]; 564 } 565 vol[iface*2+0] = cgL->volume; 566 vol[iface*2+1] = cgR->volume; 567 ++iface; 568 } 569 if (computeGradients) { 570 ierr = VecRestoreArrayRead(locGrad,&lgrad);CHKERRQ(ierr); 571 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 572 } 573 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 574 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 575 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 576 fgeom.v0 = centroid; 577 fgeom.n = normal; 578 cgeom.vol = vol; 579 /* Riemann solve */ 580 ierr = PetscFVIntegrateRHSFunction(fvm, numFaces, Nf, &fvm, 0, fgeom, cgeom, uL, uR, riemann, fluxL, fluxR, rctx);CHKERRQ(ierr); 581 /* Insert fluxes */ 582 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 583 for (face = fStart, iface = 0; face < fEnd; ++face) { 584 const PetscInt *cells; 585 PetscScalar *fL, *fR; 586 PetscInt ghost, d; 587 588 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 589 if (ghost >= 0) continue; 590 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 591 ierr = DMPlexPointGlobalRef(dm, cells[0], f, &fL);CHKERRQ(ierr); 592 ierr = DMPlexPointGlobalRef(dm, cells[1], f, &fR);CHKERRQ(ierr); 593 for (d = 0; d < pdim; ++d) { 594 if (fL) fL[d] -= fluxL[iface*pdim+d]; 595 if (fR) fR[d] += fluxR[iface*pdim+d]; 596 } 597 ++iface; 598 } 599 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 600 ierr = PetscFree7(centroid,normal,vol,uL,uR,fluxL,fluxR);CHKERRQ(ierr); 601 PetscFunctionReturn(0); 602 } 603 604 #undef __FUNCT__ 605 #define __FUNCT__ "DMPlexTSSetRHSFunctionLocal" 606 /*@C 607 DMPlexTSSetRHSFunctionLocal - set a local residual evaluation function 608 609 Logically Collective 610 611 Input Arguments: 612 + dm - DM to associate callback with 613 . func - the RHS evaluation function 614 - ctx - an optional user context 615 616 Level: beginner 617 618 .seealso: DMTSSetRHSFunctionLocal() 619 @*/ 620 PetscErrorCode DMPlexTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal time, Vec X, Vec F, void *ctx), void *ctx) 621 { 622 DMTS dmts; 623 PetscFV fvm; 624 PetscInt Nf; 625 PetscBool computeGradients; 626 PetscErrorCode ierr; 627 628 PetscFunctionBegin; 629 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 630 ierr = DMGetDMTSWrite(dm, &dmts);CHKERRQ(ierr); 631 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 632 ierr = DMGetField(dm, 0, (PetscObject *) &fvm);CHKERRQ(ierr); 633 ierr = DMPlexTSSetupGeometry(dm, fvm, dmts);CHKERRQ(ierr); 634 ierr = PetscFVGetComputeGradients(fvm, &computeGradients);CHKERRQ(ierr); 635 if (computeGradients) {ierr = DMPlexTSSetupGradient(dm, fvm, dmts);CHKERRQ(ierr);} 636 ierr = DMTSSetRHSFunctionLocal(dm, func, ctx);CHKERRQ(ierr); 637 PetscFunctionReturn(0); 638 } 639