---

GenAMRSolverImplem.H

Go to the documentation of this file.

/* _______              __
  / ___/ /  ___  __ _  / /  ___
 / /__/ _ \/ _ \/  ' \/ _ \/ _ \
 \___/_//_/\___/_/_/_/_.__/\___/
*/
//
// This software is copyright (C) by the Lawrence Berkeley
// National Laboratory.  Permission is granted to reproduce
// this software for non-commercial purposes provided that
// this notice is left intact.
//
// It is acknowledged that the U.S. Government has rights to
// this software under Contract DE-AC03-765F00098 between
// the U.S.  Department of Energy and the University of
// California.
//
// This software is provided as a professional and academic
// contribution for joint exchange. Thus it is experimental,
// is provided ``as is'', with no warranties of any kind
// whatsoever, no support, no promise of updates, or printed
// documentation. By using this software, you acknowledge
// that the Lawrence Berkeley National Laboratory and
// Regents of the University of California shall have no
// liability with respect to the infringement of other
// copyrights by any part of this software.
//

#ifndef _GENAMRSOLVERIMPLEM_H_
#define _GENAMRSOLVERIMPLEM_H_

#include "LayoutIterator.H"

/*****************/
/*****************/
template <class T>
GenAMRSolver<T>::GenAMRSolver():m_amrmgLevel()
{
  setDefaultValues();
}

/*****************/
/*****************/
template <class T>
GenAMRSolver<T>::GenAMRSolver(const Vector<DisjointBoxLayout>& a_gridsLevel,
                              const Vector<Box>&               a_domainLevel,
                              const Vector<Real>&              a_dxLevel,
                              const Vector<int>&               a_refRatio,
                              int                              a_numLevels,
                              int                              a_lBase,
                              const GenAMRLevelMGOp<T>* const  a_opin,
                              int                              a_ncomp)
:m_amrmgLevel()
{
  setDefaultValues();

  Vector<ProblemDomain> physdomains(a_domainLevel.size());

  for (int lev = 0; lev < physdomains.size(); lev++)
  {
    physdomains[lev] = ProblemDomain(a_domainLevel[lev]);
  }

  define(a_gridsLevel, physdomains, a_dxLevel,
         a_refRatio, a_numLevels, a_lBase, a_opin, a_ncomp);
}

/*****************/
/*****************/
template <class T>
GenAMRSolver<T>::GenAMRSolver(const Vector<DisjointBoxLayout>& a_gridsLevel,
                              const Vector<ProblemDomain>&     a_domainLevel,
                              const Vector<Real>&              a_dxLevel,
                              const Vector<int>&               a_refRatio,
                              int                              a_numLevels,
                              int                              a_lBase,
                              const GenAMRLevelMGOp<T>* const  a_opin,
                              int                              a_ncomp)
:m_amrmgLevel()
{
  setDefaultValues();

  define(a_gridsLevel, a_domainLevel, a_dxLevel,
         a_refRatio, a_numLevels, a_lBase, a_opin, a_ncomp);
}

/*****************/
/*****************/
template <class T>
GenAMRSolver<T>::~GenAMRSolver()
{
  clear();
}

/*****************/
/*****************/
template <class T>
bool GenAMRSolver<T>::isDefined() const
{
  return m_isDefined;
}

/*****************/
/*****************/
template <class T>
void GenAMRSolver<T>::define(const Vector<DisjointBoxLayout>& a_gridsLevel,
                             const Vector<Box>&               a_domainLevel,
                             const Vector<Real>&              a_dxLevel,
                             const Vector<int>&               a_refRatio,
                             int                              a_numLevels,
                             int                              a_lBase,
                             const GenAMRLevelMGOp<T>* const  a_opin,
                             int                              a_ncomp)
{
  setDefaultValues();

  Vector<ProblemDomain> physdomains(a_domainLevel.size());

  for (int lev = 0; lev < physdomains.size(); lev++)
  {
    physdomains[lev] = ProblemDomain(a_domainLevel[lev]);
  }

  define(a_gridsLevel, physdomains, a_dxLevel, a_refRatio, a_numLevels,
         a_lBase, a_opin, a_ncomp);
}

