[petsc-users] flux vector

[Matthew Knepley]

On Sat, Jun 12, 2021 at 5:25 AM Lawrence Mitchell <wence at gmx.li> wrote:

>
> > On 11 Jun 2021, at 13:19, Matthew Knepley <knepley at gmail.com> wrote:
> >
> > 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:
>
> Ah, I missed that the mesh is already overlapped.
>
> > 1) DMClone(dm, &fdm);
> >
> > 2) DMGetLocalSection(fdm, &ls);
> >
> > 3) Mark each face with 2 neighbors that is not a ghost
> >
> >   DMGetLabel(dm, "ghost", &ghostLabel);
>
> This label is only created if one calls DMPlexConstructGhostCells,


Shoot, that is the wrong workflow.


> I don't understand this code really, but it seems like if I call
> DMPlexDistribute(dm, overlap=1, &newdm) and then DMPlexConstructGhostCells,
> the "depth 2" faces/etc... will be marked with the ghost label, but the
> "depth 1" faces will not be.
>

Yes, I only marked those because those are the only ones you might
calculate fluxes on.


> That said, I agree with Matt's algorithm, and if the ghost label isn't
> available or set up right you can replace
>
> DMLabelGetValue(ghostLabel, f, &ghost)
>
> by checking if the point f is in the ilocal array of the DM's pointSF, I
> think.
>

It is a little more complicated. The check for "ghost" is that you have 2
cells in the support, but they are both not owned, meaning both are in the
SF.

  Thanks,

    Matt


> >   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.
>
>
> 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/20210612/4c6d83bc/attachment.html>