]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MGTransferGlobalCoarsening: initialize with MGConstrainedDoFs 15752/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 15 Jul 2023 12:34:58 +0000 (14:34 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 29 Jul 2023 16:04:47 +0000 (18:04 +0200)
include/deal.II/multigrid/mg_constrained_dofs.h
include/deal.II/multigrid/mg_transfer_global_coarsening.h
tests/multigrid-global-coarsening/multigrid_a_03.cc [new file with mode: 0644]
tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output [new file with mode: 0644]

index 9d0f2511d73fc877ef765f348778af7d7af0acd4..82f2d2711ba835862fe480b71585e52891131b1c 100644 (file)
@@ -225,6 +225,27 @@ public:
   const AffineConstraints<double> &
   get_user_constraint_matrix(const unsigned int level) const;
 
+  /**
+   * Merge selected constraints of a specifiedlevel into a given single
+   * AffineConstraints object.
+   *
+   * @param constraints AffineConstraints object to be filled.
+   * @param level Refinement to be considered.
+   * @param add_boundary_indices Add boundary indices.
+   * @param add_refinement_edge_indices Add refinement-edge indices.
+   * @param add_level_constraints Add level constraints including the one passed
+   *   during initialize() and periodicy constraints.
+   * @param add_user_constraints Add user constraints.
+   */
+  template <typename Number>
+  void
+  merge_constraints(AffineConstraints<Number> &constraints,
+                    const unsigned int         level,
+                    const bool                 add_boundary_indices,
+                    const bool                 add_refinement_edge_indices,
+                    const bool                 add_level_constraints,
+                    const bool                 add_user_constraints) const;
+
 private:
   /**
    * The indices of boundary dofs for each level.
@@ -565,6 +586,57 @@ MGConstrainedDoFs::get_user_constraint_matrix(const unsigned int level) const
 }
 
 
+template <typename Number>
+inline void
+MGConstrainedDoFs::merge_constraints(AffineConstraints<Number> &constraints,
+                                     const unsigned int         level,
+                                     const bool add_boundary_indices,
+                                     const bool add_refinement_edge_indices,
+                                     const bool add_level_constraints,
+                                     const bool add_user_constraints) const
+{
+  constraints.clear();
+
+  // determine local lines
+  IndexSet index_set(this->get_refinement_edge_indices(level).size());
+
+  if (add_boundary_indices && this->have_boundary_indices())
+    index_set.add_indices(this->get_boundary_indices(level));
+
+  if (add_refinement_edge_indices)
+    index_set.add_indices(this->get_refinement_edge_indices(level));
+
+  if (add_level_constraints)
+    index_set.add_indices(this->get_level_constraints(level).get_local_lines());
+
+  if (add_user_constraints)
+    index_set.add_indices(
+      this->get_user_constraint_matrix(level).get_local_lines());
+
+  constraints.reinit(index_set);
+
+  // merge constraints
+  if (add_boundary_indices && this->have_boundary_indices())
+    constraints.add_lines(this->get_boundary_indices(level));
+
+  if (add_refinement_edge_indices)
+    constraints.add_lines(this->get_refinement_edge_indices(level));
+
+  if (add_level_constraints)
+    constraints.merge(this->get_level_constraints(level),
+                      AffineConstraints<Number>::left_object_wins,
+                      true);
+
+  if (add_user_constraints)
+    constraints.merge(this->get_user_constraint_matrix(level),
+                      AffineConstraints<Number>::left_object_wins,
+                      true);
+
+  // finalize setup
+  constraints.close();
+}
+
+
 
 DEAL_II_NAMESPACE_CLOSE
 
index 0ef32204b449059acbd6a7863b6ad91e0bc535b8..3a72e69c9c73d1cab58a576a14751de9dec5e57a 100644 (file)
@@ -806,7 +806,23 @@ private:
 /**
  * Implementation of the MGTransferBase. In contrast to
  * other multigrid transfer operators, the user can provide separate
- * transfer operators of type MGTwoLevelTransfer between each level.
+ * transfer operators of type MGTwoLevelTransfer between each level. The
+ * sequence of functions calls for setup is:
+ * @code
+ * MGTransferGlobalCoarsening mg_transfer;
+ * mg_transfer.intitialize_two_level_transfers(two_level_transfers);
+ * mg_transfer.build(partitioners);
+ * @endcode
+ *
+ * Alternatively, this class can also be set up as in the case of
+ * MGTransferMatrixFree:
+ * @code
+ * MGTransferGlobalCoarsening mg_transfer;
+ * mg_transfer.initialize_constraints(mg_constrained_dofs);
+ * mg_transfer.build(dof_handler, partitioners);
+ * @endcode
+ * However, this is way to set up is currently only working for globally
+ * refined meshes.
  *
  * This class currently only works for the tensor-product finite elements
  * FE_Q and FE_DGQ and simplex elements FE_SimplexP and FE_SimplexDGP as well as
@@ -824,6 +840,16 @@ public:
    */
   using Number = typename VectorType::value_type;
 
+  /**
+   * Default constructor.
+   */
+  MGTransferGlobalCoarsening() = default;
+
+  /**
+   * @name Global coarsening.
+   */
+  /** @{ */
+
   /**
    * Constructor taking a collection of transfer operators (with the coarsest
    * level kept empty in @p transfer) and an optional function that initializes the
@@ -839,6 +865,14 @@ public:
     const std::function<void(const unsigned int, VectorType &)>
       &initialize_dof_vector = {});
 
+  /**
+   * Set two-level transfers.
+   */
+  template <typename MGTwoLevelTransferObject>
+  void
+  intitialize_two_level_transfers(
+    const MGLevelObject<MGTwoLevelTransferObject> &transfer);
+
   /**
    * Similar function to MGTransferMatrixFree::build() with the difference that
    * the information for the prolongation for each level has already been built.
@@ -859,6 +893,57 @@ public:
   build(const std::function<void(const unsigned int, VectorType &)>
           &initialize_dof_vector);
 
+  /** @} */
+
+  /**
+   * @name Local smoothing.
+   */
+  /** @{ */
+
+  /**
+   * Constructor with constraints. Equivalent to the default constructor
+   * followed by initialize_constraints().
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  MGTransferGlobalCoarsening(const MGConstrainedDoFs &mg_constrained_dofs);
+
+  /**
+   * Initialize the constraints to be used in build().
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  void
+  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
+
+  /**
+   * Actually build the information for the prolongation for each level.
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  void
+  build(const DoFHandler<dim> &dof_handler,
+        const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+          &external_partitioners = {});
+
+  /**
+   * Same as above but taking a lambda for initializing vector instead of
+   * partitioners.
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  void
+  build(const DoFHandler<dim> &dof_handler,
+        const std::function<void(const unsigned int, VectorType &)>
+          &initialize_dof_vector);
+
+  /** @} */
+
+  /**
+   * @name Tranfer functions.
+   */
+  /** @{ */
+
   /**
    * Perform prolongation.
    */
@@ -935,6 +1020,13 @@ public:
                     MGLevelObject<VectorType> &dst,
                     const InVector &           src) const;
 