/*****************/
/*****************/
template <class T>
void GenAMRSolver<T>::define(const Vector<DisjointBoxLayout>& a_gridsLevel,
                             const Vector<ProblemDomain>&     a_domainLevel,
                             const Vector<Real>&              a_dxLevel,
                             const Vector<int>&               a_refRatio,
                             int                              a_numLevels,
                             int                              a_lBase,
                             const GenAMRLevelMGOp<T>* const  a_opin,
                             int                              a_ncomp)
{
  clear();

  m_isDefined = true;

  m_gridsLevel = a_gridsLevel;
  m_domainLevel = a_domainLevel;
  m_dxLevel = a_dxLevel;
  m_refRatio = a_refRatio;

  m_numLevels = a_numLevels;
  m_finestLevel = a_numLevels - 1;

  m_lBase = a_lBase;
  m_ncomp = a_ncomp;

  // define amrlevelmgs
  m_amrmgLevel.resize(m_numLevels, NULL);

  // (dfm -- 12/10/01) if lBase > 1, then we only need to define levels
  // from lBase-1 (needed for c-f boundary conditions)
  int startLev = a_lBase;

  // if (startLev>0) startLev--;
  for (int ilev = startLev; ilev < m_numLevels; ilev++)
  {
    m_amrmgLevel[ilev] = new GenAMRLevelMG<T>(this, ilev, a_opin, m_ncomp);
  }

  assert(m_lBase >= 0);

  // define levelsolver stuff
  const DisjointBoxLayout& levsolvGrids = m_gridsLevel[m_lBase];
  const ProblemDomain&     levsolvDom   = m_domainLevel[m_lBase];
  const Real&              levsolvDx    = m_dxLevel[m_lBase];
  const DisjointBoxLayout* levsolvBase = NULL;

  // int levsolvRef = 2;
  int levsolvRef = -1;

  if (m_lBase > 0)
  {
    levsolvBase = &m_gridsLevel[m_lBase-1];
    levsolvRef = m_refRatio[m_lBase-1];
  }

  m_levelSolver.define(levsolvGrids, levsolvBase,
                       levsolvDom,   levsolvDx,
                       levsolvRef, a_opin, m_ncomp);

  m_levelSolver.setMaxIter(m_numVCyclesBottom);
}

template <class T>
void GenAMRSolver<T>::setNumSmoothUp(int a_numSmoothUp)
{
  assert(a_numSmoothUp >= 0);

  m_numSmoothUp = a_numSmoothUp;

  for (int ilev = 0; ilev < m_numLevels; ilev++)
  {
    m_amrmgLevel[ilev]->setnumSmoothUp(m_numSmoothUp);
  }
}

template <class T>
void GenAMRSolver<T>::setNumSmoothDown(int a_numSmoothDown)
{
  assert(a_numSmoothDown >= 0);

  m_numSmoothDown = a_numSmoothDown;

  for (int ilev = 0; ilev < m_numLevels; ilev++)
  {
    m_amrmgLevel[ilev]->setnumSmoothDown(m_numSmoothDown);
  }
}

template <class T>
void GenAMRSolver<T>::setTolerance(Real a_tolerance)
{
  assert(a_tolerance >= 0);

  m_tolerance = a_tolerance;
  m_errorTolerance = a_tolerance;
}

template <class T>
void GenAMRSolver<T>::setOperatorTolerance(Real a_tolerance)
{
  assert(a_tolerance >= 0);

  m_operatorTolerance = a_tolerance;
}

template <class T>
void GenAMRSolver<T>::setMaxIter(int a_maxIter)
{
  assert(a_maxIter >= 0);

  m_maxIter = a_maxIter;
}

