mg_transfer.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 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.
//
// ---------------------------------------------------------------------

#ifndef dealii_mg_transfer_h
#define dealii_mg_transfer_h

#include <deal.II/base/config.h>

#include <deal.II/base/mg_level_object.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/vector_memory.h>

#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>

#include <memory>


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  template <typename VectorType>
  struct MatrixSelector
  {
    using Sparsity = ::SparsityPattern;
    using Matrix   = ::SparseMatrix<typename VectorType::value_type>;

    static const bool requires_distributed_sparsity_pattern = false;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &                   matrix,
           Sparsity &                 sparsity,
           int                        level,
           const SparsityPatternType &sp,
           const DoFHandlerType &)
    {
      sparsity.copy_from(sp);
      (void)level;
      matrix.reinit(sparsity);
    }
  };

#ifdef DEAL_II_WITH_TRILINOS
  template <typename Number>
  struct MatrixSelector<LinearAlgebra::distributed::Vector<Number>>
  {
    using Sparsity = ::TrilinosWrappers::SparsityPattern;
    using Matrix   = ::TrilinosWrappers::SparseMatrix;

    static const bool requires_distributed_sparsity_pattern = false;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &matrix,
           Sparsity &,
           int                        level,
           const SparsityPatternType &sp,
           DoFHandlerType &           dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,
                                    DoFHandlerType::space_dimension>
        *dist_tria = dynamic_cast<
          const parallel::Triangulation<DoFHandlerType::dimension,
                                        DoFHandlerType::space_dimension> *>(
          &(dh.get_triangulation()));
      MPI_Comm communicator =
        dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF;

      matrix.reinit(dh.locally_owned_mg_dofs(level + 1),
                    dh.locally_owned_mg_dofs(level),
                    sp,
                    communicator,
                    true);
    }
  };

  template <>
  struct MatrixSelector<::TrilinosWrappers::MPI::Vector>
  {
    using Sparsity = ::TrilinosWrappers::SparsityPattern;
    using Matrix   = ::TrilinosWrappers::SparseMatrix;

    static const bool requires_distributed_sparsity_pattern = false;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &matrix,
           Sparsity &,
           int                        level,
           const SparsityPatternType &sp,
           DoFHandlerType &           dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,
                                    DoFHandlerType::space_dimension>
        *dist_tria = dynamic_cast<
          const parallel::Triangulation<DoFHandlerType::dimension,
                                        DoFHandlerType::space_dimension> *>(
          &(dh.get_triangulation()));
      MPI_Comm communicator =
        dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF;
      matrix.reinit(dh.locally_owned_mg_dofs(level + 1),
                    dh.locally_owned_mg_dofs(level),
                    sp,
                    communicator,
                    true);
    }
  };

#  ifdef DEAL_II_WITH_MPI
#    ifdef DEAL_II_TRILINOS_WITH_TPETRA
  template <typename Number>
  struct MatrixSelector<::LinearAlgebra::TpetraWrappers::Vector<Number>>
  {
    using Sparsity = ::TrilinosWrappers::SparsityPattern;
    using Matrix   = ::TrilinosWrappers::SparseMatrix;

    static const bool requires_distributed_sparsity_pattern = false;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &matrix,
           Sparsity &,
           int                        level,
           const SparsityPatternType &sp,
           DoFHandlerType &           dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,
                                    DoFHandlerType::space_dimension>
        *dist_tria = dynamic_cast<
          const parallel::Triangulation<DoFHandlerType::dimension,
                                        DoFHandlerType::space_dimension> *>(
          &(dh.get_triangulation()));
      MPI_Comm communicator =
        dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF;
      matrix.reinit(dh.locally_owned_mg_dofs(level + 1),
                    dh.locally_owned_mg_dofs(level),
                    sp,
                    communicator,
                    true);
    }
  };
