[petsc-users] DMPlex distribution and loops over cells/faces for finite volumes

[Matthew Knepley]

On Thu, Nov 9, 2017 at 3:08 AM, Matteo Semplice <matteo.semplice at unito.it>
wrote:

> Hi.
> I am using a DMPLex to store a grid (at present 2d symplicial but will
> need to work more in general), but after distributing it, I am finding it
> hard to organize my loops over the cells/faces.
>
> First things first: the mesh distribution. I do
>
>   DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);
>   DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE);
>

This gives you adjacent cells, but not their boundary. Usually what you
want for FV.


>   DMPlexDistribute(dm, 1, NULL, &dmDist);
>
> The 1 in DMPlexDistribute is correct to give me 1 layer of ghost cells,
> right?
>

Yes.


> Using 2 instead of 1 should add more internal ghosts, right? Actually,
> this errors out with petsc3.7.7 from debian/stable... If you deem it worth,
> I will test again with petsc3.8.
>

It should not error. If possible, send me the mesh please. I have some
tests where 2 is successful.


> Secondly, looping over cells. I do
>
>   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
>   DMPlexGetHybridBounds(dm, &cPhysical, NULL, NULL, NULL);
>   for (int c = cStart; c < cEnd; ++c) {
>       if ((cPhysical>=0)&&(c>=cPhysical))
>         {code for ghost cells outside physical domain boundary}
>       else if ( ??? )
>         {code for ghost cells}
>

I am not sure why you want to draw this distinction, but you can check to
see if the cell is present in the pointSF. If it is, then its not owned by
the process.
Here is the kind of FV loop I do (for example in TS ex11):


https://bitbucket.org/petsc/petsc/src/f9fed41fd6c28d171f9c644f97b15591be9df7d6/src/snes/utils/dmplexsnes.c?at=master&fileviewer=file-view-default#dmplexsnes.c-1634


>       else
>         {code for cells owned by this process}
>   }
>
>   What should I use as a test for ghost cells outside processor boundary?
> I have seen that (on a 2d symplicial mesh) reading the ghost label with
> DMGetLabelValue(dm, "ghost", c, &ghost) gives the value -1 for inner cells
> and domain boundary ghosts but 2 for processor boundary ghosts. Is this the
> correct test? What is the meaning of ghost=2? In general, should I test for
> any value >=0?
>

Yes, see the link I gave.


> Lastly, since I will do finite volumes, I'd like to loop over the faces to
> compute fluxes.
> DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)
> gives me the start/endpoints for faces, but they include also faces
> between ghost cells that are therefore external to the cells owned by the
> processor. I see that for faces on processor boundary, ghost=-1 on one
> process and ghost=2 on the other process, which I like. However ghost=-1
> also on faces between ghosts cells and therefore to compute fluxes I should
> test for the ghost value of the face and for the ghost value of the cells
> on both sides... Is there a better way to loop over faces?
>

See the link.

 Thanks,

   Matt


> Thanks in advance for any suggestion.
>
>     Matteo
>
>


-- 
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.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171109/e3a7ab41/attachment.html>