自学教程：C++ CH_assert函数代码示例

// Generic function to reduce the minloc,maxloc,avg of a Real over all procs onto rank0int reduce_avg_min_max_loc(Real value, Real& avg, Real& min, Real& max, int& minloc, int& maxloc){  int eek=0;#ifdef CH_MPI  struct  {    double val;    int rank;  } in, out;  in.val = value;  in.rank = procID();  Real sum;  eek = MPI_Reduce(&value, &sum, 1, MPI_CH_REAL, MPI_SUM, uniqueProc(SerialTask::compute), Chombo_MPI::comm);  CH_assert(eek == MPI_SUCCESS);  avg=sum/(Real)numProc();  eek = MPI_Reduce(&in, &out, 1, MPI_DOUBLE_INT, MPI_MINLOC, uniqueProc(SerialTask::compute), Chombo_MPI::comm);  CH_assert(eek == MPI_SUCCESS);  min = out.val;  minloc = out.rank;  eek = MPI_Reduce(&in, &out, 1, MPI_DOUBLE_INT, MPI_MAXLOC, uniqueProc(SerialTask::compute), Chombo_MPI::comm);  CH_assert(eek == MPI_SUCCESS);  max = out.val;  maxloc = out.rank;#else  avg=value;  min=value;  max=value;  minloc=0;  maxloc=0;#endif  return eek;}

