mg_transfer_global_coarsening.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 2023 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_mg_transfer_global_coarsening_h
#define dealii_mg_transfer_global_coarsening_h
 
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/mpi_remote_point_evaluation.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include <deal.II/matrix_free/constraint_info.h>
#include <deal.II/matrix_free/shape_info.h>
 
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
 
#include <deal.II/non_matching/mapping_info.h>
 
 
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
namespace internal
{
  class MGTwoLevelTransferImplementation;
}
 
namespace RepartitioningPolicyTools
{
  template <int dim, int spacedim>
  class Base;
}
#endif
 
 
 
namespace MGTransferGlobalCoarseningTools
{
  enum class PolynomialCoarseningSequenceType
  {
    bisect,
    decrease_by_one,
    go_to_one
  };
 
  unsigned int
  create_next_polynomial_coarsening_degree(
    const unsigned int                      degree,
    const PolynomialCoarseningSequenceType &p_sequence);
 
  std::vector<unsigned int>
  create_polynomial_coarsening_sequence(
    const unsigned int                      max_degree,
    const PolynomialCoarseningSequenceType &p_sequence);
 
  template <int dim, int spacedim>
  std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
  create_geometric_coarsening_sequence(
    const Triangulation<dim, spacedim> &tria);
 
  template <int dim, int spacedim>
  std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
  create_geometric_coarsening_sequence(
    Triangulation<dim, spacedim> &                        tria,
    const RepartitioningPolicyTools::Base<dim, spacedim> &policy,
    const bool preserve_fine_triangulation,
    const bool repartition_fine_triangulation);
 
  template <int dim, int spacedim>
  std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
  create_geometric_coarsening_sequence(
    const Triangulation<dim, spacedim> &                  tria,
    const RepartitioningPolicyTools::Base<dim, spacedim> &policy,
    const bool repartition_fine_triangulation = false);
 
} // namespace MGTransferGlobalCoarseningTools
 
 
template <typename VectorType>
class MGTwoLevelTransferBase : public Subscriptor
{
public:
  virtual void
  prolongate_and_add(VectorType &dst, const VectorType &src) const = 0;
 
  virtual void
  restrict_and_add(VectorType &dst, const VectorType &src) const = 0;
 
  virtual void
  interpolate(VectorType &dst, const VectorType &src) const = 0;
 
  virtual void
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_fine) = 0;
 
  virtual std::size_t
  memory_consumption() const = 0;
};
 
 
template <typename Number>
class MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>
  : public Subscriptor
{
public:
  using VectorType = LinearAlgebra::distributed::Vector<Number>;
 
  virtual void
  prolongate_and_add(VectorType &dst, const VectorType &src) const;
 
  virtual void
  restrict_and_add(VectorType &dst, const VectorType &src) const;
 
  virtual void
  interpolate(VectorType &dst, const VectorType &src) const = 0;
 
  virtual void
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_fine) = 0;
 
  virtual std::size_t
  memory_consumption() const = 0;
 
protected:
  virtual void
  prolongate_and_add_internal(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const = 0;
 
  virtual void
  restrict_and_add_internal(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const = 0;
 
  void
  update_ghost_values(
    const LinearAlgebra::distributed::Vector<Number> &vec) const;
 
  void
  compress(LinearAlgebra::distributed::Vector<Number> &vec,
           const VectorOperation::values               op) const;
 
  void
  zero_out_ghost_values(
    const LinearAlgebra::distributed::Vector<Number> &vec) const;
 
  template <int dim, std::size_t width>
  void
  internal_enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine,
    internal::MatrixFreeFunctions::ConstraintInfo<
      dim,
      VectorizedArray<Number, width>> &constraint_info_coarse,
    std::vector<unsigned int> &        dof_indices_fine);
 
  bool fine_element_is_continuous;
 
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_coarse;
 
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_fine;
 
  mutable LinearAlgebra::distributed::Vector<Number> vec_fine;
 
  mutable LinearAlgebra::distributed::Vector<Number> vec_coarse;
 
  std::shared_ptr<const Utilities::MPI::Partitioner>
    partitioner_coarse_embedded;
 
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_fine_embedded;
 
  mutable AlignedVector<Number> buffer_coarse_embedded;
 
  mutable AlignedVector<Number> buffer_fine_embedded;
};
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransfer : public MGTwoLevelTransferBase<VectorType>
{
public:
  void
  prolongate_and_add(VectorType &dst, const VectorType &src) const override;
 
  void
  restrict_and_add(VectorType &dst, const VectorType &src) const override;
 
  void
  interpolate(VectorType &dst, const VectorType &src) const override;
 
  void
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
    override;
 
  std::size_t
  memory_consumption() const override;
};
 
 
 