template <class T>
void GenAMRSolver<T>::setMinIter(int a_minIter)
{
  assert(a_minIter >= 0);

  m_minIter = a_minIter;
}

template <class T>
void GenAMRSolver<T>::setNumVCyclesBottom(int a_numVCyclesBottom)
{
  assert(a_numVCyclesBottom >= 0);

  m_numVCyclesBottom = a_numVCyclesBottom;
  m_levelSolver.setMaxIter(m_numVCyclesBottom);
}

template <class T>
void GenAMRSolver<T>::solveAMR(Vector<T *>&       a_phiLevel,
                               const Vector<T *>& a_rhsLevel)
{
  assert(isDefined());
  assert(a_phiLevel.size() > m_finestLevel);
  assert(a_rhsLevel.size() > m_finestLevel);

  int ncomp = a_phiLevel[m_lBase]->nComp();

  if (ncomp > m_ncomp)
  {
    MayDay::Warning("GenAMRSolver::solveAMR -- phi.nComp > solver.ncomp");
    ncomp = m_ncomp;
  }

  Vector<T *> residVect(m_finestLevel+1,NULL);
  Vector<T *> corrVect(m_finestLevel+1,NULL);

  // allocate residual and correction arrays
  // initial guess at the solution is zero
  for (int ilev = m_lBase; ilev <= m_finestLevel; ilev++)
  {
    const DisjointBoxLayout& levelGrids = a_phiLevel[ilev]->getBoxes();

    residVect[ilev] = new T(levelGrids, ncomp);
    corrVect[ilev]  = new T(levelGrids, ncomp, IntVect::Unit);

    T& phiLev = *a_phiLevel[ilev];
    T& corrLev = *corrVect[ilev];

    GenAMRLevelMGOp<T>& levelopCur = *(m_amrmgLevel[ilev]->m_levelopPtr);

    levelopCur.setToZero(phiLev);
    levelopCur.setToZero(corrLev);
  }

  // if m_lBase > 0, need to also create lBase-1 level correction for
  // cf BC's (set to 0)
  if (m_lBase > 0)
  {
    const DisjointBoxLayout& crseGrids = a_phiLevel[m_lBase-1]->getBoxes();

    corrVect[m_lBase-1] = new T(crseGrids, ncomp,
                                                   IntVect::Unit);

    // set crse correction to 0 as well
    T& crseCorr = *corrVect[m_lBase-1];

    GenAMRLevelMGOp<T>& crseLevelop = *(m_amrmgLevel[m_lBase-1]->m_levelopPtr);
    crseLevelop.setToZero(crseCorr);
  }

  Vector<Real> initRes = computeResidualNorm(residVect, a_phiLevel,
                                             a_rhsLevel, 0);

  Vector<Real> currentRes = initRes;
  Vector<Real> oldRes = initRes;

  if (m_verbose)
  {
    pout() << "GenAMRSolver, lBase = " << m_lBase  << endl;
    for (int comp = 0; comp < ncomp; comp++)
    {
      pout() << "initial max(residual[" << comp << "]) = "
             << initRes[comp] << endl;
    }
  }

  // if initRes = 0, don't bother solving...
  bool done = true;

  for (int comp = 0; comp < ncomp; comp++)
  {
    if (initRes[comp] != 0)
    {
      done = false;
    }
  }

  if (done)
  {
    if (m_verbose)
    {
      pout() << "GenAMRSolver::solveAMR--initial residual == 0\n returning... \n";
    }
    return;
  }

  int iter = 0;

  while (!done)
  {
    AMRVCycleMG(corrVect, residVect);

    // this is a clumsy way to do this (increment
    // solution with correction, reset corr to 0,
    // then recompute residual), but it should at least
    // do the right thing
    for (int ilev = m_lBase; ilev <= m_finestLevel; ilev++)
    {
      T& levelCorr = *corrVect[ilev];
      T& levelPhi = *a_phiLevel[ilev];

      GenAMRLevelMGOp<T>& levelopCur = *(m_amrmgLevel[ilev]->m_levelopPtr);

      levelopCur.sum(levelPhi,levelCorr);
      levelopCur.setToZero(levelCorr);
    }

    oldRes = currentRes;
    currentRes = computeResidualNorm(residVect, a_phiLevel,
                                     a_rhsLevel, 0);

    iter++;

    done = true;

    for (int comp = 0; comp < ncomp; comp++)
    {
      // DFM (3/24/2000) second test is in case where convergence
      // hangs due to machine precision (or solvability) issues...
      // done = done && (iter > m_maxIter);

      done = done &&
           ((iter > m_maxIter)
         || ((iter > m_minIter) && (currentRes[comp] >
                                    oldRes[comp]*(1.0-m_operatorTolerance)))
         || (currentRes[comp] <= m_tolerance*initRes[comp]));

      if (m_verbose)
      {
        pout() << "iteration #  " << iter
               << ": Max(res[" << comp << "]) = "
               << currentRes[comp] << endl;
      }
    }
  } // end while (!done)

  for (int comp = 0; comp < ncomp; comp++)
  {
    if (currentRes[comp] <= m_tolerance*initRes[comp])
    {
      if (m_verbose)
      {
        pout() << "GenAMRSolver: " << iter << " iterations, final max(res["
               << comp << "]) = "
               << currentRes[comp] << endl;
      }
    }
    else if (iter > m_maxIter)
    {
      if (m_verbose)
      {
        pout() << "GenAMRSolver: reached maximum number of "
               << iter << " iterations, final max(res)[" << comp
               << "] = " << currentRes[comp] << endl;
      }
    }
    else if (currentRes[comp] > oldRes[comp]*(1.0-m_operatorTolerance))
    {
      if (m_verbose)
      {
        pout() << "GenAMRSolver: reached solver hang point; " << iter
               << " iterations, final max(res[" << comp
               << "]) = " << currentRes[comp] << endl;
      }
    }
    else
    {
      if (m_verbose)
      {
        pout() << "GenAMRSolver NOT CONVERGED! - final max(res) = "
               << currentRes[comp] << " after " << iter << " iterations" << endl;
      }
    }
  }

  // clean up storage
  for (int ilev = 0; ilev <= m_finestLevel; ilev++)
  {
    if (residVect[ilev] != NULL)
    {
      delete residVect[ilev];
      residVect[ilev] = NULL;
    }

    if (corrVect[ilev] != NULL)
    {
      delete corrVect[ilev];
      corrVect[ilev] = NULL;
    }
  }
}

