[petsc-users] flux vector

[Matthew Knepley]

On Fri, Jun 11, 2021 at 7:38 AM Lawrence Mitchell <wence at gmx.li> wrote:

>
>
> > On 11 Jun 2021, at 11:29, Mark Adams <mfadams at lbl.gov> wrote:
> >
> > This is a Matt question, but You can set a "Boundary" label (
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexMarkBoundaryFaces.html)
> and that will label boundaries.
> > Next you want the dual of that ...
>

Before doing what Lawrence suggests, I would like to talk this through. I
don't think we should need communication. So,

  Goal: Section with dofs on each face that has support 2 (separates 2
cells)

Okay, suppose we loop over our local faces and get the support size of
each. Boundary faces are disqualified, which is good,
but if a face is on the boundary with another process, it will also have
support size 1. Thus, let us distribute with overlap = 1.
Now if we loop over that same face set, we will get the "right answer" for
the support size.

However, using overlap = 1 put in a bunch of new faces. We do not care
about the ones on the process boundary. They will be
handled by the other process. We do care about faces between two ghost
cells, since they will be a false positive. Luckily, these
are labeled by the "ghost" label. So I think this is our algorithm:

1) DMClone(dm, &fdm);

2) DMGetLocalSection(fdm, &ls);

3) Mark each face with 2 neighbors that is not a ghost

  DMGetLabel(dm, "ghost", &ghostLabel);
  DMGetHeightStratum(fdm, 1, &fStart, &fEnd);
  for (f = fStartl f < fEnd; ++f) {
    PetscInt  supportSize, ghost;

    DMGetSupportSize(fdm, f, &supportSize);
    DMLabelGetValue(ghostLabel, f, &ghost);
    if (ghost < 0 && supportSize == 2) PetscSectionSetDof(ls, f, 1);
  }

4) The global section will be created automatically when we ask for it

  DMGetGlobalSection(fdm, &s);

Let me know if this does not work.

  Thanks,

     Matt


> Note that if you want this to be the _global_ boundary you must call it
> before mesh distribution.
>
> You need to reduce the support size over the point SF, where you would
> make a contribution to the count only if the cell is in the owned part of
> your local cells.
>
> So (pseudocode):
>
> sf = DMGetPointSF(dm);
>
> // since you're only going to reduce over facet points you could create an
> embedded sf here, but meh.
>
> nroots, nleaves, ilocal, iremote = SFGetGraph(sf);
>
> cstart, cend = DMPlexGetChart(dm);
>
> count = PetscCalloc(nroots); // nroots should be == cend - cstart
>
> fstart, fend = DMPlexGetDepthStratum(dm, 1); // codim 1 entities
>
> for f = fstart; f < fend; f++ {
>   supportsize = DMPlexGetSupportSize(dm, f);
>   support = DMPlexGetSupport(dm, f);
>   for c in support {
>      if (c is owned) { // c not in ilocal
>          count[f]++;
>      }
>   }
> }
>
> gcount = PetscMalloc(nroots);
> PetscSFReduceBegin(sf, count, gcount, MPI_SUM);
> PetscSFReduceEnd(sf, count, gcount, MPI_SUM);
> PetscSFBcastBegin(sf, gcount, count, MPI_REPLACE);
> PetscSFBcastEnd(sf, gcount, count, MPI_REPLACE);
> interior_faces = DMLabelCreate("interior_faces")
> for f = fstart; f < fend; f++ {
>    if count[f] == 2 {
>        DMLabelSetValue(interior_faces, f, 1);
>    }
> }
>
>
> (Obviously) untested.
>
> Lawrence
>
>

-- 
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/20210611/0cfe1c56/attachment.html>