---

GenAMRLevelMGImplem.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 _GENAMRLEVELMGIMPLEM_H
#define _GENAMRLEVELMGIMPLEM_H

#include "SPMD.H"
#include "LayoutIterator.H"

#include "GenAMRSolver.H"

template <class T>
GenAMRLevelMG<T>::GenAMRLevelMG()
{
  setDefaultValues();
}

template <class T>
GenAMRLevelMG<T>::GenAMRLevelMG(const GenAMRSolver<T>* const    a_parent,
                                int                             a_level,
                                const GenAMRLevelMGOp<T>* const a_opin,
                                int                             a_ncomp)
{
  setDefaultValues();

  define(a_parent, a_level, a_opin, a_ncomp);
}

template <class T>
GenAMRLevelMG<T>::~GenAMRLevelMG()
{
  clearMemory();
}

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

template <class T>
void GenAMRLevelMG<T>::define(const GenAMRSolver<T>* const    a_parent,
                              int                             a_level,
                              const GenAMRLevelMGOp<T>* const a_opin,
                              int                             a_ncomp)
{
  clearMemory();

  assert(a_parent != NULL);
  assert(a_parent->isDefined());
  assert(a_level >= 0);
  assert(a_level <= a_parent->m_finestLevel);

  m_isDefined = true;
  m_parent = a_parent;
  m_level = a_level;

  int nCoarserLevels;
  const DisjointBoxLayout& grids = a_parent->m_gridsLevel[a_level];
  const DisjointBoxLayout* baseBAPtr = NULL;

  // put reasonable value into refinement ratio
  int reftoCoarse = 2;

  const ProblemDomain& domain = a_parent->m_domainLevel[a_level];

  if (m_level > 0)
  {
    const DisjointBoxLayout& baseGrids = a_parent->m_gridsLevel[a_level-1];

    baseBAPtr = &baseGrids;
    reftoCoarse = a_parent->m_refRatio[a_level-1];
    m_levfluxreg.define(grids, baseGrids, domain, reftoCoarse, a_ncomp);

    // nCoarserLevels = log2(reftoCoarse) - 1
    // bug fixed by petermc, 15 Apr 2002
    nCoarserLevels = -1;

    for (int interRatio = reftoCoarse; interRatio >= 2; interRatio /= 2)
    {
      nCoarserLevels++;
    }
  }
  else
  {
    nCoarserLevels = 0;
  }

  Real dx = a_parent->m_dxLevel[a_level];

  m_levelMG.define(grids, baseBAPtr, dx,
                   reftoCoarse, domain,
                   nCoarserLevels, a_opin,
                   a_ncomp);

  m_mginterp.define(grids, a_ncomp, reftoCoarse, domain);

  if (m_level > 0)
  {
    coarsen(m_coarsenedGrids, grids, reftoCoarse);
#if FIXIT
    m_resC.define(m_coarsenedGrids, a_ncomp,IntVect::TheZeroVector());
    m_residualCopier.define(m_coarsenedGrids, a_parent->m_gridsLevel[a_level-1]);
    m_averageOp.define(grids, m_coarsenedGrids, a_ncomp, reftoCoarse);
#endif
  }
  if (m_level < a_parent->m_finestLevel)
  {
    m_fineExchangeCopier.define( a_parent->m_gridsLevel[a_level+1],
                                 a_parent->m_gridsLevel[a_level+1],
                                 IntVect::TheUnitVector());
  }

  m_levelopPtr = a_opin->newOp();
#if FIXIT
  m_levelopPtr->define(grids, baseBAPtr, dx, reftoCoarse,
                       domain, false, a_ncomp);
#endif

  // residual,lofphi has no ghost cells, corr and phi have 1 ghost cell
  m_levelopPtr->defineData(m_lofPhi);
  m_levelopPtr->defineData(m_resid);

  m_levelopPtr->defineOpData(m_corr);
  m_levelopPtr->defineOpData(m_dcorr);
  m_levelopPtr->defineOpData(m_phiSave);

  // ParmParse pp("amrlevelmg");
  // m_arrayViewVerbose = pp.contains("verbose");
}

template <class T>
void GenAMRLevelMG<T>::setnumSmoothUp(int a_numSmoothUp)
{
  assert(isDefined());

  m_levelMG.setnumSmoothUp(a_numSmoothUp);
}

template <class T>
void GenAMRLevelMG<T>::setnumSmoothDown(int a_numSmoothDown)
{
  assert(isDefined());

  m_levelMG.setnumSmoothDown(a_numSmoothDown);
}

