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 <[email protected]>님이 작성: > On Fri, Apr 29, 2022 at 8:27 AM Mike Michell <[email protected]> > 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 <[email protected]> >>> 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 <[email protected]> >>>>> 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 <[email protected]> >>>>>>> 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 <[email protected]> >>>>>>>>> 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 < >>>>>>>>>>> [email protected]> 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/> >
