Regarding distributed meshes with p4est and trilinios

[Korinna Rosin]

Dear all,

I am currently working on an adaptive algorithm in parallel. To this end, I use p4est and trilinos in a parallel::distributed manner.
A DWR-error-estimator and the global refinement produce reasonable things.

I think that I found a wrong assertion and a bug, but maybe I am mistaken and did not use the given methods correctly. I would really appreciate it, if someone could look into this and either confirm my suspicion or prove me wrong.

1) Maybe I found a wrong assertion in the FETools.
    I have trouble with FETools while working in the debug mode. I cannot use the provided FETools::interpolation for the higher order reconstruction, which is typical used in the context of the DWR-method. If I am not mistaken, we have to use a ghosted vector for the InputVector and a non-ghosted one for the OutputVector. But I cannot access the locally_owned_elements of the ghosted vector, which I created based on locally_relevant_dofs.
   This is a part of the error, which I receive:
       An error occurred in line <1099> of file </home/krosin/SWD/deal_ii/dealii-8.5.1/include/deal.II/lac/trilinos_vector_base.h> in function
           dealii::IndexSet dealii::TrilinosWrappers::VectorBase::locally_owned_elements() const
       The violated condition was:
           owned_elements.size()==size()
       Additional information:
          The locally owned elements have not been properly initialized! This happens for example if this object has been initialized with exactly one overlapping IndexSet.
   
   In the release mode (and after deleting the lines 92,94-97 in the file dealii-8.5.1/source/fe/fe_tools_interpolate.cc) everything runs nicely.
   The content of those lines were:
        (92)   const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs();
        (93)   const IndexSet &dof2_local_dofs = dof2.locally_owned_dofs();
        (94)   const IndexSet u1_elements = u1.locally_owned_elements();
        (94)   Assert(u1_elements == dof1_local_dofs,
                 ExcMessage("The provided vector and DoF handler should have the same"
                          " index sets."));
       (94)   Assert(u2_elements == dof2_local_dofs,
                 ExcMessage("The provided vector and DoF handler should have the same"
                          " index sets."));
  
   Hence I suspect that the assertion in fe_tools_interpolate.cc, which is also used in e.g. interpolation_difference is simply not correct for the parallel::distributed case.
   I have also created a mini example (no. 1). Therein, I included a suggestion for a different assertion.


2) I have also a much worse problem right now, for which I have not found a satisfying solution yet.
   It seems to me that I cannot use the MeshSmooting methods for the parallel::distributed::Triangulation. In particular the flag patch_level_1 causes trouble.
   See the mini example (no. 2).

   This is the conclusion of my colleague Dustin Kumor and me:
   This patch_level_1-property is only satisfied on the locally owned cells, which is all I need for my computations, but the MeshSmoothing presumes that this property is satisfied on the complete triangulation (on one processor including artificial cells).

   However I do have a rather dirty hack, where I completly ignore everything related to coarsening flags in STEP 5 (ensure patch level 1) of prepare_coarsening_and_refinement() in tria.cc .
   Dustin and I believe a small change of the following lines in file tria.cc (version 8.5.1) solves the problem with refinement flags:
        (146)       if (cell->child(i)->active())
    (12997)                if (cell->child(0)->has_children() == true)
   to
        (146)       if (cell->child(i)->active() && cell->child(i)->is_locally_owned())
    (12997)                if (cell->child(0)->has_children() == true || !cell->child(0)->is_locally_owned())


   By the way: Is it possible to consider only the refinement flags and ignore all coarsening?


I am sorry for the lengthy text. I hope that I provided all relevant information and I am looking forward to your input.


Best,
Korinna


PS:
Used versions are:

gcc               4.8.8
openmpi        1.10
trilinos         12.4
p4est            1.1
deal.II           8.5.1 and 9.0.0-pre

[Timo Heister]

Korinna,

> 1) Maybe I found a wrong assertion in the FETools.
> I have trouble with FETools while working in the debug mode. I cannot
> use the provided FETools::interpolation for the higher order reconstruction,

It might be that we haven't been testing this function enough, so a
bug like this is not unlikely.

> The violated condition was:
> owned_elements.size()==size()
> Additional information:
> The locally owned elements have not been properly initialized!
> This happens for example if this object has been initialized with exactly
> one overlapping IndexSet.

Can you try initializing your ghosted vector with
vec.reinit(locally_owned_index_set, locally_relevant_index_set, mpi_com);
instead of only specifying the ghosted set? Maybe that fixes the exception.

> 2) I have also a much worse problem right now, for which I have not found a
> satisfying solution yet.
> It seems to me that I cannot use the MeshSmooting methods for the
> parallel::distributed::Triangulation. In particular the flag patch_level_1
> causes trouble.

Yes, mesh smoothing in parallel is a problem at the moment. It is
basically not working at all. Sorry. I was hoping to have an
implementation for that soon (we did some progress towards it last
month). See https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_pull_1711&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=fx5zPVPj5YHLwSgXHBKNn1I9OvTImEnlVgDaDmSH-Z0&s=xMjuyhQn_tZbaM_K32G208xiQJ7EjMmKAmSRfDoy_4Y&e=
https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_issues_1655&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=fx5zPVPj5YHLwSgXHBKNn1I9OvTImEnlVgDaDmSH-Z0&s=HmN6gJQSWE1WxMwrHI5csg3_TaC8JnCXyOa6GAnnxGQ&e= and issues linked in
there.
One workaround for now might be to switch to a shared::Triangulation.

