[petsc-users] Periodic boundary conditions in DMPlex

[Matthew Knepley]

On Mon, Oct 18, 2021 at 3:27 AM Pierre Seize <pierre.seize at onera.fr> wrote:

> On 15/10/21 22:30, Matthew Knepley wrote:
>
> On Fri, Oct 15, 2021 at 10:16 AM Pierre Seize <pierre.seize at onera.fr>
> wrote:
>
>> I read everything again, I think I did not understand you at first. The
>> first solution is to modify the DAG, so that the rightmost cell is linked
>> to the leftmost face, right ? To do that, do I have to manually edit the
>> DAG (the mesh is read from a file) ?
>>
> Yes, the DAG would be modified if you want it for some particular mesh
> that we cannot read automatically. For example, we can read periodic GMsh
> meshes.
>
>> If so, the mesh connectivity is like the one of a torus, then how does it
>> work with the cells/faces coordinates ?
>>
> You let the coordinate field be in a DG space, so that it can have jumps.
>
> I'm not sure I fully understand this, what I'll do is that I will
> experiment with a periodic GMSH mesh.
>

We normally assume that the coordinate field on a mesh is continuous, which
is why we associate it with vertices. However
you could define the field on cells, representing it by discontinuous
polynomials of degree 1.

> Now I think the second method may be more straightforward. What's the idea
>> ? Get the mapping with DMGetLocalToGlobalMapping, then create the mapping
>> corresponding to the periodicity with ISLocalToGlobalMappingCreate, and
>> finally ISLocalToGlobalMappingConcatenate ? I'm not sure this is the way,
>> and I did not find something like DMSetLocalToGlobalMapping to restore the
>> modified mapping.
>>
> It is more complicated. We make the LocalToGlobalMap by looking at the
> PetscSection (essentially if gives function space information)  and
> deciding which unknowns are removed from the global space.
> You would need to decide that unknowns constrained by periodicity are not
> present in the global space. Actually, this is not hard. You just mark them
> as constrained in the PetscSection, and all the layout
> functions will function correctly. However, then the LocalToGlobalMap will
> not be exactly right because the constrained unknowns will not be filled in
> (just like Dirichlet conditions). You would augment the
> map so that it fills those in by looking up their periodic counterparts.
> Jed has argued for this type of periodicity.
>
> I also don't understand. One one side of my mesh I have :   | ghost1 |
> cell1 | cell2 | ...  and on the other | cell_n-1 | cell_n | ghost_n |. Are
> not the ghosts (from DMPlexConstructGhostCells) already constrained ?
> I experimented on that too, I did:
>
> DMGetSectionSF
> PetscSFGetGraph
> augment the graph by adding the local ghost cells to the leaves, and the
> correct remote "true" cells to the roots
> PetscSFSetGraph
>
> and it seems to do what I want. Is this what you meant ? Is this a correct
> way to use the PETSc objects ? Or is this just hacky and I'm lucky it works
> ?
>

Yes, this will do what you want for field values, but not for coordinates.
This is exactly what you would get if you just connected the topology.

  Thanks,

      Matt


