[petsc-users] Periodic BC in petsc DMPlex / firedrake
Stefano Zampini
stefano.zampini at gmail.com
Thu Oct 25 09:47:43 CDT 2018
Opened the PR
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 <knepley at gmail.com>
ha scritto:
> Good catch Stefano.
>
> Matt
>
> On Thu, Oct 25, 2018 at 9:36 AM Stefano Zampini <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 <
>> 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 <
>>> knepley at gmail.com> ha scritto:
>>>
>>>> On Wed, Oct 24, 2018 at 11:36 AM Maximilian Hartig <
>>>> imilian.hartig at gmail.com> wrote:
>>>>
>>>>>
>>>>>
>>>>> On 24. Oct 2018, at 12:49, Matthew Knepley <knepley at gmail.com> wrote:
>>>>>
>>>>> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell <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 <
>>>>>> 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 ) 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/
>>>>>
>>>>>
>>>>> 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
>>>>
>>>> https://www.cse.buffalo.edu/~knepley/
>>>> <http://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
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
--
Stefano
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181025/1c48678d/attachment-0001.html>
More information about the petsc-users
mailing list