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

[Matthew Knepley]

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.

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