[petsc-users] Parallel DMPlex

Matthew Knepley knepley at gmail.com
Mon Oct 16 08:11:58 CDT 2023 

On Mon, Oct 16, 2023 at 6:54 AM erdemguer <erdemguer at proton.me> wrote:

> Hey again.
>
> This code outputs for example:
>
> After Distribution Rank: 1, cStart: 0, cEndInterior: 14, cEnd: 24
> After Distribution Rank: 0, cStart: 0, cEndInterior: 13, cEnd: 27
> [0] m: 39 n: 39
> [1] m: 42 n: 42
>
> Shouldn't it be 39 x 81 and 42 x 72 because of the overlapping cells on
> processor boundaries?
>

Here is my output

master *:~/Downloads/tmp/Guer$ /PETSc3/petsc/apple/bin/mpiexec -n 2 ./ex1
-malloc_debug 0 -dm_refine 1
Before Distribution Rank: 1, cStart: 0, cEndInterior: 0, cEnd: 0
Before Distribution Rank: 0, cStart: 0, cEndInterior: 32, cEnd: 32
After Distribution Rank: 1, cStart: 0, cEndInterior: 16, cEnd: 24
After Distribution Rank: 0, cStart: 0, cEndInterior: 16, cEnd: 24
[0] m: 48 n: 48
[1] m: 48 n: 48

The mesh is 4x4 and also split into two triangles, so 32 triangles. Then we
split it and have 8 overlap cells on each side. You can get quads using

master *:~/Downloads/tmp/Guer$ /PETSc3/petsc/apple/bin/mpiexec -n 2 ./ex1
-malloc_debug 0 -dm_plex_simplex 0 -dm_refine 1 -dm_view
Before Distribution Rank: 1, cStart: 0, cEndInterior: 0, cEnd: 0
Before Distribution Rank: 0, cStart: 0, cEndInterior: 16, cEnd: 16
After Distribution Rank: 1, cStart: 0, cEndInterior: 8, cEnd: 12
After Distribution Rank: 0, cStart: 0, cEndInterior: 8, cEnd: 12
[0] m: 24 n: 24
[1] m: 24 n: 24

It is the same 4x4 mesh, but now with quads.

  Thanks,

     Matt

