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