MultiGrid.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

// BVS, June 18, 2003

// We can assume that template class T has null construction.

#ifndef _MULTIGRID_H_
#define _MULTIGRID_H_

#include "CH_Timer.H"
#include "LinearSolver.H"
#include "ProblemDomain.H"
#include "REAL.H"
#include "Box.H"
#include "parstream.H"
#include "NamespaceHeader.H"


template <class T>
class MGLevelOp : public LinearOp<T>
{
public:

  virtual ~MGLevelOp(){;}


  virtual void createCoarser(T& a_coarse, const T& a_fine, bool ghosted) = 0;


  virtual void relax(T& a_correction, const T& a_residual, int a_iterations) = 0 ;



  virtual void restrictResidual(T& a_resCoarse, T& a_phiFine, const T& a_rhsFine) = 0;


  virtual void prolongIncrement(T& a_phiThisLevel, const T& a_correctCoarse) = 0;

};


template <class T>
class MGLevelOpFactory
{
public:
  virtual ~MGLevelOpFactory(){;}


  virtual MGLevelOp<T>* MGnewOp(const ProblemDomain& a_FineindexSpace,
                                int                  a_depth,
                                bool                 a_homoOnly = true) = 0;
};



template <class T>
class MultiGrid
{
public:
  MultiGrid();
  virtual ~MultiGrid();


  virtual void define(MGLevelOpFactory<T>& a_factory,
                      LinearSolver<T>* a_bottomSolver,
                      const ProblemDomain& a_domain,
                      int  a_maxDepth = -1);


  virtual void solve(T& a_phi, const T& a_rhs, Real a_tolerance, int a_maxIterations, int verbosity= 0);


  virtual void oneCycle(T& a_e, const T& a_res);


  // used by AMRMultiGrid--not a part of the public API
  void init(const T& a_correction, const T& a_residual);

  //internal function
  void cycle(int a_depth, T& a_correction, const T& a_residual);

  //internal function
  void clear();

  // used by AMRMultiGrid--not a part of the public API
  void setBottomSolver(LinearSolver<T>* a_bottomSolver);


  int  m_depth, m_defaultDepth, m_pre, m_post, m_bottom, m_cycle, m_numMG;
  bool m_homogeneous;
  LinearSolver<T>*  m_bottomSolver;


  Vector< MGLevelOp<T>* > getAllOperators();

protected:
  bool m_defined;

  int m_bottomCells;

  ProblemDomain m_topLevelDomain;
  Vector<MGLevelOp<T>*> m_op;
  Vector< T* >   m_residual;
  Vector< T* >   m_correction;
private:
  MultiGrid(const MultiGrid<T>& a_opin)
  {
    MayDay::Error("invalid operator");
  }

  void operator=(const MultiGrid<T>& a_opin)
  {
    MayDay::Error("invalid operator");
  }
};

// *******************************************************
// MultiGrid Implementation
// *******************************************************

template <class T>
Vector< MGLevelOp<T>* >
MultiGrid<T>::getAllOperators()
{
  Vector<MGLevelOp<T>*> retval = m_op;
  return retval;
}

template <class T>
MultiGrid<T>::MultiGrid():
  m_pre(3), m_post(3), m_bottom(0), m_cycle(1), m_numMG(1), m_homogeneous(true),
  m_defined(false){;}

template <class T>
MultiGrid<T>::~MultiGrid()
{
  clear();
}

template <class T>
void MultiGrid<T>::define(MGLevelOpFactory<T>& a_factory,
                          LinearSolver<T>*     a_bottomSolver,
                          const ProblemDomain& a_domain,
                          int                  a_maxDepth)
{

  CH_TIME("MultiGrid::define");
  m_depth = 0;
  m_bottomSolver = a_bottomSolver;
  m_residual.resize(m_depth);   // the zero'th entry is left null constructed
  m_correction.resize(m_depth); // the zero'th entry is left null constructed
  m_op.resize(m_depth);
  m_topLevelDomain = a_domain;

  m_bottomCells = a_domain.domainBox().numPts();

  MGLevelOp<T>* nextOp =  a_factory.MGnewOp(a_domain, 0, true);
  while(nextOp != NULL)
    {
      m_residual.push_back(new T());
      m_correction.push_back(new T());
      m_op.push_back(nextOp);
      m_depth++;
      if((m_depth < a_maxDepth) || (a_maxDepth < 0))
        {
          nextOp = a_factory.MGnewOp(a_domain , m_depth, true);
          if(nextOp != NULL)
            for(int i=0; i<CH_SPACEDIM; ++i) m_bottomCells /= 2;
        }
      else
        {
          nextOp = NULL;
        }
    }
  m_defaultDepth = m_depth;
  m_bottomSolver->define(m_op[m_depth-1], true);
  m_defined = true;
}

template <class T>
void MultiGrid<T>::setBottomSolver(LinearSolver<T>* a_bottomSolver)
{
  m_bottomSolver = a_bottomSolver;
  m_bottomSolver->define(m_op[m_depth-1], true);
}

