multigrid.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2016 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__multigrid_h
#define dealii__multigrid_h


#include <deal.II/base/config.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/distributed/tria.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN


template <typename VectorType>
class Multigrid : public Subscriptor
{
public:
  enum Cycle
  {
    v_cycle,
    w_cycle,
    f_cycle
  };

  typedef VectorType vector_type;
  typedef const VectorType const_vector_type;

  template <int dim>
  Multigrid(const DoFHandler<dim>              &mg_dof_handler,
            const MGMatrixBase<VectorType>     &matrix,
            const MGCoarseGridBase<VectorType> &coarse,
            const MGTransferBase<VectorType>   &transfer,
            const MGSmootherBase<VectorType>   &pre_smooth,
            const MGSmootherBase<VectorType>   &post_smooth,
            const unsigned int                 minlevel = 0,
            const unsigned int                 maxlevel = numbers::invalid_unsigned_int,
            Cycle                              cycle = v_cycle) DEAL_II_DEPRECATED;

  Multigrid(const MGMatrixBase<VectorType>     &matrix,
            const MGCoarseGridBase<VectorType> &coarse,
            const MGTransferBase<VectorType>   &transfer,
            const MGSmootherBase<VectorType>   &pre_smooth,
            const MGSmootherBase<VectorType>   &post_smooth,
            const unsigned int                 minlevel = 0,
            const unsigned int                 maxlevel = numbers::invalid_unsigned_int,
            Cycle                              cycle = v_cycle);

  void reinit (const unsigned int minlevel,
               const unsigned int maxlevel);

  void cycle ();

  void vcycle ();

  void set_edge_matrices (const MGMatrixBase<VectorType> &edge_out,
                          const MGMatrixBase<VectorType> &edge_in);

  void set_edge_flux_matrices (const MGMatrixBase<VectorType> &edge_down,
                               const MGMatrixBase<VectorType> &edge_up);

  unsigned int get_maxlevel() const;

  unsigned int get_minlevel() const;

  void set_maxlevel (const unsigned int);

  void set_minlevel (const unsigned int level,
                     bool relative = false);

  void set_cycle(Cycle);

  void set_debug (const unsigned int);

private:

  void level_v_step (const unsigned int level);

  void level_step (const unsigned int level, Cycle cycle);

  Cycle cycle_type;

  unsigned int minlevel;

  unsigned int maxlevel;

public:
  MGLevelObject<VectorType> defect;

  MGLevelObject<VectorType> solution;

private:
  MGLevelObject<VectorType> t;

  MGLevelObject<VectorType> defect2;


  SmartPointer<const MGMatrixBase<VectorType>,Multigrid<VectorType> > matrix;

  SmartPointer<const MGCoarseGridBase<VectorType>,Multigrid<VectorType> > coarse;

  SmartPointer<const MGTransferBase<VectorType>,Multigrid<VectorType> > transfer;

  SmartPointer<const MGSmootherBase<VectorType>,Multigrid<VectorType> > pre_smooth;

  SmartPointer<const MGSmootherBase<VectorType>,Multigrid<VectorType> > post_smooth;

  SmartPointer<const MGMatrixBase<VectorType> > edge_out;

  SmartPointer<const MGMatrixBase<VectorType> > edge_in;

  SmartPointer<const MGMatrixBase<VectorType>,Multigrid<VectorType> > edge_down;

  SmartPointer<const MGMatrixBase<VectorType>,Multigrid<VectorType> > edge_up;

  unsigned int debug;

  template<int dim, class OtherVectorType, class TRANSFER> friend class PreconditionMG;
};


template<int dim, typename VectorType, class TRANSFER>
class PreconditionMG : public Subscriptor
{
public:
  PreconditionMG(const DoFHandler<dim> &dof_handler,
                 Multigrid<VectorType> &mg,
                 const TRANSFER        &transfer);

  bool empty () const;

  template<class OtherVectorType>
  void vmult (OtherVectorType       &dst,
              const OtherVectorType &src) const;

  template<class OtherVectorType>
  void vmult_add (OtherVectorType       &dst,
                  const OtherVectorType &src) const;

  template<class OtherVectorType>
  void Tvmult (OtherVectorType       &dst,
               const OtherVectorType &src) const;

  template<class OtherVectorType>
  void Tvmult_add (OtherVectorType       &dst,
                   const OtherVectorType &src) const;

  IndexSet locally_owned_range_indices() const;

  IndexSet locally_owned_domain_indices() const;

  MPI_Comm get_mpi_communicator() const;

private:
  SmartPointer<const DoFHandler<dim>,PreconditionMG<dim,VectorType,TRANSFER> > dof_handler;

  SmartPointer<Multigrid<VectorType>,PreconditionMG<dim,VectorType,TRANSFER> > multigrid;

  SmartPointer<const TRANSFER,PreconditionMG<dim,VectorType,TRANSFER> > transfer;
};

#ifndef DOXYGEN
/* --------------------------- inline functions --------------------- */


