]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add test for p-multigrid 11433/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 29 Dec 2020 07:50:32 +0000 (08:50 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 29 Dec 2020 07:50:32 +0000 (08:50 +0100)
include/deal.II/matrix_free/tools.h
tests/multigrid-global-coarsening/multigrid_p_01.cc [new file with mode: 0644]
tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output [new file with mode: 0644]
tests/multigrid-global-coarsening/multigrid_util.h [new file with mode: 0644]

index d0b7f9d435a8ebdccb695c4a8682d705800c56a7..99498540c7597a29bd475185a74c9558673469f5 100644 (file)
@@ -91,7 +91,7 @@ namespace MatrixFreeTools
                                                n_components,
                                                Number,
                                                VectorizedArrayType> &) const,
-    CLASS *            owning_class,
+    const CLASS *      owning_class,
     const unsigned int dof_no                   = 0,
     const unsigned int quad_no                  = 0,
     const unsigned int first_selected_component = 0);
@@ -609,7 +609,7 @@ namespace MatrixFreeTools
                                                n_components,
                                                Number,
                                                VectorizedArrayType> &) const,
-    CLASS *            owning_class,
+    const CLASS *      owning_class,
     const unsigned int dof_no,
     const unsigned int quad_no,
     const unsigned int first_selected_component)
diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.cc b/tests/multigrid-global-coarsening/multigrid_p_01.cc
new file mode 100644 (file)
index 0000000..623d35b
--- /dev/null
@@ -0,0 +1,161 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+/**
+ * Test p-multigrid for a uniformly refined mesh both for simplex and
+ * hypercube mesh.
+ */
+
+#include "multigrid_util.h"
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int n_refinements,
+     const unsigned int fe_degree_fine,
+     const bool         do_simplex_mesh)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  Triangulation<dim> tria;
+  if (do_simplex_mesh)
+    GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2);
+  else
+    GridGenerator::subdivided_hyper_cube(tria, 2);
+  tria.refine_global(n_refinements);
+
+  const auto level_degrees = MGTransferGlobalCoarseningTools::create_p_sequence(
+    fe_degree_fine,
+    MGTransferGlobalCoarseningTools::PolynomialSequenceType::bisect);
+
+  const unsigned int min_level = 0;
+  const unsigned int max_level = level_degrees.size() - 1;
+
+  MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level, tria);
+  MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
+  MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> transfers(min_level,
+                                                               max_level);
+  MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
+
+  std::unique_ptr<Mapping<dim>> mapping_;
+
+  // set up levels
+  for (auto l = min_level; l <= max_level; ++l)
+    {
+      auto &dof_handler = dof_handlers[l];
+      auto &constraint  = constraints[l];
+      auto &op          = operators[l];
+
+      std::unique_ptr<FiniteElement<dim>> fe;
+      std::unique_ptr<Quadrature<dim>>    quad;
+      std::unique_ptr<Mapping<dim>>       mapping;
+
+      if (do_simplex_mesh)
+        {
+          fe   = std::make_unique<Simplex::FE_P<dim>>(level_degrees[l]);
+          quad = std::make_unique<Simplex::QGauss<dim>>(level_degrees[l] + 1);
+          mapping = std::make_unique<MappingFE<dim>>(Simplex::FE_P<dim>(1));
+        }
+      else
+        {
+          fe      = std::make_unique<FE_Q<dim>>(level_degrees[l]);
+          quad    = std::make_unique<QGauss<dim>>(level_degrees[l] + 1);
+          mapping = std::make_unique<MappingFE<dim>>(FE_Q<dim>(1));
+        }
+
+      if (l == max_level)
+        mapping_ = mapping->clone();
+
+      // set up dofhandler
+      dof_handler.distribute_dofs(*fe);
+
+      // set up constraints
+      IndexSet locally_relevant_dofs;
+      DoFTools::extract_locally_relevant_dofs(dof_handler,
+                                              locally_relevant_dofs);
+      constraint.reinit(locally_relevant_dofs);
+      VectorTools::interpolate_boundary_values(
+        *mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraint);
+      constraint.close();
+
+      // set up operator
+      op.reinit(*mapping, dof_handler, *quad, constraint);
+    }
+
+  // set up transfer operator
+  for (unsigned int l = min_level; l < max_level; ++l)
+    transfers[l + 1].reinit_polynomial_transfer(dof_handlers[l + 1],
+                                                dof_handlers[l],
+                                                constraints[l + 1],
+                                                constraints[l]);
+
+  MGTransferGlobalCoarsening<Operator<dim, Number>, VectorType> transfer(
+    operators, transfers);
+
+  GMGParameters mg_data; // TODO
+
+  VectorType dst, src;
+  operators[max_level].initialize_dof_vector(dst);
+  operators[max_level].initialize_dof_vector(src);
+
+  operators[max_level].rhs(src);
+
+  ReductionControl solver_control(
+    mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false);
+
+  mg_solve(solver_control,
+           dst,
+           src,
+           mg_data,
+           dof_handlers[max_level],
+           operators[max_level],
+           operators,
+           transfer);
+
+  deallog << dim << " " << fe_degree_fine << " " << n_refinements << " "
+          << (do_simplex_mesh ? "tri " : "quad") << " "
+          << solver_control.last_step() << std::endl;
+
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler(dof_handlers[max_level]);
+  data_out.add_data_vector(dst, "solution");
+  data_out.build_patches(*mapping_, 2);
+
+  static unsigned int counter = 0;
+  std::ofstream       output("test." + std::to_string(dim) + "." +
+                       std::to_string(counter++) + ".vtk");
+  data_out.write_vtk(output);
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  deallog.precision(8);
+
+  for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+    for (unsigned int degree = 2; degree <= 4; ++degree)
+      test<2>(n_refinements, degree, false /*quadrilateral*/);
+
+  return 0; // TODO: enable simplex test once MGTwoLevelTransfer works for
+            // simplex meshes
+
+  for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+    for (unsigned int degree = 2; degree <= 2; ++degree)
+      test<2>(n_refinements, degree, true /*triangle*/);
+}
diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..3031424
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL:0::2 2 2 quad 3
+DEAL:0::2 3 2 quad 3
+DEAL:0::2 4 2 quad 3
+DEAL:0::2 2 3 quad 3
+DEAL:0::2 3 3 quad 3
+DEAL:0::2 4 3 quad 3
+DEAL:0::2 2 4 quad 3
+DEAL:0::2 3 4 quad 3
+DEAL:0::2 4 4 quad 3
diff --git a/tests/multigrid-global-coarsening/multigrid_util.h b/tests/multigrid-global-coarsening/multigrid_util.h
new file mode 100644 (file)
index 0000000..e3176b6
--- /dev/null
@@ -0,0 +1,465 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_tests_multigrid_util_h
+#define dealii_tests_multigrid_util_h
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_constrained_dofs.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_global_coarsening.h>
+#include <deal.II/multigrid/multigrid.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+namespace MGTransferGlobalCoarseningTools
+{
+  enum class PolynomialSequenceType
+  {
+    bisect,
+    decrease_by_one,
+    go_to_one
+  };
+
+  unsigned int
+  generate_level_degree(const unsigned int            previous_fe_degree,
+                        const PolynomialSequenceType &p_sequence)
+  {
+    switch (p_sequence)
+      {
+        case PolynomialSequenceType::bisect:
+          return std::max(previous_fe_degree / 2, 1u);
+        case PolynomialSequenceType::decrease_by_one:
+          return std::max(previous_fe_degree - 1, 1u);
+        case PolynomialSequenceType::go_to_one:
+          return 1u;
+        default:
+          Assert(false, StandardExceptions::ExcNotImplemented());
+          return 1u;
+      }
+  }
+
+  std::vector<unsigned int>
+  create_p_sequence(const unsigned int            degree,
+                    const PolynomialSequenceType &p_sequence)
+  {
+    std::vector<unsigned int> degrees;
+    degrees.push_back(degree);
+
+    unsigned int previous_fe_degree = degree;
+    while (previous_fe_degree > 1)
+      {
+        const unsigned int level_degree =
+          generate_level_degree(previous_fe_degree, p_sequence);
+
+        degrees.push_back(level_degree);
+        previous_fe_degree = level_degree;
+      }
+
+    std::reverse(degrees.begin(), degrees.end());
+
+    return degrees;
+  }
+
+} // namespace MGTransferGlobalCoarseningTools
+
+template <int dim_, typename Number>
+class Operator : public Subscriptor
+{
+public:
+  using value_type = Number;
+  using number     = Number;
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  static const int dim = dim_;
+
+  using FECellIntegrator = FEEvaluation<dim, -1, 0, 1, number>;
+
+  void
+  reinit(const Mapping<dim> &             mapping,
+         const DoFHandler<dim> &          dof_handler,
+         const Quadrature<dim> &          quad,
+         const AffineConstraints<number> &constraints)
+  {
+    // Clear internal data structures (if operator is reused).
+    this->system_matrix.clear();
+
+    // Copy the constrains, since they might be needed for computation of the
+    // system matrix later on.
+    this->constraints.copy_from(constraints);
+
+    // Set up MatrixFree. At the quadrature points, we only need to evaluate
+    // the gradient of the solution and test with the gradient of the shape
+    // functions so that we only need to set the flag `update_gradients`.
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.mapping_update_flags = update_gradients;
+
+    matrix_free.reinit(mapping, dof_handler, constraints, quad, data);
+  }
+
+  virtual types::global_dof_index
+  m() const
+  {
+    return matrix_free.get_dof_handler().n_dofs();
+  }
+
+  Number
+  el(unsigned int, unsigned int) const
+  {
+    Assert(false, ExcNotImplemented());
+    return 0;
+  }
+
+  virtual void
+  initialize_dof_vector(VectorType &vec) const
+  {
+    matrix_free.initialize_dof_vector(vec);
+  }
+
+  virtual void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    this->matrix_free.cell_loop(
+      &Operator::do_cell_integral_range, this, dst, src, true);
+  }
+
+  void
+  Tvmult(VectorType &dst, const VectorType &src) const
+  {
+    vmult(dst, src);
+  }
+
+  void
+  compute_inverse_diagonal(VectorType &diagonal) const
+  {
+    // compute diagonal
+    MatrixFreeTools::compute_diagonal(matrix_free,
+                                      diagonal,
+                                      &Operator::do_cell_integral_local,
+                                      this);
+
+    // and invert it
+    for (auto &i : diagonal)
+      i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0;
+  }
+
+  virtual const TrilinosWrappers::SparseMatrix &
+  get_system_matrix() const
+  {
+    // Check if matrix has already been set up.
+    if (system_matrix.m() == 0 && system_matrix.n() == 0)
+      {
+        // Set up sparsity pattern of system matrix.
+        const auto &dof_handler = this->matrix_free.get_dof_handler();
+
+        TrilinosWrappers::SparsityPattern dsp(
+          dof_handler.locally_owned_dofs(),
+          matrix_free.get_task_info().communicator);
+
+        DoFTools::make_sparsity_pattern(dof_handler, dsp, this->constraints);
+
+        dsp.compress();
+        system_matrix.reinit(dsp);
+
+        // Assemble system matrix.
+        MatrixFreeTools::compute_matrix(matrix_free,
+                                        constraints,
+                                        system_matrix,
+                                        &Operator::do_cell_integral_local,
+                                        this);
+      }
+
+    return this->system_matrix;
+  }
+
+  void
+  rhs(VectorType &b) const
+  {
+    const int dummy = 0;
+
+    matrix_free.template cell_loop<VectorType, int>(
+      [](const auto &matrix_free, auto &dst, const auto &, const auto cells) {
+        FECellIntegrator phi(matrix_free, cells);
+        for (unsigned int cell = cells.first; cell < cells.second; ++cell)
+          {
+            phi.reinit(cell);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_value(1.0, q);
+
+            phi.integrate_scatter(EvaluationFlags::values, dst);
+          }
+      },
+      b,
+      dummy,
+      true);
+  }
+
+private:
+  void
+  do_cell_integral_local(FECellIntegrator &integrator) const
+  {
+    integrator.evaluate(EvaluationFlags::gradients);
+
+    for (unsigned int q = 0; q < integrator.n_q_points; ++q)
+      integrator.submit_gradient(integrator.get_gradient(q), q);
+
+    integrator.integrate(EvaluationFlags::gradients);
+  }
+
+  void
+  do_cell_integral_global(FECellIntegrator &integrator,
+                          VectorType &      dst,
+                          const VectorType &src) const
+  {
+    integrator.gather_evaluate(src, EvaluationFlags::gradients);
+
+    for (unsigned int q = 0; q < integrator.n_q_points; ++q)
+      integrator.submit_gradient(integrator.get_gradient(q), q);
+
+    integrator.integrate_scatter(EvaluationFlags::gradients, dst);
+  }
+
+  void
+  do_cell_integral_range(
+    const MatrixFree<dim, number> &              matrix_free,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &range) const
+  {
+    FECellIntegrator integrator(matrix_free, range);
+
+    for (unsigned cell = range.first; cell < range.second; ++cell)
+      {
+        integrator.reinit(cell);
+
+        do_cell_integral_global(integrator, dst, src);
+      }
+  }
+
+  MatrixFree<dim, number> matrix_free;
+
+  AffineConstraints<number> constraints;
+
+  mutable TrilinosWrappers::SparseMatrix system_matrix;
+};
+
+struct GMGParameters
+{
+  struct CoarseSolverParameters
+  {
+    std::string  type            = "cg_with_amg"; // "cg";
+    unsigned int maxiter         = 10000;
+    double       abstol          = 1e-20;
+    double       reltol          = 1e-4;
+    unsigned int smoother_sweeps = 1;
+    unsigned int n_cycles        = 1;
+    std::string  smoother_type   = "ILU";
+  };
+
+  struct SmootherParameters
+  {
+    std::string  type                = "chebyshev";
+    double       smoothing_range     = 20;
+    unsigned int degree              = 5;
+    unsigned int eig_cg_n_iterations = 20;
+  };
+
+  SmootherParameters     smoother;
+  CoarseSolverParameters coarse_solver;
+
+  unsigned int maxiter = 10000;
+  double       abstol  = 1e-20;
+  double       reltol  = 1e-4;
+};
+
+template <typename VectorType,
+          int dim,
+          typename SystemMatrixType,
+          typename LevelMatrixType,
+          typename MGTransferType>
+static void
+mg_solve(SolverControl &                       solver_control,
+         VectorType &                          dst,
+         const VectorType &                    src,
+         const GMGParameters &                 mg_data,
+         const DoFHandler<dim> &               dof,
+         const SystemMatrixType &              fine_matrix,
+         const MGLevelObject<LevelMatrixType> &mg_matrices,
+         const MGTransferType &                mg_transfer)
+{
+  AssertThrow(mg_data.smoother.type == "chebyshev", ExcNotImplemented());
+
+  const unsigned int min_level = mg_matrices.min_level();
+  const unsigned int max_level = mg_matrices.max_level();
+
+  using Number                     = typename VectorType::value_type;
+  using SmootherPreconditionerType = DiagonalMatrix<VectorType>;
+  using SmootherType               = PreconditionChebyshev<LevelMatrixType,
+                                             VectorType,
+                                             SmootherPreconditionerType>;
+  using PreconditionerType = PreconditionMG<dim, VectorType, MGTransferType>;
+
+  // Initialize level operators.
+  mg::Matrix<VectorType> mg_matrix(mg_matrices);
+
+  // Initialize smoothers.
+  MGLevelObject<typename SmootherType::AdditionalData> smoother_data(min_level,
+                                                                     max_level);
+
+  for (unsigned int level = min_level; level <= max_level; level++)
+    {
+      smoother_data[level].preconditioner =
+        std::make_shared<SmootherPreconditionerType>();
+      mg_matrices[level].compute_inverse_diagonal(
+        smoother_data[level].preconditioner->get_vector());
+      smoother_data[level].smoothing_range = mg_data.smoother.smoothing_range;
+      smoother_data[level].degree          = mg_data.smoother.degree;
+      smoother_data[level].eig_cg_n_iterations =
+        mg_data.smoother.eig_cg_n_iterations;
+    }
+
+  MGSmootherPrecondition<LevelMatrixType, SmootherType, VectorType> mg_smoother;
+  mg_smoother.initialize(mg_matrices, smoother_data);
+
+  // Initialize coarse-grid solver.
+  ReductionControl     coarse_grid_solver_control(mg_data.coarse_solver.maxiter,
+                                              mg_data.coarse_solver.abstol,
+                                              mg_data.coarse_solver.reltol,
+                                              false,
+                                              false);
+  SolverCG<VectorType> coarse_grid_solver(coarse_grid_solver_control);
+
+  PreconditionIdentity precondition_identity;
+  PreconditionChebyshev<LevelMatrixType, VectorType, DiagonalMatrix<VectorType>>
+    precondition_chebyshev;
+
+#ifdef DEAL_II_WITH_TRILINOS
+  TrilinosWrappers::PreconditionAMG precondition_amg;
+#endif
+
+  std::unique_ptr<MGCoarseGridBase<VectorType>> mg_coarse;
+
+  if (mg_data.coarse_solver.type == "cg")
+    {
+      // CG with identity matrix as preconditioner
+
+      mg_coarse =
+        std::make_unique<MGCoarseGridIterativeSolver<VectorType,
+                                                     SolverCG<VectorType>,
+                                                     LevelMatrixType,
+                                                     PreconditionIdentity>>(
+          coarse_grid_solver, mg_matrices[min_level], precondition_identity);
+    }
+  else if (mg_data.coarse_solver.type == "cg_with_chebyshev")
+    {
+      // CG with Chebyshev as preconditioner
+
+      typename SmootherType::AdditionalData smoother_data;
+
+      smoother_data.preconditioner =
+        std::make_shared<DiagonalMatrix<VectorType>>();
+      mg_matrices[min_level].compute_inverse_diagonal(
+        smoother_data.preconditioner->get_vector());
+      smoother_data.smoothing_range     = mg_data.smoother.smoothing_range;
+      smoother_data.degree              = mg_data.smoother.degree;
+      smoother_data.eig_cg_n_iterations = mg_data.smoother.eig_cg_n_iterations;
+
+      precondition_chebyshev.initialize(mg_matrices[min_level], smoother_data);
+
+      mg_coarse = std::make_unique<
+        MGCoarseGridIterativeSolver<VectorType,
+                                    SolverCG<VectorType>,
+                                    LevelMatrixType,
+                                    decltype(precondition_chebyshev)>>(
+        coarse_grid_solver, mg_matrices[min_level], precondition_chebyshev);
+    }
+  else if (mg_data.coarse_solver.type == "cg_with_amg")
+    {
+      // CG with AMG as preconditioner
+
+#ifdef DEAL_II_WITH_TRILINOS
+      TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
+      amg_data.smoother_sweeps = mg_data.coarse_solver.smoother_sweeps;
+      amg_data.n_cycles        = mg_data.coarse_solver.n_cycles;
+      amg_data.smoother_type   = mg_data.coarse_solver.smoother_type.c_str();
+
+      // CG with AMG as preconditioner
+      precondition_amg.initialize(mg_matrices[min_level].get_system_matrix(),
+                                  amg_data);
+
+      mg_coarse = std::make_unique<
+        MGCoarseGridIterativeSolver<VectorType,
+                                    SolverCG<VectorType>,
+                                    LevelMatrixType,
+                                    decltype(precondition_amg)>>(
+        coarse_grid_solver, mg_matrices[min_level], precondition_amg);
+#else
+      AssertThrow(false, ExcNotImplemented());
+#endif
+    }
+  else
+    {
+      AssertThrow(false, ExcNotImplemented());
+    }
+
+  // Create multigrid object.
+  Multigrid<VectorType> mg(
+    mg_matrix, *mg_coarse, mg_transfer, mg_smoother, mg_smoother);
+
+  // Convert it to a preconditioner.
+  PreconditionerType preconditioner(dof, mg, mg_transfer);
+
+  // Finally, solve.
+  SolverCG<VectorType>(solver_control)
+    .solve(fine_matrix, dst, src, preconditioner);
+}
+
+#endif

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.