---

GenLevelSolverImplem.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 _GENLEVELSOLVERIMPLEM_H_
#define _GENLEVELSOLVERIMPLEM_H_

#include "IntVectSet.H"
#include "parstream.H"
#include "LayoutIterator.H"

// using std::cout;
// using std::endl;

// default constructor
template <class T>
GenLevelSolver<T>::GenLevelSolver()
{
  setDefaultValues();
}

// complete constructor
template <class T>
GenLevelSolver<T>::GenLevelSolver(const DisjointBoxLayout&     a_grids,
                                  const DisjointBoxLayout*     a_baseGrids,
                                  const Box&                   a_domain,
                                  Real                         a_dxLevel,
                                  int                          a_nRefCrse,
                                  const GenLevelMGOp<T>* const a_opin,
                                  int                          a_ncomp)
{
  setDefaultValues();

  ProblemDomain physdomain(a_domain);
  define(a_grids,  a_baseGrids, physdomain, a_dxLevel,
         a_nRefCrse, a_opin, a_ncomp);
}

// complete constructor
template <class T>
GenLevelSolver<T>::GenLevelSolver(const DisjointBoxLayout&     a_grids,
                                  const DisjointBoxLayout*     a_baseGrids,
                                  const ProblemDomain&         a_domain,
                                  Real                         a_dxLevel,
                                  int                          a_nRefCrse,
                                  const GenLevelMGOp<T>* const a_opin,
                                  int                          a_ncomp)
{
  setDefaultValues();

  define(a_grids,  a_baseGrids, a_domain, a_dxLevel,
         a_nRefCrse, a_opin, a_ncomp);
}

// destructor
template <class T>
GenLevelSolver<T>::~GenLevelSolver()
{
  clearMemory();
}

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

template <class T>
void GenLevelSolver<T>::define(const DisjointBoxLayout&     a_grids,
                               const DisjointBoxLayout*     a_baseGrids,
                               const Box&                   a_domain,
                               Real                         a_dxLevel,
                               int                          a_nRefCrse,
                               const GenLevelMGOp<T>* const a_opin,
                               int                          a_ncomp)
{
  ProblemDomain physdomain(a_domain);
  define(a_grids, a_baseGrids, a_domain, a_dxLevel, a_nRefCrse,
         a_opin, a_ncomp);
}

template <class T>
void GenLevelSolver<T>::define(const DisjointBoxLayout&     a_grids,
                               const DisjointBoxLayout*     a_baseGrids,
                               const ProblemDomain&         a_domain,
                               Real                         a_dxLevel,
                               int                          a_nRefCrse,
                               const GenLevelMGOp<T>* const a_opin,
                               int                          a_ncomp)
{
  clearMemory();

  m_grids = a_grids;
  m_domain = a_domain;
  m_dxLevel = a_dxLevel;
  m_nRefCrse = a_nRefCrse;

  m_resid.define(m_grids, a_ncomp, IntVect::TheZeroVector());
  m_scratch.define(m_grids, a_ncomp, IntVect::TheZeroVector());
  m_corr.define(m_grids, a_ncomp, IntVect::TheUnitVector());

  m_isDefined = true;

  m_levelOpPtr = a_opin->newOp();
#if FIXIT
  m_levelOpPtr->define(m_grids, a_baseGrids, m_dxLevel,
                       m_nRefCrse, m_domain, false, a_ncomp);
#endif

  // coarsen boxes once to start, because coarsest level boxes
  // must also be coarsenable (for proper coarse-fine interpolation)
  // however, if grid merely covers entire physical domain, then
  // want to coarsen all the way down
  int nCoarserLevels = -1;

  // (DFM 10/21/02) don't bother doing _any_ of this if m_grids
  // is an empty set
  if (m_grids.size() > 0)
  {
    // why do it this way rather than as a dense IVS?
    IntVectSet ivsCheck(m_domain.domainBox());
    LayoutIterator lit = m_grids.layoutIterator();

    for (lit.reset(); lit.ok(); ++lit)
    {
      ivsCheck -= m_grids.get(lit());
    }

    if (ivsCheck.isEmpty())
    {
      nCoarserLevels += 1;
    }

    int refrattot = 2;
    bool getOut = false;
    while(!getOut)
    {
      bool addone = true;
      for (lit.reset(); lit.ok(); ++lit)
      {
        Box blocal = m_grids.get(lit());

        if (refine(coarsen(blocal,refrattot),refrattot) != blocal)
        {
          addone = false;
        }

        if (!addone)
        {
          break;
        }
      }

      if (addone)
      {
        nCoarserLevels += 1;
        refrattot *= 2;
      }
      else
      {
        getOut = true;
      }
    }
  }

  nCoarserLevels = Max(nCoarserLevels, 0);

  m_levelMG.define(m_grids, a_baseGrids, m_dxLevel,
                   m_nRefCrse, m_domain, nCoarserLevels, a_opin, a_ncomp);
}

