[petsc-users] Parallel DMPlex
Matthew Knepley
knepley at gmail.com
Mon Oct 16 08:26:14 CDT 2023
On Mon, Oct 16, 2023 at 9:22 AM erdemguer <erdemguer at proton.me> wrote:
> Thank you for your responses many times. Looks like I'm missing something,
> sorry for my confusion, but let's take processor 0 on your first output.
> cEndInterior: 16 and cEnd: 24.
> I'm calculating jacobian for cell=14, dof=0 (row = 42) and cell=18, dof=2
> (col = 56) have influence on it. (Cell 18 is on processor boundary)
> Shouldn't I have to write values on the (42,56)?
>
Imagine you are me getting this mail. When I mail you, I show you _exactly_
what I ran and which command line options I used. You do not. I provide
you all the output. You do not. You can see that someone would only be
guessing when replying to this email. Also note that you have two dofs per
cell, so the cell numbers are not the row numbers for the Jacobian. Please
send something reproducible when you want help on running.
Thanks,
Matt
> Thanks,
> Guer
>
> Sent with Proton Mail <https://proton.me/> secure email.
>
> ------- Original Message -------
> On Monday, October 16th, 2023 at 4:11 PM, Matthew Knepley <
> knepley at gmail.com> wrote:
>
> 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, §ion));
>>>> 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/>
>
>
>
--
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/5781ce23/attachment-0001.html>
More information about the petsc-users
mailing list