1 #include <petsc-private/dmpleximpl.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 #undef __FUNCT__ 7 #define __FUNCT__ "DMPlexTSGetGeometryFVM" 8 /*@ 9 DMPlexTSGetGeometryFVM - Return precomputed geometric data 10 11 Input Parameter: 12 . dm - The DM 13 14 Output Parameters: 15 + facegeom - The values precomputed from face geometry 16 . cellgeom - The values precomputed from cell geometry 17 - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell 18 19 Level: developer 20 21 .seealso: DMPlexTSSetRHSFunctionLocal() 22 @*/ 23 PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius) 24 { 25 DMTS dmts; 26 PetscObject obj; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 31 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 32 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", &obj);CHKERRQ(ierr); 33 if (!obj) { 34 Vec cellgeom, facegeom; 35 36 ierr = DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);CHKERRQ(ierr); 37 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject) facegeom);CHKERRQ(ierr); 38 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject) cellgeom);CHKERRQ(ierr); 39 ierr = VecDestroy(&facegeom);CHKERRQ(ierr); 40 ierr = VecDestroy(&cellgeom);CHKERRQ(ierr); 41 } 42 if (facegeom) {PetscValidPointer(facegeom, 2); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject *) facegeom);CHKERRQ(ierr);} 43 if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject *) cellgeom);CHKERRQ(ierr);} 44 if (minRadius) {ierr = DMPlexGetMinRadius(dm, minRadius);CHKERRQ(ierr);} 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "BuildGradientReconstruction" 50 static PetscErrorCode BuildGradientReconstruction(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 51 { 52 DMLabel ghostLabel; 53 PetscScalar *dx, *grad, **gref; 54 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 59 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 60 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 61 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 62 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 63 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 64 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 65 for (c = cStart; c < cEndInterior; c++) { 66 const PetscInt *faces; 67 PetscInt numFaces, usedFaces, f, d; 68 const PetscFVCellGeom *cg; 69 PetscBool boundary; 70 PetscInt ghost; 71 72 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 73 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 74 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 75 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 76 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 77 const PetscFVCellGeom *cg1; 78 PetscFVFaceGeom *fg; 79 const PetscInt *fcells; 80 PetscInt ncell, side; 81 82 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 83 ierr = DMPlexIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 84 if ((ghost >= 0) || boundary) continue; 85 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 86 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 87 ncell = fcells[!side]; /* the neighbor */ 88 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 89 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 90 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 91 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 92 } 93 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 94 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 95 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 96 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 97 ierr = DMPlexIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 98 if ((ghost >= 0) || boundary) continue; 99 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 100 ++usedFaces; 101 } 102 } 103 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNCT__ 108 #define __FUNCT__ "DMPlexTSSetupGradientFVM_Internal" 109 static PetscErrorCode DMPlexTSSetupGradientFVM_Internal(DM dm, PetscFV fvm, DM *dmGrad) 110 { 111 DM dmFace, dmCell; 112 Vec facegeom, cellgeom; 113 PetscScalar *fgeom, *cgeom; 114 PetscSection sectionGrad; 115 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 120 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 121 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 122 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 123 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 124 ierr = DMPlexTSGetGeometryFVM(dm, &facegeom, &cellgeom, NULL);CHKERRQ(ierr); 125 ierr = VecGetDM(facegeom, &dmFace);CHKERRQ(ierr); 126 ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr); 127 ierr = VecGetArray(facegeom, &fgeom);CHKERRQ(ierr); 128 ierr = VecGetArray(cellgeom, &cgeom);CHKERRQ(ierr); 129 ierr = BuildGradientReconstruction(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 130 ierr = VecRestoreArray(facegeom, &fgeom);CHKERRQ(ierr); 131 ierr = VecRestoreArray(cellgeom, &cgeom);CHKERRQ(ierr); 132 /* Create storage for gradients */ 133 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 134 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 135 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 136 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 137 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 138 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 139 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "DMPlexTSGetGradientDM" 145 /*@C 146 DMPlexTSGetGradientDM - Return gradient data layout 147 148 Input Parameters: 149 + dm - The DM 150 - fv - The PetscFV 151 152 Output Parameter: 153 . dmGrad - The layout for gradient values 154 155 Level: developer 156 157 .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal() 158 @*/ 159 PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad) 160 { 161 DMTS dmts; 162 PetscObject obj; 163 PetscBool computeGradients; 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 168 PetscValidHeaderSpecific(fv,PETSCFV_CLASSID,2); 169 PetscValidPointer(dmGrad,3); 170 ierr = PetscFVGetComputeGradients(fv, &computeGradients);CHKERRQ(ierr); 171 if (!computeGradients) {*dmGrad = NULL; PetscFunctionReturn(0);} 172 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 173 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);CHKERRQ(ierr); 174 if (!obj) { 175 DM dmGrad; 176 177 ierr = DMPlexTSSetupGradientFVM_Internal(dm, fv, &dmGrad);CHKERRQ(ierr); 178 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);CHKERRQ(ierr); 179 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 180 } 181 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 #undef __FUNCT__ 186 #define __FUNCT__ "DMPlexInsertBoundaryValuesFVM_Static" 187 static PetscErrorCode DMPlexInsertBoundaryValuesFVM_Static(DM dm, PetscFV fvm, PetscReal time, Vec locX, Vec Grad) 188 { 189 Vec faceGeometry, cellGeometry; 190 DM dmFace, dmCell, dmGrad; 191 const PetscScalar *facegeom, *cellgeom, *grad; 192 PetscScalar *x, *fx; 193 PetscInt numBd, b, dim, pdim, fStart, fEnd; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 198 ierr = DMPlexTSGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 199 ierr = DMPlexTSGetGradientDM(dm, fvm, &dmGrad);CHKERRQ(ierr); 200 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 201 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 202 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 203 if (Grad) { 204 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 205 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 206 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 207 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 208 } 209 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 210 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 211 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 212 for (b = 0; b < numBd; ++b) { 213 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*); 214 DMLabel label; 215 const char *labelname; 216 const PetscInt *ids; 217 PetscInt numids, i; 218 void *ctx; 219 220 ierr = DMPlexGetBoundary(dm, b, NULL, NULL, &labelname, NULL, (void (**)()) &func, &numids, &ids, &ctx);CHKERRQ(ierr); 221 ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr); 222 for (i = 0; i < numids; ++i) { 223 IS faceIS; 224 const PetscInt *faces; 225 PetscInt numFaces, f; 226 227 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 228 if (!faceIS) continue; /* No points with that id on this process */ 229 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 230 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 231 for (f = 0; f < numFaces; ++f) { 232 const PetscInt face = faces[f], *cells; 233 const PetscFVFaceGeom *fg; 234 235 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 236 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 237 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 238 if (Grad) { 239 const PetscFVCellGeom *cg; 240 const PetscScalar *cx, *cgrad; 241 PetscScalar *xG; 242 PetscReal dx[3]; 243 PetscInt d; 244 245 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 246 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 247 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 248 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 249 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 250 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 251 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 252 } else { 253 const PetscScalar *xI; 254 PetscScalar *xG; 255 256 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 257 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 258 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 259 } 260 } 261 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 262 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 263 } 264 } 265 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 266 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 267 if (Grad) { 268 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 269 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "DMPlexTSComputeRHSFunctionFVM" 276 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *ctx) 277 { 278 PetscDS prob; 279 void (*riemann)(const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[], void *); 280 void *rctx; 281 PetscFV fvm; 282 PetscLimiter lim; 283 Vec faceGeometry, cellGeometry; 284 Vec Grad = NULL, locGrad; 285 DM dmFace, dmCell, dmGrad; 286 DMLabel ghostLabel; 287 PetscFVFaceGeom *fgeom; 288 const PetscScalar *facegeom, *cellgeom, *x, *lgrad; 289 PetscScalar *grad, *f, *uL, *uR, *fluxL, *fluxR; 290 PetscReal *vol, *cellPhi; 291 PetscInt Nf, dim, pdim, fStart, fEnd, numFaces = 0, face, iface, cell, cStart, cEnd, cEndInterior; 292 PetscErrorCode ierr; 293 294 PetscFunctionBegin; 295 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 296 PetscValidHeaderSpecific(locX,VEC_CLASSID,3); 297 PetscValidHeaderSpecific(F,VEC_CLASSID,5); 298 /* 1: Get sizes from dm and dmAux */ 299 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 300 ierr = PetscDSGetRiemannSolver(prob, 0, &riemann);CHKERRQ(ierr); 301 ierr = PetscDSGetContext(prob, 0, &rctx);CHKERRQ(ierr); 302 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 303 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 304 ierr = DMGetField(dm, 0, (PetscObject *) &fvm);CHKERRQ(ierr); 305 ierr = PetscFVGetLimiter(fvm, &lim);CHKERRQ(ierr); 306 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 307 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 308 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 309 /* 2: Get geometric data */ 310 ierr = DMPlexTSGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 311 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 312 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 313 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 314 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 315 /* 3: Get gradient data */ 316 ierr = DMPlexTSGetGradientDM(dm, fvm, &dmGrad);CHKERRQ(ierr); 317 if (dmGrad) { 318 ierr = DMGetGlobalVector(dmGrad, &Grad);CHKERRQ(ierr); 319 ierr = VecZeroEntries(Grad);CHKERRQ(ierr); 320 ierr = VecGetArray(Grad, &grad);CHKERRQ(ierr); 321 } 322 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 323 /* Count faces and reconstruct gradients */ 324 for (face = fStart; face < fEnd; ++face) { 325 const PetscInt *cells; 326 const PetscFVFaceGeom *fg; 327 const PetscScalar *cx[2]; 328 PetscScalar *cgrad[2]; 329 PetscBool boundary; 330 PetscInt ghost, c, pd, d; 331 332 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 333 if (ghost >= 0) continue; 334 ++numFaces; 335 if (!dmGrad) continue; 336 ierr = DMPlexIsBoundaryPoint(dm, face, &boundary);CHKERRQ(ierr); 337 if (boundary) continue; 338 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 339 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 340 for (c = 0; c < 2; ++c) { 341 ierr = DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);CHKERRQ(ierr); 342 ierr = DMPlexPointGlobalRef(dmGrad, cells[c], grad, &cgrad[c]);CHKERRQ(ierr); 343 } 344 for (pd = 0; pd < pdim; ++pd) { 345 PetscScalar delta = cx[1][pd] - cx[0][pd]; 346 347 for (d = 0; d < dim; ++d) { 348 if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta; 349 if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta; 350 } 351 } 352 } 353 /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */ 354 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 355 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 356 ierr = DMGetWorkArray(dm, pdim, PETSC_REAL, &cellPhi);CHKERRQ(ierr); 357 for (cell = dmGrad && lim ? cStart : cEnd; cell < cEndInterior; ++cell) { 358 const PetscInt *faces; 359 const PetscScalar *cx; 360 const PetscFVCellGeom *cg; 361 PetscScalar *cgrad; 362 PetscInt coneSize, f, pd, d; 363 364 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 365 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 366 ierr = DMPlexPointLocalRead(dm, cell, x, &cx);CHKERRQ(ierr); 367 ierr = DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);CHKERRQ(ierr); 368 ierr = DMPlexPointGlobalRef(dmGrad, cell, grad, &cgrad);CHKERRQ(ierr); 369 if (!cgrad) continue; /* Unowned overlap cell, we do not compute */ 370 /* Limiter will be minimum value over all neighbors */ 371 for (d = 0; d < pdim; ++d) cellPhi[d] = PETSC_MAX_REAL; 372 for (f = 0; f < coneSize; ++f) { 373 const PetscScalar *ncx; 374 const PetscFVCellGeom *ncg; 375 const PetscInt *fcells; 376 PetscInt face = faces[f], ncell, ghost; 377 PetscReal v[3]; 378 PetscBool boundary; 379 380 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 381 ierr = DMPlexIsBoundaryPoint(dm, face, &boundary);CHKERRQ(ierr); 382 if ((ghost >= 0) || boundary) continue; 383 ierr = DMPlexGetSupport(dm, face, &fcells);CHKERRQ(ierr); 384 ncell = cell == fcells[0] ? fcells[1] : fcells[0]; 385 ierr = DMPlexPointLocalRead(dm, ncell, x, &ncx);CHKERRQ(ierr); 386 ierr = DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);CHKERRQ(ierr); 387 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v); 388 for (d = 0; d < pdim; ++d) { 389 /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */ 390 PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / DMPlex_DotD_Internal(dim, &cgrad[d*dim], v); 391 392 ierr = PetscLimiterLimit(lim, flim, &phi);CHKERRQ(ierr); 393 cellPhi[d] = PetscMin(cellPhi[d], phi); 394 } 395 } 396 /* Apply limiter to gradient */ 397 for (pd = 0; pd < pdim; ++pd) 398 /* Scalar limiter applied to each component separately */ 399 for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd]; 400 } 401 ierr = DMRestoreWorkArray(dm, pdim, PETSC_REAL, &cellPhi);CHKERRQ(ierr); 402 /* Handle boundary values */ 403 ierr = DMPlexInsertBoundaryValuesFVM_Static(dm, fvm, time, locX, Grad);CHKERRQ(ierr); 404 if (dmGrad) { 405 ierr = VecRestoreArray(Grad, &grad);CHKERRQ(ierr); 406 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 407 ierr = DMGlobalToLocalBegin(dmGrad, Grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 408 ierr = DMGlobalToLocalEnd(dmGrad, Grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 409 ierr = DMRestoreGlobalVector(dmGrad, &Grad);CHKERRQ(ierr); 410 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 411 } 412 ierr = PetscMalloc6(numFaces,&fgeom,numFaces*2,&vol,numFaces*pdim,&uL,numFaces*pdim,&uR,numFaces*pdim,&fluxL,numFaces*pdim,&fluxR);CHKERRQ(ierr); 413 /* Read out values */ 414 for (face = fStart, iface = 0; face < fEnd; ++face) { 415 const PetscInt *cells; 416 const PetscFVFaceGeom *fg; 417 const PetscFVCellGeom *cgL, *cgR; 418 const PetscScalar *xL, *xR, *gL, *gR; 419 PetscInt ghost, d; 420 421 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 422 if (ghost >= 0) continue; 423 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 424 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 425 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 426 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 427 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xL);CHKERRQ(ierr); 428 ierr = DMPlexPointLocalRead(dm, cells[1], x, &xR);CHKERRQ(ierr); 429 if (dmGrad) { 430 PetscReal dxL[3], dxR[3]; 431 432 ierr = DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);CHKERRQ(ierr); 433 ierr = DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);CHKERRQ(ierr); 434 DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL); 435 DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR); 436 for (d = 0; d < pdim; ++d) { 437 uL[iface*pdim+d] = xL[d] + DMPlex_DotD_Internal(dim, &gL[d*dim], dxL); 438 uR[iface*pdim+d] = xR[d] + DMPlex_DotD_Internal(dim, &gR[d*dim], dxR); 439 } 440 } else { 441 for (d = 0; d < pdim; ++d) { 442 uL[iface*pdim+d] = xL[d]; 443 uR[iface*pdim+d] = xR[d]; 444 } 445 } 446 for (d = 0; d < dim; ++d) { 447 fgeom[iface].centroid[d] = fg->centroid[d]; 448 fgeom[iface].normal[d] = fg->normal[d]; 449 } 450 vol[iface*2+0] = cgL->volume; 451 vol[iface*2+1] = cgR->volume; 452 ++iface; 453 } 454 if (dmGrad) { 455 ierr = VecRestoreArrayRead(locGrad,&lgrad);CHKERRQ(ierr); 456 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 457 } 458 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 459 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 460 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 461 /* Riemann solve */ 462 ierr = PetscFVIntegrateRHSFunction(fvm, numFaces, Nf, &fvm, 0, fgeom, vol, uL, uR, riemann, fluxL, fluxR, rctx);CHKERRQ(ierr); 463 /* Insert fluxes */ 464 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 465 for (face = fStart, iface = 0; face < fEnd; ++face) { 466 const PetscInt *cells; 467 PetscScalar *fL, *fR; 468 PetscInt ghost, d; 469 470 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 471 if (ghost >= 0) continue; 472 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 473 ierr = DMPlexPointGlobalRef(dm, cells[0], f, &fL);CHKERRQ(ierr); 474 ierr = DMPlexPointGlobalRef(dm, cells[1], f, &fR);CHKERRQ(ierr); 475 for (d = 0; d < pdim; ++d) { 476 if (fL) fL[d] -= fluxL[iface*pdim+d]; 477 if (fR) fR[d] += fluxR[iface*pdim+d]; 478 } 479 ++iface; 480 } 481 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 482 ierr = PetscFree6(fgeom,vol,uL,uR,fluxL,fluxR);CHKERRQ(ierr); 483 PetscFunctionReturn(0); 484 } 485 486 #undef __FUNCT__ 487 #define __FUNCT__ "DMPlexTSComputeIFunctionFEM" 488 /*@ 489 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 490 491 Input Parameters: 492 + dm - The mesh 493 . t - The time 494 . X - Local solution 495 . X_t - Local solution time derivative, or NULL 496 - user - The user context 497 498 Output Parameter: 499 . F - Local output vector 500 501 Level: developer 502 503 .seealso: DMPlexComputeJacobianActionFEM() 504 @*/ 505 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, void *user) 506 { 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 ierr = DMPlexComputeResidualFEM_Internal(dm, X, X_t, F, user);CHKERRQ(ierr); 511 PetscFunctionReturn(0); 512 } 513