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