+  /** @} */
+
+  /**
+   * @name Utility functions.
+   */
+  /** @{ */
+
   /**
    * Return the memory consumption of the allocated memory in this class.
    *
@@ -956,13 +1048,47 @@ public:
   unsigned int
   max_level() const;
 
+  /** @} */
+
 private:
+  /**
+   * Initial internal transfer operator.
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  void
+  intitialize_internal_transfer(
+    const DoFHandler<dim> &                      dof_handler,
+    const SmartPointer<const MGConstrainedDoFs> &mg_constrained_dofs);
+
+  /**
+   * Set references to two-level transfer operators to be used.
+   */
+  template <typename MGTwoLevelTransferObject>
+  void
+  intitialize_transfer_references(
+    const MGLevelObject<MGTwoLevelTransferObject> &transfer);
+
   /**
    * Function to initialize internal level vectors.
    */
   void
   initialize_dof_vector(const unsigned int level, VectorType &vector) const;
 
+  /**
+   * MGConstrainedDoFs passed during build().
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  SmartPointer<const MGConstrainedDoFs> mg_constrained_dofs;
+
+  /**
+   * Internal transfer operator.
+   *
+   * @note See also MGTransferMatrixFree.
+   */
+  MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> internal_transfer;
+
   /**
    * Collection of the two-level transfer operators.
    */