P.S. It looks like I should use PetscFV or something like that at the first
> place. At first I thought, "I will just use SNES, I will compute only
> residual and jacobian on cells so why do bother with PetscFV?" So
>
> Thanks,
> E.
> Sent with Proton Mail <https://proton.me/> secure email.
>
> ------- Original Message -------
> On Friday, October 13th, 2023 at 3:00 PM, Matthew Knepley <
> knepley at gmail.com> wrote:
>
> On Fri, Oct 13, 2023 at 7:26 AM erdemguer <erdemguer at proton.me> wrote:
>
>> Hi, unfortunately it's me again.
>>
>> I have some weird troubles with creating matrix with DMPlex. Actually I
>> might not need to create matrix explicitly, but SNESSolve crashes at there
>> too. So, I updated the code you provided. When I tried to use
>> DMCreateMatrix() at first, I got an error "Unknown discretization type
>> for field 0" at first I applied DMSetLocalSection() and this error is gone.
>> But this time when I run the code with multiple processors, sometimes I got
>> an output like:
>>
>
> Some setup was out of order so the section size on proc1 was 0, and I was
> not good about checking this.
> I have fixed it and attached.
>
> Thanks,
>
> Matt
>
> Before Distribution Rank: 0, cStart: 0, cEndInterior: 27, cEnd: 27
>> Before Distribution Rank: 1, cStart: 0, cEndInterior: 0, cEnd: 0
>> [1] ghost cell 14
>> [1] ghost cell 15
>> [1] ghost cell 16
>> [1] ghost cell 17
>> [1] ghost cell 18
>> [1] ghost cell 19
>> [1] ghost cell 20
>> [1] ghost cell 21
>> [1] ghost cell 22
>> After Distribution Rank: 1, cStart: 0, cEndInterior: 14, cEnd: 23
>> [0] ghost cell 13
>> [0] ghost cell 14
>> [0] ghost cell 15
>> [0] ghost cell 16
>> [0] ghost cell 17
>> [0] ghost cell 18
>> [0] ghost cell 19
>> [0] ghost cell 20
>> [0] ghost cell 21
>> [0] ghost cell 22
>> [0] ghost cell 23
>> After Distribution Rank: 0, cStart: 0, cEndInterior: 13, cEnd: 24
>> Fatal error in internal_Waitall: Unknown error class, error stack:
>> internal_Waitall(82)......................: MPI_Waitall(count=1,
>> array_of_requests=0xaaaaf5f72264, array_of_statuses=0x1) failed
>> MPIR_Waitall(1099)........................:
>> MPIR_Waitall_impl(1011)...................:
>> MPIR_Waitall_state(976)...................:
>> MPIDI_CH3i_Progress_wait(187).............: an error occurred while
>> handling an event returned by MPIDI_CH3I_Sock_Wait()
>> MPIDI_CH3I_Progress_handle_sock_event(411):
>> ReadMoreData(744).........................: ch3|sock|immedread
>> 0xffff8851c5c0 0xaaaaf5e81cd0 0xaaaaf5e8a880
>> MPIDI_CH3I_Sock_readv(2553)...............: the supplied buffer contains
>> invalid memory (set=0,sock=1,errno=14:Bad address)
>>
>> Sometimes the error message isn't appearing but for example I'm trying to
>> print size of the matrix but it isn't working.
>> If necessary, my Configure options --download-mpich --download-hwloc
>> --download-pastix --download-hypre --download-ml --download-ctetgen
>> --download-triangle --download-exodusii --download-netcdf --download-zlib
>> --download-pnetcdf --download-ptscotch --download-hdf5 --with-cc=clang-16
>> --with-cxx=clang++-16 COPTFLAGS="-g -O2" CXXOPTFLAGS="-g -O2" FOPTFLAGS="-g
>> -O2" --with-debugging=1
>>
>> Version: Petsc Release Version 3.20.0
>>
>> Thank you,
>> Guer
>>
>> Sent with Proton Mail <https://proton.me/> secure email.
>>
>> ------- Original Message -------
>> On Thursday, October 12th, 2023 at 12:59 AM, erdemguer <
>> erdemguer at proton.me> wrote:
>>
>> Thank you! That's exactly what I need.
>>
>> Sent with Proton Mail <https://proton.me/> secure email.
>>
>> ------- Original Message -------
>> On Wednesday, October 11th, 2023 at 4:17 PM, Matthew Knepley <
>> knepley at gmail.com> wrote:
>>
>> On Wed, Oct 11, 2023 at 4:42 AM erdemguer <erdemguer at proton.me> wrote:
>>
>>> Hi again,
>>>
>>
>> I see the problem. FV ghosts mean extra boundary cells added in FV
>> methods using DMPlexCreateGhostCells() in order to impose boundary
>> conditions. They are not the "ghost" cells for overlapping parallel
>> decompositions. I have changed your code to give you what you want. It is
>> attached.
>>
>> Thanks,
>>
>> Matt
>>
>>> Here is my code:
>>> #include <petsc.h>
>>> static char help[] = "dmplex";
>>>
>>> int main(int argc, char **argv)
>>> {
>>> PetscCall(PetscInitialize(&argc, &argv, NULL, help));
>>> DM dm, dm_dist;
>>> PetscSection section;
>>> PetscInt cStart, cEndInterior, cEnd, rank;
>>> PetscInt nc[3] = {3, 3, 3};
>>> PetscReal upper[3] = {1, 1, 1};
>>>
>>> PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
>>>
>>> DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 3, PETSC_FALSE, nc, NULL, upper,
>>> NULL, PETSC_TRUE, &dm);
>>> DMViewFromOptions(dm, NULL, "-dm1_view");
>>> PetscCall(DMSetFromOptions(dm));
>>> DMViewFromOptions(dm, NULL, "-dm2_view");
>>>
>>> PetscCall(DMPlexGetDepthStratum(dm, 3, &cStart, &cEnd));
>>> DMPlexComputeCellTypes(dm);
>>> PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_INTERIOR_GHOST,
>>> &cEndInterior, NULL));
>>> PetscPrintf(PETSC_COMM_SELF, "Before Distribution Rank: %d, cStart: %d,
>>> cEndInterior: %d, cEnd: %d\n", rank, cStart,
>>> cEndInterior, cEnd);
>>>
>>> PetscInt nField = 1, nDof = 3, field = 0;
>>> PetscCall(PetscSectionCreate(PETSC_COMM_WORLD, &section));
>>> PetscSectionSetNumFields(section, nField);
>>> PetscCall(PetscSectionSetChart(section, cStart, cEnd));
>>> for (PetscInt p = cStart; p < cEnd; p++)
>>> {
>>> PetscCall(PetscSectionSetFieldDof(section, p, field, nDof));
>>> PetscCall(PetscSectionSetDof(section, p, nDof));
>>> }
>>>
>>> PetscCall(PetscSectionSetUp(section));
>>>
>>> DMSetLocalSection(dm, section);
>>> DMViewFromOptions(dm, NULL, "-dm3_view");
>>>
>>> DMSetAdjacency(dm, field, PETSC_TRUE, PETSC_TRUE);
>>> DMViewFromOptions(dm, NULL, "-dm4_view");
>>> PetscCall(DMPlexDistribute(dm, 1, NULL, &dm_dist));
>>> if (dm_dist)
>>> {
>>> DMDestroy(&dm);
>>> dm = dm_dist;
>>> }
>>> DMViewFromOptions(dm, NULL, "-dm5_view");
>>> PetscCall(DMPlexGetDepthStratum(dm, 3, &cStart, &cEnd));
>>> DMPlexComputeCellTypes(dm);
>>> PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST,
>>> &cEndInterior, NULL));
>>> PetscPrintf(PETSC_COMM_SELF, "After Distribution Rank: %d, cStart: %d,
>>> cEndInterior: %d, cEnd: %d\n", rank, cStart,
>>> cEndInterior, cEnd);
>>>
>>> DMDestroy(&dm);
>>> PetscCall(PetscFinalize());
>>> }
>>>
>>> This codes output is currently (on 2 processors) is:
>>> Before Distribution Rank: 1, cStart: 0, cEndInterior: -1, cEnd: 14
>>> Before Distribution Rank: 0, cStart: 0, cEndInterior: -1, cEnd: 13
>>> After Distribution Rank: 0, cStart: 0, cEndInterior: -1, cEnd: 27
>>> After Distribution Rank: 1, cStart: 0, cEndInterior: -1, cEnd: 24
>>>
>>> DMView outputs:
>>> dm1_view (after creation):
>>> DM Object: 2 MPI processes
>>> type: plex
>>> DM_0x84000004_0 in 3 dimensions:
>>> Number of 0-cells per rank: 64 0
>>> Number of 1-cells per rank: 144 0
>>> Number of 2-cells per rank: 108 0
>>> Number of 3-cells per rank: 27 0
>>> Labels:
>>> marker: 1 strata with value/size (1 (218))
>>> Face Sets: 6 strata with value/size (6 (9), 5 (9), 3 (9), 4 (9), 1 (9),
>>> 2 (9))
>>> depth: 4 strata with value/size (0 (64), 1 (144), 2 (108), 3 (27))
>>> celltype: 4 strata with value/size (7 (27), 0 (64), 4 (108), 1 (144))
>>>
>>> dm2_view (after setfromoptions):
>>> DM Object: 2 MPI processes
>>> type: plex
>>> DM_0x84000004_0 in 3 dimensions:
>>> Number of 0-cells per rank: 40 46
>>> Number of 1-cells per rank: 83 95
>>> Number of 2-cells per rank: 57 64
>>> Number of 3-cells per rank: 13 14
>>> Labels:
>>> depth: 4 strata with value/size (0 (40), 1 (83), 2 (57), 3 (13))
>>> marker: 1 strata with value/size (1 (109))
>>> Face Sets: 5 strata with value/size (1 (6), 2 (1), 3 (7), 5 (5), 6 (4))
>>> celltype: 4 strata with value/size (0 (40), 1 (83), 4 (57), 7 (13))
>>>
>>> dm3_view (after setting local section):
>>> DM Object: 2 MPI processes
>>> type: plex
>>> DM_0x84000004_0 in 3 dimensions:
>>> Number of 0-cells per rank: 40 46
>>> Number of 1-cells per rank: 83 95
>>> Number of 2-cells per rank: 57 64
>>> Number of 3-cells per rank: 13 14
>>> Labels:
>>> depth: 4 strata with value/size (0 (40), 1 (83), 2 (57), 3 (13))
>>> marker: 1 strata with value/size (1 (109))
>>> Face Sets: 5 strata with value/size (1 (6), 2 (1), 3 (7), 5 (5), 6 (4))
>>> celltype: 4 strata with value/size (0 (40), 1 (83), 4 (57), 7 (13))
>>> Field Field_0:
>>> adjacency FEM
>>>
>>> dm4_view (after setting adjacency):
>>> DM Object: 2 MPI processes
>>> type: plex
>>> DM_0x84000004_0 in 3 dimensions:
>>> Number of 0-cells per rank: 40 46
>>> Number of 1-cells per rank: 83 95
>>> Number of 2-cells per rank: 57 64
>>> Number of 3-cells per rank: 13 14
>>> Labels:
>>> depth: 4 strata with value/size (0 (40), 1 (83), 2 (57), 3 (13))
>>> marker: 1 strata with value/size (1 (109))
>>> Face Sets: 5 strata with value/size (1 (6), 2 (1), 3 (7), 5 (5), 6 (4))
>>> celltype: 4 strata with value/size (0 (40), 1 (83), 4 (57), 7 (13))
>>> Field Field_0:
>>> adjacency FVM++
>>>
>>> dm5_view (after distribution):
>>> DM Object: Parallel Mesh 2 MPI processes
>>> type: plex
>>> Parallel Mesh in 3 dimensions:
>>> Number of 0-cells per rank: 64 60
>>> Number of 1-cells per rank: 144 133
>>> Number of 2-cells per rank: 108 98
>>> Number of 3-cells per rank: 27 24
>>> Labels:
>>> depth: 4 strata with value/size (0 (64), 1 (144), 2 (108), 3 (27))
>>> marker: 1 strata with value/size (1 (218))
>>> Face Sets: 6 strata with value/size (1 (9), 2 (9), 3 (9), 4 (9), 5 (9),
>>> 6 (9))
>>> celltype: 4 strata with value/size (0 (64), 1 (144), 4 (108), 7 (27))
>>> Field Field_0:
>>> adjacency FVM++
>>>
>>> Thanks,
>>> Guer.
>>> Sent with Proton Mail <https://proton.me/> secure email.
>>>
>>> ------- Original Message -------
>>> On Wednesday, October 11th, 2023 at 3:33 AM, Matthew Knepley <
>>> knepley at gmail.com> wrote:
>>>
>>> On Tue, Oct 10, 2023 at 7:01 PM erdemguer <erdemguer at proton.me> wrote:
>>>
>>>>
>>>> Hi,
>>>> Sorry for my late response. I tried with your suggestions and I think I
>>>> made a progress. But I still got issues. Let me explain my latest mesh
>>>> routine:
>>>>
>>>>
>>>>    1. DMPlexCreateBoxMesh
>>>>    2. DMSetFromOptions
>>>>    3. PetscSectionCreate
>>>>    4. PetscSectionSetNumFields
>>>>    5. PetscSectionSetFieldDof
>>>>    6. PetscSectionSetDof
>>>>    7. PetscSectionSetUp
>>>>    8. DMSetLocalSection
>>>>    9. DMSetAdjacency
>>>>    10. DMPlexDistribute
>>>>
>>>>
>>>> It's still not working but it's promising, if I call
>>>> DMPlexGetDepthStratum for cells, I can see that after distribution
>>>> processors have more cells.
>>>>
>>>
>>> Please send the output of DMPlexView() for each incarnation of the mesh.
>>> What I do is put
>>>
>>> DMViewFromOptions(dm, NULL, "-dm1_view")
>>>
>>>
>>> with a different string after each call.
>>>
>>>> But I couldn't figure out how to decide where the ghost/processor
>>>> boundary cells start.
>>>>
>>>
>>> Please send the actual code because the above is not specific enough.
>>> For example, you will not have
>>> "ghost cells" unless you partition with overlap. This is because by
>>> default cells are the partitioned quantity,
>>> so each process gets a unique set.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>> In older mails I saw there is a function DMPlexGetHybridBounds but I
>>>> think that function is deprecated. I tried to use,
>>>> DMPlexGetCellTypeStratum as in ts/tutorials/ex11_sa.c but I'm getting
>>>> -1 as cEndInterior before and after distribution. I tried it for DM_POLYTOPE_FV_GHOST,
>>>> DM_POLYTOPE_INTERIOR_GHOST polytope types. I also tried calling
>>>> DMPlexComputeCellTypes before DMPlexGetCellTypeStratum but nothing changed.
>>>> I think I can calculate the ghost cell indices using cStart/cEnd before &
>>>> after distribution but I think there is a better way I'm currently missing.
>>>>
>>>> Thanks again,
>>>> Guer.
>>>>
>>>> ------- Original Message -------
>>>> On Thursday, September 28th, 2023 at 10:42 PM, Matthew Knepley <
>>>> knepley at gmail.com> wrote:
>>>>
>>>> On Thu, Sep 28, 2023 at 3:38 PM erdemguer via petsc-users <
>>>> petsc-users at mcs.anl.gov> wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> I am currently using DMPlex in my code. It runs serially at the
>>>>> moment, but I'm interested in adding parallel options. Here is my workflow:
>>>>>
>>>>> Create a DMPlex mesh from GMSH.
>>>>> Reorder it with DMPlexPermute.
>>>>> Create necessary pre-processing arrays related to the mesh/problem.
>>>>> Create field(s) with multi-dofs.
>>>>> Create residual vectors.
>>>>> Define a function to calculate the residual for each cell and, use
>>>>> SNES.
>>>>> As you can see, I'm not using FV or FE structures (most examples do).
>>>>> Now, I'm trying to implement this in parallel using a similar approach.
>>>>> However, I'm struggling to understand how to create corresponding vectors
>>>>> and how to obtain index sets for each processor. Is there a tutorial or
>>>>> paper that covers this topic?
>>>>>
>>>>
>>>> The intention was that there is enough information in the manual to do
>>>> this.
>>>>
>>>> Using PetscFE/PetscFV is not required. However, I strongly encourage
>>>> you to use PetscSection. Without this, it would be incredibly hard to do
>>>> what you want. Once the DM has a Section, it can do things like
>>>> automatically create vectors and matrices for you. It can redistribute
>>>> them, subset them, etc. The Section describes how dofs are assigned to
>>>> pieces of the mesh (mesh points). This is in the manual, and there are a
>>>> few examples that do it by hand.
>>>>
>>>> So I suggest changing your code to use PetscSection, and then letting
>>>> us know if things still do not work.
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>>> Thank you.
>>>>> Guer.
>>>>>
>>>>> Sent with Proton Mail <https://proton.me/> secure email.
>>>>>
>>>>
>>>>
>>>> --
>>>> 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/20231016/7953816e/attachment-0001.html>

More information about the petsc-users mailing list