1 #include <petscdmplex.h> /*I "petscdmplex.h" I*/ 2 #include <petscfv.h> 3 4 struct _FVPhysics { 5 void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx); 6 }; 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "TSComputeRHSFunctionLocalDMPlex" 10 PetscErrorCode TSComputeRHSFunctionLocalDMPlex(DM dm, PetscReal time, Vec locX, Vec F, void *ctx) 11 { 12 void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx) = ((struct _FVPhysics *) ctx)->riemann; 13 PetscFV fvm; 14 Vec faceGeometry, cellGeometry; 15 DM dmFace, dmCell; 16 DMLabel ghostLabel; 17 PetscCellGeometry fgeom, cgeom; 18 const PetscScalar *facegeom, *cellgeom, *x; 19 PetscScalar *f, *uL, *uR, *fluxL, *fluxR; 20 PetscReal *centroid, *normal, *vol; 21 PetscInt Nf, dim, pdim, fStart, fEnd, numFaces = 0, face, iface; 22 PetscErrorCode ierr; 23 24 PetscFunctionBeginUser; 25 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 26 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 27 ierr = DMGetField(dm, 0, (PetscObject *) &fvm);CHKERRQ(ierr); 28 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 29 ierr = DMPlexInsertBoundaryValuesFVM(dm, time, locX);CHKERRQ(ierr); 30 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 31 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 32 /* TODO Pull this geometry calculation into the library */ 33 ierr = PetscObjectQuery((PetscObject) dm, "FaceGeometry", (PetscObject *) &faceGeometry);CHKERRQ(ierr); 34 ierr = PetscObjectQuery((PetscObject) dm, "CellGeometry", (PetscObject *) &cellGeometry);CHKERRQ(ierr); 35 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 36 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 37 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 38 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 39 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 40 for (face = fStart; face < fEnd; ++face) { 41 PetscInt ghost; 42 43 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 44 if (ghost >= 0) continue; 45 ++numFaces; 46 } 47 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); 48 for (face = fStart, iface = 0; face < fEnd; ++face) { 49 const PetscInt *cells; 50 const FaceGeom *fg; 51 const CellGeom *cgL, *cgR; 52 const PetscScalar *xL, *xR; 53 PetscInt ghost, d; 54 55 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 56 if (ghost >= 0) continue; 57 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 58 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 59 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 60 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 61 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xL);CHKERRQ(ierr); 62 ierr = DMPlexPointLocalRead(dm, cells[1], x, &xR);CHKERRQ(ierr); 63 for (d = 0; d < pdim; ++d) { 64 uL[iface*pdim+d] = xL[d]; 65 uR[iface*pdim+d] = xR[d]; 66 } 67 for (d = 0; d < dim; ++d) { 68 centroid[iface*dim+d] = fg->centroid[d]; 69 normal[iface*dim+d] = fg->normal[d]; 70 } 71 vol[iface*2+0] = cgL->volume; 72 vol[iface*2+1] = cgR->volume; 73 ++iface; 74 } 75 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 76 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 77 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 78 fgeom.v0 = centroid; 79 fgeom.n = normal; 80 cgeom.vol = vol; 81 ierr = PetscFVIntegrateRHSFunction(fvm, numFaces, Nf, &fvm, 0, fgeom, cgeom, uL, uR, riemann, fluxL, fluxR, ctx);CHKERRQ(ierr); 82 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 83 for (face = fStart, iface = 0; face < fEnd; ++face) { 84 const PetscInt *cells; 85 PetscScalar *fL, *fR; 86 PetscInt ghost, d; 87 88 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 89 if (ghost >= 0) continue; 90 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 91 ierr = DMPlexPointGlobalRef(dm, cells[0], f, &fL);CHKERRQ(ierr); 92 ierr = DMPlexPointGlobalRef(dm, cells[1], f, &fR);CHKERRQ(ierr); 93 for (d = 0; d < pdim; ++d) { 94 if (fL) fL[d] -= fluxL[iface*pdim+d]; 95 if (fR) fR[d] += fluxR[iface*pdim+d]; 96 } 97 ++iface; 98 } 99 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 100 ierr = PetscFree7(centroid,normal,vol,uL,uR,fluxL,fluxR);CHKERRQ(ierr); 101 PetscFunctionReturn(0); 102 } 103