template <typename VectorType>
template <int dim>
Multigrid<VectorType>::Multigrid (const DoFHandler<dim>              &mg_dof_handler,
                                  const MGMatrixBase<VectorType>     &matrix,
                                  const MGCoarseGridBase<VectorType> &coarse,
                                  const MGTransferBase<VectorType>   &transfer,
                                  const MGSmootherBase<VectorType>   &pre_smooth,
                                  const MGSmootherBase<VectorType>   &post_smooth,
                                  const unsigned int                 min_level,
                                  const unsigned int                 max_level,
                                  Cycle                              cycle)
  :
  cycle_type(cycle),
  minlevel(min_level),
  matrix(&matrix, typeid(*this).name()),
  coarse(&coarse, typeid(*this).name()),
  transfer(&transfer, typeid(*this).name()),
  pre_smooth(&pre_smooth, typeid(*this).name()),
  post_smooth(&post_smooth, typeid(*this).name()),
  edge_down(0, typeid(*this).name()),
  edge_up(0, typeid(*this).name()),
  debug(0)
{
  const unsigned int dof_handler_max_level
    = mg_dof_handler.get_triangulation().n_global_levels()-1;
  if (max_level == numbers::invalid_unsigned_int)
    maxlevel = dof_handler_max_level;
  else
    maxlevel = max_level;

  reinit(minlevel, maxlevel);
}



template <typename VectorType>
Multigrid<VectorType>::Multigrid (const MGMatrixBase<VectorType>        &matrix,
                                  const MGCoarseGridBase<VectorType>    &coarse,
                                  const MGTransferBase<VectorType>      &transfer,
                                  const MGSmootherBase<VectorType>      &pre_smooth,
                                  const MGSmootherBase<VectorType>      &post_smooth,
                                  const unsigned int                    min_level,
                                  const unsigned int                    max_level,
                                  Cycle                                 cycle)
  :
  cycle_type(cycle),
  matrix(&matrix, typeid(*this).name()),
  coarse(&coarse, typeid(*this).name()),
  transfer(&transfer, typeid(*this).name()),
  pre_smooth(&pre_smooth, typeid(*this).name()),
  post_smooth(&post_smooth, typeid(*this).name()),
  edge_out(0, typeid(*this).name()),
  edge_in(0, typeid(*this).name()),
  edge_down(0, typeid(*this).name()),
  edge_up(0, typeid(*this).name()),
  debug(0)
{
  if (max_level == numbers::invalid_unsigned_int)
    maxlevel = matrix.get_maxlevel();
  else
    maxlevel = max_level;
  reinit(min_level, maxlevel);
}



template <typename VectorType>
inline
unsigned int
Multigrid<VectorType>::get_maxlevel () const
{
  return maxlevel;
}



template <typename VectorType>
inline
unsigned int
Multigrid<VectorType>::get_minlevel () const
{
  return minlevel;
}


/* --------------------------- inline functions --------------------- */


template<int dim, typename VectorType, class TRANSFER>
PreconditionMG<dim, VectorType, TRANSFER>
::PreconditionMG(const DoFHandler<dim>  &dof_handler,
                 Multigrid<VectorType>  &mg,
                 const TRANSFER         &transfer)
  :
  dof_handler(&dof_handler),
  multigrid(&mg),
  transfer(&transfer)
{}

template<int dim, typename VectorType, class TRANSFER>
inline bool
PreconditionMG<dim, VectorType, TRANSFER>::empty () const
{
  return false;
}

template<int dim, typename VectorType, class TRANSFER>
template<class OtherVectorType>
void
PreconditionMG<dim, VectorType, TRANSFER>::vmult
(OtherVectorType       &dst,
 const OtherVectorType &src) const
{
  transfer->copy_to_mg(*dof_handler,
                       multigrid->defect,
                       src);
  multigrid->cycle();

  transfer->copy_from_mg(*dof_handler,
                         dst,
                         multigrid->solution);
}


template<int dim, typename VectorType, class TRANSFER>
IndexSet
PreconditionMG<dim, VectorType, TRANSFER>::locally_owned_range_indices() const
{
  return dof_handler->locally_owned_dofs();
}


template<int dim, typename VectorType, class TRANSFER>
IndexSet
PreconditionMG<dim, VectorType, TRANSFER>::locally_owned_domain_indices() const
{
  return dof_handler->locally_owned_dofs();
}


template<int dim, typename VectorType, class TRANSFER>
MPI_Comm
PreconditionMG<dim, VectorType, TRANSFER>::get_mpi_communicator() const
{
  // currently parallel GMG works with distributed Triangulation only,
  // so it should be a safe bet to use it to query MPI communicator:
  const Triangulation<dim> &tria = dof_handler->get_triangulation();
  const parallel::distributed::Triangulation<dim> *ptria = dynamic_cast<const parallel::distributed::Triangulation<dim> *>(&tria);
  Assert (ptria != NULL, ExcInternalError());
  return ptria->get_communicator ();
}


template<int dim, typename VectorType, class TRANSFER>
template<class OtherVectorType>
void
PreconditionMG<dim, VectorType, TRANSFER>::vmult_add
(OtherVectorType       &dst,
 const OtherVectorType &src) const
{
  transfer->copy_to_mg(*dof_handler,
                       multigrid->defect,
                       src);
  multigrid->cycle();
  transfer->copy_from_mg_add(*dof_handler,
                             dst,
                             multigrid->solution);
}


template<int dim, typename VectorType, class TRANSFER>
template<class OtherVectorType>
void
PreconditionMG<dim, VectorType, TRANSFER>::Tvmult
(OtherVectorType &,
 const OtherVectorType &) const
{
  Assert(false, ExcNotImplemented());
}


template<int dim, typename VectorType, class TRANSFER>
template<class OtherVectorType>
void
PreconditionMG<dim, VectorType, TRANSFER>::Tvmult_add
(OtherVectorType &,
 const OtherVectorType &) const
{
  Assert(false, ExcNotImplemented());
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif