]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix DoFInfo::assign_ghosts for hp in parallel 11060/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 20 Oct 2020 11:15:16 +0000 (13:15 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 21 Oct 2020 04:44:55 +0000 (06:44 +0200)
include/deal.II/matrix_free/dof_info.templates.h
tests/matrix_free/matrix_vector_hp_parallel_01.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=1.with_trilinos=true.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=5.with_trilinos=true.output [new file with mode: 0644]

index f8d115f6d89406ae08c5325ccb666605cad395e3..c94c42819ca0b2516cb7bdce4fd9c8fe1b96122f 100644 (file)
@@ -473,7 +473,7 @@ namespace internal
                         (cell_active_fe_index.size() == 0 ||
                          dofs_per_cell.size() == 1) ?
                           0 :
-                          cell_active_fe_index[i];
+                          cell_active_fe_index[boundary_cells[i]];
                       AssertIndexRange(fe_index, dofs_per_cell.size());
                       const unsigned int *row_end =
                         data_ptr + dofs_per_cell[fe_index];
diff --git a/tests/matrix_free/matrix_vector_hp_parallel_01.cc b/tests/matrix_free/matrix_vector_hp_parallel_01.cc
new file mode 100644 (file)
index 0000000..53f81eb
--- /dev/null
@@ -0,0 +1,293 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Similar to matrix_vector_hp.cc but in parallel.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim, typename Number>
+class MatrixFreeTestHP
+{
+public:
+  MatrixFreeTestHP(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  void
+  local_apply(const MatrixFree<dim, Number> &                   data,
+              LinearAlgebra::distributed::Vector<Number> &      dst,
+              const LinearAlgebra::distributed::Vector<Number> &src,
+              const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    // ask MatrixFree for cell_range for different orders
+    std::pair<unsigned int, unsigned int> subrange_deg =
+      data.create_cell_subrange_hp(cell_range, 1);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 1, LinearAlgebra::distributed::Vector<Number>, 2>(
+        data, dst, src, subrange_deg);
+    subrange_deg = data.create_cell_subrange_hp(cell_range, 2);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 2, LinearAlgebra::distributed::Vector<Number>, 3>(
+        data, dst, src, subrange_deg);
+    subrange_deg = data.create_cell_subrange_hp(cell_range, 3);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 3, LinearAlgebra::distributed::Vector<Number>, 4>(
+        data, dst, src, subrange_deg);
+    subrange_deg = data.create_cell_subrange_hp(cell_range, 4);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 4, LinearAlgebra::distributed::Vector<Number>, 5>(
+        data, dst, src, subrange_deg);
+    subrange_deg = data.create_cell_subrange_hp(cell_range, 5);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 5, LinearAlgebra::distributed::Vector<Number>, 6>(
+        data, dst, src, subrange_deg);
+    subrange_deg = data.create_cell_subrange_hp(cell_range, 6);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 6, LinearAlgebra::distributed::Vector<Number>, 7>(
+        data, dst, src, subrange_deg);
+    subrange_deg = data.create_cell_subrange_hp(cell_range, 7);
+    if (subrange_deg.second > subrange_deg.first)
+      helmholtz_operator<dim, 7, LinearAlgebra::distributed::Vector<Number>, 8>(
+        data, dst, src, subrange_deg);
+  }
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    dst = 0;
+    data.cell_loop(&MatrixFreeTestHP<dim, Number>::local_apply, this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  if (fe_degree > 1)
+    return;
+
+  typedef double                            number;
+  const SphericalManifold<dim>              manifold;
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_ball(tria);
+  typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
+                                                    endc = tria.end();
+  for (; cell != endc; ++cell)
+    for (const unsigned int f : GeometryInfo<dim>::face_indices())
+      if (cell->at_boundary(f))
+        cell->face(f)->set_all_manifold_ids(0);
+  tria.set_manifold(0, manifold);
+  tria.refine_global(1);
+
+  // refine a few cells
+  for (unsigned int i = 0; i < 11 - 3 * dim; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator cell =
+                                                          tria.begin_active(),
+                                                        endc = tria.end();
+      unsigned int counter                                   = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (cell->is_locally_owned() && (counter % (7 - i) == 0))
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  const unsigned int max_degree = 9 - 2 * dim;
+
+  hp::FECollection<dim> fe_collection;
+  hp::QCollection<dim>  quadrature_collection;
+  hp::QCollection<1>    quadrature_collection_mf;
+
+  for (unsigned int deg = 1; deg <= max_degree; ++deg)
+    {
+      fe_collection.push_back(FE_Q<dim>(QGaussLobatto<1>(deg + 1)));
+      quadrature_collection.push_back(QGauss<dim>(deg + 1));
+      quadrature_collection_mf.push_back(QGauss<1>(deg + 1));
+    }
+
+  hp::DoFHandler<dim> dof(tria);
+  // set the active FE index in a random order
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator cell =
+                                                         dof.begin_active(),
+                                                       endc = dof.end();
+    for (; cell != endc; ++cell)
+      {
+        if (cell->is_locally_owned() == false)
+          continue;
+
+        const unsigned int fe_index = Testing::rand() % max_degree;
+        cell->set_active_fe_index(fe_index);
+      }
+  }
+
+  // setup DoFs
+  dof.distribute_dofs(fe_collection);
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values(dof,
+                                           0,
+                                           Functions::ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close();
+
+  TrilinosWrappers::SparsityPattern dsp(dof.locally_owned_dofs(),
+                                        MPI_COMM_WORLD);
+  DoFTools::make_sparsity_pattern(dof, dsp, constraints, false);
+  dsp.compress();
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+  system_matrix.reinit(dsp);
+
+  // std::cout << "Number of cells: " <<
+  // dof.get_triangulation().n_active_cells() << std::endl; std::cout << "Number
+  // of degrees of freedom: " << dof.n_dofs() << std::endl; std::cout << "Number
+  // of constraints: " << constraints.n_constraints() << std::endl;
+
+  // set up MatrixFree
+  MatrixFree<dim, number>                          mf_data;
+  typename MatrixFree<dim, number>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+  mf_data.reinit(dof, constraints, quadrature_collection_mf, data);
+  MatrixFreeTestHP<dim, number> mf(mf_data);
+
+  // assemble sparse matrix with (\nabla v,
+  // \nabla u) + (v, 10 * u)
+  {
+    hp::FEValues<dim>                    hp_fe_values(fe_collection,
+                                   quadrature_collection,
+                                   update_values | update_gradients |
+                                     update_JxW_values);
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+
+    typename hp::DoFHandler<dim>::active_cell_iterator cell =
+                                                         dof.begin_active(),
+                                                       endc = dof.end();
+    for (; cell != endc; ++cell)
+      {
+        if (cell->is_locally_owned() == false)
+          continue;
+
+        const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+        cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+        cell_matrix = 0;
+        hp_fe_values.reinit(cell);
+        const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+
+        for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+             ++q_point)
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                cell_matrix(i, j) += ((fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) +
+                                       10. * fe_values.shape_value(i, q_point) *
+                                         fe_values.shape_value(j, q_point)) *
+                                      fe_values.JxW(q_point));
+            }
+        local_dof_indices.resize(dofs_per_cell);
+        cell->get_dof_indices(local_dof_indices);
+
+        constraints.distribute_local_to_global(cell_matrix,
+                                               local_dof_indices,
+                                               system_matrix);
+      }
+  }
+  system_matrix.compress(VectorOperation::values::add);
+
+  // fill a right hand side vector with random
+  // numbers in unconstrained degrees of freedom
+  LinearAlgebra::distributed::Vector<double> src;
+  mf_data.initialize_dof_vector(src);
+
+
+  LinearAlgebra::distributed::Vector<double> result_spmv(src), result_mf(src);
+
+  for (unsigned int i = 0; i < dof.n_locally_owned_dofs(); ++i)
+    {
+      if (constraints.is_constrained(
+            src.get_partitioner()->local_to_global(i)) == false)
+        src.local_element(i) = random_value<double>();
+    }
+
+  // now perform matrix-vector product and check
+  // its correctness
+  system_matrix.vmult(result_spmv, src);
+  mf.vmult(result_mf, src);
+
+  result_mf -= result_spmv;
+  const double diff_norm = result_mf.linfty_norm();
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+  mpi_initlog();
+
+  {
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    test<3, 2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=1.with_trilinos=true.output b/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=1.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..7252835
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=5.with_trilinos=true.output b/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=5.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..7252835
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Norm of difference: 0
+DEAL:3d::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.