> Pierre
>
>
> To me, the first kind is much more straightforward, but maybe this is
> because I find the topology code more clear.
>
>   Thanks,
>
>       Matt
>
>> Pierre
>>
>> On 15/10/21 15:33, Pierre Seize wrote:
>>
>> When I first tried to handle the periodicity, I found the
>> DMPlexCreateBoxMesh function (I cannot find the cylinder one).
>>
>> From reading the sources, I understand that we do some work either in
>> DMPlexCreateCubeMesh_Internal or with DMSetPeriodicity.
>>
>> I tried to use DMSetPeriodicity before, for example with a 2x2 box on
>> length 10. I did something like:
>> const PetscReal maxCell[] = {2, 2};
>> const PetscReal L[] = {10, 10};
>> const DMBoundaryType bd[] = {DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC};
>> DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bd);
>> // or:
>> DMSetPeriodicity(dm, PETSC_TRUE, NULL, L, bd);
>>
>> but it did not work:
>> VecSet(X, 1);
>> DMGetLocalVector(dm, &locX);
>> VecZeroEntries(locX);
>> DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
>> DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
>> VecView(locX, PETSC_VIEWER_STDOUT_WORLD);
>>
>> but the ghost cells values are all 0, only the real cells are 1. So I
>> guess DMSetPeriodicity alone is not sufficient to handle the
>> periodicity. Is there a way to do what I want ? That is set up my DMPlex in
>> a way that  DMGlobalToLocalBegin/DMGlobalToLocalEnd do exchange values
>> between procs AND exchange the periodic values?
>>
>>
>> Thanks for the help
>>
>>
>> Pierre
>>
>> On 15/10/21 14:03, Matthew Knepley wrote:
>>
>> On Fri, Oct 15, 2021 at 7:31 AM Pierre Seize <pierre.seize at onera.fr>
>> wrote:
>>
>>> It makes sense, thank you. In fact, both ways seems better than my way.
>>> The first one looks the most straightforward. Unfortunately I do not know
>>> how to implement either of them. Could you please direct me to the
>>> corresponding PETSc functions ?
>>>
>> The first way is implemented for example in DMPlexCreateBoxMesh() and
>> DMPlexCreateCylinderMesh(). The second is not implemented since
>> there did not seem to be a general way to do it. I would help if you
>> wanted to try coding it up.
>>
>>   Thanks,
>>
>>     Matt
>>
>>> Pierre
>>>
>>> On 15/10/21 13:25, Matthew Knepley wrote:
>>>
>>> On Fri, Oct 15, 2021 at 7:08 AM Pierre Seize <pierre.seize at onera.fr>
>>> wrote:
>>>
>>>> Hi,
>>>>
>>>> I'm writing a code using PETSc to solve NS equations with FV on an
>>>> unstructured mesh. Therefore I use DMPlex.
>>>>
>>>> Regarding periodicity, I manage to implement it this way:
>>>>
>>>>    - for each couple of boundaries that is linked with periodicity, I
>>>> create a buffer vector with an ISLocalToGlobalMapping
>>>>
>>>>    - then, when I need to fill the ghost cells corresponding to the
>>>> periodicity, the i "true" cell of the local vector fills the buffer
>>>> vector on location i with VecSetValuesBlockedLocal, then
>>>> VecAssemblyBegin/VecAssemblyEnd ensure each value is send to the
>>>> correct
>>>> location thanks to the mapping, then the i "ghost" cell of the local
>>>> vector reads the vector on location i to get it's value.
>>>>
>>>>
>>>> It works, but it seems to me there is a better way, with maybe PetscSF,
>>>> VecScatter, or something I don't know yet. Does anyone have any advice ?
>>>>
>>>
>>> There are at least two other ways to handle this. First, the method that
>>> is advocated in
>>> Plex is to actually make a periodic geometry, meaning connect the cells
>>> that are meant
>>> to be connected. Then, if you partition with overlap = 1,
>>> PetscGlobalToLocal() will fill in
>>> these cell values automatically.
>>>
>>> Second, you could use a non-periodic geometry, but alter the
>>> LocalToGlobal map such
>>> that the cells gets filled in anyway. Many codes use this scheme and it
>>> is straightforward
>>> with Plex just by augmenting the map it makes automatically.
>>>
>>> Does this make sense?
>>>
>>>   Thanks,
>>>
>>>      Matt
>>>
>>>
>>>> Pierre Seize
>>>>
>>> --
>>> 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/%7Eknepley/>
>>>
>>>
>>>
>>
>> --
>> 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/%7Eknepley/>
>>
>>
>>
>>
>
> --
> 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/%7Eknepley/>
>
>
>

-- 
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/20211018/ed679571/attachment.html>