template <class T>
void GenAMRSolver<T>::AMRVCycleMG(Vector<T *>&       a_corrLevel,
                                  const Vector<T *>& a_residLevel)
{
  assert(isDefined());
  assert(a_corrLevel.size() > m_finestLevel);
  assert(a_residLevel.size() > m_finestLevel);

  // this is a kluge to make this work in residual-correction form
  // copy a_residLevel into each level's m_resid
  Interval compInterval = a_residLevel[m_lBase]->interval();
  for (int ilev = m_lBase; ilev <= m_finestLevel; ilev++)
  {
    T& levelResid = m_amrmgLevel[ilev]->m_resid;

    a_residLevel[ilev]->copyTo(compInterval, levelResid, compInterval);
  }

  // sweep down vcycle
  // for (int ilev = m_finestLevel; ilev >= m_lBase; ilev--)
  //{
  //  // m_amrmgLevel[ilev]->setArrayViewVerbose(true);
  //}

  // shouldn't need this anymore
  // m_amrmgLevel[m_finestLevel]->computeAMRResidual(a_phiLevel,a_rhsLevel);

  // for (int ilev = m_finestLevel; ilev >= m_lBase; ilev--)
  //{
  //  // m_amrmgLevel[ilev]->setArrayViewVerbose(true);
  //}

  for (int ilev = m_finestLevel; ilev > m_lBase; ilev--)
  {
    m_amrmgLevel[ilev]->downSweep(a_corrLevel,a_residLevel);
  }

  // solve at level lBase
  T & bottomCorr = m_amrmgLevel[m_lBase]->m_corr;
  const T & bottomRes = m_amrmgLevel[m_lBase]->m_resid;

  // this only does m_numVCyclesBottom iterations
  if (m_finestLevel == m_lBase)
  {
    // if only one level is being solved, it makes sense to
    // do bottom smoothing on resid/phi directly
    m_levelSolver.levelSolveH(*a_corrLevel[m_lBase],
                              *a_residLevel[m_lBase]);
  }
  else
  {
    m_levelSolver.levelSolveH(bottomCorr, bottomRes);

    // add correction to lBase's phi
    T & bottomPhi = *a_corrLevel[m_lBase];
    GenAMRLevelMGOp<T>& bottomLevelop = *(m_amrmgLevel[m_lBase]->m_levelopPtr);

    bottomLevelop.sum(bottomPhi,bottomCorr);
  }

  // sweep back up vcycle
  for (int ilev = m_lBase+1; ilev <= m_finestLevel; ilev++)
  {
    m_amrmgLevel[ilev]->upSweep(a_corrLevel,a_residLevel);
  }
}

