[petsc-users] Periodic BC in petsc DMPlex / firedrake

Thu Oct 25 09:50:53 CDT 2018

Did you forget to attach "The nice picture attached" 
Stephen 

> From: "Stefano Zampini" <stefano.zampini at gmail.com>
> To: "petsc-maint" <knepley at gmail.com>
> Cc: "imilian hartig" <imilian.hartig at gmail.com>, "PETSc users list"
> <petsc-users at mcs.anl.gov>
> Sent: Thursday, October 25, 2018 4:47:43 PM
> Subject: Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

> Opened the PR [
> https://bitbucket.org/petsc/petsc/pull-requests/1203/fix-dump-vtk-field-with-periodic-meshes/diff
> |
> https://bitbucket.org/petsc/petsc/pull-requests/1203/fix-dump-vtk-field-with-periodic-meshes/diff
> ]

> Il giorno gio 25 ott 2018 alle ore 17:44 Matthew Knepley < [
> mailto:knepley at gmail.com | knepley at gmail.com ] > ha scritto:

>> Good catch Stefano.
>> Matt

>> On Thu, Oct 25, 2018 at 9:36 AM Stefano Zampini < [
>> mailto:stefano.zampini at gmail.com | stefano.zampini at gmail.com ] > wrote:

>>> Maybe this is a fix
>>> diff --git a/src/dm/impls/plex/plexvtu.c b/src/dm/impls/plex/plexvtu.c
>>> index acdea12c2f..1a8bbada6a 100644
>>> --- a/src/dm/impls/plex/plexvtu.c
>>> +++ b/src/dm/impls/plex/plexvtu.c
>>> @@ -465,10 +465,11 @@ PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer
>>> viewer)
>>> if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
>>> PetscScalar *xpoint;
>>> - ierr = DMPlexPointLocalRead(dm,v,x,&xpoint);CHKERRQ(ierr);
>>> + ierr = DMPlexPointLocalRead(dm,closure[v],x,&xpoint);CHKERRQ(ierr);
>>> y[cnt + off++] = xpoint[i];
>>> }
>>> }
>>> + cnt += off;
>>> ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize,
>>> &closure);CHKERRQ(ierr);
>>> }
>>> }

>>> Max, does this fix your problem? If you confirm, I'll fix this in the maint
>>> branch

>>> If I run the below command line with the patch and with snes tutorials ex12 I
>>> get the nice picture attached
>>> $ ./ex12 -quiet -run_type test -interpolate 1 -bc_type dirichlet
>>> -petscspace_degree 1 -vec_view vtk:test.vtu:vtk_vtu -f
>>> ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh
>>> -dm_plex_gmsh_periodic -x_periodicity periodic -y_periodicity periodic
>>> -dm_refine 4

>>> Il giorno gio 25 ott 2018 alle ore 15:11 Stefano Zampini < [
>>> mailto:stefano.zampini at gmail.com | stefano.zampini at gmail.com ] > ha scritto:

>>>> Matt,

>>>> you can reproduce it via
>>>> $ valgrind ./ex12 -quiet -run_type test -interpolate 1 -bc_type dirichlet
>>>> -petscspace_degree 1 -vec_view vtk:test.vtu:vtk_vtu -f
>>>> ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh
>>>> -dm_plex_gmsh_periodic

>>>> Long time ago I added support for viewing meshes with periodic vertices in the
>>>> VTK_VTU viewer, but I did not fix the part that writes fields

>>>> Il giorno mer 24 ott 2018 alle ore 21:04 Matthew Knepley < [
>>>> mailto:knepley at gmail.com | knepley at gmail.com ] > ha scritto:

>>>>> On Wed, Oct 24, 2018 at 11:36 AM Maximilian Hartig < [
>>>>> mailto:imilian.hartig at gmail.com | imilian.hartig at gmail.com ] > wrote:

>>>>>>> On 24. Oct 2018, at 12:49, Matthew Knepley < [ mailto:knepley at gmail.com |
>>>>>>> knepley at gmail.com ] > wrote:

>>>>>>> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell < [ mailto:wence at gmx.li |
>>>>>>> wence at gmx.li ] > wrote:

>>>>>>>> Hi Max,

