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

[Thibault Bridel-Bertomeu]

Hi Matthew,

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

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 ?
Also, does it still create a « ghost » label even if it’s empty ? Do you
think that will hold for later uses in the GetFaceFields and similar
routines ?
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.

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210518/584af4b8/attachment-0001.html>