]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Test case
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 18 Sep 2019 20:46:02 +0000 (22:46 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 12 Nov 2019 11:07:27 +0000 (12:07 +0100)
tests/matrix_free/parallel_multigrid_mf_02.cc [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_mf_02.with_mpi=true.with_p4est=true.with_lapack=true.mpirun=3.output [new file with mode: 0644]

diff --git a/tests/matrix_free/parallel_multigrid_mf_02.cc b/tests/matrix_free/parallel_multigrid_mf_02.cc
new file mode 100644 (file)
index 0000000..8577651
--- /dev/null
@@ -0,0 +1,277 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+// same test as parallel_multigrid_mf but using the Laplace operator from
+// MatrixFreeOperators::Laplace and handing the vector partitioning of those
+// operators to the MG transfer
+
+#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/matrix_free/operators.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 <typename MatrixType, typename Number>
+class MGCoarseIterative
+  : public MGCoarseGridBase<LinearAlgebra::distributed::Vector<Number>>
+{
+public:
+  MGCoarseIterative()
+  {}
+
+  void
+  initialize(const MatrixType &matrix)
+  {
+    coarse_matrix = &matrix;
+  }
+
+  virtual void
+  operator()(const unsigned int                                level,
+             LinearAlgebra::distributed::Vector<Number> &      dst,
+             const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    ReductionControl solver_control(1e4, 1e-50, 1e-10);
+    SolverCG<LinearAlgebra::distributed::Vector<Number>> solver_coarse(
+      solver_control);
+    solver_coarse.solve(*coarse_matrix, dst, src, PreconditionIdentity());
+  }
+
+  const MatrixType *coarse_matrix;
+};
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, 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;
+
+  MappingQ<dim> mapping(fe_degree + 1);
+  MatrixFreeOperators::LaplaceOperator<
+    dim,
+    fe_degree,
+    n_q_points_1d,
+    1,
+    LinearAlgebra::distributed::Vector<double>>
+    fine_matrix;
+
+  // Dirichlet BC
+  Functions::ZeroFunction<dim>                        zero_function;
+  std::map<types::boundary_id, const Function<dim> *> dirichlet_boundary;
+  dirichlet_boundary[0] = &zero_function;
+
+  // fine-level constraints
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs);
+  AffineConstraints<double> constraints;
+  constraints.reinit(locally_relevant_dofs);
+  VectorTools::interpolate_boundary_values(dof,
+                                           dirichlet_boundary,
+                                           constraints);
+  constraints.close();
+
+  std::shared_ptr<MatrixFree<dim, double>> fine_level_data(
+    new MatrixFree<dim, double>());
+
+  typename MatrixFree<dim, double>::AdditionalData fine_level_additional_data;
+  fine_level_data->reinit(mapping,
+                          dof,
+                          constraints,
+                          QGauss<1>(n_q_points_1d),
+                          fine_level_additional_data);
+
+  fine_matrix.initialize(fine_level_data);
+
+  LinearAlgebra::distributed::Vector<double> in, sol;
+  fine_matrix.initialize_dof_vector(in);
+  fine_matrix.initialize_dof_vector(sol);
+
+  // set constant rhs vector
+  in = 1.;
+
+  // set up multigrid in analogy to step-37
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof, dirichlet_boundary);
+
+  typedef MatrixFreeOperators::LaplaceOperator<
+    dim,
+    fe_degree,
+    n_q_points_1d,
+    1,
+    LinearAlgebra::distributed::Vector<number>>
+    LevelMatrixType;
+
+  MGLevelObject<LevelMatrixType>                          mg_matrices;
+  MGLevelObject<std::shared_ptr<MatrixFree<dim, number>>> mg_level_data;
+  mg_matrices.resize(0, dof.get_triangulation().n_global_levels() - 1);
+  mg_level_data.resize(0, dof.get_triangulation().n_global_levels() - 1);
+  for (unsigned int level = 0;
+       level < dof.get_triangulation().n_global_levels();
+       ++level)
+    {
+      typename MatrixFree<dim, number>::AdditionalData mg_additional_data;
+      mg_additional_data.tasks_parallel_scheme =
+        MatrixFree<dim, number>::AdditionalData::none;
+      mg_additional_data.tasks_block_size = 3;
+      mg_additional_data.level_mg_handler = level;
+
+      AffineConstraints<double> level_constraints;
+      IndexSet                  relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof, level, relevant_dofs);
+      level_constraints.reinit(relevant_dofs);
+      level_constraints.add_lines(
+        mg_constrained_dofs.get_boundary_indices(level));
+      level_constraints.close();
+
+      mg_level_data[level].reset(new MatrixFree<dim, number>());
+      mg_level_data[level]->reinit(mapping,
+                                   dof,
+                                   level_constraints,
+                                   QGauss<1>(n_q_points_1d),
+                                   mg_additional_data);
+      mg_matrices[level].initialize(mg_level_data[level],
+                                    mg_constrained_dofs,
+                                    level);
+      mg_matrices[level].compute_diagonal();
+    }
+
+  MGTransferMatrixFree<dim, number>                               mg_transfer;
+  std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>> partitioners(
+    dof.get_triangulation().n_global_levels());
+  for (unsigned int level = 0; level < partitioners.size(); ++level)
+    partitioners[level] =
+      mg_level_data[level]->get_dof_info().vector_partitioner;
+  mg_transfer.build(dof, partitioners);
+
+  MGCoarseIterative<LevelMatrixType, number> mg_coarse;
+  mg_coarse.initialize(mg_matrices[0]);
+
+  typedef PreconditionChebyshev<LevelMatrixType,
+                                LinearAlgebra::distributed::Vector<number>>
+    SMOOTHER;
+  MGSmootherPrecondition<LevelMatrixType,
+                         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;
+      smoother_data[level].preconditioner =
+        mg_matrices[level].get_matrix_diagonal_inverse();
+    }
+
+  // temporarily disable deallog for the setup of the preconditioner that
+  // involves a CG solver for eigenvalue estimation
+  mg_smoother.initialize(mg_matrices, smoother_data);
+
+  mg::Matrix<LinearAlgebra::distributed::Vector<number>> mg_matrix(mg_matrices);
+
+  Multigrid<LinearAlgebra::distributed::Vector<number>> mg(
+    dof, mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
+  PreconditionMG<dim,
+                 LinearAlgebra::distributed::Vector<number>,
+                 MGTransferMatrixFree<dim, number>>
+    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, int fe_degree, typename number>
+void
+test()
+{
+  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, fe_degree, fe_degree + 1, number>(dof);
+    }
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+
+  {
+    test<2, 1, double>();
+    test<2, 2, float>();
+
+    test<3, 1, double>();
+    test<3, 2, float>();
+  }
+}
diff --git a/tests/matrix_free/parallel_multigrid_mf_02.with_mpi=true.with_p4est=true.with_lapack=true.mpirun=3.output b/tests/matrix_free/parallel_multigrid_mf_02.with_mpi=true.with_p4est=true.with_lapack=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..604fbf3
--- /dev/null
@@ -0,0 +1,33 @@
+
+DEAL::Testing FE_Q<2>(1)
+DEAL::Number of degrees of freedom: 81
+DEAL:cg::Starting value 9.00000
+DEAL:cg::Convergence step 4 value 2.07070e-08
+DEAL::Testing FE_Q<2>(1)
+DEAL::Number of degrees of freedom: 289
+DEAL:cg::Starting value 17.0000
+DEAL:cg::Convergence step 4 value 2.77789e-08
+DEAL::Testing FE_Q<2>(2)
+DEAL::Number of degrees of freedom: 289
+DEAL:cg::Starting value 17.0000
+DEAL:cg::Convergence step 4 value 3.93059e-08
+DEAL::Testing FE_Q<2>(2)
+DEAL::Number of degrees of freedom: 1089
+DEAL:cg::Starting value 33.0000
+DEAL:cg::Convergence step 4 value 8.11904e-08
+DEAL::Testing FE_Q<3>(1)
+DEAL::Number of degrees of freedom: 125
+DEAL:cg::Starting value 11.1803
+DEAL:cg::Convergence step 4 value 3.45976e-07
+DEAL::Testing FE_Q<3>(1)
+DEAL::Number of degrees of freedom: 729
+DEAL:cg::Starting value 27.0000
+DEAL:cg::Convergence step 4 value 1.07425e-07
+DEAL::Testing FE_Q<3>(2)
+DEAL::Number of degrees of freedom: 729
+DEAL:cg::Starting value 27.0000
+DEAL:cg::Convergence step 4 value 2.46943e-06
+DEAL::Testing FE_Q<3>(2)
+DEAL::Number of degrees of freedom: 4913
+DEAL:cg::Starting value 70.0928
+DEAL:cg::Convergence step 4 value 4.57526e-06

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.