Shared Triangulation and Locally Owned Dofs issue

[GivAlz]

Hi everyone!
I believe there is some problem at least one of "locally_owned_dofs" and "map_dofs_to_support_points" on shared triangulation. The following test is a modified half of step-60, where I do the following:

• build a shared triangulation
• deform it
• ask for support points
• ask for the indices of locally relevant support points
• check if the support points with such indexes are really numbers, or contain -nan (I check only the first coordinate, but it is enough to trigger the problem).

There might be something wrong in my lengthy procedure, but with 1 process everything is fine (also by printing and checking manually the points), while with multiple cores it sometimes works, sometimes doesn't (depending on the deformation function and number of cores used)...

I am sorry if for the 250 lines; I have been trying to reduce it, but I couldn't understand how to produce a mesh_configuration_function without the parameters...(and, being a part of step-60, on the first run will fail and generate a .prm file )

#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>

#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/fe/mapping_fe_field.h>

#include <deal.II/dofs/dof_tools.h>
#include <deal.II/base/parsed_function.h>

#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/shared_tria.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/non_matching/coupling.h>

#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/linear_operator_tools.h>

#include <deal.II/base/mpi.h>

#include <math.h>
#include <iostream>
#include <fstream>

namespace Step60
{
  using namespace dealii;

  template <int dim, int spacedim = dim>
  class DistributedLagrangeProblem
  {
  public:
    class Parameters : public ParameterAcceptor
    {
    public:
      Parameters();

      unsigned int initial_mesh_refinement = 8;

      unsigned int embedding_space_finite_element_degree = 1;

      unsigned int mesh_space_finite_element_degree = 1;

      unsigned int mesh_configuration_finite_element_degree = 1;

      bool initialized = false;
    };

    DistributedLagrangeProblem(const Parameters &parameters);

    void run();

  private:
    // Object containing the actual parameters
    const Parameters &parameters;

    void setup_grids_and_dofs();

    void setup_mesh_dofs();

    std::unique_ptr< parallel::shared::Triangulation<dim,spacedim> > mesh_grid;
    std::unique_ptr<FiniteElement<dim, spacedim>> mesh_fe;
    std::unique_ptr<DoFHandler<dim, spacedim>>    mesh_dh;

    std::unique_ptr<FiniteElement<dim, spacedim>> mesh_configuration_fe;
    std::unique_ptr<DoFHandler<dim, spacedim>>    mesh_configuration_dh;
    Vector<double>                                mesh_configuration;

    ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
      mesh_configuration_function;

    std::unique_ptr<Mapping<dim, spacedim>> mesh_mapping;

  };


  template <int dim, int spacedim>
  DistributedLagrangeProblem<dim, spacedim>::Parameters::Parameters()
    : ParameterAcceptor("/Distributed Lagrange<" +
                        Utilities::int_to_string(dim) + "," +
                        Utilities::int_to_string(spacedim) + ">/")
  {
    add_parameter("Initial embedded space refinement",
                  initial_mesh_refinement);

    add_parameter("Embedded space finite element degree",
                  mesh_space_finite_element_degree);

    add_parameter("Embedded configuration finite element degree",
                  mesh_configuration_finite_element_degree);

    parse_parameters_call_back.connect([&]() -> void { initialized = true; });
  }

  template <int dim, int spacedim>
  DistributedLagrangeProblem<dim, spacedim>::DistributedLagrangeProblem(
    const Parameters &parameters)
    : parameters(parameters)
    , mesh_configuration_function("Embedded configuration", spacedim)
  {
    mesh_configuration_function.declare_parameters_call_back.connect(
      []() -> void {
        ParameterAcceptor::prm.set("Function constants", "R=.3, Cx=.4, Cy=.4");


        ParameterAcceptor::prm.set("Function expression",
                                   "R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy");
      });
  }