template <class T>
void GenAMRLevelMG<T>::applyAMROperator(Vector<T* >& a_phiLevel,
                                        T&           a_LOfPhi)
{
  assert(isDefined());

  T* crsePhi = NULL;

  // operator on this level including interpolated bc's from Level-1
  if (m_level > 0)
  {
    crsePhi = a_phiLevel[m_level-1];
  }

  m_levelopPtr->setToZero(a_LOfPhi);

  T& phi = *a_phiLevel[m_level];
  m_levelopPtr->applyOpI(phi, crsePhi, a_LOfPhi);

  // now reflux to enforce flux-matching from finer levels...
  if (m_level < m_parent->m_finestLevel)
  {
    reflux(a_phiLevel, a_LOfPhi);
  }
}

template <class T>
void GenAMRLevelMG<T>::applyAMROperatorHphys(Vector<T* >& a_phiLevel,
                                             T&           a_LOfPhi)
{
  assert(isDefined());

  m_levelopPtr->setToZero(a_LOfPhi);

  T* crsePhi = NULL;

  // operator on this level including interpolated bc's from Level-1
  if (m_level > 0)
  {
    crsePhi = a_phiLevel[m_level-1];
  }

  T& phi = *a_phiLevel[m_level];
  m_levelopPtr->applyOpIcfHphys(phi, crsePhi, a_LOfPhi);

  // now reflux to enforce flux-matching from finer levels...
  if (m_level < m_parent->m_finestLevel)
  {
    reflux(a_phiLevel, a_LOfPhi);
  }
}

template <class T>
void GenAMRLevelMG<T>::computeAMRResidual(Vector<T* >&       a_phiLevel,
                                          const Vector<T* >& a_rhsLevel)
{
  assert(isDefined());

  T& rhs = *a_rhsLevel[m_level];

  applyAMROperator(a_phiLevel,m_resid);
  m_levelopPtr->negDiff(m_resid,rhs);
}

template <class T>
void GenAMRLevelMG<T>::computeAMRResidual(T&                 a_resid,
                                          Vector<T* >&       a_phiLevel,
                                          const Vector<T* >& a_rhsLevel)
{
  assert(isDefined());

  T& rhs = *a_rhsLevel[m_level];

  applyAMROperator(a_phiLevel, a_resid);
  m_levelopPtr->diff(a_resid,rhs);
  m_levelopPtr->negate(a_resid);
}

template <class T>
void GenAMRLevelMG<T>::computeAMRResidualHphys(Vector<T* >&       a_phiLevel,
                                               const Vector<T* >& a_rhsLevel)
{
  assert(isDefined());

  T& rhs = *a_rhsLevel[m_level];

  applyAMROperatorHphys(a_phiLevel, m_resid);
  m_levelopPtr->diff(m_resid,rhs);
  m_levelopPtr->negate(m_resid);
}

template <class T>
void GenAMRLevelMG<T>::upSweep(Vector<T *>&       a_phiLevel,
                               const Vector<T *>& a_rhsLevel)
{
  assert(isDefined());

  T& phi = *a_phiLevel[m_level];
  const T& rhs = *a_rhsLevel[m_level];

  assert(phi.nComp() == rhs.nComp());

  // Interpolate corrections from next coarser level to the correction
  // field at this level.
  if (m_level > 0)
  {
    // Initialize correction at the next coarser level.
    // this uses the two-level only operator (Lnf)
    GenAMRLevelMG& nextCoarser = *m_parent->m_amrmgLevel[m_level-1];
    m_mginterp.interpToFine(m_corr, nextCoarser.m_corr);

    m_levelopPtr->applyOpIcfHphys(m_corr, &nextCoarser.m_corr, m_lofPhi);
  }

  // Finish resid calc,
  // initialize the correction to the correction,
  // and apply smoother.
  m_levelopPtr->diff(m_resid,m_lofPhi);

  m_levelopPtr->setToZero(m_dcorr);
  smooth(m_dcorr,m_resid);

  // Increment correction.
  // Increment saved value of phi with corr, and overwrite phi.
  m_levelopPtr->sum(m_corr,m_dcorr);
  m_levelopPtr->sum(m_phiSave,m_corr);

  m_levelopPtr->copy(phi,m_phiSave);
}

