Hi,
I am writing a simple code to solve linear advection equation using a first
order cell centered finite volume scheme on 2D unstructured girds using
DMPlex. In order to set values to a Global vector (for example, to set the
initial value of the solution vector), I am looping over the cells owned by
a process (including partition ghost cells) and checking the "vtk" label
for each cell, assigned by the DMPlexConstructGhostCells() function, to
prevent it from writing to ghost cells. If I don't do this check, the code
gives a segmentation fault, which, as far as I understand, is caused by
trying to write into the ghost points which do not exist on the Global
Vector. Following is the function that I have written to set the initial
condition:
PetscErrorCode SetIC(DM dm, Vec U)
{
PetscErrorCode ierr;
PetscScalar *u;
PetscFunctionBegin;
PetscInt c, cStart, cEnd; // cells
PetscReal area, centroid[3], normal[3]; // geometric data
// get cell stratum owned by processor
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
// get array for U
ierr = VecGetArray(U, &u);
// 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)
{
PetscReal X[2]; // cell centroid
ierr = DMPlexComputeCellGeometryFVM(dm, c, &area, centroid, normal); CHKERRQ
(ierr);
X[0] = centroid[0]; X[1] = centroid[1];
u[c] = initial_condition(X);
}
}
ierr = VecRestoreArray(U, &u);
PetscFunctionReturn(0);
}
This gives me the desired output, but, I wanted to ask if there is a better
and more efficient way to achieve this, i.e. to write to a global vector,
than checking the label for each cell. I am also attaching a sample code
which sets the initial condition and writes a corresponding vtk file.
Kindly correct me if I am wrong in my understanding and give your
suggestions to improve upon this.
Regards,
Shashwat
static char help[] = "Linear advection using first order FVM\n";
#include <petsc.h>
// initial condition
PetscReal initial_condition(const PetscReal *X)
{
PetscReal x = X[0];
PetscReal y = X[1];
PetscReal r2 = PetscPowReal(x-0.5,2) + PetscPowReal(y,2);
return PetscExpReal(-50.0*r2);
}
// set initial condition for the given dm and vector
PetscErrorCode SetIC(DM dm, Vec U)
{
PetscErrorCode ierr;
PetscScalar *u;
PetscFunctionBegin;
PetscInt c, cStart, cEnd; // cells
PetscReal area, centroid[3], normal[3]; // geometric data
// get cell stratum owned by processor
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
// get array for U
ierr = VecGetArray(U, &u);
// 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)
{
PetscReal X[2]; // cell centroid
ierr = DMPlexComputeCellGeometryFVM(dm, c, &area, centroid, normal); CHKERRQ(ierr);
X[0] = centroid[0]; X[1] = centroid[1];
u[c] = initial_condition(X);
}
}
ierr = VecRestoreArray(U, &u);
PetscFunctionReturn(0);
}
static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
{
PetscFunctionBegin;
PetscErrorCode ierr;
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, filename);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
int main(int argc, char **argv)
{
PetscErrorCode ierr;
PetscViewer viewer;
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];
PetscInt bcField[1]; // number of BCs for each field
IS bcPointsIS[1]; // contains boundary faces
PetscBool interpolate = PETSC_TRUE, status;
Vec U; // solution vector
PetscMPIInt rank, size;
char filename[PETSC_MAX_PATH_LEN];
ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);
// get mesh filename
ierr = PetscOptionsGetString(NULL, NULL, "-mesh", filename, PETSC_MAX_PATH_LEN, &status);
if(status) // gmsh file provided by user
{
char file[PETSC_MAX_PATH_LEN];
ierr = PetscStrcpy(file, filename); CHKERRQ(ierr);
ierr = PetscSNPrintf(filename, sizeof filename,"./%s", file); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Reading gmsh %s ... ", file); CHKERRQ(ierr);
ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, filename, interpolate, &dm); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); CHKERRQ(ierr);
}
else
{
// ask user for the mesh file
ierr = PetscPrintf(PETSC_COMM_WORLD, "Please provide a gmsh file with using: -mesh <filename>\n"); CHKERRQ(ierr);
return 0;
}
// distribute mesh over processes;
ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist); CHKERRQ(ierr);
if(dmDist)
{
ierr = DMDestroy(&dm); CHKERRQ(ierr);
dm = dmDist;
}
// print mesh information
ierr = PetscPrintf(PETSC_COMM_WORLD, "overlap: %d, "
"distributed among %d processors\n",
overlap, size); CHKERRQ(ierr);
// 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;
}
// create scalar field for solution u
numComp[0] = 1;
for(i=0; i<numFields*(dim+1); ++i) numDof[i] = 0;
// define u at cell centers
numDof[0*(dim+1)+dim] = 1;
// set BC
numBC = 1;
bcField[0] = 0; // Dirichlet BC on u at boundary
ierr = DMGetStratumIS(dm, "Face Sets", 1, &bcPointsIS[0]); CHKERRQ(ierr);
// create a PetscSection with this layout
ierr = DMSetNumFields(dm, numFields); CHKERRQ(ierr);
ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC,
bcField, NULL, bcPointsIS, NULL, &sec); CHKERRQ(ierr);
ierr = ISDestroy(&bcPointsIS[0]); CHKERRQ(ierr);
ierr = PetscSectionSetFieldName(sec, 0, "u"); CHKERRQ(ierr);
// set the section for dm
ierr = DMSetLocalSection(dm, sec); CHKERRQ(ierr);
ierr = DMSetAdjacency(dm, 0, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr);
ierr = DMSetUp(dm); CHKERRQ(ierr);
ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);
// get global vector for solution
ierr = DMCreateGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = SetIC(dm, U); CHKERRQ(ierr);
// write initial solution to vtk
ierr = OutputVTK(dm, "sol-000.vtu", &viewer); 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);
ierr = DMRestoreGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr);
ierr = DMDestroy(&dm); CHKERRQ(ierr);
ierr = PetscFinalize(); CHKERRQ(ierr);
return ierr;
}