#    endif

  template <>
  struct MatrixSelector<::LinearAlgebra::EpetraWrappers::Vector>
  {
    using Sparsity = ::TrilinosWrappers::SparsityPattern;
    using Matrix   = ::TrilinosWrappers::SparseMatrix;

    static const bool requires_distributed_sparsity_pattern = false;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &matrix,
           Sparsity &,
           int                        level,
           const SparsityPatternType &sp,
           DoFHandlerType &           dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,
                                    DoFHandlerType::space_dimension>
        *dist_tria = dynamic_cast<
          const parallel::Triangulation<DoFHandlerType::dimension,
                                        DoFHandlerType::space_dimension> *>(
          &(dh.get_triangulation()));
      MPI_Comm communicator =
        dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF;
      matrix.reinit(dh.locally_owned_mg_dofs(level + 1),
                    dh.locally_owned_mg_dofs(level),
                    sp,
                    communicator,
                    true);
    }
  };
#  endif

#else
  // ! DEAL_II_WITH_TRILINOS
  template <typename Number>
  struct MatrixSelector<LinearAlgebra::distributed::Vector<Number>>
  {
    using Sparsity = ::SparsityPattern;
    using Matrix   = ::SparseMatrix<Number>;

    static const bool requires_distributed_sparsity_pattern = false;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &,
           Sparsity &,
           int,
           const SparsityPatternType &,
           const DoFHandlerType &)
    {
      AssertThrow(
        false,
        ExcNotImplemented(
          "ERROR: MGTransferPrebuilt with LinearAlgebra::distributed::Vector currently "
          "needs deal.II to be configured with Trilinos."));
    }
  };

#endif

#ifdef DEAL_II_WITH_PETSC
  template <>
  struct MatrixSelector<::PETScWrappers::MPI::Vector>
  {
    using Sparsity = ::DynamicSparsityPattern;
    using Matrix   = ::PETScWrappers::MPI::SparseMatrix;

    static const bool requires_distributed_sparsity_pattern = true;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void
    reinit(Matrix &matrix,
           Sparsity &,
           int                        level,
           const SparsityPatternType &sp,
           const DoFHandlerType &     dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,
                                    DoFHandlerType::space_dimension>
        *dist_tria = dynamic_cast<
          const parallel::Triangulation<DoFHandlerType::dimension,
                                        DoFHandlerType::space_dimension> *>(
          &(dh.get_triangulation()));
      MPI_Comm communicator =
        dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF;
      // Reinit PETSc matrix
      matrix.reinit(dh.locally_owned_mg_dofs(level + 1),
                    dh.locally_owned_mg_dofs(level),
                    sp,
                    communicator);
    }
  };
#endif
} // namespace internal

/*
 * MGTransferBase is defined in mg_base.h
 */




template <typename VectorType>
class MGLevelGlobalTransfer : public MGTransferBase<VectorType>
{
public:
  void
  clear();

  template <int dim, class InVector, int spacedim>
  void
  copy_to_mg(const DoFHandler<dim, spacedim> &mg_dof,
             MGLevelObject<VectorType> &      dst,
             const InVector &                 src) const;

  template <int dim, class OutVector, int spacedim>
  void
  copy_from_mg(const DoFHandler<dim, spacedim> &mg_dof,
               OutVector &                      dst,
               const MGLevelObject<VectorType> &src) const;

  template <int dim, class OutVector, int spacedim>
  void
  copy_from_mg_add(const DoFHandler<dim, spacedim> &mg_dof,
                   OutVector &                      dst,
                   const MGLevelObject<VectorType> &src) const;

  void
  set_component_to_block_map(const std::vector<unsigned int> &map);

  std::size_t
  memory_consumption() const;

  void
  print_indices(std::ostream &os) const;

protected:
  template <int dim, int spacedim>
  void
  fill_and_communicate_copy_indices(const DoFHandler<dim, spacedim> &mg_dof);

  std::vector<types::global_dof_index> sizes;

  std::vector<
    std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
    copy_indices;

  std::vector<
    std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
    copy_indices_global_mine;

  std::vector<
    std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
    copy_indices_level_mine;

  bool perform_plain_copy;

  std::vector<unsigned int> component_to_block_map;

  SmartPointer<const MGConstrainedDoFs, MGLevelGlobalTransfer<VectorType>>
    mg_constrained_dofs;
};



