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