template <class T>
Vector<Real> GenAMRSolver<T>::computeResidualNorm(Vector<T *>&       a_phiLevel,
                                                  const Vector<T *>& a_rhsLevel,
                                                  int                a_normType)
{
  assert(isDefined());
  assert(a_phiLevel.size() > m_finestLevel);
  assert(a_rhsLevel.size() > m_finestLevel);
  assert((a_normType >= 0) && (a_normType <= 2));

  int ncomp = a_phiLevel[m_lBase]->nComp();

  if (ncomp > m_ncomp)
  {
    ncomp = m_ncomp;
  }

  Vector<Real> normTot(ncomp,0);
  Vector<Real> normLevel(ncomp,0.0);

  for (int ilev = m_finestLevel; ilev >= m_lBase; ilev--)
  {
    m_amrmgLevel[ilev]->setArrayViewVerbose(false);
    m_amrmgLevel[ilev]->computeAMRResidual(a_phiLevel, a_rhsLevel);

    normLevel = m_amrmgLevel[ilev]->computeResidualNorm(a_normType);

    for (int comp = 0; comp < ncomp; comp++)
    {
      if (a_normType == 0)
      {
        normTot[comp] = Max(normTot[comp], normLevel[comp]);
      }
      else
      {
        normTot[comp] += normLevel[comp];
      }
    }
  } // end loop over levels

  if (a_normType == 2)
  {
    for (int comp = 0; comp < ncomp; comp++)
    {
      normTot[comp] = sqrt(normTot[comp]);
    }
  }

  return normTot;
}

template <class T>
Vector<Real> GenAMRSolver<T>::computeResidualNorm(Vector<T *>&       a_residLevel,
                                                  Vector<T *>&       a_phiLevel,
                                                  const Vector<T *>& a_rhsLevel,
                                                  int                a_normType)
{
  assert(isDefined());
  assert(a_phiLevel.size() > m_finestLevel);
  assert(a_rhsLevel.size() > m_finestLevel);
  assert((a_normType >= 0) && (a_normType <= 2));

  int ncomp = a_phiLevel[m_lBase]->nComp();
  if (ncomp > m_ncomp)
  {
    ncomp = m_ncomp;
  }

  Vector<Real> normTot(ncomp,0);

  Vector<Real> normLevel(ncomp,0.0);

  for (int ilev = m_finestLevel; ilev >= m_lBase; ilev--)
  {
    T& levelResid = *a_residLevel[ilev];
    m_amrmgLevel[ilev]->computeAMRResidual(levelResid, a_phiLevel, a_rhsLevel);

    normLevel = m_amrmgLevel[ilev]->computeNorm(levelResid,a_normType);

    for (int comp = 0; comp < ncomp; comp++)
    {
      if (a_normType == 0)
      {
        normTot[comp] = Max(normTot[comp], normLevel[comp]);
      }
      else
      {
        normTot[comp] += normLevel[comp];
      }
    }
  } // end loop over levels

  if (a_normType == 2)
  {
    for (int comp = 0; comp < ncomp; comp++)
    {
      normTot[comp] = sqrt(normTot[comp]);
    }
  }

  return normTot;
}