template <int dim, typename Number>
class MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>
  : public MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>
{
  using VectorizedArrayType = VectorizedArray<Number>;
 
public:
  void
  reinit_geometric_transfer(
    const DoFHandler<dim> &          dof_handler_fine,
    const DoFHandler<dim> &          dof_handler_coarse,
    const AffineConstraints<Number> &constraint_fine =
      AffineConstraints<Number>(),
    const AffineConstraints<Number> &constraint_coarse =
      AffineConstraints<Number>(),
    const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
    const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  void
  reinit_polynomial_transfer(
    const DoFHandler<dim> &          dof_handler_fine,
    const DoFHandler<dim> &          dof_handler_coarse,
    const AffineConstraints<Number> &constraint_fine =
      AffineConstraints<Number>(),
    const AffineConstraints<Number> &constraint_coarse =
      AffineConstraints<Number>(),
    const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
    const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  void
  reinit(const DoFHandler<dim> &          dof_handler_fine,
         const DoFHandler<dim> &          dof_handler_coarse,
         const AffineConstraints<Number> &constraint_fine =
           AffineConstraints<Number>(),
         const AffineConstraints<Number> &constraint_coarse =
           AffineConstraints<Number>(),
         const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
         const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  static bool
  fast_polynomial_transfer_supported(const unsigned int fe_degree_fine,
                                     const unsigned int fe_degree_coarse);
 
  void
  interpolate(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const override;
 
  void
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
    override;
 
  std::size_t
  memory_consumption() const override;
 
protected:
  void
  prolongate_and_add_internal(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const override;
 
  void
  restrict_and_add_internal(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const override;
 
private:
  struct MGTransferScheme
  {
    unsigned int n_coarse_cells;
 
    unsigned int n_dofs_per_cell_coarse;
 
    unsigned int n_dofs_per_cell_fine;
 
    unsigned int degree_coarse;
 
    unsigned int degree_fine;
 
    AlignedVector<VectorizedArrayType> prolongation_matrix;
 
    AlignedVector<VectorizedArrayType> prolongation_matrix_1d;
 
    AlignedVector<VectorizedArrayType> restriction_matrix;
 
    AlignedVector<VectorizedArrayType> restriction_matrix_1d;
 
    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
      shape_info_coarse;
  };
 
  std::vector<MGTransferScheme> schemes;
 
  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArrayType>
    constraint_info_coarse;
 
  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArrayType>
    constraint_info_fine;
 
  std::vector<Number> weights; // TODO: vectorize
 
  AlignedVector<VectorizedArrayType> weights_compressed;
 
  unsigned int n_components;
 
  friend class internal::MGTwoLevelTransferImplementation;
};
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransferNonNested : public MGTwoLevelTransferBase<VectorType>
{
public:
  void
  prolongate_and_add(VectorType &dst, const VectorType &src) const override;
 
  void
  restrict_and_add(VectorType &dst, const VectorType &src) const override;
 
  void
  interpolate(VectorType &dst, const VectorType &src) const override;
 
  std::size_t
  memory_consumption() const override;
};
 
 
 
template <int dim, typename Number>
class MGTwoLevelTransferNonNested<dim,
                                  LinearAlgebra::distributed::Vector<Number>>
  : public MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>
{
private:
  using VectorizedArrayType = VectorizedArray<Number, 1>;
 
public:
  void
  reinit(const DoFHandler<dim> &          dof_handler_fine,
         const DoFHandler<dim> &          dof_handler_coarse,
         const Mapping<dim> &             mapping_fine,
         const Mapping<dim> &             mapping_coarse,
         const AffineConstraints<Number> &constraint_fine =
           AffineConstraints<Number>(),
         const AffineConstraints<Number> &constraint_coarse =
           AffineConstraints<Number>());
 
  void
  interpolate(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const override;
 
  void
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
    override;
 
  std::size_t
  memory_consumption() const override;
 
protected:
  void
  prolongate_and_add_internal(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const override;
 
  void
  restrict_and_add_internal(
    LinearAlgebra::distributed::Vector<Number> &      dst,
    const LinearAlgebra::distributed::Vector<Number> &src) const override;
 
private:
  Utilities::MPI::RemotePointEvaluation<dim> rpe;
 
  std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> mapping_info;
 
  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArrayType>
    constraint_info;
 
  std::unique_ptr<FiniteElement<dim>> fe_coarse;
 
  std::vector<unsigned int> level_dof_indices_fine;
 
  std::vector<unsigned int> level_dof_indices_fine_ptrs;
};
 
 
 
template <int dim, typename VectorType>
class MGTransferGlobalCoarsening : public ::MGTransferBase<VectorType>
{
public:
  using Number = typename VectorType::value_type;
 
  template <typename MGTwoLevelTransferObject>
  MGTransferGlobalCoarsening(
    const MGLevelObject<MGTwoLevelTransferObject> &transfer,
    const std::function<void(const unsigned int, VectorType &)>
      &initialize_dof_vector = {});
 
  void
  build(const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
          &external_partitioners = {});
 
  void
  build(const std::function<void(const unsigned int, VectorType &)>
          &initialize_dof_vector);
 
  void
  prolongate(const unsigned int to_level,
             VectorType &       dst,
             const VectorType & src) const override;
 
  void
  prolongate_and_add(const unsigned int to_level,
                     VectorType &       dst,
                     const VectorType & src) const override;
 
  virtual void
  restrict_and_add(const unsigned int from_level,
                   VectorType &       dst,
                   const VectorType & src) const override;
 
  template <class InVector, int spacedim>
  void
  copy_to_mg(const DoFHandler<dim, spacedim> &dof_handler,
             MGLevelObject<VectorType> &      dst,
             const InVector &                 src) const;
 
  template <class OutVector, int spacedim>
  void
  copy_from_mg(const DoFHandler<dim, spacedim> &dof_handler,
               OutVector &                      dst,
               const MGLevelObject<VectorType> &src) const;
 
  template <class InVector>
  void
  interpolate_to_mg(MGLevelObject<VectorType> &dst, const InVector &src) const;
 
  template <class InVector, int spacedim>
  void
  interpolate_to_mg(const DoFHandler<dim, spacedim> &dof_handler,
                    MGLevelObject<VectorType> &      dst,
                    const InVector &                 src) const;
 
  std::size_t
  memory_consumption() const;
 
  unsigned int
  min_level() const;
 
  unsigned int
  max_level() const;
 
private:
  void
  initialize_dof_vector(const unsigned int level, VectorType &vector) const;
 
  MGLevelObject<SmartPointer<MGTwoLevelTransferBase<VectorType>>> transfer;
 
  std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
    external_partitioners;
};
 
 
 
template <int dim, typename VectorType>
class MGTransferBlockGlobalCoarsening
  : public MGTransferBlockMatrixFreeBase<
      dim,
      typename VectorType::value_type,
      MGTransferGlobalCoarsening<dim, VectorType>>
{
public:
  MGTransferBlockGlobalCoarsening(
    const MGTransferGlobalCoarsening<dim, VectorType> &transfer_operator);
 
protected:
  const MGTransferGlobalCoarsening<dim, VectorType> &
  get_matrix_free_transfer(const unsigned int b) const override;
 
private:
  const MGTransferGlobalCoarsening<dim, VectorType> &transfer_operator;
};
 
 
 
#ifndef DOXYGEN
 
/* ----------------------- Inline functions --------------------------------- */
 
 
 
template <int dim, typename VectorType>
template <typename MGTwoLevelTransferObject>
MGTransferGlobalCoarsening<dim, VectorType>::MGTransferGlobalCoarsening(
  const MGLevelObject<MGTwoLevelTransferObject> &transfer,
  const std::function<void(const unsigned int, VectorType &)>
    &initialize_dof_vector)
{
  const unsigned int min_level = transfer.min_level();
  const unsigned int max_level = transfer.max_level();
 
  this->transfer.resize(min_level, max_level);
 
  for (unsigned int l = min_level; l <= max_level; ++l)
    this->transfer[l] = &const_cast<MGTwoLevelTransferBase<VectorType> &>(
      static_cast<const MGTwoLevelTransferBase<VectorType> &>(
        Utilities::get_underlying_value(transfer[l])));
 
  this->build(initialize_dof_vector);
}
 
 
 
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::build(
  const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
    &external_partitioners)
{
  this->external_partitioners = external_partitioners;
 
  if (this->external_partitioners.size() > 0)
    {
      const unsigned int min_level = transfer.min_level();
      const unsigned int max_level = transfer.max_level();
 
      AssertDimension(this->external_partitioners.size(), transfer.n_levels());
 
      for (unsigned int l = min_level + 1; l <= max_level; ++l)
        transfer[l]->enable_inplace_operations_if_possible(
          this->external_partitioners[l - 1 - min_level],
          this->external_partitioners[l - min_level]);
    }
}
 
 
 
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::build(
  const std::function<void(const unsigned int, VectorType &)>
    &initialize_dof_vector)
{
  if (initialize_dof_vector)
    {
      const unsigned int min_level = transfer.min_level();
      const unsigned int max_level = transfer.max_level();
      const unsigned int n_levels  = transfer.n_levels();
 
      std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
        external_partitioners(n_levels);
 
      for (unsigned int l = min_level; l <= max_level; ++l)
        {
          LinearAlgebra::distributed::Vector<typename VectorType::value_type>
            vector;
          initialize_dof_vector(l, vector);
          external_partitioners[l - min_level] = vector.get_partitioner();
        }
 
      this->build(external_partitioners);
    }
  else
    {
      this->build();
    }
}
 
 
 
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::initialize_dof_vector(
  const unsigned int level,
  VectorType &       vec) const
{
  AssertDimension(transfer.n_levels(), external_partitioners.size());
  AssertIndexRange(level, transfer.max_level() + 1);
 
  const auto &external_partitioner =
    external_partitioners[level - transfer.min_level()];
 
  if (vec.get_partitioner().get() != external_partitioner.get())
    vec.reinit(external_partitioner);
}
 
 
 
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::prolongate(
  const unsigned int to_level,
  VectorType &       dst,
  const VectorType & src) const
{
  dst = 0;
  prolongate_and_add(to_level, dst, src);
}
 
 
 
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::prolongate_and_add(
  const unsigned int to_level,
  VectorType &       dst,
  const VectorType & src) const
{
  this->transfer[to_level]->prolongate_and_add(dst, src);
}
 
 
 
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::restrict_and_add(
  const unsigned int from_level,
  VectorType &       dst,
  const VectorType & src) const
{
  this->transfer[from_level]->restrict_and_add(dst, src);
}
 
 
 
template <int dim, typename VectorType>
template <class InVector, int spacedim>
void
MGTransferGlobalCoarsening<dim, VectorType>::copy_to_mg(
  const DoFHandler<dim, spacedim> &dof_handler,
  MGLevelObject<VectorType> &      dst,
  const InVector &                 src) const
{
  (void)dof_handler;
 
  for (unsigned int level = dst.min_level(); level <= dst.max_level(); ++level)
    {
      initialize_dof_vector(level, dst[level]);
 
      if (level == dst.max_level())
        dst[level].copy_locally_owned_data_from(src);
      else
        dst[level] = 0.0;
    }
}
 
 
 
template <int dim, typename VectorType>
template <class OutVector, int spacedim>
void
MGTransferGlobalCoarsening<dim, VectorType>::copy_from_mg(
  const DoFHandler<dim, spacedim> &dof_handler,
  OutVector &                      dst,
  const MGLevelObject<VectorType> &src) const
{
  (void)dof_handler;
 
  dst.copy_locally_owned_data_from(src[src.max_level()]);
}
 
 
 
template <int dim, typename VectorType>
template <class InVector>
void
MGTransferGlobalCoarsening<dim, VectorType>::interpolate_to_mg(
  MGLevelObject<VectorType> &dst,
  const InVector &           src) const
{
  const unsigned int min_level = transfer.min_level();
  const unsigned int max_level = transfer.max_level();
 
  AssertDimension(min_level, dst.min_level());
  AssertDimension(max_level, dst.max_level());
 
  for (unsigned int level = min_level; level <= max_level; ++level)
    initialize_dof_vector(level, dst[level]);
 
  dst[transfer.max_level()].copy_locally_owned_data_from(src);
 
  for (unsigned int l = max_level; l > min_level; --l)
    this->transfer[l]->interpolate(dst[l - 1], dst[l]);
}
 
 
 
template <int dim, typename VectorType>
template <class InVector, int spacedim>
void
MGTransferGlobalCoarsening<dim, VectorType>::interpolate_to_mg(
  const DoFHandler<dim, spacedim> &dof_handler,
  MGLevelObject<VectorType> &      dst,
  const InVector &                 src) const
{
  (void)dof_handler;
 
  this->interpolate_to_mg(dst, src);
}
 
 
 
template <int dim, typename VectorType>
std::size_t
MGTransferGlobalCoarsening<dim, VectorType>::memory_consumption() const
{
  std::size_t size = 0;
 
  const unsigned int min_level = transfer.min_level();
  const unsigned int max_level = transfer.max_level();
 
  for (unsigned int l = min_level + 1; l <= max_level; ++l)
    size += this->transfer[l]->memory_consumption();
 
  return size;
}
 
 
 
template <int dim, typename VectorType>
inline unsigned int
MGTransferGlobalCoarsening<dim, VectorType>::min_level() const
{
  return transfer.min_level();
}
 
 
 
template <int dim, typename VectorType>
inline unsigned int
MGTransferGlobalCoarsening<dim, VectorType>::max_level() const
{
  return transfer.max_level();
}
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif