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