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), 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
static char help[] = "testing periodic bc\n";
#include <petsc.h>
#define nvar 2
// set initial condition for the given dm and vector
PetscErrorCode SetIC(DM dm, Vec U)
{
PetscErrorCode ierr;
PetscScalar *uArr, *u;
PetscFunctionBegin;
PetscInt c, cStart, cEnd; // cells
// get cell stratum owned by processor
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
// get array for U
ierr = VecGetArray(U, &uArr); CHKERRQ(ierr);
// loop over cells and assign values
for(c=cStart; c<cEnd; ++c)
{
PetscInt label;
ierr = DMGetLabelValue(dm, "vtk", c, &label); CHKERRQ(ierr);
// write into Global vector if the cell is a real cell
if(label == 1)
{
ierr = DMPlexPointLocalRef(dm, c, uArr, &u); CHKERRQ(ierr);
for(PetscInt i=0; i<nvar; ++i)
u[i] = i+1;
}
}
ierr = VecRestoreArray(U, &uArr);
PetscFunctionReturn(0);
}
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, status;
Vec U; // solution
PetscMPIInt rank, size;
PetscScalar lower[3], upper[3];
PetscInt cells[2];
ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);
// set domain properties
cells[0] = 4; cells[1] = 4;
lower[0] = 0; lower[1] = 0; lower[2] = 0;
upper[0] = 1; upper[1] = 1; upper[2] = 0;
// get number of cells in x and y direction
ierr = PetscOptionsGetInt(NULL, NULL, "-nx", &cells[0], &status); CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL, NULL, "-ny", &cells[1], &status); CHKERRQ(ierr);
// create the mesh
ierr = PetscPrintf(PETSC_COMM_WORLD, "Generating mesh ... "); CHKERRQ(ierr);
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE,
cells, lower, upper, NULL, interpolate, &dm); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); 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, PETSCPARTITIONERPARMETIS); 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;
}
// construct ghost cells
PetscInt nGhost; // number of ghost cells
DM dmG; // DM with ghost cells
ierr = DMPlexConstructGhostCells(dm, NULL, &nGhost, &dmG); CHKERRQ(ierr);
if(dmG)
{
ierr = DMDestroy(&dm); CHKERRQ(ierr);
dm = dmG;
}
ierr = DMSetNumFields(dm, numFields); CHKERRQ(ierr);
ierr = DMSetAdjacency(dm, 0, useCone, useClosure); 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);
ierr = SetIC(dm, U); CHKERRQ(ierr);
// update ghost cell values
ierr = VecGhostUpdateBegin(U, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(U, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
// get local vector and view it
Vec Ul;
ierr = DMGetLocalVector(dm, &Ul); CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Ul); CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Ul); CHKERRQ(ierr);
ierr = VecView(Ul, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
PetscScalar *Uarr;
ierr = VecGetArray(Ul, &Uarr); CHKERRQ(ierr);
PetscInt g, gStart, gEnd;
ierr = DMPlexGetGhostCellStratum(dm, &gStart, &gEnd); CHKERRQ(ierr);
for(g=gStart; g<gEnd; ++g)
{
PetscScalar *u;
ierr = DMPlexPointLocalRead(dm, g, Uarr, &u); CHKERRQ(ierr);
for(PetscInt i=0; i<nvar; ++i)
PetscPrintf(PETSC_COMM_WORLD, "u[%d] = %f\n", i, u[i]); CHKERRQ(ierr);
}
// cleanup
ierr = VecRestoreArray(U, &Uarr); CHKERRQ(ierr);
ierr = DMRestoreLocalVector(dm, &Ul); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr);
ierr = DMDestroy(&dm); CHKERRQ(ierr);
ierr = PetscFinalize(); CHKERRQ(ierr);
return ierr;
}
