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:
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 0xaaaaf5e8a880MPIDI_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 <erdemg...@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
> <knep...@gmail.com> wrote:
>
>> On Wed, Oct 11, 2023 at 4:42 AM erdemguer <erdemg...@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: 27After
>>> 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
>>> <knep...@gmail.com> wrote:
>>>
>>>> On Tue, Oct 10, 2023 at 7:01 PM erdemguer <erdemg...@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:
>>>>>
>>>>> - DMPlexCreateBoxMesh
>>>>>
>>>>> - DMSetFromOptions
>>>>> - PetscSectionCreate
>>>>> - PetscSectionSetNumFields
>>>>> - PetscSectionSetFieldDof
>>>>>
>>>>> - PetscSectionSetDof
>>>>>
>>>>> - PetscSectionSetUp
>>>>> - DMSetLocalSection
>>>>> - DMSetAdjacency
>>>>> - 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,
>>>>> DMPlexGetCellTypeStratumas 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
>>>>> <knep...@gmail.com> wrote:
>>>>>
>>>>>> On Thu, Sep 28, 2023 at 3:38 PM erdemguer via petsc-users
>>>>>> <petsc-users@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/)#include <petsc.h>
static char help[] = "dmplex";
// Ghost cells are in the point SF
static PetscErrorCode GetLeafCells(DM dm, PetscInt *cEndInterior)
{
PetscSF sf;
const PetscInt *leaves;
PetscInt Nl, cStart, cEnd;
PetscMPIInt rank;
PetscFunctionBeginUser;
PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
*cEndInterior = cEnd;
PetscCall(DMGetPointSF(dm, &sf));
PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
for (PetscInt l = 0; l < Nl; ++l)
{
const PetscInt leaf = leaves[l];
if (leaf >= cStart && leaf < cEnd)
{
*cEndInterior = PetscMin(leaf, *cEndInterior);
PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] ghost cell %"
PetscInt_FMT "\n", rank, leaf));
}
}
PetscFunctionReturn(0);
}
int main(int argc, char **argv)
{
Mat J;
DM dm, dm_dist;
PetscSection section;
PetscInt cStart, cEndInterior, cEnd, rank, m, n;
PetscInt nField = 1, nDof = 3, field = 0;
PetscCall(PetscInitialize(&argc, &argv, NULL, help));
PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
PetscCall(DMSetType(dm, DMPLEX));
PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
PetscCall(DMSetFromOptions(dm));
PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
PetscCall(GetLeafCells(dm, &cEndInterior));
PetscCall(PetscPrintf(PETSC_COMM_SELF, "Before Distribution Rank: %d,
cStart: %d, cEndInterior: %d, cEnd: %d\n",
rank, cStart, cEndInterior, cEnd));
PetscCall(PetscSectionCreate(PETSC_COMM_WORLD, §ion));
PetscCall(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));
PetscCall(DMSetLocalSection(dm, section));
PetscCall(DMSetAdjacency(dm, field, PETSC_TRUE, PETSC_FALSE));
PetscCall(DMPlexDistribute(dm, 1, NULL, &dm_dist));
if (dm_dist)
{
PetscCall(DMDestroy(&dm));
dm = dm_dist;
}
PetscCall(DMSetLocalSection(dm, section));
PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
PetscCall(GetLeafCells(dm, &cEndInterior));
PetscCall(PetscPrintf(PETSC_COMM_SELF, "After Distribution Rank: %d,
cStart: %d, cEndInterior: %d, cEnd: %d\n",
rank, cStart, cEndInterior, cEnd));
PetscCall(DMCreateMatrix(dm, &J));
PetscCall(MatGetSize(J, &m, &n));
PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] m: %d n: %d\n", rank, m, n));
PetscCall(DMDestroy(&dm));
PetscCall(MatDestroy(&J));
PetscCall(PetscSectionDestroy(§ion));
PetscCall(PetscFinalize());
return 0;
}
/*TEST
test:
args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3
TEST*/
