[petsc-users] DMPlex box mesh periodicity bug (?)

[Matthew Knepley]

On Mon, Aug 2, 2021 at 4:23 PM Thibault Bridel-Bertomeu <
thibault.bridelbertomeu at gmail.com> wrote:

> Hi Matt,
>
> Thank you for taking the time to take a look  !
>
> 1) Yea ... I know most visualization softwares and most storage formats do
> not like periodicity (a.k.a. infinity) ... but fo CFD-research it's a great
> tool so we keep trying ... x) The thing that surprised me is that my mwe
> works perfectly well in 2D, and not in 3D : that's why I ventured to call
> it a bug !
>
> 2) Hmm I tried adding a "DMLocalizeCoordinates" after the
> "DMPlexConstructGhostCells" but it does not change anything, the result is
> still exactly the same as I showed above.
>

Hmm, I will try it. It should error telling you that VTK cannot handle
meshes with localized coordinates, so something is going wrong there.


> As for your code, it compiles fine and I can execute it as well, thanks !
> However ...
> 1) when I run it in 2D and ask for the HDF5 output, and then run the
> petsc_gen_xdmf.py script, I get the attached result (sol_2D.png) : it seems
> either Paraview cannot handle what's in the HDF file or the HDF file does
> not contain something consistent,
>

This is exactly what I get. Something is wrong with the periodic cut when I
have a double point. If you use

  -dm_plex_box_bd periodic,none

it looks fine. This is the bug with periodic visualization that I was
talking about.


> 2) when I run it in 3D and ask for the HDF5 output + run petsc_gen_xdmf.py
> I can see ... nothing : paraview does not show anything when I try to open
> the XDMF file, although it seems syntaxically correct (the HDF file as well)
>

Yes, Paraview misunderstands the mesh connections, so I must be telling it
something wrong with the periodic cut. You can see that each periodic
dimension
vanishes from the display, but all the questions I ask it in tests are
correct.


> 3) when I run it in 3D and ask for VTU output, I get exactly the same
> thing as what I got with my F90 program (see sol_3D.png attached) - my
> command line is ./thibault -dm_plex_dim 3 -dm_plex_simplex 0
> -dm_plex_box_faces 16,16,16 -dm_plex_box_lower -5,-5,-5 -dm_plex_box_upper
> 5,5,5 -dm_plex_box_bd periodic,periodic,periodic
> -dm_plex_create_fv_ghost_cells -dm_plex_periodic_cut -vec_view vtk:sol.vtu
>

Yes, VTU output does not work at all with periodicity. It was only really
intended to work with DMDA. I do everything with HDF5 now since it can
handle
multiple meshes, etc.

I think the right thing to do is fix the "periodic cut" support for
multiply periodic things, and for 3D.

  Thanks,

     Matt


> Does your piece of code yields different results on your end ?
>
> (P.S. I am using the main branch, commit id ae6adb75dd).
>
> Thank you very much for your support !!
>
> Thibault
>
>
>
> Le lun. 2 août 2021 à 15:48, Matthew Knepley <knepley at gmail.com> a écrit :
>
>> On Sat, Jul 31, 2021 at 6:01 AM Thibault Bridel-Bertomeu <
>> thibault.bridelbertomeu at gmail.com> wrote:
>>
>>> Dear all,
>>>
>>> I have noticed what I think is a bug with a 3D DMPlex box mesh with
>>> periodic boundaries.
>>> When I project a function onto it, it behaves as if the last row of
>>> cells in X and in Y direction do not have the right coordinates.
>>>
>>> I attach to this email a minimal example that reproduces the bug (files
>>> mwe_periodic_3d.F90, wrapper_petsc.c, wrapper_petsc.h90, makefile), as well
>>> as the output of this code (initmesh.vtu for an output of the DM,
>>> solution.vtu for an output of the data projected onto the mesh). There is
>>> also a screenshot of what's going on.
>>> If one considers the function I project onto the mesh, what should
>>> normally happen is that there is a "hole" is the density field around the
>>> x=0, y=0 region, the rest being equal to one.
>>>
>>> I hope it is just a mishandling from my end !!
>>>
>>
>> Hi Thibault,
>>
>> 1) The main thing happening here is that visualization has some problems
>> for completely periodic things. It is on my list to fix, but below other
>> things.
>>
>> 2) Second, it looks like you need to localize coordinates again after
>> creating ghost cells. Everything is in the right order if you use the
>> command line to create
>>     the mesh. I have attached a C example where I do this. In your code,
>> I think just adding Localize again will work.
>>
>> In my example, you can see 2D as well as 3D, and change from non-periodic
>> to periodic to see what is happening.
>>
>>   Thanks,
>>
>>      Matt
>>
>>
>>> Thank you in advance for your help,
>>>
>>> Cheers,
>>>
>>> Thibault
>>>
>>>
>>
>> --
>> 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/20210802/1fe325aa/attachment.html>