template <class T>
void MultiGrid<T>::init(const T& a_e, const T& a_residual)
{
  CH_TIME("MutliGrid::init");
  if (m_depth > 1)
    {
      m_op[0]->createCoarser(*(m_residual[1]), a_residual, false);
      m_op[0]->createCoarser(*(m_correction[1]), a_e, true);
    }

  for (int i = 2; i < m_depth; i++)
    {
      m_op[i-1]->createCoarser(*(m_residual[i]), *(m_residual[i-1]), false);
      m_op[i-1]->createCoarser(*(m_correction[i]), *(m_correction[i-1]), true);
    }
}

template <class T>
void MultiGrid<T>::solve(T& a_phi, const T& a_rhs, Real a_tolerance, int a_maxIterations, int verbosity)
{
  CH_TIME("MultiGrid::solve");
  this->init(a_phi, a_rhs);

  T correction, residual;
  m_op[0]->create(correction, a_phi);
  m_op[0]->create(residual, a_rhs);
  m_op[0]->setToZero(a_phi);
  m_op[0]->residual(residual, a_phi, a_rhs, false);

  Real errorno = m_op[0]->norm(residual, 0);
  if(verbosity > 2)
    {
      pout() << "multigrid::solve initial residual = " << errorno << std::endl;
    }
  Real compval = a_tolerance*errorno;
  Real epsilon = 1.0e-16;
  compval = Max(compval, epsilon);
  Real error = errorno;
  int iter = 0;
  while((error > compval) && (error > a_tolerance*errorno) && (iter < a_maxIterations))
    {
      m_op[0]->setToZero(correction);
      m_op[0]->residual(residual, a_phi, a_rhs, false);
      error = m_op[0]->norm(residual, 0);
      if(verbosity > 3)
        {
          pout() << "multigrid::solve iter = " << iter <<  ",  residual = " << error << std::endl;
        }

      this->cycle(0, correction, residual);
      m_op[0]->incr(a_phi, correction, 1.0);

      iter++;
    }
  if(verbosity > 2)
    {
      pout() << "multigrid::solve final residual = " << error << std::endl;
    }

  m_op[0]->clear(correction);
  m_op[0]->clear(residual);
}

template <class T>
void MultiGrid<T>::oneCycle(T& a_e, const T& a_residual)
{
  //this->init(a_e, a_residual); oneCycle is only used by AMRMultigrid
  if (m_homogeneous)
    {
      cycle(0, a_e, a_residual);
    }
  else
    {
      T correction, residual;
      m_op[0]->create(correction, a_e);
      m_op[0]->create(residual, a_residual);

      m_op[0]->residual(residual, a_e, a_residual);

      m_op[0]->setToZero(correction);
      this->cycle(0, correction, residual);

      m_op[0]->incr(a_e, correction, 1.0);
      m_op[0]->clear(correction);
      m_op[0]->clear(residual);
    }
}

template <class T>
void MultiGrid<T>::cycle(int depth, T& correction, const T& residual)
{
  // currently I can't drop a timer in this function because it is recursive.  doh

  if(depth == m_depth-1)
    {
      if(m_bottomCells == 1)
        {
          m_op[depth  ]->relax(correction, residual, 1);
        }
      else
        {
          m_op[depth  ]->relax(correction, residual, m_bottom);
          m_bottomSolver->solve(correction, residual);
        }
    }
  else
    {
      int cycles = m_cycle;
      if( cycles < 0 ) 
        {
          cycles = -cycles;
          // F cycle multigrid
          m_op[depth  ]->restrictResidual(*(m_residual[depth+1]), correction, residual);
          m_op[depth  ]->setToZero(*(m_correction[depth+1]));

          // recursive call
          cycle(depth+1, *(m_correction[depth+1]), *(m_residual[depth+1]));

          m_op[depth  ]->prolongIncrement(correction, *(m_correction[depth+1]));

          m_op[depth  ]->relax(correction, residual, m_pre);

          for(int img = 0; img < cycles; img++)
            {
              m_op[depth  ]->restrictResidual(*(m_residual[depth+1]), correction, residual);
              m_op[depth  ]->setToZero(*(m_correction[depth+1]));

              m_cycle = 1;          // hack to get a V-cycle
              cycle(depth+1, *(m_correction[depth+1]), *(m_residual[depth+1]));
              m_cycle = -cycles;
              //
              m_op[depth  ]->prolongIncrement(correction, *(m_correction[depth+1]));
            }

          m_op[depth  ]->relax(correction, residual, m_post);
        }
      else 
        {
          m_op[depth  ]->relax( correction, residual, m_pre );

          m_op[depth  ]->restrictResidual(*(m_residual[depth+1]), correction, residual);
          m_op[depth  ]->setToZero(*(m_correction[depth+1]));

          for(int img = 0; img < cycles; img++)
            {
              cycle(depth+1, *(m_correction[depth+1]), *(m_residual[depth+1]));
              if(depth+1 == m_depth-1) break; // coarsest grid - 1 iteration
            }
          m_op[depth  ]->prolongIncrement(correction, *(m_correction[depth+1]));

          m_op[depth  ]->relax(correction, residual, m_post);
        }
    }
}

template <class T>
void MultiGrid<T>::clear()
{
  for(int i=0; i<m_op.size(); ++i)
    {
      delete m_op[i];
      delete m_correction[i];
      delete m_residual[i];
    }
  m_op.resize(0);
  m_residual.resize(0);
  m_correction.resize(0);
  m_defined = false;
}

#include "NamespaceFooter.H"
#endif /*_MULTIGRID_H_*/

---