---

GenLevelMGImplem.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 _GENLEVELMGIMPLEM_H_
#define _GENLEVELMGIMPLEM_H_

#include "GenLevelMGF_F.H"

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

// this calls the levelmg big define function
template <class T>
GenLevelMG<T>::GenLevelMG(const DisjointBoxLayout&     a_ba,
                          const DisjointBoxLayout*     a_baseBaPtr,
                          Real                         a_dxLevel,
                          int                          a_refRatio,
                          const Box&                   a_domain,
                          int                          a_nCoarserLevels,
                          const GenLevelMGOp<T>* const a_opin,
                          int                          a_ncomp)
{
  setDefaultValues();

  ProblemDomain physdomain(a_domain);
  define(a_ba,a_baseBaPtr, a_dxLevel, a_refRatio, physdomain,
         a_nCoarserLevels, a_opin, a_ncomp);
}

// this calls the levelmg big define function
template <class T>
GenLevelMG<T>::GenLevelMG(const DisjointBoxLayout&     a_ba,
                          const DisjointBoxLayout*     a_baseBaPtr,
                          Real                         a_dxLevel,
                          int                          a_refRatio,
                          const ProblemDomain&         a_domain,
                          int                          a_nCoarserLevels,
                          const GenLevelMGOp<T>* const a_opin,
                          int                          a_ncomp)
{
  setDefaultValues();

  define(a_ba,a_baseBaPtr, a_dxLevel, a_refRatio, a_domain,
         a_nCoarserLevels, a_opin, a_ncomp);
}

// --------------------------------------------------------------
// constructor for coarsened version of levelMG
// refcoarse is the ratio from coarsened thing to this thing
// --------------------------------------------------------------
template <class T>
GenLevelMG<T>::GenLevelMG(const GenLevelMG&      a_L,
                          int                    a_refCoarse,
                          const GenLevelMGOp<T>* a_opin)
{
  // set pointers to NULL so define does not try to delete them
  setDefaultValues();

  define(a_L, a_refCoarse, a_opin);
}

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

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

// corresponding define function
template <class T>
void GenLevelMG<T>::define(const DisjointBoxLayout&     a_ba,
                           const DisjointBoxLayout*     a_baseBaPtr,
                           Real                         a_dxLevel,
                           int                          a_refRatio,
                           const  Box&                  a_domain,
                           int                          a_nCoarserLevels,
                           const GenLevelMGOp<T>* const a_opin,
                           int                          a_ncomp)
{
  ProblemDomain physdomain(a_domain);
  define(a_ba, a_baseBaPtr, a_dxLevel, a_refRatio, a_domain,
         a_nCoarserLevels, a_opin, a_ncomp);
}

// corresponding define function
template <class T>
void GenLevelMG<T>::define(const DisjointBoxLayout&     a_ba,
                           const DisjointBoxLayout*     a_baseBaPtr,
                           Real                         a_dxLevel,
                           int                          a_refRatio,
                           const  ProblemDomain&        a_domain,
                           int                          a_nCoarserLevels,
                           const GenLevelMGOp<T>* const a_opin,
                           int                          a_ncomp)
{
  assert(a_opin != NULL);
  assert(a_nCoarserLevels >= 0);
  assert(!a_domain.isEmpty());
  assert(a_ba.checkPeriodic(a_domain));

  // either there are no coarser levels
  // or we have a legitimate ref ratio
  assert((a_refRatio > 0) || (a_baseBaPtr == NULL));
  assert(a_dxLevel > 0);

  clearMemory();

  m_ba = a_ba;
  m_baseBaPtr = a_baseBaPtr;

  m_dxLevel = a_dxLevel;
  m_domain  = a_domain;

  m_nCoarserLevels = a_nCoarserLevels;

  m_isDefined = true;
  m_levelopPtr = a_opin->newOp();
  // no need to do inhomgeneous CFInterp here
  // all levelmg operations (smooth, etc) use
  // only homogeneous cf interpolation
#if FIXIT
  bool homogeneousOnly = true;
  m_levelopPtr->define(m_ba,  m_baseBaPtr,
                       m_dxLevel,  a_refRatio,
                       m_domain, homogeneousOnly, a_ncomp);
#endif

  m_levelopPtr->defineData(m_resid);
  // recursively define MG hierarchy
  if (m_nCoarserLevels != 0)
  {
    // this constructor is called at the finest mg level.
    // this means that it might need an averaging operator
    // but it does not need an interpolation operator
    m_refToCoar = 2;
    if (m_baCoarsened.isClosed())
    {
      m_baCoarsened = DisjointBoxLayout();
    }
    coarsen(m_baCoarsened, m_ba, m_refToCoar);
#if FIXIT      
    m_crseCorr.define(m_baCoarsened,a_ncomp,IntVect::TheUnitVector());
    m_crseResid.define(m_baCoarsened,a_ncomp,IntVect::TheZeroVector());
    m_averageOp.define(m_ba, m_baCoarsened, a_ncomp, m_refToCoar);
#endif

    m_lCoarsePtr = new GenLevelMG<T>(*this, m_refToCoar, a_opin);
    if (m_lCoarsePtr == NULL)
    {
      MayDay::Error("Out of Memory in GenLevelMG::define");
    }
  }
}

// --------------------------------------------------------------
// constructor for coarsened version of levelMG
// refcoarse is the ratio from coarsened thing to this thing
// --------------------------------------------------------------
template <class T>
void GenLevelMG<T>::define(const GenLevelMG&            a_L,
                           int                          a_refCoarse,
                           const GenLevelMGOp<T>* const a_opin)
{
  assert(a_refCoarse == 2);
  assert(a_opin != NULL);

  clearMemory();

  m_isDefined = true;

  m_ba = a_L.m_baCoarsened;

  Real dxcoarse = a_refCoarse*a_L.m_dxLevel;
  int ncomp = a_L.m_resid.nComp();

  ProblemDomain dprobc = coarsen(a_L.m_domain, a_refCoarse);
  m_nCoarserLevels = a_L.m_nCoarserLevels - 1;

  assert(m_nCoarserLevels >= 0);

  m_baseBaPtr = a_L.m_baseBaPtr;
  m_dxLevel = dxcoarse;
  m_domain = dprobc;
  m_resid.define(m_ba, ncomp, IntVect::TheZeroVector());
  m_levelopPtr = a_L.m_levelopPtr->newOp();

#if FIXIT
  m_levelopPtr->define(a_L.m_levelopPtr, a_refCoarse);
#endif

  // recursively define MG hierarchy
  if (m_nCoarserLevels != 0)
  {
    m_refToCoar = 2;
    coarsen(m_baCoarsened, m_ba, m_refToCoar);
#if FIXIT
    m_crseCorr.define(m_baCoarsened,ncomp,IntVect::TheUnitVector());
    m_crseResid.define(m_baCoarsened,ncomp,IntVect::TheZeroVector());
    m_averageOp.define(m_ba, m_baCoarsened, ncomp, m_refToCoar);
#endif

    m_lCoarsePtr = new GenLevelMG(*this, m_refToCoar, a_opin);
    if (m_lCoarsePtr == NULL)
    {
      MayDay::Error("Out of Memory in GenLevelMG::define");
    }
  }
  else
  {
    m_lCoarsePtr = NULL;
  }
}

template <class T>
void GenLevelMG<T>::mgRelax(T&       a_soln,
                            const T& a_rhs,
                            bool     a_bottomSolveFlag)
{
  assert(isDefined());

  int ncomp = a_rhs.nComp();

  assert(a_soln.nComp() == ncomp);
  assert(ncomp == m_resid.nComp());
  assert(a_soln.ghostVect() >= IntVect::TheUnitVector());

  // now, either coarsen or call bottom solver
  if (m_lCoarsePtr != NULL)
  {
    // first, smooth on this level
    for (int iter = 0; iter < m_numSmoothDown; iter++)
    {
      m_levelopPtr->smooth(a_soln,a_rhs);
    }

    GenLevelMGOp<T>& crseLevelop = *(m_lCoarsePtr->m_levelopPtr);
    crseLevelop.setToZero(m_crseCorr);

    // define residual on this grid
    // using homogeneous cf bcs
    m_levelopPtr->applyOpH(a_soln,m_resid);
    m_levelopPtr->diff(m_resid,a_rhs);
    m_levelopPtr->negate(m_resid);

    // coarsen residual and physical boundary condition
#if FIXIT
    assert(m_averageOp.isDefined());

    m_averageOp.averageToCoarse(m_crseResid, m_resid);
#endif

    // relax recursively on coarsened level
    m_lCoarsePtr->mgRelax(m_crseCorr, m_crseResid, a_bottomSolveFlag);

    // interpolate correction back up and modify solution
    crseCorrect(a_soln, m_crseCorr, m_refToCoar);

    // smooth again on this level
    for (int iter = 0; iter < m_numSmoothUp; iter++)
    {
      m_levelopPtr->smooth(a_soln,a_rhs);
    }

  }
  else if (m_domain.domainBox().numPts() == 1)
  {
    // grid is 1x1, just need to smooth
    // first, smooth on this level
    for (int iter = 0; iter < m_numBottomGSRB; iter++)
    {
      m_levelopPtr->smooth(a_soln,a_rhs);
    }
  }
  else if (a_bottomSolveFlag)
  {
    // first, smooth on this level
    for (int iter = 0; iter < m_numBottomGSRB; iter++)
    {
      m_levelopPtr->smooth(a_soln,a_rhs);
    }
    // use bottom smoother
    m_levelopPtr->bottomSmoother(a_soln,a_rhs);
  }
  else
  {
    // smooth on this level
    // do numSmoothDown iterations
    for (int iter = 0; iter < m_numSmoothDown; iter++)
    {
      m_levelopPtr->smooth(a_soln,a_rhs);
    }
  }
}

template <class T>
void GenLevelMG<T>::setnumBottomGSRB(int a_numBottomGSRB)
{
  m_numBottomGSRB = a_numBottomGSRB;
}

template <class T>
void GenLevelMG<T>::setnumSmoothUp(int a_numSmoothUp)
{
  m_numSmoothUp = a_numSmoothUp;
}

template <class T>
void GenLevelMG<T>::setnumSmoothDown(int a_numSmoothDown)
{
  m_numSmoothDown = a_numSmoothDown;
}

template <class T>
GenLevelMGOp<T>* GenLevelMG<T>::levelOpPtr()
{
  return m_levelopPtr;
}

template <class T>
GenLevelMG<T>* GenLevelMG<T>::lCoarsePtr()
{
  return m_lCoarsePtr;
}

template <class T>
void GenLevelMG<T>::setDefaultValues()
{
  // set defaults
  m_numBottomGSRB = 16;
  m_numSmoothUp = 4;
  m_numSmoothDown = 4;

  m_lCoarsePtr = NULL;
  m_levelopPtr = NULL;

  m_dxLevel = -1;

  m_isDefined = false;
}

template <class T>
void GenLevelMG<T>::clearMemory()
{
  if (m_lCoarsePtr != NULL)
  {
    delete m_lCoarsePtr;
    m_lCoarsePtr = NULL;
  }

  if (m_levelopPtr != NULL)
  {
    delete m_levelopPtr;
    m_levelopPtr = NULL;
  }

  m_isDefined = false;
}

template <class T>
void GenLevelMG<T>::crseCorrect(T&       a_fine,
                                const T& a_crse,
                                int      a_refRat)
{
#if FIXIT
  DataIterator crseIt = a_crse.dataIterator();
  const DisjointBoxLayout& crseBoxes = a_crse.getBoxes();
  const DisjointBoxLayout& fineBoxes = a_fine.getBoxes();

  for (crseIt.reset(); crseIt.ok(); ++crseIt)
  {
    const Box& crseBox = crseBoxes[crseIt()];

    const T& crseFab = a_crse[crseIt()];
    T& fineFab = a_fine[crseIt()];

    Box fineBox = fineBoxes[crseIt()];

    fineBox.coarsen(a_refRat);
    fineBox &= crseBox;

    // c set refinement box
    Box nrefbox(IntVect::TheZeroVector(),
                (a_refRat-1)*IntVect::TheUnitVector());

    if (!fineBox.isEmpty())
    {
      FORT_GENINTERPMG(CHF_FRA(fineFab),
                       CHF_FRA(crseFab),
                       CHF_BOX(fineBox),
                       CHF_CONST_INT(a_refRat),
                       CHF_BOX(nrefbox));
     }
  }
#endif
}

#endif

---