lib/src/AMRTools/QuadCFInterp.cpp
/***********************/
//get extended phi (lives next to interpivs)
/***********************/
get extended phi (lives next to interpivs)
get extended phi ( lives next to interpivs )
void
QuadCFInterp::getPhiStar(BaseFab<Real> & a_phistar,
                         const BaseFab<Real> & a_phic,
                         const QuadCFStencil& a_qcfs,
                         const Side::LoHiSide a_hiorlo,
                         const int a_idir,
                         const Interval& a_variables) const
{
  CH_TIMERS("QuadCFInterp::getPhiStar");
  CH_TIMER("QuadCFInterp::slopes", tslopes);
  CH_TIMER("QuadCFInterp::notPacked", tnp);
  CH_TIMER("QuadCFInterp::preamble", tpreamble);
  CH_assert(isDefined());
  CH_assert(a_qcfs.isDefined());
#if (CH_SPACEDIM > 1)
  Real dxf = m_dxFine;
  Real dxc = m_refRatio*dxf;
#endif
#if (CH_SPACEDIM == 3)
  int you1, you2;
  if (a_idir == 0)
    {
      you1 = 1;
      you2 = 2;
    }
  else if (a_idir == 1)
    {
      you1 = 0;
      you2 = 2;
    }
  else
    {
      you1 = 0;
      you2 = 1;
    }
#else 
  int you1;
  if (a_idir == 0)
    {
      you1 = 1;
    }
  else
    {
      you1 = 0;
    }
#endif
  if (!a_qcfs.isEmpty())
    {
      CH_START(tpreamble);
      CH_assert(m_level > 0);
      const IntVectSet& interp_ivs = a_qcfs.getFineIVS();
      const IntVectSet& coarsl_ivs = a_qcfs.getCoarIVS();
      if (!coarsl_ivs.isDense())
        {
          MayDay::Error("What the hell?? TreeIntVectSet ???");
        }
      if (!interp_ivs.isDense())
        {
          MayDay::Error("What the hell?? TreeIntVectSet ???");
        }
      Box coarinterpbox = coarsl_ivs.minBox();
      int ncomp = a_phic.nComp();
      CH_assert(ncomp == m_nComp);
      CH_assert(a_phic.box().contains((coarinterpbox)));
      int ihilo = sign(a_hiorlo);
      Box phistarbox = interp_ivs.minBox();
      phistarbox.shift(a_idir, ihilo);
      a_phistar.define(phistarbox, ncomp);
      CH_STOP(tpreamble);
      for (int ivar = a_variables.begin(); ivar <= a_variables.end(); ivar++)
        {
          CH_START(tslopes);
          BaseFab<Real>  coarslope(coarinterpbox, SpaceDim);
          BaseFab<Real>  coarcurva(coarinterpbox, SpaceDim);
#if (CH_SPACEDIM == 3)
          BaseFab<Real>  coarmixed(coarinterpbox, 1);
#endif
          IVSIterator coar_ivsit(coarsl_ivs);
          for (coar_ivsit.begin(); coar_ivsit.ok(); ++coar_ivsit)
            {
#if (CH_SPACEDIM > 1)
              const IntVect& coariv = coar_ivsit();
              coarslope(coariv, you1) =
                a_qcfs.computeFirstDerivative (a_phic, you1, ivar, coariv, dxc);
              coarcurva(coariv, you1) =
                a_qcfs.computeSecondDerivative(a_phic, you1, ivar, coariv, dxc);
#endif
#if (CH_SPACEDIM == 3)
              coarslope(coariv, you2) =
                a_qcfs.computeFirstDerivative (a_phic, you2, ivar, coariv, dxc);
              coarcurva(coariv, you2) =
                a_qcfs.computeSecondDerivative(a_phic, you2, ivar, coariv, dxc);
              coarmixed(coariv) =
                a_qcfs.computeMixedDerivative(a_phic, ivar, coariv, dxc);
#endif
            } 
          CH_STOP(tslopes);
          if (a_qcfs.finePacked() && CH_SPACEDIM==3)
            {
              const IntVect& iv = phistarbox.smallEnd();
              IntVect civ(iv);
              civ.coarsen(m_refRatio);
              Box region = a_qcfs.packedBox();
#if (CH_SPACEDIM == 3)
              FORT_PHISTAR(CHF_FRA_SHIFT(a_phistar, iv),
                           CHF_BOX_SHIFT(region, iv),
                           CHF_CONST_FRA_SHIFT(a_phic, civ),
                           CHF_FRA_SHIFT(coarslope, civ),
                           CHF_FRA_SHIFT(coarcurva, civ),
                           CHF_FRA_SHIFT(coarmixed, civ),
                           CHF_CONST_REAL(dxf),
                           CHF_CONST_INT(ivar),
                           CHF_CONST_INT(a_idir),
                           CHF_CONST_INT(ihilo),
                           CHF_CONST_INT(m_refRatio));
#endif
            }
          else
            {
              CH_START(tnp);
              IntVect ivf, ivc, ivstar;
#if (CH_SPACEDIM > 1)
              int jf, jc;
              Real xf, xc, x1;
#endif
              Real pc, update1=0, update2=0, update3=0;
              IVSIterator fine_ivsit(interp_ivs);
              for (fine_ivsit.begin(); fine_ivsit.ok(); ++fine_ivsit)
                {
                  ivf = fine_ivsit();
                  ivc = coarsen(ivf, m_refRatio);
                  ivstar = ivf;
                  ivstar.shift(a_idir, ihilo);
                  pc = a_phic(ivc,ivar);
#if (CH_SPACEDIM > 1)
                  jf = ivf[you1];
                  jc = ivc[you1];
                  xf = (jf+0.5)*dxf;
                  xc = (jc+0.5)*dxc;
                  x1 = xf-xc;
                  update1= x1*coarslope(ivc, you1) + 0.5*x1*x1*coarcurva(ivc, you1);
#endif
#if (CH_SPACEDIM==3)
                  Real x2;
                  jf = ivf[you2];
                  jc = ivc[you2];
                  xf = (jf+0.5)*dxf;
                  xc = (jc+0.5)*dxc;
                  x2 = xf-xc;
                  update2 =  x2*coarslope(ivc, you2) + 0.5*x2*x2*coarcurva(ivc, you2);
                  update3 =  x1*x2*coarmixed(ivc);
#endif
                  a_phistar(ivstar, ivar) = pc+update1+update2+update3;
                } 
              CH_STOP(tnp);
            } 
        }
    } 
}
void QuadCFInterp :: getPhiStar ( BaseFab < Real > & a_phistar , const BaseFab < Real > & a_phic , const QuadCFStencil & a_qcfs , const Side :: LoHiSide a_hiorlo , const int a_idir , const Interval & a_variables ) const { CH_TIMERS ( " " ) ; CH_TIMER ( " " , tslopes ) ; CH_TIMER ( " " , tnp ) ; CH_TIMER ( " " , tpreamble ) ; CH_assert ( isDefined ( ) ) ; CH_assert ( a_qcfs . isDefined ( ) ) ; #if ( CH_SPACEDIM > 1 ) Real dxf = m_dxFine ; Real dxc = m_refRatio * dxf ; #endif #if ( CH_SPACEDIM == 3 ) int you1 , you2 ; if ( a_idir == 0 ) { you1 = 1 ; you2 = 2 ; } else if ( a_idir == 1 ) { you1 = 0 ; you2 = 2 ; } else { you1 = 0 ; you2 = 1 ; } #else int you1 ; if ( a_idir == 0 ) { you1 = 1 ; } else { you1 = 0 ; } #endif if ( ! a_qcfs . isEmpty ( ) ) { CH_START ( tpreamble ) ; CH_assert ( m_level > 0 ) ; const IntVectSet & interp_ivs = a_qcfs . getFineIVS ( ) ; const IntVectSet & coarsl_ivs = a_qcfs . getCoarIVS ( ) ; if ( ! coarsl_ivs . isDense ( ) ) { MayDay :: Error ( " " ) ; } if ( ! interp_ivs . isDense ( ) ) { MayDay :: Error ( " " ) ; } Box coarinterpbox = coarsl_ivs . minBox ( ) ; int ncomp = a_phic . nComp ( ) ; CH_assert ( ncomp == m_nComp ) ; CH_assert ( a_phic . box ( ) . contains ( ( coarinterpbox ) ) ) ; int ihilo = sign ( a_hiorlo ) ; Box phistarbox = interp_ivs . minBox ( ) ; phistarbox . shift ( a_idir , ihilo ) ; a_phistar . define ( phistarbox , ncomp ) ; CH_STOP ( tpreamble ) ; for ( int ivar = a_variables . begin ( ) ; ivar <= a_variables . end ( ) ; ivar ++ ) { CH_START ( tslopes ) ; BaseFab < Real > coarslope ( coarinterpbox , SpaceDim ) ; BaseFab < Real > coarcurva ( coarinterpbox , SpaceDim ) ; #if ( CH_SPACEDIM == 3 ) BaseFab < Real > coarmixed ( coarinterpbox , 1 ) ; #endif IVSIterator coar_ivsit ( coarsl_ivs ) ; for ( coar_ivsit . begin ( ) ; coar_ivsit . ok ( ) ; ++ coar_ivsit ) { #if ( CH_SPACEDIM > 1 ) const IntVect & coariv = coar_ivsit ( ) ; coarslope ( coariv , you1 ) = a_qcfs . computeFirstDerivative ( a_phic , you1 , ivar , coariv , dxc ) ; coarcurva ( coariv , you1 ) = a_qcfs . computeSecondDerivative ( a_phic , you1 , ivar , coariv , dxc ) ; #endif #if ( CH_SPACEDIM == 3 ) coarslope ( coariv , you2 ) = a_qcfs . computeFirstDerivative ( a_phic , you2 , ivar , coariv , dxc ) ; coarcurva ( coariv , you2 ) = a_qcfs . computeSecondDerivative ( a_phic , you2 , ivar , coariv , dxc ) ; coarmixed ( coariv ) = a_qcfs . computeMixedDerivative ( a_phic , ivar , coariv , dxc ) ; #endif } CH_STOP ( tslopes ) ; if ( a_qcfs . finePacked ( ) && CH_SPACEDIM == 3 ) { const IntVect & iv = phistarbox . smallEnd ( ) ; IntVect civ ( iv ) ; civ . coarsen ( m_refRatio ) ; Box region = a_qcfs . packedBox ( ) ; #if ( CH_SPACEDIM == 3 ) FORT_PHISTAR ( CHF_FRA_SHIFT ( a_phistar , iv ) , CHF_BOX_SHIFT ( region , iv ) , CHF_CONST_FRA_SHIFT ( a_phic , civ ) , CHF_FRA_SHIFT ( coarslope , civ ) , CHF_FRA_SHIFT ( coarcurva , civ ) , CHF_FRA_SHIFT ( coarmixed , civ ) , CHF_CONST_REAL ( dxf ) , CHF_CONST_INT ( ivar ) , CHF_CONST_INT ( a_idir ) , CHF_CONST_INT ( ihilo ) , CHF_CONST_INT ( m_refRatio ) ) ; #endif } else { CH_START ( tnp ) ; IntVect ivf , ivc , ivstar ; #if ( CH_SPACEDIM > 1 ) int jf , jc ; Real xf , xc , x1 ; #endif Real pc , update1 = 0 , update2 = 0 , update3 = 0 ; IVSIterator fine_ivsit ( interp_ivs ) ; for ( fine_ivsit . begin ( ) ; fine_ivsit . ok ( ) ; ++ fine_ivsit ) { ivf = fine_ivsit ( ) ; ivc = coarsen ( ivf , m_refRatio ) ; ivstar = ivf ; ivstar . shift ( a_idir , ihilo ) ; pc = a_phic ( ivc , ivar ) ; #if ( CH_SPACEDIM > 1 ) jf = ivf [ you1 ] ; jc = ivc [ you1 ] ; xf = ( jf + 0.5 ) * dxf ; xc = ( jc + 0.5 ) * dxc ; x1 = xf - xc ; update1 = x1 * coarslope ( ivc , you1 ) + 0.5 * x1 * x1 * coarcurva ( ivc , you1 ) ; #endif #if ( CH_SPACEDIM == 3 ) Real x2 ; jf = ivf [ you2 ] ; jc = ivc [ you2 ] ; xf = ( jf + 0.5 ) * dxf ; xc = ( jc + 0.5 ) * dxc ; x2 = xf - xc ; update2 = x2 * coarslope ( ivc , you2 ) + 0.5 * x2 * x2 * coarcurva ( ivc , you2 ) ; update3 = x1 * x2 * coarmixed ( ivc ) ; #endif a_phistar ( ivstar , ivar ) = pc + update1 + update2 + update3 ; } CH_STOP ( tnp ) ; } } } }
get extended phi (lives next to interpivs)
get extended phi ( lives next to interpivs )
//CH_TIMER("QuadCFInterp::computeFirstDerivative",  t1st); //CH_TIMER("QuadCFInterp::computesecondDerivative", t2nd); //CH_TIMER("QuadCFInterp::computemixedDerivative",  tmixed); // if we think of a_idir as the "me" direction, then // the other directions can be "you1" and "you2" // (CH_SPACEDIM == 2) //if cfsten is empty, nothing to interpolate. // allocate phistar here //phi = phino + slope*x + half*x*x*curvature // coarslope.setVal(0.); //first find extended value phistar. get slopes of coarse solution // this isn't relevant for 1D // coarslope(coariv, a_idir) = 0.0; // coarcurva(coariv, a_idir) = 0.0; //end loop over coarse intvects // ifdef is here to prevent unused variable wasrnings in 1D // for 1D, none of this is necessary -- just copy // coarse value into phiStar //add in mixed derivative component //end loop over fine intvects // end if for not packed optimization //end loop over variables //end if (level>0 && !hocfs.isempty())