[petsc-users] Parallel DMPlex

Fri Oct 13 07:00:01 CDT 2023

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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231013/dc82fc2a/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex1.c
Type: application/octet-stream
Size: 3127 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231013/dc82fc2a/attachment-0001.obj>