C++ CH_assert函数代码示例

// Generic function to reduce the minloc,maxloc,avg of a Real over all procs onto rank0
int 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;
}

void
EBCellFAB::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 levelGSRB
void 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 fluxes
void 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
    }
}

// ---------------------------------------------------------
void
NodeQCFI::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];
}

void
NoFlowAdvectBC::
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);
                }
            }
        }
    }
}

void
EBMGInterp::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);
    }
}