[petsc-users] DMPlex from periodic Gmsh and ghost cells
Matthew Knepley
knepley at gmail.com
Wed May 19 06:32:50 CDT 2021
On Wed, May 19, 2021 at 4:36 AM Thibault Bridel-Bertomeu <
thibault.bridelbertomeu at gmail.com> wrote:
> Hello Matt,
>
> Le mar. 18 mai 2021 à 20:31, Matthew Knepley <knepley at gmail.com> a écrit :
>
>> 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.
>>
>
> OK thanks ! I had to write the F90 wrapper for DMLabelCreate but aside
> from that, everything indeed goes through fine with the ConstructGhostCells.
>
>
>>
>>
>>> 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.
>>
>
> Actually, there is an issue with the corner cells. My mesh is periodic in
> both directions and it turns out the corner cells accumulate way too much
> numerical flux and the computation crashes after 4 iterations or so - it
> happens both with quads and tris. Whether I use my version of the RHS
> computing routine or the DMPlexTSComputeRHSFunctionFVM leads to the same
> crash.
> I'll keep investigating, but if you have any idea what could go wrong ...
>
Hmm, I have not run a fully periodic case. I guess I could set one up.
However, it must be the case that we are conservative, since we first find
a flux on every edge, and then update cells. So the corner cells should not
be able to accumulate more than the total mass/4. I would first check at
each step that total mass is conserved. If not, we can track that down.
Thanks,
Matt
> Thanks !!
>
> Thibault
>
>
>>
>> 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/>
>>
>
--
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/20210519/21d2c3de/attachment-0001.html>
More information about the petsc-users
mailing list