[petsc-users] DMPlex from periodic Gmsh and ghost cells

Tue May 18 13:30:58 CDT 2021

On Tue, May 18, 2021 at 12:16 PM Thibault Bridel-Bertomeu <
thibault.bridelbertomeu at gmail.com> wrote:

> Hi Matthew,
>
> Thank you very much for your quick answer, as always !
>

Cool.

> Le mar. 18 mai 2021 à 17:46, Matthew Knepley <knepley at gmail.com> a écrit :
>
>> On Tue, May 18, 2021 at 5:19 AM Thibault Bridel-Bertomeu <
>> thibault.bridelbertomeu at gmail.com> wrote:
>>
>>> Dear all,
>>>
>>> I am playing around with creating DMPlex from a periodic Gmsh mesh
>>> (still in the same finite volume code that solves the Euler equations)
>>> because I want to run some classical periodic test cases from the
>>> literature (like convecting a vortex and so on).
>>> I got from the mailing list and the examples that I gotta use
>>> -dm_plex_gmsh_periodic for the reader and the DM creator to do their job
>>> properly and read the $Periodic section from the Gmsh mesh. That works
>>> fine, I can write the DM in HDF5 format afterwards and the periodicity is
>>> clearly there (i can see in Paraview that opposite sides are "attached" to
>>> each other).
>>> Now, in the context of the finite volume solver, I wrote my one
>>> right-hand side term computing routine, and I rely heavily on routines such
>>> as DMPlexGetFaceFields. Those routines, in turn rely on the ghost cells and
>>> ... I am having an issue with DMPlexConstructGhostCells. When I use
>>> -dm_plex_gmsh_periodic, DMPlexCreateFromFile works fine,
>>> DMSetBasicAdjacency(dm, true, false) works fine, DMPlexDistribute works
>>> fine, but when I get to DMPlexConstructGhostCells, I get
>>>
>>>  "DM has boundary face 2048 with 2 support cells"
>>>
>>> I understand the DM begin periodic, even the "boundary" faces have two
>>> neighbors because of the topological folding, so ... yeah, the ghost cells
>>> cannot really exist. But then how can I proceed ? Because I need that
>>> "ghost" label for all the routines I call from the RHS computing routine ...
>>>
>>
>> Thanks for the example. I have run this myself.  It looks to me like the
>> mesh you are providing. Has no boundary. It is doubly periodic (a torus).
>> When I tell ConstructGhostCells() to
>> create a label for the boundary, rather than using "Face Sets", it
>> creates an empty label because no faces have only 1 neighbor.
>>
>
> Yes indeed, up is down and right is left for this mesh, the idea being of
> emulating an infinite medium.
> So you manage to get through ConstructGhostCells without a segfault ? What
> did you specify instead of PETSC_NULL_CHARACTER ?
>

The NULL tells it to use the label "Face Sets", but that label coming from
Gmsh marks edges that are not actually boundaries. I asked PETSc
to create a boundary label named "marker" using
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexMarkBoundaryFaces.html
It was empty, but everything goes through fine.

> Also, does it still create a « ghost » label even if it’s empty ?
>

Yes.

> Do you think that will hold for later uses in the GetFaceFields and
> similar routines ?
>

Yes.

> I read through the source and it seems it won’t like it much if there are
> no ghost at all, but since I did not pass the ConstructGhostCell step I’m
> not sure.
>

It should work fine, you just will not need boundary conditions.

  Thanks,

    Matt

