[petsc-users] Setting values to a Global Vector using Labels

[Shashwat Tiwari]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200502/808a9269/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.c
Type: text/x-csrc
Size: 5707 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200502/808a9269/attachment.bin>