// returns GenLevelSolver to a basically undefined state
template <class T>
void GenLevelSolver<T>::clearMemory()
{
  if (m_levelOpPtr != NULL)
  {
    delete m_levelOpPtr;
    m_levelOpPtr = NULL;
  }
}

template <class T>
void GenLevelSolver<T>::setDefaultValues()
{
  m_levelOpPtr = NULL;
  m_isDefined = false;
  // set this to an invalid value
  m_nRefCrse = -1;
  m_dxLevel = -1.0;

  m_maxIter = 33;
  m_minIter = 4;

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

  m_verbose = false;
}

// solve on just this level, homogeneous bcs
template <class T>
void GenLevelSolver<T>::levelSolveH(T&       a_phi,
                                    const T& a_rhs,
                                    bool     a_initializePhiToZero)
{
  assert(isDefined());
  assert(a_phi.nComp() == a_rhs.nComp());
  assert(a_phi.nComp() == m_resid.nComp());

  assert(m_levelOpPtr->isDefined());

  if (a_initializePhiToZero)
  {
    // initial guess of phi is zero
    m_levelOpPtr->setToZero(a_phi);
  }

  m_levelOpPtr->applyOpH(a_phi, m_resid);
  m_levelOpPtr->diff(m_resid,a_rhs);
  m_levelOpPtr->negate(m_resid);

  m_levelOpPtr->setToZero(m_corr);

  Vector<Real> initRes(m_resid.nComp());
  Vector<Real> oldRes(m_resid.nComp());

  bool done = true;

  for (int comp = 0; comp < m_resid.nComp(); comp++)
  {
#if FIXIT
    Interval thisInterval(comp,comp);
    initRes[comp] = norm(m_resid, thisInterval, 0);
    // also initialize oldRes to be initRes;
    oldRes[comp] = initRes[comp];
#else
    initRes[comp] = 0.0;
#endif

// if initRes = 0, don't bother solving
    if (initRes[comp] > 0.01*m_tolerance)
    {
      done = false;
    }
  }

  if (done)
  {
    return;
  }

  int iter = 0;
  while (!done)
  {
    iter++;

    // compute modified residual LofPhi = res - LNf(corr)
    // use homogeneous form of ApplyOp because want to use
    // homogeneous CF BC's for correction
    m_levelMG.mgRelax(m_corr, m_resid, true);

    m_levelOpPtr->applyOpH(m_corr, m_scratch);
    m_levelOpPtr->diff(m_scratch,m_resid);
    m_levelOpPtr->negate(m_scratch);

    done = true;
    for (int comp = 0; comp < m_scratch.nComp(); comp++)
    {
#if FIXIT
      Interval thisInterval(comp,comp);
      Real currentRes = norm(m_scratch, thisInterval, 0);
#else
      Real currentRes = 0.0;
#endif

      if (m_verbose)
      {
        pout() << "iteration #= " << iter
               << ": Max(res) = " << currentRes
               << endl;
      }

      // DFM (9/20/02) third test is to catch case where
      // convergence hangs due to machine precision (or
      // solvability) issues.  This follows the approach
      // used in AMRSolver
      done = done &&
        ((iter > m_maxIter)
      || (currentRes <= m_tolerance*initRes[comp])
      || ((iter < m_minIter) && ((currentRes/oldRes[comp]) > 1.0-m_operatorTolerance)));

      // save currentRes into oldRes
      oldRes[comp] = currentRes;
    } // end loop over components
  } // end solver iterations

  // update phi
  m_levelOpPtr->sum(a_phi,m_corr);
}