voidEBCellFAB::setInvalidData(const Real& a_val,                          const int& a_comp){  CH_assert(a_comp >= 0);  CH_assert(a_comp < m_nComp);  if (m_ebisBox.isAllRegular())    {      return;    }  else if (m_ebisBox.isAllCovered())    {      m_regFAB.setVal(a_val, a_comp);    }  else    {      for (BoxIterator bit(m_region); bit.ok(); ++bit)        {          const IntVect& iv = bit();          if (m_ebisBox.isCovered(iv))            {              m_regFAB(iv, a_comp) = a_val;            }        }      //also set the multivalued cells      for (IVSIterator ivsit(getMultiCells());          ivsit.ok(); ++ivsit)        {          const IntVect& iv = ivsit();          m_regFAB(iv, a_comp) = a_val;        }      //also set the multivalued cells    }}

// this preconditioner first initializes phihat to (IA)phihat = rhshat// (diagonization of L -- A is the matrix version of L)// then smooths with a couple of passes of levelGSRBvoid VCAMRPoissonOp2::preCond(LevelData<FArrayBox>&       a_phi,                              const LevelData<FArrayBox>& a_rhs){  CH_TIME("VCAMRPoissonOp2::preCond");  // diagonal term of this operator in:  //  //       alpha * a(i)  //     + beta  * sum_over_dir (b(i-1/2*e_dir) + b(i+1/2*e_dir)) / (dx*dx)  //  // The inverse of this is our initial multiplier.  int ncomp = a_phi.nComp();  CH_assert(m_lambda.isDefined());  CH_assert(a_rhs.nComp()    == ncomp);  CH_assert(m_bCoef->nComp() == ncomp);  // Recompute the relaxation coefficient if needed.  resetLambda();  // don't need to use a Copier -- plain copy will do  DataIterator dit = a_phi.dataIterator();  for (dit.begin(); dit.ok(); ++dit)    {      // also need to average and sum face-centered bCoefs to cell-centers      Box gridBox = a_rhs[dit].box();      // approximate inverse      a_phi[dit].copy(a_rhs[dit]);      a_phi[dit].mult(m_lambda[dit], gridBox, 0, 0, ncomp);    }  relax(a_phi, a_rhs, 2);}

void DenseIntVectSet::coarsen(int iref){  if (iref == 1) return;  CH_assert(iref >= 1);  // int refinements = iref/2;  CH_assert((iref/2)*2 == iref); // check iref for power of 2  Box newDomain(m_domain);  newDomain.coarsen(iref);  DenseIntVectSet newSet(newDomain, false);  BoxIterator bit(m_domain);  int count=0;  for (bit.begin(); bit.ok(); ++bit, ++count)    {      if (m_bits[count])        {          IntVect iv(bit());          iv.coarsen(iref);          long index = newDomain.index(iv);          newSet.m_bits.setTrue(index);        }    }  *this = newSet;}

// Set boundary slopes://   The boundary slopes in a_dW are already set to one sided difference//   approximations.  If this function doesn't change them they will be//   used for the slopes at the boundaries.void ExplosionIBC::setBdrySlopes(FArrayBox&       a_dW,                                 const FArrayBox& a_W,                                 const int&       a_dir,                                 const Real&      a_time){  CH_assert(m_isFortranCommonSet == true);  CH_assert(m_isDefined == true);  // In periodic case, this doesn't do anything  if (!m_domain.isPeriodic(a_dir))  {    Box loBox,hiBox,centerBox,domain;    int hasLo,hasHi;    Box slopeBox = a_dW.box();    slopeBox.grow(a_dir,1);    // Generate the domain boundary boxes, loBox and hiBox, if there are    // domain boundarys there    loHiCenter(loBox,hasLo,hiBox,hasHi,centerBox,domain,               slopeBox,m_domain,a_dir);    // Set the boundary slopes if necessary    if ((hasLo != 0) || (hasHi != 0))    {      FORT_SLOPEBCSF(CHF_FRA(a_dW),                     CHF_CONST_FRA(a_W),                     CHF_CONST_INT(a_dir),                     CHF_BOX(loBox),                     CHF_CONST_INT(hasLo),                     CHF_BOX(hiBox),                     CHF_CONST_INT(hasHi));    }  }}

void EBPoissonOp::relax(LevelData<EBCellFAB>&       a_e,      const LevelData<EBCellFAB>& a_residual,      int                         a_iterations){  CH_TIME("EBPoissonOp::relax");  CH_assert(a_e.ghostVect() == m_ghostCellsPhi);  CH_assert(a_residual.ghostVect() == m_ghostCellsRHS);  CH_assert(a_e.nComp() == 1);  CH_assert(a_residual.nComp() == 1);  if (m_relaxType == 0)    {      for (int i = 0; i < a_iterations; i++)        {          levelJacobi(a_e,a_residual);        }    }  else  if (m_relaxType == 1)    {      for (int i = 0; i < a_iterations; i++)        {          levelMulticolorGS(a_e,a_residual);        }    }  else    {      MayDay::Error("EBPoissonOp::relax - invalid relaxation type");    }}

void PatchGodunov::updateState(FArrayBox&       a_U,                               FluxBox&         a_F,                               Real&            a_maxWaveSpeed,                               const FArrayBox& a_S,                               const Real&      a_dt,                               const Box&       a_box){  CH_assert(isDefined());  CH_assert(a_box == m_currentBox);  int numPrim = m_gdnvPhysics->numPrimitives();  int numFlux = m_gdnvPhysics->numFluxes();  FluxBox whalf(a_box,numPrim);  whalf.setVal(0.0);  a_F.resize(a_box,numFlux);  a_F.setVal(0.0);  computeWHalf(whalf, a_U, a_S, a_dt, a_box);  FArrayBox dU(a_U.box(),a_U.nComp());  computeUpdate(dU, a_F, a_U, whalf, a_dt, a_box);  a_U += dU;  // Get and return the maximum wave speed on this patch/grid  a_maxWaveSpeed = m_gdnvPhysics->getMaxWaveSpeed(a_U, m_currentBox);}

// Set boundary fluxesvoid RampIBC::primBC(FArrayBox&            a_WGdnv,                     const FArrayBox&      a_Wextrap,                     const FArrayBox&      a_W,                     const int&            a_dir,                     const Side::LoHiSide& a_side,                     const Real&           a_time){  CH_assert(m_isFortranCommonSet == true);  CH_assert(m_isDefined == true);  // Neither the x or y direction can be periodic  if ((a_dir == 0 || a_dir == 1) && m_domain.isPeriodic(a_dir))    {      MayDay::Error("RampIBC::primBC: Neither the x or y boundaries can be periodic");    }  Box boundaryBox;  getBoundaryFaces(boundaryBox, a_WGdnv.box(), a_dir, a_side);  if (! boundaryBox.isEmpty() )    {      // Set the boundary fluxes      int lohisign = sign(a_side);      FORT_RAMPBCF(CHF_FRA(a_WGdnv),                   CHF_CONST_FRA(a_Wextrap),                   CHF_CONST_FRA(a_W),                   CHF_CONST_REAL(a_time),                   CHF_CONST_INT(lohisign),                   CHF_CONST_REAL(m_dx),                   CHF_CONST_INT(a_dir),                   CHF_BOX(boundaryBox));    }}

NewPoissonOp* NewPoissonOpFactory::MGnewOp(const ProblemDomain& a_FineindexSpace,                                           int   a_depth,                                           bool  a_homoOnly){  CH_assert(a_depth >= 0 );  CH_assert(m_bc != NULL);  NewPoissonOp* newOp = new NewPoissonOp();  RealVect dx = m_dx;  ProblemDomain domain = a_FineindexSpace;  for (int i=0; i<a_depth; i++)    {      Box d = domain.domainBox();      d.coarsen(8);      d.refine(8);      if (domain.domainBox() == d)        {          dx*=2;          domain.coarsen(2);        }      else        {          return NULL;        }    }  newOp->define(dx, domain, m_bc);  return newOp;}

// -----------------------------------------------------------------------------// Interpolate ghosts at CF interface using zeros on coarser grids.// This version acts on a set of IntVects.// -----------------------------------------------------------------------------void interpOnIVSHomo (LevelData<FArrayBox>& a_phif,                      const DataIndex&      a_index,                      const int             a_dir,                      const Side::LoHiSide  a_side,                      const IntVectSet&     a_interpIVS,                      const Real            a_fineDxDir,                      const Real            a_crseDxDir){    CH_TIME("interpOnIVSHomo");    // Sanity checks    CH_assert((a_dir >= 0) && (a_dir < SpaceDim));    CH_assert(a_phif.ghostVect()[a_dir] >= 1);    IVSIterator fine_ivsit(a_interpIVS);    FArrayBox& a_phi = a_phif[a_index];    const int isign = sign(a_side);    if (a_phi.box().size(a_dir) == 3) {        // Linear interpolation of fine ghosts assuming        // all zeros on coarser level.        // we are in a 1-wide box        IntVect iv;        Real pa;        Real factor = 1.0 - 2.0 * a_fineDxDir / (a_fineDxDir + a_crseDxDir);        for (fine_ivsit.begin(); fine_ivsit.ok(); ++fine_ivsit) {            iv = fine_ivsit();            iv[a_dir] -= isign;            // Use linear interpolation            for (int ivar = 0; ivar < a_phif.nComp(); ivar++) {                pa = a_phi(iv, ivar);                a_phi(fine_ivsit(), ivar) = factor * pa;            }        }    } else {        // Quadratic interpolation of fine ghosts assuming        // all zeros on coarser level.        // Symbolic reduced version of CF quadratic stencil        Real pa, pb;        Real c1 = 2.0*(a_crseDxDir-a_fineDxDir)/(a_crseDxDir+    a_fineDxDir); //first inside point        Real c2 =    -(a_crseDxDir-a_fineDxDir)/(a_crseDxDir+3.0*a_fineDxDir); // next point inward        IntVect ivf;        for (fine_ivsit.begin(); fine_ivsit.ok(); ++fine_ivsit) {            ivf = fine_ivsit();            // Use quadratic interpolation            for (int ivar = 0; ivar < a_phif.nComp(); ++ivar) {                ivf[a_dir]-=2*isign;                pa = a_phi(ivf, ivar);                ivf[a_dir]+=isign;                pb = a_phi(ivf, ivar);                ivf[a_dir]+=isign;                a_phi(fine_ivsit(), ivar) = c1*pb + c2*pa;            } //end loop over components        } //end loop over fine intvects    }}

// ---------------------------------------------------------voidNodeQCFI::coarseFineInterp(LevelData<NodeFArrayBox>& a_phiFine,                           const LevelData<NodeFArrayBox>& a_phiCoarse,                           bool a_inhomogeneous){  CH_assert(isDefined());  CH_assert(a_phiFine.nComp() == m_ncomp);  CH_assert(a_phiCoarse.nComp() == m_ncomp);  if (m_coarsenings == 1)    {      m_qcfi2[0]->coarseFineInterp(a_phiFine, a_phiCoarse);    }  else // m_coarsenings >= 2    {      // if m_coarsenings == 2:      // m_qcfi2[1] = new NodeQuadCFInterp2(a_grids, true);      // m_qcfi2[0] = new NodeQuadCFInterp2(a_grids.coarsen(2), false);      // m_inter[0] = new LevelData(a_grids.coarsen(2));      // m_qcfi2[0] interpolates m_inter[0] from a_phiCoarse on all of it.      // m_qcfi2[1] interpolates a_phiFine from m_inter[0] on interface only.      m_qcfi2[0]->coarseFineInterp(*m_inter[0], a_phiCoarse);      m_inter[0]->exchange(m_inter[0]->interval());      // Set domain and dx for coarsest level refined by 2.      ProblemDomain domain(m_domainPenultimate);      Real dx = m_dxPenultimate;      bool homogeneous = !a_inhomogeneous;      const DisjointBoxLayout& dbl0 = m_inter[0]->disjointBoxLayout();      for (DataIterator dit = dbl0.dataIterator(); dit.ok(); ++dit)        {          m_bc((*m_inter[0])[dit()],  dbl0.get(dit()), domain,  dx, homogeneous);        }      for (int interlev = 1; interlev < m_coarsenings - 1; interlev++)        {          m_qcfi2[interlev]->coarseFineInterp(*m_inter[interlev],                                              *m_inter[interlev-1]);          m_inter[interlev]->exchange(m_inter[interlev]->interval());          domain.refine(2);          dx *= 0.5;          const DisjointBoxLayout& dbl = m_inter[interlev]->disjointBoxLayout();          for (DataIterator dit = dbl.dataIterator(); dit.ok(); ++dit)            {              m_bc((*m_inter[interlev])[dit()],  dbl.get(dit()), domain,  dx, homogeneous);            }        }      m_qcfi2[m_coarsenings-1]->coarseFineInterp(a_phiFine,                                                 *m_inter[m_coarsenings-2]);      // Don't need to set boundary conditions for a_phiFine      // because the calling function will do that.    }}

IntVectSet& IntVectSet::grow(int idir, int igrow){  CH_assert(idir >= 0);  CH_assert(idir < SpaceDim);  if (m_isdense) m_dense.grow(idir, igrow);  else          m_ivs.grow(idir, igrow);  return *this;}

// ---------------------------------------------------------FArrayBox&FluxBox::getFlux(const int dir){  CH_assert(m_nvar >0);  CH_assert(dir < SpaceDim);  CH_assert(m_fluxes[dir] != NULL);  return *m_fluxes[dir];}

voidNoFlowAdvectBC::fluxBC(EBFluxFAB&            a_primGdnv,       const EBCellFAB&      a_primCenter,       const EBCellFAB&      a_primExtrap,       const Side::LoHiSide& a_side,       const Real&           a_time,       const EBISBox&        a_ebisBox,       const DataIndex&      a_dit,       const Box&            a_box,       const Box&            a_faceBox,       const int&            a_dir){  CH_assert(m_isDefined);  Box FBox = a_faceBox;  Box cellBox = FBox;  CH_assert(a_primGdnv[a_dir].nComp()==1);  // Determine which side and thus shifting directions  int isign = sign(a_side);  cellBox.shiftHalf(a_dir,isign);  // Is there a domain boundary next to this grid  if (!m_domain.contains(cellBox))    {      cellBox &= m_domain;      // Find the strip of cells next to the domain boundary      Box bndryBox = adjCellBox(cellBox, a_dir, a_side, 1);      // Shift things to all line up correctly      bndryBox.shift(a_dir,-isign);      IntVectSet ivs(bndryBox);      for (VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)        {          const VolIndex& vof = vofit();          Vector<FaceIndex> bndryFaces = a_ebisBox.getFaces(vof, a_dir, a_side);          for (int iface= 0; iface < bndryFaces.size(); iface++)            {              const FaceIndex& face = bndryFaces[iface];              //set all fluxes to zero then fix momentum flux              //solid wall              if (a_dir == m_velComp)                {                  a_primGdnv[a_dir](face, 0) = 0.0;                }              else                {                  a_primGdnv[a_dir](face, 0) = a_primExtrap(vof, 0);                }            }        }    }}

voidEBMGInterp::pwlInterp(LevelData<EBCellFAB>&       a_fineData,                      const LevelData<EBCellFAB>& a_coarData,                      const Interval&             a_variables){  CH_TIME("EBMGInterp::pwlInterp");  CH_assert(a_fineData.ghostVect() == m_ghost);  CH_assert(a_coarData.ghostVect() == m_ghost);  CH_assert(m_doLinear); //otherwise stencils have not been defined  if (m_layoutChanged)    {      if (m_coarsenable)        {          CH_TIME("EBMGInterp::pwlInterp::coarsenable");          EBCellFactory ebcellfact(m_buffEBISL);          LevelData<EBCellFAB> coarsenedFineData(m_buffGrids, m_nComp, m_ghost, ebcellfact);          a_coarData.copyTo(a_variables, coarsenedFineData, a_variables);          fillGhostCellsPWC(coarsenedFineData, m_buffEBISL, m_coarDomain);          for (DataIterator dit = m_fineGrids.dataIterator(); dit.ok(); ++dit)            {              //does incrementonly = true              pwlInterpFAB(a_fineData[dit()],                           m_buffGrids[dit()],                           coarsenedFineData[dit()],                           dit(),                           a_variables);            }        }      else        {          CH_TIME("EBMGInterp::pwlInterp::uncoarsenable");          EBCellFactory ebcellfact(m_buffEBISL);          fillGhostCellsPWC((LevelData<EBCellFAB>&)a_coarData, m_coarEBISL, m_coarDomain);          LevelData<EBCellFAB> refinedCoarseData(m_buffGrids, m_nComp, m_ghost, ebcellfact);          for (DataIterator dit = m_coarGrids.dataIterator(); dit.ok(); ++dit)            {              refinedCoarseData[dit()].setVal(0.);              pwlInterpFAB(refinedCoarseData[dit()],                           m_coarGrids[dit()],                           a_coarData[dit()],                           dit(),                           a_variables);            }          EBAddOp op;          refinedCoarseData.copyTo(a_variables, a_fineData, a_variables, m_copierRCtoF, op);        }    }  else    {      pwcInterpMG(a_fineData, a_coarData, a_variables);    }}

// Compute the speed of soundvoid PolytropicPhysics::soundSpeed(FArrayBox& a_speed,                                   const FArrayBox& a_U,                                   const Box&       a_box){    CH_assert(isDefined());    CH_assert(a_U.contains(a_box));    CH_assert(a_speed.contains(a_box));    FORT_SOUNDSPEEDF(CHF_CONST_FRA1(a_speed, 0),                     CHF_CONST_FRA(a_U),                     CHF_BOX(a_box));}

void EBPoissonOp::GSColorAllRegular(LevelData<EBCellFAB>&        a_phi,                  const LevelData<EBCellFAB>&  a_rhs,                  const IntVect&               a_color,                  const Real&                  a_weight,                  const bool&                  a_homogeneousPhysBC){  CH_TIME("EBPoissonOp::GSColorAllRegular");  CH_assert(a_rhs.ghostVect() == m_ghostCellsRHS);  CH_assert(a_phi.ghostVect() == m_ghostCellsPhi);  int nComps = a_phi.nComp();  for (DataIterator dit = a_phi.dataIterator(); dit.ok(); ++dit)    {      Box dblBox(m_eblg.getDBL().get(dit()));      BaseFab<Real>& phiFAB       = (a_phi[dit()] ).getSingleValuedFAB();      const BaseFab<Real>& rhsFAB = (a_rhs[dit()] ).getSingleValuedFAB();      Box loBox[SpaceDim],hiBox[SpaceDim];      int hasLo[SpaceDim],hasHi[SpaceDim];      applyDomainFlux(loBox, hiBox, hasLo, hasHi,                      dblBox, nComps, phiFAB,                      a_homogeneousPhysBC, dit(),m_beta);      IntVect loIV = dblBox.smallEnd();      IntVect hiIV = dblBox.bigEnd();      for (int idir = 0; idir < SpaceDim; idir++)        {          if (loIV[idir] % 2 != a_color[idir])            {              loIV[idir]++;            }        }      if (loIV <= hiIV)        {          Box coloredBox(loIV, hiIV);          for (int comp=0; comp<a_phi.nComp(); comp++)            {              FORT_DOALLREGULARMULTICOLOR(CHF_FRA1(phiFAB,comp),                                          CHF_CONST_FRA1(rhsFAB,comp),                                          CHF_CONST_REAL(a_weight),                                          CHF_CONST_REAL(m_alpha),                                          CHF_CONST_REAL(m_beta),                                          CHF_CONST_REALVECT(m_dx),                                          CHF_BOX(coloredBox));            }        }    }}

voidEBPlanarShockIBC::initialize(LevelData<EBCellFAB>& a_conState,           const EBISLayout& a_ebisl) const{ CH_assert(m_isDefined); CH_assert(m_isFortranCommonSet);  // Iterator of all grids in this level  for(DataIterator dit = a_conState.dataIterator(); dit.ok(); ++dit)    {      const EBISBox& ebisBox = a_ebisl[dit()];      // Storage for current grid      EBCellFAB& conFAB = a_conState[dit()];      BaseFab<Real>& regConFAB = conFAB.getSingleValuedFAB();      // Box of current grid      Box uBox = regConFAB.box();      uBox &= m_domain;      // Set up initial condition in this grid      /**/      FORT_PLANARSHOCKINIT(CHF_CONST_FRA(regConFAB),                           CHF_CONST_REAL(m_dx),                           CHF_BOX(uBox));      /**/      //now for the multivalued cells.  Since it is pointwise,      //the regular calc is correct for all single-valued cells.      IntVectSet ivs = ebisBox.getMultiCells(uBox);      for(VoFIterator vofit(ivs, ebisBox.getEBGraph()); vofit.ok(); ++vofit)        {          const VolIndex& vof = vofit();          const IntVect& iv = vof.gridIndex();          RealVect momentum;          Real energy, density;          /**/          FORT_POINTPLANARSHOCKINIT(CHF_REAL(density),                                    CHF_REALVECT(momentum),                                    CHF_REAL(energy),                                    CHF_CONST_INTVECT(iv),                                    CHF_CONST_REAL(m_dx));          /**/          conFAB(vof, CRHO) = density;          conFAB(vof, CENG) = energy;          for(int idir = 0; idir < SpaceDim; idir++)            {              conFAB(vof, CMOMX+idir) = momentum[idir];            }        }//end loop over multivalued cells    } //end loop over boxes}

// Define new AMR levelvoid AMRLevelPluto::define(AMRLevel*            a_coarserLevelPtr,                           const ProblemDomain& a_problemDomain,                           int                  a_level,                           int                  a_refRatio){  if (s_verbosity >= 3)  {    pout() << "AMRLevelPluto::define " << a_level << endl;  }  // Call inherited define  AMRLevel::define(a_coarserLevelPtr,                   a_problemDomain,                   a_level,                   a_refRatio);  // Get setup information from the next coarser level  if (a_coarserLevelPtr != NULL)  {    AMRLevelPluto* amrGodPtr = dynamic_cast<AMRLevelPluto*>(a_coarserLevelPtr);    if (amrGodPtr != NULL)    {      m_cfl = amrGodPtr->m_cfl;      m_domainLength = amrGodPtr->m_domainLength;      m_refineThresh = amrGodPtr->m_refineThresh;      m_tagBufferSize = amrGodPtr->m_tagBufferSize;    }    else    {      MayDay::Error("AMRLevelPluto::define: a_coarserLevelPtr is not castable to AMRLevelPluto*");    }  }  // Compute the grid spacing  m_dx = m_domainLength / (a_problemDomain.domainBox().bigEnd(0)-a_problemDomain.domainBox().smallEnd(0)+1.);  // Nominally, one layer of ghost cells is maintained permanently and  // individual computations may create local data with more  m_numGhost = 1;  CH_assert(m_patchPluto != NULL);  CH_assert(isDefined());  m_patchPluto->define(m_problem_domain,m_dx,m_level,m_numGhost);  // Get additional information from the patch integrator  m_numStates      = m_patchPluto->numConserved();  m_ConsStateNames = m_patchPluto->ConsStateNames();  m_PrimStateNames = m_patchPluto->PrimStateNames();}

void LinElastPhysics::consToPrim(FArrayBox&       a_W,    const FArrayBox& a_U,    const Box&       a_box){    //JK Our primitives are our conserved variables    //JK pout() << "LinElastPhysics::consToPrim" << endl;    CH_assert(isDefined());    CH_assert(a_U.box().contains(a_box));    CH_assert(a_W.box().contains(a_box));    FORT_COPYF(CHF_FRA(a_W),        CHF_CONST_FRA(a_U),        CHF_BOX(a_box));}

void EBLevelTransport::divergeF(LevelData<EBCellFAB>&         a_divergeF,         LevelData<BaseIVFAB<Real> >&  a_massDiff,         EBFluxRegister&               a_fineFluxRegister,         EBFluxRegister&               a_coarFluxRegister,         LevelData<EBCellFAB>&         a_consState,   //not really changed         LevelData<EBCellFAB>&         a_normalVel,         const LevelData<EBFluxFAB>&   a_advVel,         const LayoutData< Vector <BaseIVFAB<Real> * > >& a_coveredAdvVelMinu,         const LayoutData< Vector <BaseIVFAB<Real> * > >& a_coveredAdvVelPlus,         const LevelData<EBCellFAB>&   a_source,         const LevelData<EBCellFAB>&   a_consStateCoarseOld,         const LevelData<EBCellFAB>&   a_consStateCoarseNew,         const LevelData<EBCellFAB>&   a_normalVelCoarseOld,         const LevelData<EBCellFAB>&   a_normalVelCoarseNew,         const Real&                   a_time,         const Real&                   a_coarTimeOld,         const Real&                   a_coarTimeNew,         const Real&                   a_dt){  Interval consInterv(0, m_nCons-1);  Interval fluxInterv(0, m_nFlux-1);  CH_assert(isDefined());  CH_assert(a_consState.disjointBoxLayout() == m_thisGrids);  //fill a_consState with data interpolated between constate and coarser data  fillCons(a_consState, a_normalVel, a_consStateCoarseOld, a_consStateCoarseNew, a_normalVelCoarseOld, a_normalVelCoarseNew, a_time, a_coarTimeOld, a_coarTimeNew);   // clear flux registers with fine level  //remember this level is the coarse level to the fine FR  if(m_hasFiner)    {      a_fineFluxRegister.setToZero();    }  //compute flattening coefficients. this saves a ghost cell. woo hoo.  computeFlattening(a_consState, a_time, a_dt);  //this includes copying flux into flux interpolant and updating  //regular grids and incrementing flux registers.  doRegularUpdate(a_divergeF, a_consState, a_normalVel, a_advVel, a_coveredAdvVelMinu, a_coveredAdvVelPlus, a_fineFluxRegister, a_coarFluxRegister, a_source, a_time, a_dt);  //this does irregular update and deals with flux registers.  //also computes the mass increment  doIrregularUpdate(a_divergeF, a_consState, a_fineFluxRegister, a_coarFluxRegister, a_massDiff, a_time, a_dt);}

void EBPoissonOp::residual(LevelData<EBCellFAB>&       a_residual,         const LevelData<EBCellFAB>& a_phi,         const LevelData<EBCellFAB>& a_rhs,         bool                        a_homogeneousPhysBC){  CH_TIME("EBPoissonOp::residual");  //this is a multigrid operator so only homogeneous CF BC  //and null coar level  CH_assert(a_residual.ghostVect() == m_ghostCellsRHS);  CH_assert(a_phi.ghostVect() == m_ghostCellsPhi);  DataIterator dit = a_residual.dataIterator();  applyOp(a_residual,a_phi, a_homogeneousPhysBC, dit );  axby(a_residual,a_residual,a_rhs,-1.0, 1.0);}

voidcompareError(const LevelData<EBCellFAB>& a_errorFine,             const EBISLayout& a_ebislFine,             const DisjointBoxLayout& a_gridsFine,             const Box& a_domainFine,             const LevelData<EBCellFAB>& a_errorCoar,             const EBISLayout& a_ebislCoar,             const DisjointBoxLayout& a_gridsCoar,             const Box& a_domainCoar){    CH_assert(a_errorFine.nComp() == 1);    CH_assert(a_errorCoar.nComp() == 1);    pout() << "==============================================" << endl;    EBNormType::NormMode normtype = EBNormType::OverBoth;    pout() << endl << "Using all uncovered cells." << endl  ;    for (int inorm = 0; inorm <= 2; inorm++)    {        if (inorm == 0)        {            pout() << endl << "Using max norm." << endl;        }        else        {            pout() << endl << "Using L-" << inorm << "norm." << endl;        }        int comp = 0;        Real coarnorm = EBArith::norm(a_errorCoar,                                      a_gridsCoar, a_ebislCoar,                                      comp, inorm, normtype);        Real finenorm = EBArith::norm(a_errorFine,                                      a_gridsFine, a_ebislFine,                                      comp, inorm, normtype);        pout() << "Coarse Error Norm = " << coarnorm << endl;        pout() << "Fine   Error Norm = " << finenorm << endl;        if ((Abs(finenorm) > 1.0e-12) && (Abs(coarnorm) > 1.0e-12))        {            Real order = log(Abs(coarnorm/finenorm))/log(2.0);            pout() << "Order of scheme = " << order << endl;        }    }    pout() << "==============================================" << endl ;;}