On Tue, May 18, 2021 at 5:19 AM Thibault Bridel-Bertomeu < [email protected]> 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/>