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

[Matthew Knepley]

On Sat, May 2, 2020 at 3:36 AM Shashwat Tiwari <shaswat121994 at gmail.com>
wrote:

> 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.
>

1) This will work, but I think the intent is to use the "ghost" label to
identify ghost cells, It has the same information as "vtk" since you do not
want to output those
     cells for visualization either, but it might be more clear in the code.

2) If I were doing this many times, I would create a new IS with the cells
I was looping over, namely [cStart, cEnd) - {ghost cells}. That way you
have only a lookup rather
    than a search. However, if you do it only once or twice, I don't think
there is much advantage.

  Thanks,

    Matt


> 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200502/17b05ff6/attachment-0001.html>