template <typename Number>
class MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>
  : public MGTransferBase<LinearAlgebra::distributed::Vector<Number>>
{
public:
  void
  clear();

  template <int dim, typename Number2, int spacedim>
  void
  copy_to_mg(const DoFHandler<dim, spacedim> &                          mg_dof,
             MGLevelObject<LinearAlgebra::distributed::Vector<Number>> &dst,
             const LinearAlgebra::distributed::Vector<Number2> &src) const;

  template <int dim, typename Number2, int spacedim>
  void
  copy_from_mg(
    const DoFHandler<dim, spacedim> &                                mg_dof,
    LinearAlgebra::distributed::Vector<Number2> &                    dst,
    const MGLevelObject<LinearAlgebra::distributed::Vector<Number>> &src) const;

  template <int dim, typename Number2, int spacedim>
  void
  copy_from_mg_add(
    const DoFHandler<dim, spacedim> &                                mg_dof,
    LinearAlgebra::distributed::Vector<Number2> &                    dst,
    const MGLevelObject<LinearAlgebra::distributed::Vector<Number>> &src) const;

  void
  set_component_to_block_map(const std::vector<unsigned int> &map);

  std::size_t
  memory_consumption() const;

  void
  print_indices(std::ostream &os) const;

protected:
  template <int dim, typename Number2, int spacedim>
  void
  copy_to_mg(const DoFHandler<dim, spacedim> &                          mg_dof,
             MGLevelObject<LinearAlgebra::distributed::Vector<Number>> &dst,
             const LinearAlgebra::distributed::Vector<Number2> &        src,
             const bool solution_transfer) const;

  template <int dim, int spacedim>
  void
  fill_and_communicate_copy_indices(const DoFHandler<dim, spacedim> &mg_dof);

  std::vector<types::global_dof_index> sizes;

  std::vector<std::vector<std::pair<unsigned int, unsigned int>>> copy_indices;


  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
    solution_copy_indices;

  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
    copy_indices_global_mine;

  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
    solution_copy_indices_global_mine;

  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
    copy_indices_level_mine;

  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
    solution_copy_indices_level_mine;

  bool perform_plain_copy;

  bool perform_renumbered_plain_copy;

  std::vector<unsigned int> component_to_block_map;

  SmartPointer<
    const MGConstrainedDoFs,
    MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>>
    mg_constrained_dofs;

  mutable LinearAlgebra::distributed::Vector<Number> ghosted_global_vector;

  mutable LinearAlgebra::distributed::Vector<Number>
    solution_ghosted_global_vector;

  mutable MGLevelObject<LinearAlgebra::distributed::Vector<Number>>
    ghosted_level_vector;

  mutable MGLevelObject<LinearAlgebra::distributed::Vector<Number>>
    solution_ghosted_level_vector;
};



template <typename VectorType>
class MGTransferPrebuilt : public MGLevelGlobalTransfer<VectorType>
{
public:
  MGTransferPrebuilt() = default;

  MGTransferPrebuilt(const MGConstrainedDoFs &mg_constrained_dofs);

  DEAL_II_DEPRECATED
  MGTransferPrebuilt(const AffineConstraints<double> &constraints,
                     const MGConstrainedDoFs &        mg_constrained_dofs);

  virtual ~MGTransferPrebuilt() override = default;

  void
  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);

  DEAL_II_DEPRECATED
  void
  initialize_constraints(const AffineConstraints<double> &constraints,
                         const MGConstrainedDoFs &        mg_constrained_dofs);

  void
  clear();

  template <int dim, int spacedim>
  void
  build_matrices(const DoFHandler<dim, spacedim> &mg_dof);

  virtual void
  prolongate(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;

  DeclException0(ExcNoProlongation);

  DeclException0(ExcMatricesNotBuilt);

  std::size_t
  memory_consumption() const;

  void
  print_matrices(std::ostream &os) const;

private:
  std::vector<
    std::shared_ptr<typename internal::MatrixSelector<VectorType>::Sparsity>>
    prolongation_sparsities;

  std::vector<
    std::shared_ptr<typename internal::MatrixSelector<VectorType>::Matrix>>
    prolongation_matrices;

  std::vector<std::vector<bool>> interface_dofs;
};


DEAL_II_NAMESPACE_CLOSE

#endif