]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New test case
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 12 Oct 2021 13:09:53 +0000 (15:09 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 19 Oct 2021 12:35:46 +0000 (14:35 +0200)
tests/matrix_free/parallel_multigrid_interleave.cc [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output [new file with mode: 0644]

diff --git a/tests/matrix_free/parallel_multigrid_interleave.cc b/tests/matrix_free/parallel_multigrid_interleave.cc
new file mode 100644 (file)
index 0000000..20b59fc
--- /dev/null
@@ -0,0 +1,473 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2021 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 test as parallel_multigrid_mf, but using the functionality of
+// PreconditionChebyshev to embed vector updates into the matrix-free loops of
+// a suitable operator class.
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_matrix.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/multigrid.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+template <int dim, typename number = double>
+class LaplaceOperator : public Subscriptor
+{
+public:
+  using value_type = number;
+
+  LaplaceOperator()
+    : n_calls_vmult(0)
+  {}
+
+  ~LaplaceOperator()
+  {
+    if (n_calls_vmult > 0)
+      deallog << "Number of calls to special vmult for Operator of size " << m()
+              << ": " << n_calls_vmult << std::endl;
+  }
+
+  void
+  initialize(const Mapping<dim> &                mapping,
+             const DoFHandler<dim> &             dof_handler,
+             const std::set<types::boundary_id> &dirichlet_boundaries,
+             const unsigned int level = numbers::invalid_unsigned_int)
+  {
+    n_calls_vmult = 0;
+    const QGauss<1> quad(dof_handler.get_fe().degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData addit_data;
+    addit_data.tasks_parallel_scheme =
+      MatrixFree<dim, number>::AdditionalData::none;
+    addit_data.mg_level = level;
+
+    // extract the constraints due to Dirichlet boundary conditions
+    AffineConstraints<double>                           constraints;
+    Functions::ZeroFunction<dim>                        zero;
+    std::map<types::boundary_id, const Function<dim> *> functions;
+    for (std::set<types::boundary_id>::const_iterator it =
+           dirichlet_boundaries.begin();
+         it != dirichlet_boundaries.end();
+         ++it)
+      functions[*it] = &zero;
+    if (level == numbers::invalid_unsigned_int)
+      VectorTools::interpolate_boundary_values(dof_handler,
+                                               functions,
+                                               constraints);
+    else
+      {
+        std::vector<types::global_dof_index>    local_dofs;
+        typename DoFHandler<dim>::cell_iterator cell = dof_handler.begin(level),
+                                                endc = dof_handler.end(level);
+        for (; cell != endc; ++cell)
+          {
+            if (dof_handler.get_triangulation().locally_owned_subdomain() !=
+                  numbers::invalid_subdomain_id &&
+                cell->level_subdomain_id() == numbers::artificial_subdomain_id)
+              continue;
+            const FiniteElement<dim> &fe = cell->get_fe();
+            local_dofs.resize(fe.dofs_per_face);
+
+            for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
+              if (cell->at_boundary(face_no) == true)
+                {
+                  const typename DoFHandler<dim>::face_iterator face =
+                    cell->face(face_no);
+                  const types::boundary_id bi = face->boundary_id();
+                  if (functions.find(bi) != functions.end())
+                    {
+                      face->get_mg_dof_indices(level, local_dofs);
+                      for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+                        constraints.add_line(local_dofs[i]);
+                    }
+                }
+          }
+      }
+    constraints.close();
+
+    data.reinit(mapping, dof_handler, constraints, quad, addit_data);
+
+    compute_inverse_diagonal();
+  }
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<number> &      dst,
+        const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    data.cell_loop(&LaplaceOperator::local_apply, this, dst, src, true);
+    for (unsigned int i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<number> &      dst,
+        const LinearAlgebra::distributed::Vector<number> &src,
+        const std::function<void(const unsigned int, const unsigned int)>
+          &operation_before_loop,
+        const std::function<void(const unsigned int, const unsigned int)>
+          &operation_after_loop) const
+  {
+    ++n_calls_vmult;
+    data.cell_loop(&LaplaceOperator::local_apply,
+                   this,
+                   dst,
+                   src,
+                   operation_before_loop,
+                   operation_after_loop);
+    for (unsigned int i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  Tvmult(LinearAlgebra::distributed::Vector<number> &      dst,
+         const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    vmult(dst, src);
+  }
+
+  void
+  vmult_add(LinearAlgebra::distributed::Vector<number> &      dst,
+            const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    data.cell_loop(&LaplaceOperator::local_apply, this, dst, src);
+    for (unsigned int i : data.get_constrained_dofs())
+      dst.local_element(i) += src.local_element(i);
+  }
+
+  void
+  Tvmult_add(LinearAlgebra::distributed::Vector<number> &      dst,
+             const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    vmult_add(dst, src);
+  }
+
+  types::global_dof_index
+  m() const
+  {
+    return data.get_vector_partitioner()->size();
+  }
+
+  types::global_dof_index
+  n() const
+  {
+    return data.get_vector_partitioner()->size();
+  }
+
+  number
+  el(const unsigned int row, const unsigned int col) const
+  {
+    AssertThrow(false,
+                ExcMessage("Matrix-free does not allow for entry access"));
+    return number();
+  }
+
+  void
+  initialize_dof_vector(
+    LinearAlgebra::distributed::Vector<number> &vector) const
+  {
+    if (!vector.partitioners_are_compatible(
+          *data.get_dof_info(0).vector_partitioner))
+      data.initialize_dof_vector(vector);
+    Assert(vector.partitioners_are_globally_compatible(
+             *data.get_dof_info(0).vector_partitioner),
+           ExcInternalError());
+  }
+
+  const LinearAlgebra::distributed::Vector<number> &
+  get_matrix_diagonal_inverse() const
+  {
+    Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
+    return inverse_diagonal_entries;
+  }
+
+
+private:
+  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
+  {
+    FEEvaluation<dim, -1, 0, 1, number> phi(data);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values(src);
+        phi.evaluate(EvaluationFlags::gradients);
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          phi.submit_gradient(phi.get_gradient(q), q);
+        phi.integrate(EvaluationFlags::gradients);
+        phi.distribute_local_to_global(dst);
+      }
+  }
+
+  void
+  compute_inverse_diagonal()
+  {
+    data.initialize_dof_vector(inverse_diagonal_entries);
+    unsigned int dummy = 0;
+    data.cell_loop(&LaplaceOperator::local_diagonal_cell,
+                   this,
+                   inverse_diagonal_entries,
+                   dummy);
+
+    for (unsigned int i = 0; i < inverse_diagonal_entries.local_size(); ++i)
+      if (std::abs(inverse_diagonal_entries.local_element(i)) > 1e-10)
+        inverse_diagonal_entries.local_element(i) =
+          1. / inverse_diagonal_entries.local_element(i);
+      else
+        inverse_diagonal_entries.local_element(i) = 1.;
+  }
+
+  void
+  local_diagonal_cell(
+    const MatrixFree<dim, number> &             data,
+    LinearAlgebra::distributed::Vector<number> &dst,
+    const unsigned int &,
+    const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, -1, 0, 1, number> phi(data);
+
+    AlignedVector<VectorizedArray<number>> local_diagonal_vector(
+      phi.dofs_per_cell);
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        phi.reinit(cell);
+
+        for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < phi.dofs_per_cell; ++j)
+              phi.begin_dof_values()[j] = VectorizedArray<number>();
+            phi.begin_dof_values()[i] = 1.;
+            phi.evaluate(EvaluationFlags::gradients);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_gradient(phi.get_gradient(q), q);
+            phi.integrate(EvaluationFlags::gradients);
+            local_diagonal_vector[i] = phi.begin_dof_values()[i];
+          }
+        for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+          phi.begin_dof_values()[i] = local_diagonal_vector[i];
+        phi.distribute_local_to_global(dst);
+      }
+  }
+
+  MatrixFree<dim, number>                    data;
+  LinearAlgebra::distributed::Vector<number> inverse_diagonal_entries;
+  mutable unsigned int                       n_calls_vmult;
+};
+
+
+
+template <int dim, typename MatrixType>
+class MGTransferPrebuiltMF
+  : public MGTransferMatrixFree<dim, typename MatrixType::value_type>
+{
+public:
+  MGTransferPrebuiltMF(const MGLevelObject<MatrixType> &laplace,
+                       const MGConstrainedDoFs &        mg_constrained_dofs)
+    : MGTransferMatrixFree<dim, typename MatrixType::value_type>(
+        mg_constrained_dofs)
+    , laplace_operator(laplace){};
+
+  /**
+   * Overload copy_to_mg from MGTransferPrebuilt to get the vectors compatible
+   * with MatrixFree and bypass the crude initialization in MGTransferPrebuilt
+   */
+  template <class InVector, int spacedim>
+  void
+  copy_to_mg(
+    const DoFHandler<dim, spacedim> &mg_dof,
+    MGLevelObject<
+      LinearAlgebra::distributed::Vector<typename MatrixType::value_type>> &dst,
+    const InVector &src) const
+  {
+    for (unsigned int level = dst.min_level(); level <= dst.max_level();
+         ++level)
+      laplace_operator[level].initialize_dof_vector(dst[level]);
+    MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<
+      typename MatrixType::value_type>>::copy_to_mg(mg_dof, dst, src);
+  }
+
+private:
+  const MGLevelObject<MatrixType> &laplace_operator;
+};
+
+
+
+template <int dim, typename number>
+void
+do_test(const DoFHandler<dim> &dof)
+{
+  deallog << "Testing " << dof.get_fe().get_name();
+  deallog << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+
+  const unsigned int           fe_degree = dof.get_fe().degree;
+  MappingQ<dim>                mapping(fe_degree + 1);
+  LaplaceOperator<dim, double> fine_matrix;
+  std::set<types::boundary_id> dirichlet_boundaries;
+  dirichlet_boundaries.insert(0);
+  fine_matrix.initialize(mapping, dof, dirichlet_boundaries);
+
+  LinearAlgebra::distributed::Vector<double> in, sol;
+  fine_matrix.initialize_dof_vector(in);
+  fine_matrix.initialize_dof_vector(sol);
+
+  // set constant rhs vector, except boundary
+  in = 1.;
+  {
+    std::map<types::global_dof_index, double> boundary_values;
+    VectorTools::interpolate_boundary_values(dof,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             boundary_values);
+    for (const auto it : boundary_values)
+      if (dof.locally_owned_dofs().is_element(it.first))
+        in(it.first) = 0.;
+  }
+
+  // set up multigrid in analogy to step-37
+  using LevelMatrixType = LaplaceOperator<dim, number>;
+
+  MGLevelObject<LevelMatrixType> mg_matrices;
+  mg_matrices.resize(0, dof.get_triangulation().n_global_levels() - 1);
+  for (unsigned int level = 0;
+       level < dof.get_triangulation().n_global_levels();
+       ++level)
+    {
+      mg_matrices[level].initialize(mapping, dof, dirichlet_boundaries, level);
+    }
+
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof);
+  mg_constrained_dofs.make_zero_boundary_constraints(dof, {0});
+
+  MGTransferPrebuiltMF<dim, LevelMatrixType> mg_transfer(mg_matrices,
+                                                         mg_constrained_dofs);
+  mg_transfer.build(dof);
+
+  using SMOOTHER =
+    PreconditionChebyshev<LevelMatrixType,
+                          LinearAlgebra::distributed::Vector<number>>;
+  mg::SmootherRelaxation<SMOOTHER, LinearAlgebra::distributed::Vector<number>>
+    mg_smoother;
+
+  MGLevelObject<typename SMOOTHER::AdditionalData> smoother_data;
+  smoother_data.resize(0, dof.get_triangulation().n_global_levels() - 1);
+  for (unsigned int level = 0;
+       level < dof.get_triangulation().n_global_levels();
+       ++level)
+    {
+      smoother_data[level].smoothing_range     = 15.;
+      smoother_data[level].degree              = 5;
+      smoother_data[level].eig_cg_n_iterations = 15;
+      auto preconditioner                      = std::make_shared<
+        DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>>();
+      preconditioner->reinit(mg_matrices[level].get_matrix_diagonal_inverse());
+      smoother_data[level].preconditioner = std::move(preconditioner);
+    }
+
+  mg_smoother.initialize(mg_matrices, smoother_data);
+  MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<number>>
+    mg_coarse(mg_smoother);
+
+  mg::Matrix<LinearAlgebra::distributed::Vector<number>> mg_matrix(mg_matrices);
+
+  Multigrid<LinearAlgebra::distributed::Vector<number>> mg(
+    mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
+  PreconditionMG<dim,
+                 LinearAlgebra::distributed::Vector<number>,
+                 MGTransferPrebuiltMF<dim, LevelMatrixType>>
+    preconditioner(dof, mg, mg_transfer);
+
+  {
+    // avoid output from inner (coarse-level) solver
+    deallog.depth_file(2);
+    ReductionControl control(30, 1e-20, 1e-7);
+    SolverCG<LinearAlgebra::distributed::Vector<double>> solver(control);
+    solver.solve(fine_matrix, sol, in, preconditioner);
+  }
+}
+
+
+
+template <int dim, typename number>
+void
+test(const unsigned int fe_degree)
+{
+  for (unsigned int i = 5; i < 7; ++i)
+    {
+      parallel::distributed::Triangulation<dim> tria(
+        MPI_COMM_WORLD,
+        Triangulation<dim>::limit_level_difference_at_vertices,
+        parallel::distributed::Triangulation<
+          dim>::construct_multigrid_hierarchy);
+      GridGenerator::hyper_cube(tria);
+      tria.refine_global(i - dim);
+
+      FE_Q<dim>       fe(fe_degree);
+      DoFHandler<dim> dof(tria);
+      dof.distribute_dofs(fe);
+      dof.distribute_mg_dofs();
+
+      do_test<dim, number>(dof);
+    }
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+
+  {
+    test<2, double>(1);
+    test<2, float>(3);
+
+    test<3, double>(1);
+    test<3, float>(2);
+  }
+}
diff --git a/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..34c42cc
--- /dev/null
@@ -0,0 +1,65 @@
+
+DEAL::Testing FE_Q<2>(1)
+DEAL::Number of degrees of freedom: 81
+DEAL:cg::Starting value 7.00000
+DEAL:cg::Convergence step 3 value 1.44113e-07
+DEAL::Number of calls to special vmult for Operator of size 4: 15
+DEAL::Number of calls to special vmult for Operator of size 9: 27
+DEAL::Number of calls to special vmult for Operator of size 25: 27
+DEAL::Number of calls to special vmult for Operator of size 81: 27
+DEAL::Testing FE_Q<2>(1)
+DEAL::Number of degrees of freedom: 289
+DEAL:cg::Starting value 15.0000
+DEAL:cg::Convergence step 3 value 1.49408e-06
+DEAL::Number of calls to special vmult for Operator of size 4: 15
+DEAL::Number of calls to special vmult for Operator of size 9: 27
+DEAL::Number of calls to special vmult for Operator of size 25: 27
+DEAL::Number of calls to special vmult for Operator of size 81: 27
+DEAL::Number of calls to special vmult for Operator of size 289: 27
+DEAL::Testing FE_Q<2>(3)
+DEAL::Number of degrees of freedom: 625
+DEAL:cg::Starting value 23.0000
+DEAL:cg::Convergence step 4 value 6.30917e-08
+DEAL::Number of calls to special vmult for Operator of size 16: 20
+DEAL::Number of calls to special vmult for Operator of size 49: 36
+DEAL::Number of calls to special vmult for Operator of size 169: 36
+DEAL::Number of calls to special vmult for Operator of size 625: 36
+DEAL::Testing FE_Q<2>(3)
+DEAL::Number of degrees of freedom: 2401
+DEAL:cg::Starting value 47.0000
+DEAL:cg::Convergence step 4 value 1.99839e-07
+DEAL::Number of calls to special vmult for Operator of size 16: 20
+DEAL::Number of calls to special vmult for Operator of size 49: 36
+DEAL::Number of calls to special vmult for Operator of size 169: 36
+DEAL::Number of calls to special vmult for Operator of size 625: 36
+DEAL::Number of calls to special vmult for Operator of size 2401: 36
+DEAL::Testing FE_Q<3>(1)
+DEAL::Number of degrees of freedom: 125
+DEAL:cg::Starting value 5.19615
+DEAL:cg::Convergence step 3 value 5.16655e-09
+DEAL::Number of calls to special vmult for Operator of size 8: 15
+DEAL::Number of calls to special vmult for Operator of size 27: 27
+DEAL::Number of calls to special vmult for Operator of size 125: 27
+DEAL::Testing FE_Q<3>(1)
+DEAL::Number of degrees of freedom: 729
+DEAL:cg::Starting value 18.5203
+DEAL:cg::Convergence step 3 value 1.01773e-06
+DEAL::Number of calls to special vmult for Operator of size 8: 15
+DEAL::Number of calls to special vmult for Operator of size 27: 27
+DEAL::Number of calls to special vmult for Operator of size 125: 27
+DEAL::Number of calls to special vmult for Operator of size 729: 27
+DEAL::Testing FE_Q<3>(2)
+DEAL::Number of degrees of freedom: 729
+DEAL:cg::Starting value 18.5203
+DEAL:cg::Convergence step 4 value 5.65213e-09
+DEAL::Number of calls to special vmult for Operator of size 27: 20
+DEAL::Number of calls to special vmult for Operator of size 125: 36
+DEAL::Number of calls to special vmult for Operator of size 729: 36
+DEAL::Testing FE_Q<3>(2)
+DEAL::Number of degrees of freedom: 4913
+DEAL:cg::Starting value 58.0948
+DEAL:cg::Convergence step 4 value 3.37720e-08
+DEAL::Number of calls to special vmult for Operator of size 27: 20
+DEAL::Number of calls to special vmult for Operator of size 125: 36
+DEAL::Number of calls to special vmult for Operator of size 729: 36
+DEAL::Number of calls to special vmult for Operator of size 4913: 36

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.