Best,
Timo

--
Timo Heister
http://www.math.clemson.edu/~heister/

[Korinna Rosin]

Hi Timo,

Am Donnerstag, 7. September 2017 13:40:19 UTC+2 schrieb Timo Heister:

Korinna,

> 1) Maybe I found a wrong assertion in the FETools.
>     I have trouble with FETools while working in the debug mode. I cannot
> use the provided FETools::interpolation for the higher order reconstruction,

It might be that we haven't been testing this function enough, so a
bug like this is not unlikely.

>        The violated condition was:
>            owned_elements.size()==size()
>        Additional information:
>           The locally owned elements have not been properly initialized!
> This happens for example if this object has been initialized with exactly
> one overlapping IndexSet.

Can you try initializing your ghosted vector with
vec.reinit(locally_owned_index_set, locally_relevant_index_set, mpi_com);
instead of only specifying the ghosted set? Maybe that fixes the exception.

Thanks for great tipp. That was rather easy.

But should I use the locally_relevant_index_set as the second argument or locally_relevant_dofs.subtract_set(locally_owned_dofs)?

Both versions work for me so far.

> 2) I have also a much worse problem right now, for which I have not found a
> satisfying solution yet.
>    It seems to me that I cannot use the MeshSmooting methods for the
> parallel::distributed::Triangulation. In particular the flag patch_level_1
> causes trouble.

Yes, mesh smoothing in parallel is a problem at the moment. It is
basically not working at all. Sorry. I was hoping to have an
implementation for that soon (we did some progress towards it last
month). See https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_pull_1711&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=fx5zPVPj5YHLwSgXHBKNn1I9OvTImEnlVgDaDmSH-Z0&s=xMjuyhQn_tZbaM_K32G208xiQJ7EjMmKAmSRfDoy_4Y&e=
https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_issues_1655&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=fx5zPVPj5YHLwSgXHBKNn1I9OvTImEnlVgDaDmSH-Z0&s=HmN6gJQSWE1WxMwrHI5csg3_TaC8JnCXyOa6GAnnxGQ&e=  and issues linked in
there.
One workaround for now might be to switch to a shared::Triangulation.

Oh. Okay. Thank you for the information. Sorry for the duplicate thread. I will read those threads carefully next.

So I guess it is not possible right now to only prepare and excecute refinement without considering any coarsening flags?

Best,

Korinna

[Timo Heister]

> Thanks for great tipp. That was rather easy.
> But should I use the locally_relevant_index_set as the second argument or
> locally_relevant_dofs.subtract_set(locally_owned_dofs)?
> Both versions work for me so far.

Great! That doesn't matter (we actually subtract that set to get the
ghost values inside that function).

> Oh. Okay. Thank you for the information. Sorry for the duplicate thread. I
> will read those threads carefully next.

No problem. That is not a thing that is easy to find if you don't know
what you are looking for.

> So I guess it is not possible right now to only prepare and excecute
> refinement without considering any coarsening flags?

No, I don't think so (without modifying the source code). The flags
are set in prepare_coarsening_and_refinement().
So one hack you could do is the following
class MyTriangulation : public parallel::distributed::Triangulation
{
public:
bool prepare_coarsening_and_refinement() {
// 1. call dealii::Triangulation<dim>::prepare_coarsening_and_refinement();
// 2. remove all coarsening flags
// 3. if any flags changed, return true, else return false
}
};

and use that. Take a look at the implementation in distributed/tria.cc
for inspiration. Note that this will break if you have periodic
boundary conditions.

Let us know if this lets you run with patch_level_1

[Korinna Rosin]

Hello again,

it took me longer than expected to get my code running without any obvious errors.
Basically, I did not manage to fix my problem with the patch structure on parallel::distributed triangulations without touching the source code of deal.II.

Am Donnerstag, 7. September 2017 15:52:57 UTC+2 schrieb Timo Heister:

[...]

> So I guess it is not possible right now to only prepare and excecute
> refinement without considering any coarsening flags?

No, I don't think so (without modifying the source code). The flags
are set in prepare_coarsening_and_refinement().
So one hack you could do is the following
class MyTriangulation : public parallel::distributed::Triangulation
{
public:
bool prepare_coarsening_and_refinement() {
// 1. call dealii::Triangulation<dim>::prepare_coarsening_and_refinement();
// 2. remove all coarsening flags
// 3. if any flags changed, return true, else return false
}
};

and use that. Take a look at the implementation in distributed/tria.cc
for inspiration. Note that this will break if you have periodic
boundary conditions.

 

I have not touched the distributed/tria.cc but grid/tria.cc. Hence the periodic boundary conditions will break most definitely.

Let us know if this lets you run with patch_level_1

Now, I do have a version that seems to enforce the patch_level_1 condition on the locally owned cells. I believe that on the artificial ones the patch_level_1 condition must not be set.

For this I have been working with the grid/tria.cc (version 8.5.1) and have changed the methods cell_is_patch_level_1() and prepare_coarsening_and_refinement().

Especially Steps 3 and 5.

My modification slows down the algorithm a little bit as I cannot use the actual patch_level_1 - presumption.

But maybe someone else is interested in the code. (Which is not rigorously tested.)

Also, there should not be any conflicts with the current 9.0.0-pre version. However, I only checked the source code with a diff-tool.

Best,
Korinna


PS: I took a look at the github. Do you think this should be further discussed in Issue 1655?