>>>>>>>> (I'm cc'ing in the petsc-users mailing list which may have more advice, if you
>>>>>>>> are using PETSc you should definitely subscribe!

>>>>>>>>> On 24 Oct 2018, at 09:27, Maximilian Hartig < [ mailto:imilian.hartig at gmail.com
>>>>>>>> > | imilian.hartig at gmail.com ] > wrote:

>>>>>>>> > Hello Lawrence,

>>>>>>>>> sorry to message you out of the blue. My name is Max and I found your post on
>>>>>>>>> GitHub ( [ https://github.com/firedrakeproject/firedrake/issues/1246 |
>>>>>>>>> https://github.com/firedrakeproject/firedrake/issues/1246 ] ) on DMPlex being
>>>>>>>>> able to read periodic gmsh files. I am currently trying to do just that
>>>>>>>>> (creating a periodic DMPlex mesh with gmsh) in the context of my PhD work. So
>>>>>>>>> far I haven’ t found any documentation on the periodic BC’s with DMPlex and
>>>>>>>> > gmsh in the official petsc documentation.
>>>>>>>>> I was wondering whether you’d be so kind as to point me in a general direction
>>>>>>>>> concerning how to achieve this. You seem experienced in using petsc and I would
>>>>>>>> > greatly appreciate your help.

>>>>>>>> I think the answer is "it depends". If you're just using DMPlex directly and all
>>>>>>>> the of the functionality with PetscDS, then I /think/ that reading periodic
>>>>>>>> meshes via gmsh (assuming you're using the appropriate gmsh mesh format [v2])
>>>>>>>> "just works".

>>>>>>> There are two phases here: topological and geometric. DMPlex represents the
>>>>>>> periodic topological entity directly. For example, a circle is just a segment
>>>>>>> with one end hooked to the other. Vertices are not duplicated, or mapped to
>>>>>>> each other. This makes topology simple and easy to implement. However, then
>>>>>>> geometry is more complicated. What Plex does is allow coordinates to be
>>>>>>> represented by a discontinuous field taking values on cells, in addition to
>>>>>>> vertices. In our circle example, each cells near the cut will have 2
>>>>>>> coordinates, one for each vertex, but they will not agree across the cut. If
>>>>>>> you define a periodic domain, then Plex can construct this coordinate field
>>>>>>> automatically using DMPlexLocalize(). These DG coordinates are then used by the
>>>>>>> integration routines.

>>>>>> Ok, I think I understand the concept. DMPlex reads information about both
>>>>>> topology and coordinates from the .msh file. Creating a periodic mesh in gmsh
>>>>>> then should allow DMPlex to identify the periodic boundaries as the “cut” and
>>>>>> build the mesh topology accordingly. Coordinate information is handled
>>>>>> separately.
>>>>>> That means, as Lawrence suggested, building periodic meshes in gmsh and reading
>>>>>> them in to petsc’s DMPlex should indeed “just work”. (From the user
>>>>>> perspective). The only extra step is to call DMLocalizeCoordinates() after
>>>>>> DMPlexReadFromFile(). Sorry to reiterate, I am just trying to make sense of
>>>>>> this.

>>>>>>>> From my side, the issue is to do with mapping that coordinate field into one
>>>>>>>> that I understand (within Firedrake). You may not have this problem.

>>>>>>> Firedrake uses its own coordinate mapping and integration routines, so they must
>>>>>>> manage the second part independently. I hope to get change this slightly soon
>>>>>>> by making the Firedrake representation a DMField, so that it looks the same to
>>>>>>> Plex.

>>>>>>> Thanks,

>>>>>>> Matt

>>>>>>>> Thanks,

>>>>>>>> Lawrence
>>>>>>> --
>>>>>>> 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/ | https://www.cse.buffalo.edu/~knepley/
>>>>>>> ]

>>>>>>> To read periodic meshes from GMSH, you need to use the option
>>>>>>> -dm_plex_gmsh_periodic and DMPlexCreateFromFile

>>>>>> Ahh, thanks. I was missing the option " -dm_plex_gmsh_periodic “. But using this
>>>>>> option I now generate a segmentation fault error when calling VecView() on the
>>>>>> solution vector with vtk and hdf5 viewers. Any suggestions?

>>>>> Small example? VTK is deprecated. HDF5 should work, although it will require you
>>>>> to have proper coordinates I think. We have to
>>>>> think about what you mean. If its for a checkpoint, there is no problem, but for
>>>>> viz, those programs do not understand periodicity. Thus I embed it in a higher
>>>>> dimensional space.

>>>>> Matt

>>>>>>> See src/dm/impls/plex/examples/tests/ex1.c. An example runs

>>>>>>> $ ./ex1 -filename
>>>>>>> ${PETSC_DIR}/share/petsc/datafiles/meshes/cube_periodic_bin.msh
>>>>>>> -dm_plex_gmsh_periodic -dm_view ::ascii_info_detail -interpolate -test_shape

>>>>>>> generating periodic meshes in gmsh may be tricky, Lisandro for sure may advice.

>>>>>> Thanks,
>>>>>> Max

>>>>> --
>>>>> 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

>>>>> [ http://www.cse.buffalo.edu/~knepley/ | https://www.cse.buffalo.edu/~knepley/ ]

>>>> --
>>>> Stefano

>>> --
>>> Stefano

>> --
>> 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

>> [ http://www.cse.buffalo.edu/~knepley/ | https://www.cse.buffalo.edu/~knepley/ ]

> --
> Stefano
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181025/f682ac3c/attachment.html>