@@ -1021,6 +1147,80 @@ MGTransferGlobalCoarsening<dim, VectorType>::MGTransferGlobalCoarsening(
   const MGLevelObject<MGTwoLevelTransferObject> &transfer,
   const std::function<void(const unsigned int, VectorType &)>
     &initialize_dof_vector)
+{
+  this->intitialize_transfer_references(transfer);
+  this->build(initialize_dof_vector);
+}
+
+
+
+template <int dim, typename VectorType>
+MGTransferGlobalCoarsening<dim, VectorType>::MGTransferGlobalCoarsening(
+  const MGConstrainedDoFs &mg_constrained_dofs)
+{
+  this->initialize_constraints(mg_constrained_dofs);
+}
+
+
+
+template <int dim, typename VectorType>
+template <typename MGTwoLevelTransferObject>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::intitialize_two_level_transfers(
+  const MGLevelObject<MGTwoLevelTransferObject> &transfer)
+{
+  this->intitialize_transfer_references(transfer);
+}
+
+
+
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::initialize_constraints(
+  const MGConstrainedDoFs &mg_constrained_dofs)
+{
+  this->mg_constrained_dofs = &mg_constrained_dofs;
+}
+
+
+
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::intitialize_internal_transfer(
+  const DoFHandler<dim> &                      dof_handler,
+  const SmartPointer<const MGConstrainedDoFs> &mg_constrained_dofs)
+{
+  const unsigned int min_level = 0;
+  const unsigned int max_level =
+    dof_handler.get_triangulation().n_global_levels() - 1;
+
+  MGLevelObject<AffineConstraints<typename VectorType::value_type>> constraints(
+    min_level, max_level);
+
+  if (mg_constrained_dofs)
+    for (unsigned int l = min_level; l <= max_level; ++l)
+      mg_constrained_dofs->merge_constraints(
+        constraints[l],
+        l,
+        /*add_boundary_indices*/ true,
+        /*add_refinement_edge_indices*/ false,
+        /*add_level_constraints*/ true,
+        /*add_user_constraints*/ true);
+
+  this->internal_transfer.resize(min_level, max_level);
+
+  for (unsigned int l = min_level; l < max_level; ++l)
+    internal_transfer[l + 1].reinit_geometric_transfer(
+      dof_handler, dof_handler, constraints[l + 1], constraints[l], l + 1, l);
+}
+
+
+
+template <int dim, typename VectorType>
+template <typename MGTwoLevelTransferObject>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::intitialize_transfer_references(
+  const MGLevelObject<MGTwoLevelTransferObject> &transfer)
 {
   const unsigned int min_level = transfer.min_level();
   const unsigned int max_level = transfer.max_level();
@@ -1031,8 +1231,6 @@ MGTransferGlobalCoarsening<dim, VectorType>::MGTransferGlobalCoarsening(
     this->transfer[l] = &const_cast<MGTwoLevelTransferBase<VectorType> &>(
       static_cast<const MGTwoLevelTransferBase<VectorType> &>(
         Utilities::get_underlying_value(transfer[l])));
-
-  this->build(initialize_dof_vector);
 }
 
 
@@ -1094,6 +1292,34 @@ MGTransferGlobalCoarsening<dim, VectorType>::build(
 
 
 
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::build(
+  const DoFHandler<dim> &dof_handler,
+  const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+    &external_partitioners)
+{
+  this->intitialize_internal_transfer(dof_handler, mg_constrained_dofs);
+  this->intitialize_transfer_references(internal_transfer);
+  this->build(external_partitioners);
+}
+
+
+
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::build(
+  const DoFHandler<dim> &dof_handler,
+  const std::function<void(const unsigned int, VectorType &)>
+    &initialize_dof_vector)
+{
+  this->intitialize_internal_transfer(dof_handler, mg_constrained_dofs);
+  this->intitialize_transfer_references(internal_transfer);
+  this->build(initialize_dof_vector);
+}
+
+
+
 template <int dim, typename VectorType>
 void
 MGTransferGlobalCoarsening<dim, VectorType>::initialize_dof_vector(
diff --git a/tests/multigrid-global-coarsening/multigrid_a_03.cc b/tests/multigrid-global-coarsening/multigrid_a_03.cc
new file mode 100644 (file)
index 0000000..0f4be4b
--- /dev/null
@@ -0,0 +1,153 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 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 global-coarsening multigrid so that it works also on refinement levels.
+ */
+
+#include "multigrid_util.h"
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int n_refinements, const unsigned int fe_degree_fine)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  parallel::distributed::Triangulation<dim> tria(
+    MPI_COMM_WORLD,
+    Triangulation<dim>::MeshSmoothing::none,
+    parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::subdivided_hyper_cube(tria, 2);
+  tria.refine_global(n_refinements);
+
+  const unsigned int min_level = 0;
+  const unsigned int max_level = n_refinements;
+
+  MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
+  MGLevelObject<Operator<dim, Number>>     operators(min_level, max_level);
+
+  std::unique_ptr<FiniteElement<dim>> fe;
+  std::unique_ptr<Quadrature<dim>>    quad;
+  std::unique_ptr<Mapping<dim>>       mapping;
+
+  fe      = std::make_unique<FE_Q<dim>>(fe_degree_fine);
+  quad    = std::make_unique<QGauss<dim>>(fe_degree_fine + 1);
+  mapping = std::make_unique<MappingFE<dim>>(FE_Q<dim>(1));
+
+  // set up dofhandler
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(*fe);
+  dof_handler.distribute_mg_dofs();
+
+  const std::set<types::boundary_id> dirichlet_boundary = {0};
+  MGConstrainedDoFs                  mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof_handler);
+  mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
+                                                     dirichlet_boundary);
+
+  // set up levels
+  for (auto l = min_level; l <= max_level; ++l)
+    {
+      auto &constraint = constraints[l];
+      auto &op         = operators[l];
+
+      // set up constraints
+      IndexSet relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof_handler,
+                                                    l,
+                                                    relevant_dofs);
+      constraint.reinit(relevant_dofs);
+      constraint.add_lines(mg_constrained_dofs.get_boundary_indices(l));
+      constraint.close();
+
+      constraint.print(std::cout);
+
+      // set up operator
+      op.reinit(*mapping, dof_handler, *quad, constraint, l);
+    }
+
+  // set up transfer operator
+  MGTransferGlobalCoarsening<dim, VectorType> transfer(mg_constrained_dofs);
+
+  transfer.build(dof_handler, [&](const auto l, auto &vec) {
+    operators[l].initialize_dof_vector(vec);
+  });
+
+  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_handler,
+           operators[max_level],
+           operators,
+           transfer);
+
+  constraints[max_level].distribute(dst);
+
+  deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
+          << solver_control.last_step() << std::endl;
+
+  return;
+
+  static unsigned int counter = 0;
+
+  MGLevelObject<VectorType> results(min_level, max_level);
+
+  transfer.interpolate_to_mg(dof_handler, results, dst);
+
+  for (unsigned int l = min_level; l <= max_level; ++l)
+    {
+      DataOut<dim> data_out;
+
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(
+        results[l],
+        "solution",
+        DataOut_DoFData<dim, dim>::DataVectorType::type_dof_data);
+      data_out.build_patches(*mapping, 2);
+
+      std::ofstream output("test." + std::to_string(dim) + "." +
+                           std::to_string(counter) + "." + std::to_string(l) +
+                           ".vtk");
+      data_out.write_vtk(output);
+    }
+
+  counter++;
+}
+
+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);
+}
diff --git a/tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_a_03.mpirun=1.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..e460d32
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL:0::2 2 2 3
+DEAL:0::2 3 2 3
+DEAL:0::2 4 2 3
+DEAL:0::2 2 3 3
+DEAL:0::2 3 3 3
+DEAL:0::2 4 3 3
+DEAL:0::2 2 4 3
+DEAL:0::2 3 4 3
+DEAL:0::2 4 4 3

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.