template <class T>
void GenAMRLevelMG<T>::downSweep(Vector<T* >&       a_phiLevel,
                                 const Vector<T* >& a_rhsLevel)
{
  assert(isDefined());

  T & phi = *a_phiLevel[m_level];
  T & rhs = *a_rhsLevel[m_level];

  assert(phi.nComp() == rhs.nComp());

  m_levelopPtr->setToZero(m_corr);
  m_levelopPtr->copy(m_phiSave,phi);

  // Apply smoother.
  smooth(m_corr,m_resid);

  // Update phi.
  m_levelopPtr->sum(phi,m_corr);

  // Compute residual for the next coarser level.
  // Form residual on next coarser level.
  if (m_level > 0)
  {
    // Initialize correction at the next coarser level.
    // this has to be done before applyOpI
    // int reftoCoarse = m_parent->m_refRatio[m_level-1];
    GenAMRLevelMG& nextCoarser = *m_parent->m_amrmgLevel[m_level-1];
    GenAMRLevelMGOp<T>& coarLevelop = *(nextCoarser.m_levelopPtr);

    coarLevelop.setToZero(nextCoarser.m_corr);

    // petermc, 18 Apr 2002, replaced this with applyOpH
    // because coarser level is zero anyway.
    // m_levelopPtr->applyOpIcfHphys(m_corr, &nextCoarser.m_corr, m_lofPhi);
    m_levelopPtr->applyOpH(m_corr,m_lofPhi);
    m_levelopPtr->diff(m_lofPhi,m_resid);
    m_levelopPtr->negate(m_lofPhi);

#if FIXIT
    m_averageOp.averageToCoarse(m_resC, m_lofPhi);
#endif

    nextCoarser.computeAMRResidualHphys(a_phiLevel, a_rhsLevel);

    // overwrite residual on Level-1 with
    // coarsened residual from this level
    m_resC.copyTo(m_resC.interval(),
                  nextCoarser.m_resid,
                  nextCoarser.m_resid.interval(), m_residualCopier);
  }
}

template <class T>
Vector<Real> GenAMRLevelMG<T>::computeResidualNorm(int a_normType) const
{
  assert(isDefined());

  Vector<Real> normPeterson(m_resid.nComp(),0.0);
  normPeterson = computeNorm(m_resid, a_normType);

  return normPeterson;
}

template <class T>
GenAMRLevelMGOp<T>* GenAMRLevelMG<T>::levelOpPtr() const
{
  return m_levelopPtr;
}

// returns normType norm of mfab
template <class T>
Vector<Real> GenAMRLevelMG<T>::computeNorm(const T& a_mfinput,
                                           int      a_normType) const
{
  assert(isDefined());
  assert((a_normType >= 0) && (a_normType <= 2));

  // create temps so inputs do not get changed
  T mfab;
  mfab.define(a_mfinput, a_mfinput.interval());

  // compute the bloody norm
  int ncomp = mfab.nComp();

  Vector<Real> normTot(ncomp,0.0);

  if (m_level < m_parent->m_finestLevel)
  {
#if FIXIT
    const DisjointBoxLayout& fineGrids = m_parent->m_gridsLevel[m_level+1];
    int nrefFine =  m_parent->m_refRatio[m_level];

    DataIterator dit = mfab.dataIterator();
    for (dit.ok(); dit.ok(); ++dit)
    {
      LayoutIterator litFine = fineGrids.layoutIterator();

      for (litFine.reset(); litFine.ok(); ++litFine)
      {
        Box overlayBox = mfab[dit()].box();
        Box coarsenedGrid = coarsen(fineGrids[litFine()], nrefFine);

        overlayBox &= coarsenedGrid;

        if (!overlayBox.isEmpty())
        {
          mfab[dit()].setVal(0.0,overlayBox,0, ncomp);
        }
      }
    }
#endif
  }

  for (int comp = 0; comp < mfab.nComp(); comp++)
  {
#if FIXIT
    Interval thisInterval(comp,comp);
    // this function does the wacky boxlib-type norm
    normTot[comp] = norm(mfab, thisInterval, a_normType);
#endif

    // so we normalize it for sanity's sake
    if (a_normType != 0)
    {
      Real dx = m_parent->m_dxLevel[m_level];
      Real expon = SpaceDim/a_normType;
      Real normfac = pow(dx, expon);
      normTot[comp] *= normfac;
    }
  } // end loop over components

  return normTot;
}

template <class T>
void GenAMRLevelMG<T>::smooth(T&       a_soln,
                              const T& a_rhs)
{
  assert(isDefined());
  assert(a_soln.ghostVect() >= IntVect::TheUnitVector());
  assert(a_soln.getBoxes() == m_parent->m_gridsLevel[m_level]);
  assert(a_rhs.getBoxes() == m_parent->m_gridsLevel[m_level]);

  // no bottom solve here
  m_levelMG.mgRelax(a_soln, a_rhs, false);
}

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

template <class T>
void GenAMRLevelMG<T>::setDefaultValues()
{
  m_isDefined = false;
  m_parent = NULL;
  m_level = -1;
  m_isDefined = false;
  m_arrayViewVerbose = false;
  m_levelopPtr = NULL;
}

template <class T>
void GenAMRLevelMG<T>::reflux(Vector<T *>& a_phiLevel,
                              T&           a_LOfPhi)
{
  assert(isDefined());
  assert(m_level < m_parent->m_finestLevel);

  GenAMRLevelMG& higher = *(m_parent->m_amrmgLevel[m_level + 1]);

  assert(higher.isDefined());

#if FIXIT
  LevelFluxRegister& finelevfr = higher.m_levfluxreg;
#endif

  // finelevfr.setToZero(); called in initFRCoarse
  initFRCoarse(a_phiLevel);

  incrementFRFine(a_phiLevel);
  Real dxlevel = m_parent->m_dxLevel[m_level];

  assert(dxlevel > 0);

#if FIXIT
  Real scale = 1.0/dxlevel;
  finelevfr.reflux(a_LOfPhi, scale);
#endif
}


template <class T>
void GenAMRLevelMG<T>::initFRCoarse(Vector<T *>& a_phiLevel)
{
  assert(isDefined());
  assert(m_level < m_parent->m_finestLevel);

  GenAMRLevelMG& higher = *(m_parent->m_amrmgLevel[m_level + 1]);

  assert(higher.isDefined());

  LevelFluxRegister& finelevfr = higher.m_levfluxreg;

  // looping through every box in this grid
  // and fill fluxregisters with coarse flux
  finelevfr.setToZero();

#if FIXIT
  Interval interv(0,a_phiLevel[m_level]->nComp()-1);
  T& phic = *a_phiLevel[m_level];

  DataIterator dit = phic.dataIterator();
  for (dit.reset(); dit.ok(); ++dit)
  {
    const T& coarfab = phic[dit()];

    for (int idir = 0; idir < SpaceDim; idir++)
    {
      T coarflux;
      m_levelopPtr->getFlux(coarflux, coarfab, dit(), idir);

      Real scale = 1.0;
      finelevfr.incrementCoarse(coarflux,
                                scale,dit(),
                                interv,interv,idir);
    }
  }
#endif
}

template <class T>
void GenAMRLevelMG<T>::incrementFRFine(Vector<T *>& a_phiLevel)
{
  assert(isDefined());
  assert(m_level < m_parent->m_finestLevel);

  GenAMRLevelMG& higher = *(m_parent->m_amrmgLevel[m_level+1]);

  assert(higher.isDefined());

#if FIXIT
  LevelFluxRegister& finelevfr = higher.m_levfluxreg;
#endif

  // flux register stuff
  // looping through every box in finer grid
  // and increment fluxregisters with fine flux
#if FIXIT
  T& phif = *a_phiLevel[m_level + 1];
  T& phic = *a_phiLevel[m_level];
  Interval interv = phif.interval();

  higher.m_levelopPtr->CFInterp(phif, phic);
  phif.exchange(phif.interval(),  m_fineExchangeCopier);

  DataIterator ditf = phif.dataIterator();
  for (ditf.reset(); ditf.ok(); ++ditf)
  {
    const T& phifFab = phif[ditf()];

    for (int idir = 0; idir < SpaceDim; idir++)
    {
      SideIterator sit;
      for (sit.begin(); sit.ok(); sit.next())
      {
        Side::LoHiSide hiorlo = sit();
        Box fabbox;

        if (sit() == Side::Lo)
        {
          fabbox = adjCellLo(phifFab.box(),idir, 2);
          fabbox.shift(idir, 2);
        }
        else
        {
          fabbox = adjCellHi(phifFab.box(),idir, 2);
          fabbox.shift(idir, -2);
        }

        assert(!fabbox.isEmpty());

        T phifab(fabbox, phif.nComp());
        phifab.copy(phifFab);

        T fineflux;
        higher.m_levelopPtr->getFlux(fineflux, phifab, ditf(), idir);

        Real scale = 1.0;
        finelevfr.incrementFine(fineflux, scale, ditf(),
                                interv, interv, idir, hiorlo);
      }
    }
  }
#endif
}

template <class T>
void GenAMRLevelMG<T>::setArrayViewVerbose(bool a_verbosity)
{
  m_arrayViewVerbose = a_verbosity;
}

#endif

---