// solve on just this level, inhomogeneous bcs
template <class T>
void GenLevelSolver<T>::levelSolve(T&       a_phi,
                                   const T* a_phic,
                                   const T& a_rhs,
                                   bool     a_initializePhiToZero)
{
  assert(isDefined());
  assert(a_phi.nComp() == a_rhs.nComp());
  assert(a_phi.nComp() == m_resid.nComp());

  /* compute initial residual
     -- can't use residual() fn because
     don't want to zero out higher levels
     -- also note that we use LNf,
     because this is a LEVEL solve
  */

  assert(m_levelOpPtr->isDefined());

  if (a_initializePhiToZero)
  {
    // initial guess of phi is zero
    m_levelOpPtr->setToZero(a_phi);
  }

  m_levelOpPtr->applyOpI(a_phi, a_phic, m_resid);
  m_levelOpPtr->negDiff(m_resid,a_rhs);

  m_levelOpPtr->setToZero(m_corr);

  Vector<Real> initRes(m_resid.nComp());
  Vector<Real> oldRes(m_resid.nComp());

  bool done = true;

  for (int comp = 0; comp < m_resid.nComp(); comp++)
  {
    Interval thisInterval(comp,comp);
#if FIXIT
    initRes[comp] = norm(m_resid, thisInterval, 0);
    // initialize oldRes to be initRes as well
    oldRes[comp] = initRes[comp];
#endif

    // if initRes = 0, don't bother solving
    if (initRes[comp] > 0.01*m_tolerance)
    {
      done = false;
    }
  }

  if (done)
  {
    return;
  }

  int iter = 0;
  while (!done)
  {
    m_levelMG.mgRelax(m_corr, m_resid, true);

    m_levelOpPtr->applyOpH(m_corr, m_scratch);
    m_levelOpPtr->negDiff(m_scratch,m_resid);

    iter++;

    done = true;
    for (int comp = 0; comp < m_resid.nComp(); comp++)
    {
      Real currentRes = 0.0;
      Interval thisInterval(comp,comp);

#if FIXIT
      currentRes = norm(m_scratch, thisInterval, 0);
#endif

      if (m_verbose)
      {
        pout() << "iteration #= " << iter
               << ": Max(res) = " << currentRes
               << endl;
      }

      // DFM (9/20/02) third test is to catch case where
      // convergence hangs due to machine precision (or
      // solvability) issues.  This follows the approach
      // used in AMRSolver
      done = done &&
        ((iter > m_maxIter)
      || (currentRes <= m_tolerance*initRes[comp] )
      || ((iter > m_minIter) && ((currentRes/oldRes[comp]) > 1.0-m_operatorTolerance)));
      // save currentRes into oldRes
      oldRes[comp] = currentRes;
    } // end loop over components
  } // end solver iterations

  m_levelOpPtr->sum(a_phi,m_corr);
}

template <class T>
void GenLevelSolver<T>::setTolerance(Real a_tolerance)
{
  m_tolerance = a_tolerance;
}

// set tolerance for stopping due to hung convergence
template <class T>
void GenLevelSolver<T>::setOperatorTolerance(Real a_operatorTolerance)
{
  assert(a_tolerance >= 0.0);

  m_operatorTolerance = a_operatorTolerance;
}

template <class T>
void GenLevelSolver<T>::setVerbose(bool a_verbose)
{
  m_verbose = a_verbose;
}

template <class T>
void GenLevelSolver<T>::setMaxIter(int a_maxIter)
{
  m_maxIter = a_maxIter;
}

template <class T>
void GenLevelSolver<T>::setMinIter(int a_minIter)
{
  m_minIter = a_minIter;
}

#endif

---