[petsc-users] DMPlex Parallel Output and DMPlexCreateSection Crashes when DMPlexCreateGmshFromFile used
Mike Michell
mi.mike1021 at gmail.com
Fri Apr 29 11:40:02 CDT 2022
Thanks for helpful idea. For "median-dual" control volume, which is shaped
by connecting {cell centroids - face centroids - edge midpoints}, I
realized below:
(a) Verticies in original DMPlex (read from a mesh file) will be a list of
cells in a reconstructed median-dual DMPlex.
(b) Cells, and interpolated edges and faces in the original DMPlex will be
a list of vertices in the median-dual DMPlex.
(c) Once the median-dual DMPlex built by (a) and (b), new faces and edges
can be interpolated.
Is there a smart way to create new DMPlex by doing (a) and (b)? I was
thinking about writing the original topologies into a dat file then using
DMPlexCreateCellVertexFromFile() to create new plex, but it looks not an
efficient way. How can I say to new DMPlex to mapping all the original
vertices into the new cell list in DAG? Any recommendations?
Thanks,
Mike
2022년 4월 29일 (금) 오전 7:51, Matthew Knepley <knepley at gmail.com>님이 작성:
> On Fri, Apr 29, 2022 at 8:27 AM Mike Michell <mi.mike1021 at gmail.com>
> wrote:
>
>>
>> Thanks for the answers and I agree. Creating dual mesh when the code
>> starts and uses that dm will be the easiest way.
>> But it is confusing how to achieve that. The entire DAG of the original
>> mesh should change, except for vertices. How can PETSc rebuild DAG for
>> cells and edges (and faces for 3D)?
>>
>
> In this paper (https://arxiv.org/pdf/0908.4427.pdf), Section 3.1, we show
> that the dual is just the original DAG with the arrows
> reversed, so it is trivial to construct. There are two issues:
>
> 1) You have to compute the geometry once you have the topology, but you
> already do that
>
> 2) In parallel you will produce a mesh with cell overlap, but I think
> this is normal for dual methods. We have special methods
> to get around this (we use them in mesh partitioning), which we can
> use to parallelize this step.
>
> Thanks,
>
> Matt
>
>
>> In 2D for example, I need to reconstruct median-dual edges and cells by
>> connecting the cell-centroids and edge-centers in the original mesh.
>>
>> Thanks,
>> Mike
>>
>> On Wed, Apr 27, 2022 at 11:31 PM Mike Michell <mi.mike1021 at gmail.com>
>>> wrote:
>>>
>>>> Thanks for the answers. Major reason that not to store those dual cells
>>>> is to share the same grid file with other codes that are required for
>>>> coupled physics modeling.
>>>>
>>>
>>> Yes, I would create the dual on startup and throw away the original mesh.
>>>
>>>
>>>> As a follow-up question, there is struggling to use PetscSection. I
>>>> attached my example code.
>>>> The PETSc-provided example, "ex1f90.F90" uses DMPlexCreateSection() to
>>>> make a PetscSection object.
>>>> However, by reading .msh & creating a DMPlex, none of label is working
>>>> with that, except for 'marker' label, which should be provided through
>>>> "-dm_plex_gmsh_use_marker" option in case gmsh file is used.
>>>>
>>>
>>> I have replaced your startup code with the canonical code we use in C
>>> examples. It should be more flexible and robust.
>>> Here is the concise output for your mesh:
>>>
>>> master *:~/Downloads/tmp/Mitchell/Question_3$ ./test -dm_plex_filename
>>> new.msh -dm_view
>>> DM Object: 1 MPI processes
>>> type: plex
>>> DM_0x84000000_1 in 2 dimensions:
>>> Number of 0-cells per rank: 44
>>> Number of 1-cells per rank: 109
>>> Number of 2-cells per rank: 66
>>> Labels:
>>> celltype: 3 strata with value/size (0 (44), 3 (66), 1 (109))
>>> depth: 3 strata with value/size (0 (44), 1 (109), 2 (66))
>>> Face Sets: 4 strata with value/size (5 (5), 7 (5), 6 (5), 8 (5))
>>>
>>> The original code was built to take GMsh markers to canonical names
>>> because users wanted that. However, if you want your names preserved, you
>>> can use
>>>
>>> master *:~/Downloads/tmp/Mitchell/Question_3$ ./test -dm_plex_filename
>>> new.msh -dm_view -dm_plex_gmsh_use_regions
>>> MPI information ... 0 1
>>> DM Object: 1 MPI processes
>>> type: plex
>>> DM_0x84000000_1 in 2 dimensions:
>>> Number of 0-cells per rank: 44
>>> Number of 1-cells per rank: 109
>>> Number of 2-cells per rank: 66
>>> Labels:
>>> celltype: 3 strata with value/size (0 (44), 3 (66), 1 (109))
>>> depth: 3 strata with value/size (0 (44), 1 (109), 2 (66))
>>> inlet: 1 strata with value/size (5 (5))
>>> upper: 1 strata with value/size (7 (5))
>>> outlet: 1 strata with value/size (6 (5))
>>> lower: 1 strata with value/size (8 (5))
>>>
>>> There is also code to mark vertices, which some users do not want, that
>>> you can turn on using -dm_plex_gmsh_mark_vertices. I use both of these
>>> options with PyLith.
>>>
>>>
>>>> In my .msh file, there are several boundary conditions, and if I
>>>> declare labels using them and try to give them to DMGetStratumIS(), the
>>>> code crashed. It only works with 'marker' label. Any comments?
>>>> More importantly, trying to use DMPlexCreateSection() is quite painful.
>>>>
>>>
>>> What is not working for you?
>>>
>>>
>>>> Therefore, I tried to create PetscSection object through
>>>> PetscSectionCreate() as attached. But declared vectors are not printed out
>>>> if I use this way.
>>>>
>>>
>>> If you used the default viewer (PETSC_VIEWER_STDOUT_WORLD), you can see
>>> them printed. You are using VTK with a Section that put values on all mesh
>>> points. VTK only understands cell fields and vertex fields, so the code
>>> does not know how to output this vector. If you define values only on cells
>>> or vertices, you will see the output. I believe.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>
>>>> I guess the section object is not properly linked to dm object.
>>>> Basically, my vector objects (x, y, vel) are not seen from dm viewer &
>>>> relevant output file. Do you have any recommendations?
>>>>
>>>> Thanks,
>>>> Mike
>>>>
>>>>
>>>> On Tue, Apr 26, 2022 at 7:27 PM Mike Michell <mi.mike1021 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Thanks for the answers. One last question related to designing my
>>>>>> code.
>>>>>> What I want to do is to build a finite-volume code discretized on
>>>>>> "median dual" unstructured mesh. So basically it will use a vertex-centered
>>>>>> discretization scheme by looping over "median dual faces", which is not
>>>>>> physically existing in my mesh file.
>>>>>>
>>>>>> I can acheive this by allocating some vectors to compute required
>>>>>> geometrical metrics (e.g., face normal vectors, face area, etc. that
>>>>>> required to integrate over median dual control volume) as I did before even
>>>>>> without PETSc. However, I think this definitely cannot be an optimal way,
>>>>>> and I believe there should be much larger benefit that can be obtained, if
>>>>>> I well design my DMPlex's data structure (e.g., PetscSection).
>>>>>> Is there any example with this median dual control volume case? or
>>>>>> Any guideline will be helpful.
>>>>>> I was thinking to borrow some functions designed for finite element.
>>>>>>
>>>>>
>>>>> I have never done a dual scheme before, so let me ask some questions
>>>>> to better understand the underlying ideas.
>>>>>
>>>>> 1) Is there a reason not to store the dual directly? That seems like
>>>>> the easiest approach.
>>>>>
>>>>> 2) You can use Plex to layout auxiliary geometric data. I do it the
>>>>> following way:
>>>>>
>>>>> DMClone(dm, dmGeom);
>>>>> <Create the Section s you want>
>>>>> DMSetLocalSection(dmGeom, s);
>>>>> PetscSectionDestroy(&s);
>>>>>
>>>>> then
>>>>>
>>>>> DMCreateLocalVector(dmGeom, &lv);
>>>>>
>>>>> will give you a Vec with the local data in it that can be addressed by
>>>>> mesh point (global vectors too). Also you would be able
>>>>> to communicate this data if the mesh is redistributed, or replicate
>>>>> this data if you overlap cells in parallel.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>>
>>>>>> Thanks,
>>>>>> Mike
>>>>>>
>>>>>>
>>>>>> On Tue, Apr 26, 2022 at 4:48 PM Mike Michell <mi.mike1021 at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Below two interesting things are found:
>>>>>>>> - If DMPlexCreateGmshFromFile() is used, DMPlexDistribute() is not
>>>>>>>> done by default. I should call DMPlexDistribute() to distribute DMPlex over
>>>>>>>> the procs.
>>>>>>>>
>>>>>>>
>>>>>>> Here is the explanation. The default distribution comes from
>>>>>>> DMSetFromOptions(), which I am guessing was not called after
>>>>>>> DMPlexCreateGmshFromFile().
>>>>>>>
>>>>>>>
>>>>>>>> - If {DMCreate(), DMSetType(), DMSetFromOptions()} in serially used
>>>>>>>> with "-dm_plex_filename ./new.msh" option in command line,
>>>>>>>> DMPlexDistribute() is done by default even though DMPlex object is created
>>>>>>>> by reading the same .msh file.
>>>>>>>>
>>>>>>>> Thanks for the modification to the example ex1f90.F90, now it works
>>>>>>>> with 2 procs.
>>>>>>>> But still, if I print my output to "sol.vtk", the file has only a
>>>>>>>> part of whole mesh. However, the file is okay if I print to "sol.vtu".
>>>>>>>> From "sol.vtu" I can see the entire field with rank. Is using .vtu
>>>>>>>> format preferred by petsc?
>>>>>>>>
>>>>>>>
>>>>>>> VTK is generally for debugging, but it should work. I will take a
>>>>>>> look.
>>>>>>>
>>>>>>> VTU and HDF5 are the preferred formats.
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>> Matt
>>>>>>>
>>>>>>>
>>>>>>>>
>>>>>>>> On Tue, Apr 26, 2022 at 9:33 AM Mike Michell <mi.mike1021 at gmail.com>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Thank you for the answers.
>>>>>>>>>> For the first question, basically, I cannot
>>>>>>>>>> run "/dm/impls/plex/ex1f90.F90" example with more than 1 proc. I removed
>>>>>>>>>> DMPlexDistribute() following your comment and what I tried is:
>>>>>>>>>>
>>>>>>>>>> - no modification to ex1f90.F90 (as it is)
>>>>>>>>>> - make "ex1f90"
>>>>>>>>>> - mpirun -np 2 ./ex1f90
>>>>>>>>>>
>>>>>>>>>> It gives me "Bad termination of one of ..." for Rank 1. The code
>>>>>>>>>> runs okay with "mpirun -np 1 ./ex1f90".
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> You are correct. It evidently never worked in parallel, since
>>>>>>>>> those checks will only work in serial.
>>>>>>>>> I have fixed the code, and added a parallel test. I have attached
>>>>>>>>> the new file, but it is also in this MR:
>>>>>>>>>
>>>>>>>>> https://gitlab.com/petsc/petsc/-/merge_requests/5173
>>>>>>>>>
>>>>>>>>> Thanks,
>>>>>>>>>
>>>>>>>>> Matt
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>> Thanks,
>>>>>>>>>> Mike
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On Mon, Apr 25, 2022 at 9:41 PM Mike Michell <
>>>>>>>>>>> mi.mike1021 at gmail.com> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Dear PETSc developer team,
>>>>>>>>>>>>
>>>>>>>>>>>> I'm trying to learn DMPlex to build a parallel finite volume
>>>>>>>>>>>> code in 2D & 3D. More specifically, I want to read a grid from .msh file by
>>>>>>>>>>>> Gmsh.
>>>>>>>>>>>> For practice, I modified /dm/impls/plex/ex1f90.F90 case to read
>>>>>>>>>>>> & distribute my sample 2D grid, which is attached. I have two questions as
>>>>>>>>>>>> below:
>>>>>>>>>>>>
>>>>>>>>>>>> (1) First, if I do not use my grid, but use the default box
>>>>>>>>>>>> grid built by ex1f90.F90 and if I try "DMPlexDistribute" over mpi
>>>>>>>>>>>> processors, the output file (sol.vtk) has only some portion of the entire
>>>>>>>>>>>> mesh. How can I print out the entire thing into a single file? Is there any
>>>>>>>>>>>> example for parallel output? (Related files attached to "/Question_1/")
>>>>>>>>>>>> Paraview gives me some error messages about data size mismatching.
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> For the last release, we made parallel distribution the default.
>>>>>>>>>>> Thus, you do not need to call DMPlexDIstribute() explicitly here. Taking it
>>>>>>>>>>> out, I can run your example.
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>> (2) If I create DMPlex object through
>>>>>>>>>>>> "DMPlexCreateGmshFromFile", "DMPlexCreateSection" part is crashed. I do not
>>>>>>>>>>>> understand why my example code does not work, because the only change was
>>>>>>>>>>>> switching from "DMCreate" to "DMPlexCreateGmshFromFile" and providing
>>>>>>>>>>>> "new.msh" file. Without the PetscSection object, the code works fine. Any
>>>>>>>>>>>> comments about this? (Related files attached to "/Question_2/")
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> If I remove DMPlexDistribute() from this code, it is clear that
>>>>>>>>>>> the problem is with the "marker" label. We do not create this by default
>>>>>>>>>>> from GMsh since we assume people have defined their own labels. You can pass
>>>>>>>>>>>
>>>>>>>>>>> -dm_plex_gmsh_use_marker
>>>>>>>>>>>
>>>>>>>>>>> to your code. WHen I do this, you example runs for me.
>>>>>>>>>>>
>>>>>>>>>>> Thanks,
>>>>>>>>>>>
>>>>>>>>>>> Matt
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>> Thanks,
>>>>>>>>>>>> Mike
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> --
>>>>>>>>>>> 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/>
>>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> 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/>
>>>>>
>>>>
>>>
>>> --
>>> 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/20220429/e405c652/attachment-0001.html>
More information about the petsc-users
mailing list