template <class T>
void GenAMRSolver<T>::computeAMRResidual(Vector<T *>&       a_phiLevel,
                                         const Vector<T *>& a_rhsLevel,
                                         T&                 a_res,
                                         int                a_ilev)
{
  assert(isDefined());
  assert(a_ilev <= m_finestLevel);
  assert(a_ilev >= 0);
  assert(a_phiLevel.size() > m_finestLevel);
  assert(a_rhsLevel.size() > m_finestLevel);
  assert(a_rhsLevel[a_ilev]->getBoxes() == m_gridsLevel[a_ilev]);
  assert(a_phiLevel[a_ilev]->getBoxes() == m_gridsLevel[a_ilev]);
  assert(a_res.getBoxes() == m_gridsLevel[a_ilev]);

  // compute residual internal to amrmgLevel
  m_amrmgLevel[a_ilev]->computeAMRResidual(a_phiLevel,a_rhsLevel);
  const T & amrmgRes = m_amrmgLevel[a_ilev]->m_resid;

  assert(amrmgRes.getBoxes() == m_gridsLevel[a_ilev]);

  // copy residual into output array
  GenAMRLevelMGOp<T>& levelop = *(m_amrmgLevel[a_ilev]->m_levelopPtr);

  levelop->copy(a_res,amrmgRes);
}

template <class T>
void GenAMRSolver<T>::applyAMROperator(Vector<T *>& a_phiLevel,
                                       T&           a_LOfPhi,
                                       int          a_ilev)
{
  assert(isDefined());
  assert(a_ilev <= m_finestLevel);
  assert(a_ilev >= 0);
  assert(a_phiLevel.size() > m_finestLevel);

  m_amrmgLevel[a_ilev]->applyAMROperator(a_phiLevel,a_LOfPhi);
}

template <class T>
void GenAMRSolver<T>::applyAMROperatorHphys(Vector<T* >& a_phiLevel,
                                            T&           a_LOfPhi,
                                            int          a_ilev)
{
  assert(isDefined());
  assert(a_ilev <= m_finestLevel);
  assert(a_ilev >= 0);
  assert(a_phiLevel.size() > m_finestLevel);

  m_amrmgLevel[a_ilev]->applyAMROperatorHphys(a_phiLevel,a_LOfPhi);
}

template <class T>
void GenAMRSolver<T>::setVerbose(bool a_verbose)
{
  assert(isDefined());

  m_verbose = a_verbose;
}

/*****************/
/*****************/
template <class T>
void GenAMRSolver<T>::setDefaultValues()
{
  m_numLevels = -1;
  m_finestLevel = -1;
  m_refRatio.resize(0);
  m_numVCyclesBottom = 1;

  m_tolerance = 1.0e-10;
#ifdef CH_USE_FLOAT
  m_tolerance = 10.*sqrt(m_tolerance);
#endif
  m_errorTolerance = m_tolerance;
  m_operatorTolerance = 1.0e-4;

  m_maxIter = 42;
  m_minIter = 4;

  m_numSmoothUp = 4;
  m_numSmoothDown = 4;

  m_isDefined = false;

  m_ncomp = 1;

  m_verbose = false;
}

/*****************/
/*****************/
template <class T>
void GenAMRSolver<T>::clear()
{
  for (int ilev = 0; ilev < m_amrmgLevel.size(); ilev++)
  {
    if (m_amrmgLevel[ilev] != NULL)
    {
      delete m_amrmgLevel[ilev];
    }
  }
}

#endif

---