  template <int dim, int spacedim>
  void DistributedLagrangeProblem<dim, spacedim>::setup_grids_and_dofs()
  {
    mesh_grid = std_cxx14::make_unique<parallel::shared::Triangulation<dim, spacedim>>(MPI_COMM_WORLD);
    GridGenerator::hyper_cube(*mesh_grid);
    mesh_grid->refine_global(parameters.initial_mesh_refinement);

    mesh_configuration_fe = std_cxx14::make_unique<FESystem<dim, spacedim>>(
      FE_Q<dim, spacedim>(
        parameters.mesh_configuration_finite_element_degree),
      spacedim);

    mesh_configuration_dh =
      std_cxx14::make_unique<DoFHandler<dim, spacedim>>(*mesh_grid);

    mesh_configuration_dh->distribute_dofs(*mesh_configuration_fe);
    mesh_configuration.reinit(mesh_configuration_dh->n_dofs());

    VectorTools::interpolate(*mesh_configuration_dh,
                             mesh_configuration_function,
                             mesh_configuration);

    mesh_mapping =
        std_cxx14::make_unique<MappingFEField<dim,
                                              spacedim,
                                              Vector<double>,
                                              DoFHandler<dim, spacedim>>>(
          *mesh_configuration_dh, mesh_configuration);

    setup_mesh_dofs();

    std::vector<Point<spacedim>> support_points(mesh_dh->n_dofs());
    DoFTools::map_dofs_to_support_points(*mesh_mapping,
                                           *mesh_dh,
                                           support_points);

    // We own only a part of the support points, those with the following index
    auto owned_dofs_idx = mesh_dh->locally_owned_dofs();

    // Cheking if the local support points are fine
    for(auto idx: owned_dofs_idx)
    {
      if( isnan(support_points[idx][0]) )
      {
        std::cout << "Non-valid point found: " << support_points[idx] << std::endl;
        AssertThrow(false,
                    ExcMessage("Sorry, this time not everything is fine... "
                               "A point containing NaN as coordinates is not valid! "
                               "Something is wrong! "));
      }
    }

    deallog << "Everything seems fine" << std::endl;
  }

  template <int dim, int spacedim>
  void DistributedLagrangeProblem<dim, spacedim>::setup_mesh_dofs()
  {
    mesh_dh =
      std_cxx14::make_unique<DoFHandler<dim, spacedim>>(*mesh_grid);
    mesh_fe = std_cxx14::make_unique<FE_Q<dim, spacedim>>(
      parameters.mesh_space_finite_element_degree);
    mesh_dh->distribute_dofs(*mesh_fe);
  }

  template <int dim, int spacedim>
  void DistributedLagrangeProblem<dim, spacedim>::run()
  {
    AssertThrow(parameters.initialized, ExcNotInitialized());
    deallog.depth_console(10);

    setup_grids_and_dofs();
  }
} // namespace Step60


int main(int argc, char **argv)
{
  try
    {
      using namespace dealii;
      using namespace Step60;
      auto mpi_initialization = Utilities::MPI::MPI_InitFinalize(argc, argv, 1);

      const unsigned int dim = 1, spacedim = 2;

      DistributedLagrangeProblem<dim, spacedim>::Parameters parameters;
      DistributedLagrangeProblem<dim, spacedim>             problem(parameters);

      std::string parameter_file;
      if (argc > 1)
        parameter_file = argv[1];
      else
        parameter_file = "parameters.prm";

      ParameterAcceptor::initialize(parameter_file, "used_parameters.prm");
      problem.run();
    }
  catch (std::exception &exc)
    {
      std::cerr << std::endl
                << std::endl
                << "----------------------------------------------------"
                << std::endl;
      std::cerr << "Exception on processing: " << std::endl
                << exc.what() << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;
      return 1;
    }
  catch (...)
    {
      std::cerr << std::endl
                << std::endl
                << "----------------------------------------------------"
                << std::endl;
      std::cerr << "Unknown exception!" << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;
      return 1;
    }
  return 0;
}

P.S. So, @luca-heltai , turns out I really found a problem...

[masterleinad]

I fail to reproduce your problem with the test case above (copying the main function from step-60, constructing an MPI_InitFinalize object and replacing isnan by std::isnan) running with anything from 1 to 10 MPI processes. Can you give an example for the number of processes such that the code above fails?

[GivAlz]

Sorry, I didn't copy properly the whole code. Now I added the main function (which should be the same as the one you copied from step-60).
For me running it with -np 3 does the trick...but now I wonder if there's a problem with my configuration...

[masterleinad]

Hmm... Is still can't reproduce the problem. For this to build I still need to replace isnan by std::isnan, but that should not change anything. Removing the check and printing everything instead, this is the sorted output with distinct lines (sort -u) I get when running with 3 processes:
outpu3.txt

Maybe, it is related to partitioning. Which partitioning are you using? Can you attach your detailed.log?
The revision I am using is ff5462e I can try later with a newer version.

[masterleinad]

I can reproduce the issue with 28fa42f.

[GivAlz]

Ok thank you @masterleinad !
I believe it is related to the partitioning (grids are non-matching, so find_active_cell_around_point is getting an input which can't be handled locally)...hopefully in the next days I'll have time to start a PR (if the problem is non-matching grids, as I think).
Still it's strange that with a different build you didn't get the problem.

Anyway here is the detailed.log
detailed.log

[luca-heltai]

If you replace Triangulation with parallel::shared::Triangulation, you also have to replace compute_point_locations with distributed_compute_point_locations

[masterleinad]

This seems to be fixed by #7527. At least, I don't see any problems anymore after applying the patch.