Thank you for the help. I followed your suggestion and removed the "DMPlexConstructGhostCells" function and now I am using a gmsh generated mesh with periodicity in both directions, and distributing the mesh with 1 level of overlap. I can see that the periodic faces now have two neighbouring cells on the opposite side of the domain which I can use to get the periodic cell values using "DMGlobalToLocal" without the need to construct ghost cells. There are two main issues that I am facing now, regarding geometric information needed for implementing periodic bc.
One, I may need the centroid of periodic counterpart of a cell which is possible to get from a "ghost" cell created on the periodic face, and containing the solution data from periodic counterpart. But now, I would need to shift the coordinates of the neighbouring cell accordingly. If you see the sample output given below along with the screenshot of the sample 2D mesh that I have used (please find attached), while reconstructing the left and right solution values on face 32 for a second order scheme, I need the centroid of a ghost cell on left side of the face which is at (x, y) = (-1.25, 8.75). I wanted to know if there is some way to get these shifted "ghost" coordinates directly from the geometric information obtained using "PetscFV" class. Kindly let me know if you have some suggestions regarding this. Sample Output: face: 32, centroid = 0.000000, 8.750000, normal = 2.500000, -0.000000 support[0] = 0, cent = 1.250000, 8.750000, area = 6.250000 support[1] = 3, cent = 8.750000, 8.750000, area = 6.250000 face: 33, centroid = 1.250000, 7.500000, normal = -0.000000, -2.500000 support[0] = 0, cent = 1.250000, 8.750000, area = 6.250000 support[1] = 4, cent = 1.250000, 6.250000, area = 6.250000 Second, the face normals of the periodic faces obtained using the face geometry data from "DMPlexGetDataFVM" function seem to be in the wrong direction i.e., from the right cell "support[1]" to the left cell "support[0]", while the other faces have the opposite direction. For face: 32, left cell is 0 and the right cell is 3, but, the normal points the opposite way (cell 3 would be on the left side of cell 0 on periodic boundary), while, for face: 33 which is an interior face, the normal direction is from left to right cell. Also, I am getting incorrect normal values for some of the faces (including some interior faces). Here is a sample output of the same: face: 59, centroid = 3.750000, 2.500000, normal = 0.000000, -7.500000 support[0] = 11, cent = 8.750000, 3.750000, area = 6.250000 support[1] = 15, cent = 8.750000, 1.250000, area = 6.250000 this face has the correct direction (from 11 to 15), but the value is wrong (all the faces are 2.5 units). This problem does not occur when I ignore periodicity by setting -dm_plex_gmsh_periodic to 0. These issues occur regardless of the number of processes. I am attaching the test code that I am using along with the gmsh file for your reference. Please suggest to me, what might be the issue here. Regards, Shashwat On Mon, Jun 15, 2020 at 2:54 PM Stefano Zampini <[email protected]> wrote: > It is enough if you use DMPlexDistribute with 1 level of overlap, set the > local section for your dofs, and call DMGlobalToLocal . The local vector > will contain data for what you call “ghost” cells. I call them “not-owned’ > cells. > > On Jun 15, 2020, at 12:09 PM, Shashwat Tiwari <[email protected]> > wrote: > > The way I'm trying to implement periodic bc is, when I loop over the > boundary faces, say, on the left boundary of the domain to compute flux and > residual, I need solution values from the two cells neighbouring the face, > i.e. the left cell and the right cell (the face normal pointing from the > left cell to the right cell and the right cell being a boundary ghost cell > for boundary faces). For the boundary to be periodic, I need the value that > I get from the right cell (boundary ghost cell) to be the solution value at > its periodic counterpart, i.e. the solution value at left cell of the face > on right boundary of the domain in this case. My question is how do I > update the value at a boundary ghost cell with the value at the real cell > which is its periodic counterpart from the other side of the domain. Is > there some kind of mapping of the boundary ghost cells to their > corresponding real cells which I can use to update the solution values at > ghost cells? > > Regards, > Shashwat > > On Sun, Jun 14, 2020 at 5:11 AM Matthew Knepley <[email protected]> wrote: > >> On Fri, Jun 12, 2020 at 3:19 PM Shashwat Tiwari <[email protected]> >> wrote: >> >>> Hi, >>> I am writing a first order 2D solver for unstructured grids with >>> periodic boundaries using DMPlex. After generating the mesh, I use >>> "DMSetPeriodicity" function to set periodicity in both directions. After >>> which I partition the mesh (DMPlexDistribute), construct ghost cells >>> (DMPlexConstructGhostCells), >>> >> >> These ghost cells are for FVM boundary conditions. If you want cells to >> be shared across parallel partitions, then you want to give overlap=1 >> to DMPlexDIstribute(). Is that what you want? >> >> Thanks, >> >> Matt >> >> >>> create a section, and set some initial values in the global vector. Then >>> I use "VecGhostUpdateBegin" to start updating the boundary ghost cell >>> values, but, I get the following error in case I use multiple processors: >>> >>> [0]PETSC ERROR: Invalid argument >>> [0]PETSC ERROR: Vector is not ghosted >>> [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> >>> if I run with a single process, there is no error but the values remain >>> empty (zero) and are not updated. Kindly let me know, if I am missing some >>> crucial step before I can update the ghost values in order to implement the >>> periodic bc, or if there is any other approach to achieve it. I am >>> attaching a small code to demonstrate the issue for your reference. >>> >>> Regards, >>> Shashwat >>> >> >> >> -- >> What most experimenters take for granted before they begin their >> experiments is infinitely more interesting than any results to which their >> experiments lead. >> -- Norbert Wiener >> >> https://www.cse.buffalo.edu/~knepley/ >> <http://www.cse.buffalo.edu/~knepley/> >> > >
static char help[] = "testing periodic bc\n";
#include <petsc.h>
#define nvar 2
int main(int argc, char **argv)
{
PetscErrorCode ierr;
DM dm, dmDist = NULL;
PetscSection sec;
PetscInt dim = 2, numFields = 1, overlap = 1, numBC, i;
PetscInt numComp[1]; // number of components per field
PetscInt numDof[3];
PetscBool interpolate = PETSC_TRUE, useCone = PETSC_TRUE, useClosure = PETSC_TRUE;
Vec U; // solution
PetscMPIInt rank, size;
ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);
// read mesh from gmsh file
ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, "ug_periodic.msh", interpolate, &dm); CHKERRQ(ierr);
// set periodicity information (periodic in x and y)
const DMBoundaryType bd[] = {DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC};
ierr = DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, bd); CHKERRQ(ierr);
// set partitioner type to Parmetis
PetscPartitioner part;
ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr);
ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE); CHKERRQ(ierr);
ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr);
// distribute mesh over processes;
ierr = DMSetBasicAdjacency(dm, useCone, useClosure); CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist); CHKERRQ(ierr);
if(dmDist)
{
ierr = DMDestroy(&dm); CHKERRQ(ierr);
dm = dmDist;
}
ierr = DMSetNumFields(dm, numFields); CHKERRQ(ierr);
ierr = DMSetAdjacency(dm, 0, useCone, useClosure); CHKERRQ(ierr);
// setup FV
PetscFV fv;
ierr = PetscFVCreate(PETSC_COMM_WORLD, &fv); CHKERRQ(ierr);
ierr = PetscFVSetType(fv, PETSCFVUPWIND); CHKERRQ(ierr);
ierr = PetscFVSetNumComponents(fv, nvar); CHKERRQ(ierr);
ierr = PetscFVSetSpatialDimension(fv, 2); CHKERRQ(ierr);
ierr = PetscFVSetFromOptions(fv); CHKERRQ(ierr);
ierr = PetscFVSetUp(fv); CHKERRQ(ierr);
// create field for solution u
numComp[0] = nvar; // six components
for(i=0; i<numFields*(dim+1); ++i) numDof[i] = 0;
// define u at cell centers
numDof[0*(dim+1)+dim] = nvar;
numBC = 0;
// create a PetscSection with this layout
ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC,
NULL, NULL, NULL, NULL, &sec); CHKERRQ(ierr);
// name the field and its components
ierr = PetscSectionSetFieldName(sec, 0, "U"); CHKERRQ(ierr);
// set the section for dm
ierr = DMSetLocalSection(dm, sec); CHKERRQ(ierr);
ierr = DMSetUp(dm); CHKERRQ(ierr);
ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);
// get global vector for solution and set IC
ierr = DMCreateGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) U, "Sol"); CHKERRQ(ierr);
// output initial condition
PetscViewer viewer;
ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm), &viewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, "sol-000.vtu");CHKERRQ(ierr);
ierr = VecView(U, viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Wrote sol-000.vtu\n"); CHKERRQ(ierr);
// get geometric data
DM dmFace, dmCell, dummy;
Vec cellGeom, faceGeom;
const PetscScalar *fgeom, *cgeom; // arrays corresponding to vectors
PetscFVCellGeom *cg; // struct with cell geometry information
PetscFVFaceGeom *fg; // struct with face geometry information
ierr = DMPlexGetDataFVM(dm, fv, &cellGeom, &faceGeom, &dummy); CHKERRQ(ierr);
ierr = VecGetDM(faceGeom, &dmFace); CHKERRQ(ierr);
ierr = VecGetDM(cellGeom, &dmCell); CHKERRQ(ierr);
ierr = VecGetArrayRead(faceGeom, &fgeom); CHKERRQ(ierr);
ierr = VecGetArrayRead(cellGeom, &cgeom); CHKERRQ(ierr);
// print geometric information
PetscInt f, fStart, fEnd, ss;
const PetscInt *supp; // support
ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd); CHKERRQ(ierr);
for(f=fStart; f<fEnd; ++f)
{
ierr = DMPlexGetSupportSize(dm, f, &ss); CHKERRQ(ierr);
ierr = DMPlexGetSupport(dm, f, &supp); CHKERRQ(ierr);
// face geometry
ierr = DMPlexPointLocalRead(dmFace, f, fgeom, &fg); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "face: %d, centroid = %f, %f, normal = %f, %f\n", f, fg->centroid[0], fg->centroid[1], fg->normal[0], fg->normal[1]);
// cell geometry
for(PetscInt i=0; i<ss; ++i)
{
ierr = DMPlexPointLocalRead(dmCell, supp[i], cgeom, &cg); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "support[%d] = %d, cent = %f, %f, area = %f\n", i, supp[i], cg->centroid[0], cg->centroid[1], cg->volume);
}
PetscPrintf(PETSC_COMM_WORLD, "\n");
}
// cleanup
ierr = PetscFVDestroy(&fv); CHKERRQ(ierr);
ierr = VecRestoreArrayRead(faceGeom, &fgeom); CHKERRQ(ierr);
ierr = VecRestoreArrayRead(cellGeom, &cgeom); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr);
ierr = DMDestroy(&dm); CHKERRQ(ierr);
ierr = PetscFinalize(); CHKERRQ(ierr);
return ierr;
}
ug_periodic.geo
Description: Binary data