> Thanks !
>
> Thibault
>
>
>>   Thanks,
>>
>>      Matt
>>
>>
>>> I attach to this email the periodic Gmsh mesh I am playing with, and the
>>> routine from my code where I create the DM and so on:
>>>
>>> subroutine initmesh
>>>
>>> PetscErrorCode :: ierr
>>> DM :: dmDist, dmGhost
>>> integer :: overlap
>>> PetscViewer :: vtkViewer
>>>
>>> call PetscLogEventBegin(timer_initmesh, ierr); CHKERRA(ierr)
>>>
>>> ! Number of neighbours taken into account in MP communications(1 - Order
>>> 1; 2 - Order 2)
>>> overlap = 1
>>>
>>> call PetscPrintf(PETSC_COMM_WORLD, "Initializing mesh...\n", ierr) ;
>>> CHKERRA(ierr)
>>>
>>> ! Force DMPlex to use gmsh marker
>>> ! call PetscOptionsSetValue(PETSC_NULL_OPTIONS,
>>> "-dm_plex_gmsh_use_marker", "true", ierr); CHKERRA(ierr)
>>>
>>> ! Read mesh from file name 'meshname'
>>> call DMPlexCreateFromFile(PETSC_COMM_WORLD, meshname, PETSC_TRUE, dm,
>>> ierr); CHKERRA(ierr)
>>>
>>> ! Distribute on processors
>>> ! Start with connectivity
>>> call DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE, ierr) ;
>>> CHKERRA(ierr)
>>>
>>> ! Distribute on processors
>>> call DMPlexDistribute(dm, overlap, PETSC_NULL_SF, dmDist, ierr) ;
>>> CHKERRA(ierr)
>>>
>>> ! Security check
>>> if (dmDist /= PETSC_NULL_DM) then
>>> ! Destroy previous dm
>>> call DMDestroy(dm, ierr) ; CHKERRA(ierr)
>>> ! Replace with dmDist
>>> dm = dmDist
>>> end if
>>>
>>> ! Finalize setup of the object
>>> call DMSetFromOptions(dm, ierr) ; CHKERRA(ierr)
>>>
>>> ! Boundary condition with ghost cells
>>> call DMPlexConstructGhostCells(dm, PETSC_NULL_CHARACTER,
>>> PETSC_NULL_INTEGER, dmGhost, ierr); CHKERRA(ierr)
>>>
>>> ! Security check
>>> if (dmGhost /= PETSC_NULL_DM) then
>>> ! Destroy previous dm
>>> call DMDestroy(dm, ierr) ; CHKERRA(ierr)
>>> ! Replace with dmGhost
>>> dm = dmGhost
>>> end if
>>>
>>> if (debug) then
>>> ! Show in terminal
>>> call PetscPrintf(PETSC_COMM_WORLD, ":: [DEBUG] Visualizing DM in
>>> console ::\n", ierr); CHKERRA(ierr)
>>> call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr) ; CHKERRA(ierr)
>>> ! VTK viewer
>>> call PetscViewerCreate(PETSC_COMM_WORLD, vtkViewer, ierr) ;
>>> CHKERRA(ierr)
>>> call PetscViewerSetType(vtkViewer, PETSCVIEWERHDF5, ierr) ;
>>> CHKERRA(ierr)
>>> call PetscViewerFileSetMode(vtkViewer, FILE_MODE_WRITE, ierr) ;
>>> CHKERRA(ierr)
>>> call PetscViewerFileSetName(vtkViewer, "debug_initmesh.h5", ierr) ;
>>> CHKERRA(ierr)
>>> call DMView(dm, vtkViewer, ierr) ; CHKERRA(ierr)
>>> call PetscViewerDestroy(vtkViewer, ierr) ; CHKERRA(ierr)
>>> end if
>>>
>>> call PetscPrintf(PETSC_COMM_WORLD, "Done !\n", ierr) ; CHKERRA(ierr)
>>>
>>> call PetscLogEventEnd(timer_initmesh, ierr); CHKERRA(ierr)
>>>
>>> end subroutine initmesh
>>>
>>> Thank you very much for your help !
>>>
>>> Best regards,
>>>
>>> Thibault Bridel-Bertomeu
>>>
>>
>>
>> --
>> 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/>
>>
> --
> Thibault Bridel-Bertomeu
> —
> Eng, MSc, PhD
> Research Engineer
> CEA/CESTA
> 33114 LE BARP
> Tel.: (+33)557046924
> Mob.: (+33)611025322
> Mail: thibault.bridelbertomeu at gmail.com
>

-- 
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/20210518/2413542b/attachment.html>