C++ FArrayBox类代码示例

/* ********************************************************************* */
void PatchPluto::updateSolution(FArrayBox&       a_U,
                                FArrayBox&       a_Utmp,
                                const FArrayBox& a_dV,
                                FArrayBox&       split_tags,
                                BaseFab<unsigned char>& a_Flags,
                                FluxBox&         a_F,
                                Time_Step        *Dts,
                                const Box&       UBox, 
                                Grid *grid)
/*
 *
 *
 *
 *
 *********************************************************************** */
{
  CH_assert(isDefined());
  CH_assert(UBox == m_currentBox);

  int nv, in;
  int nxf, nyf, nzf, indf;
  int nxb, nyb, nzb;
  int *i, *j, *k;
  int ii, jj, kk;

  double ***UU[NVAR];
  double *inv_dl, dl2, cylr;
  static Data d;
 #ifdef SKIP_SPLIT_CELLS
  double ***splitcells;
 #endif
 #if (PARABOLIC_FLUX & EXPLICIT)
  static double **dcoeff;
 #endif   
  Index indx;
  static State_1D state;

  Riemann_Solver *Riemann;
  Riemann = rsolver;

#if TIME_STEPPING == RK2 
  double wflux = 0.5;
#else
  double wflux = 1.;
#endif

/* -----------------------------------------------------------------
               Check algorithm compatibilities
   ----------------------------------------------------------------- */

  if (NX1_TOT > NMAX_POINT || NX2_TOT > NMAX_POINT || NX3_TOT > NMAX_POINT){
    print ("!updateSolution (Euler): need to re-allocate matrix\n");
    QUIT_PLUTO(1);
  }

/* -----------------------------------------------------------------
                          Allocate memory
   ----------------------------------------------------------------- */

  #if GEOMETRY != CARTESIAN
   for (nv = 0; nv < NVAR; nv++) a_U.divide(a_dV,0,nv);
   #if CHOMBO_CONS_AM == YES
    #if ROTATING_FRAME == YES
     Box curBox = a_U.box();
     for(BoxIterator bit(curBox); bit.ok(); ++bit) {
       const IntVect& iv = bit();
       a_U(iv,iMPHI) /= a_dV(iv,1);
       a_U(iv,iMPHI) -= a_U(iv,RHO)*a_dV(iv,1)*g_OmegaZ;
     }
    #else
     a_U.divide(a_dV,1,iMPHI);
    #endif
   #endif
  #else
   if (g_stretch_fact != 1.) a_U /= g_stretch_fact; 
  #endif

  for (nv = 0; nv < NVAR; nv++){
    UU[nv] = ArrayMap(NX3_TOT, NX2_TOT, NX1_TOT, a_U.dataPtr(nv));
  }
  #ifdef SKIP_SPLIT_CELLS
   splitcells = ArrayBoxMap(KBEG, KEND, JBEG, JEND, IBEG, IEND,
                            split_tags.dataPtr(0));
  #endif
  #if (TIME_STEPPING == RK2)
   d.flag = ArrayCharMap(NX3_TOT, NX2_TOT, NX1_TOT,a_Flags.dataPtr(0));
  #endif
  #if RESISTIVE_MHD != NO
   if (d.J == NULL) d.J = ARRAY_4D(3,NX3_MAX, NX2_MAX, NX1_MAX, double);
  #endif
 
/* -----------------------------------------------------------
         Allocate static memory areas
   -----------------------------------------------------------  */

  if (state.flux == NULL){

    MakeState (&state);

//.........这里部分代码省略.........

void
CellConservativeLinear::interp (const FArrayBox& crse,
                                int              crse_comp,
                                FArrayBox&       fine,
                                int              fine_comp,
                                int              ncomp,
                                const Box&       fine_region,
                                const IntVect&   ratio,
                                const Geometry&  crse_geom,
                                const Geometry&  fine_geom,
                                Array<BCRec>&    bcr,
                                int              actual_comp,
                                int              actual_state)
{
    BL_PROFILE("CellConservativeLinear::interp()");
    BL_ASSERT(bcr.size() >= ncomp);

    //
    // Make box which is intersection of fine_region and domain of fine.
    //
    Box target_fine_region = fine_region & fine.box();
    //
    // crse_bx is coarsening of target_fine_region, grown by 1.
    //
    Box crse_bx = CoarseBox(target_fine_region,ratio);
    //
    // Slopes are needed only on coarsening of target_fine_region.
    //
    Box cslope_bx(crse_bx);
    cslope_bx.grow(-1);
    //
    // Make a refinement of cslope_bx
    //
    Box fine_version_of_cslope_bx = BoxLib::refine(cslope_bx,ratio);
    //
    // Get coarse and fine edge-centered volume coordinates.
    //
    Array<Real> fvc[BL_SPACEDIM];
    Array<Real> cvc[BL_SPACEDIM];
    int dir;
    for (dir = 0; dir < BL_SPACEDIM; dir++)
    {
        fine_geom.GetEdgeVolCoord(fvc[dir],fine_version_of_cslope_bx,dir);
        crse_geom.GetEdgeVolCoord(cvc[dir],crse_bx,dir);
    }
    //
    // alloc tmp space for slope calc.
    //
    // In ucc_slopes and lcc_slopes , there is a slight abuse of
    // the number of compenents argument
    // --> there is a slope for each component in each coordinate
    //     direction
    //
    FArrayBox ucc_slopes(cslope_bx,ncomp*BL_SPACEDIM);
    FArrayBox lcc_slopes(cslope_bx,ncomp*BL_SPACEDIM);
    FArrayBox slope_factors(cslope_bx,BL_SPACEDIM);

    FArrayBox  cmax(cslope_bx,ncomp);
    FArrayBox  cmin(cslope_bx,ncomp);
    FArrayBox alpha(cslope_bx,ncomp);

    Real* fdat       = fine.dataPtr(fine_comp);
    const Real* cdat = crse.dataPtr(crse_comp);
    Real* ucc_xsldat = ucc_slopes.dataPtr(0);
    Real* lcc_xsldat = lcc_slopes.dataPtr(0);
    Real* xslfac_dat = slope_factors.dataPtr(0);
#if (BL_SPACEDIM>=2)
    Real* ucc_ysldat = ucc_slopes.dataPtr(ncomp);
    Real* lcc_ysldat = lcc_slopes.dataPtr(ncomp);
    Real* yslfac_dat = slope_factors.dataPtr(1);
#endif
#if (BL_SPACEDIM==3)
    Real* ucc_zsldat = ucc_slopes.dataPtr(2*ncomp);
    Real* lcc_zsldat = lcc_slopes.dataPtr(2*ncomp);
    Real* zslfac_dat = slope_factors.dataPtr(2);
#endif

    const int* flo    = fine.loVect();
    const int* fhi    = fine.hiVect();
    const int* clo    = crse.loVect();
    const int* chi    = crse.hiVect();
    const int* fblo   = target_fine_region.loVect();
    const int* fbhi   = target_fine_region.hiVect();
    const int* csbhi  = cslope_bx.hiVect();
    const int* csblo  = cslope_bx.loVect();
    int lin_limit     = (do_linear_limiting ? 1 : 0);
    const int* cvcblo = crse_bx.loVect();
    const int* fvcblo = fine_version_of_cslope_bx.loVect();
    int slope_flag    = 1;

    int cvcbhi[BL_SPACEDIM];
    int fvcbhi[BL_SPACEDIM];

    for (dir=0; dir<BL_SPACEDIM; dir++)
    {
        cvcbhi[dir] = cvcblo[dir] + cvc[dir].size() - 1;
        fvcbhi[dir] = fvcblo[dir] + fvc[dir].size() - 1;
    }

    D_TERM(Real* voffx = new Real[fvc[0].size()]; ,

/* ********************************************************************* */
void PatchPluto::updateSolution(FArrayBox&  a_U,
                                FArrayBox&  a_Utmp,
                                FArrayBox&  split_tags,
                                BaseFab<unsigned char>& a_Flags,
                                FluxBox&    a_F,
                                Time_Step   *Dts,
                                const Box&  UBox,
                                Grid *grid)
/*
 *
 *
 *
 *
 *********************************************************************** */
{
  CH_assert(isDefined());
  CH_assert(UBox == m_currentBox);

  int nv, in;
  int nxf, nyf, nzf, indf;
  int nxb, nyb, nzb;
  int *i, *j, *k;
  int ii, jj, kk;

  double ***UU[NVAR];
 #ifdef SKIP_SPLIT_CELLS
  double ***splitcells;
 #endif
  double *inv_dl, dl2;
  static Data d;
  static Data_Arr UU_1;
  static double ***C_dt[NVAR];
  static double **u;
  #if (PARABOLIC_FLUX & EXPLICIT)
   static double **dcoeff;
  #endif   
  Index indx;
  static State_1D state;

  Riemann_Solver *Riemann;
  Riemann = rsolver;

/* -----------------------------------------------------------------
               Check algorithm compatibilities
   ----------------------------------------------------------------- */

  #if !(GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL || GEOMETRY == SPHERICAL)
   print1 ("! RK only works in cartesian or cylindrical/spherical coordinates\n");
   QUIT_PLUTO(1);
  #endif     

  if (NX1_TOT > NMAX_POINT || NX2_TOT > NMAX_POINT || NX3_TOT > NMAX_POINT){
    print ("!updateSolution (RK_MID): need to re-allocate matrix\n");
    QUIT_PLUTO(1);
  }

/* -----------------------------------------------------------------
                          Allocate memory
   ----------------------------------------------------------------- */

  for (nv = 0; nv < NVAR; nv++){
    UU[nv] = ArrayMap(NX3_TOT, NX2_TOT, NX1_TOT, a_U.dataPtr(nv));
  }

  #ifdef SKIP_SPLIT_CELLS
   splitcells = ArrayBoxMap(KBEG, KEND, JBEG, JEND, IBEG, IEND, 
                                split_tags.dataPtr(0));
  #endif

/* -----------------------------------------------------------
         Allocate static memory areas
   -----------------------------------------------------------  */

  if (state.flux == NULL){

    MakeState (&state);

    nxf = nyf = nzf = 1;
    D_EXPAND(nxf = NMAX_POINT;  ,

// Set boundary fluxes
void VelIBC::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_isDefined == true);

  // In periodic case, this doesn't do anything
  if (!m_domain.isPeriodic(a_dir))
    {
      // This needs to be fixed
      // CH_assert(m_isBCvalSet);

      int lohisign;
      Box tmp = a_WGdnv.box() & m_domain;
      Real bcVal;

      // Determine which side and thus shifting directions
      if (a_side == Side::Lo)
        {
          lohisign = -1;
          bcVal = m_bcVal[a_dir][0];
        }
      else
        {
          lohisign = 1;
          bcVal = m_bcVal[a_dir][1];
        }

      tmp.shiftHalf(a_dir,lohisign);

      // Is there a domain boundary next to this grid
      if (!m_domain.contains(tmp))
        {
          tmp &= m_domain;

          Box boundaryBox;

          // Find the strip of cells next to the domain boundary
          if (a_side == Side::Lo)
            {
              boundaryBox = bdryLo(tmp,a_dir);
            }
          else
            {
              boundaryBox = bdryHi(tmp,a_dir);
            }

          // Set the boundary fluxes
          FORT_SOLIDVELBCF(CHF_FRA(a_WGdnv),
                           CHF_CONST_FRA(a_Wextrap),
                           CHF_CONST_REAL(bcVal),
                           CHF_CONST_INT(lohisign),
                           CHF_CONST_REAL(m_dx),
                           CHF_CONST_INT(a_dir),
                           CHF_BOX(boundaryBox));
        }
    }
}

void
LevelFluxRegisterEdge::incrementFine(
                                     FArrayBox& a_fineFlux,
                                     Real a_scale,
                                     const DataIndex& a_fineDataIndex,
                                     const Interval& a_srcInterval,
                                     const Interval& a_dstInterval,
                                     int a_dir,
                                     Side::LoHiSide a_sd)
{
  CH_assert(isDefined());
  CH_assert(!a_fineFlux.box().isEmpty());
  CH_assert(a_srcInterval.size() == a_dstInterval.size());
  CH_assert(a_srcInterval.begin() >= 0);
  CH_assert(a_srcInterval.end() < a_fineFlux.nComp());
  CH_assert(a_dstInterval.begin() >= 0);
  CH_assert(a_dstInterval.end() < m_nComp);

  CH_assert(a_dir >= 0);
  CH_assert(a_dir < SpaceDim);
  CH_assert((a_sd == Side::Lo)||(a_sd == Side::Hi));
  //
  //
  //denom is the number of fine faces per coarse face
  //this is intrinsically dimension-dependent
#if (CH_SPACEDIM == 2)
  Real denom = 1;
#elif (CH_SPACEDIM == 3)
  Real denom = m_nRefine;
#else
  // This code doesn't make any sense in 1D, and hasn't been implemented
  // for DIM > 3
  Real denom = -1.0;
  MayDay::Error("LevelFluxRegisterEdge -- bad SpaceDim");
#endif

  Real scale = a_scale/denom;

  // need which fluxbox face we're doing this for
  Box thisBox = a_fineFlux.box();
  int fluxComp = -1;
  for (int sideDir=0; sideDir<SpaceDim; sideDir++)
  {
    // we do nothing in the direction normal to face
    if (sideDir != a_dir)
    {
      if (thisBox.type(sideDir) == IndexType::CELL)
      {
        fluxComp = sideDir;
      }
    }
  }
  CH_assert (fluxComp >= 0);
  int regcomp = getRegComp(a_dir, fluxComp);

  FluxBox& thisReg = m_fabFine[index(a_dir, a_sd)][a_fineDataIndex];
  FArrayBox& reg = thisReg[regcomp];

  a_fineFlux.shiftHalf(a_dir, sign(a_sd));

  // this is a way of geting a face-centered domain
  // box which we can then use to intersect with things
  // to screen out cells outside the physical domain
  // (nothing is screened out in periodic case)
  Box shiftedValidDomain = m_domainCoarse.domainBox();
  shiftedValidDomain.grow(2);
  shiftedValidDomain &= m_domainCoarse;
  shiftedValidDomain.surroundingNodes(regcomp);

  BoxIterator regIt(reg.box() & shiftedValidDomain);
  for (regIt.begin(); regIt.ok(); ++regIt)
    {
      const IntVect& coarseIndex = regIt();
      // create a cell-centered box, then shift back to face-centered
      Box box(coarseIndex, coarseIndex);
      box.shiftHalf(regcomp,-1);
      // to avoid adding in edges which do not overlie coarse-grid
      // edges, will refine only in non-fluxComp directions to
      // determine box from which to grab fluxes.
      IntVect refineVect(m_nRefine*IntVect::Unit);
      //refineVect.setVal(fluxComp,1);
      box.refine(refineVect);
      if (a_sd == Side::Lo) box.growLo(a_dir,-(m_nRefine-1));
      else                 box.growHi(a_dir,-(m_nRefine-1));
      BoxIterator fluxIt(box);
      for (fluxIt.begin(); fluxIt.ok(); ++fluxIt)
        {
          int src = a_srcInterval.begin();
          int dest =  a_dstInterval.begin();
          for ( ; src <=a_srcInterval.end(); ++src,++dest)
            reg(coarseIndex, dest) += scale*a_fineFlux(fluxIt(), src);
        }
    }
  a_fineFlux.shiftHalf(a_dir, -sign(a_sd));
}

int main( int a_argc, char* a_argv[] )
{
#ifdef CH_MPI
   // Start MPI
   MPI_Init( &a_argc, &a_argv );
   setChomboMPIErrorHandler();
#endif

   string RZ_mapping_file_name;
   string field_file_name;
   string block_name;

   if (a_argc==4) {
      RZ_mapping_file_name = a_argv[1];
      field_file_name = a_argv[2];
      block_name = a_argv[3];
   }
   else {
      cout << "Usage:  magGeomData...ex <RZ_mapping_file> [<field_file>] <blockname>," << endl;
      cout << "        where <RZ_mapping_file> defines the coordinate mapping, and" << endl;
      cout << "        where <field_file> defines the magnetic field, and" << endl;
      cout << "        where <blockname> is one of lcore, rcore, lcsol, rcsol, lsol, rsol, lpf, rpf, mcore, mcsol" << endl;
      return -1;
   }

   // Construct the field data object
   FieldData field_data(RZ_mapping_file_name, field_file_name, block_name);

   /*
     Fill an FArrayBox with physical coordinates.  In the
     new mesh generator code using the Brackbill-Saltzmann
     approach, the Euler equation will solved for these
     coordinates.  Here, however, we simply create a grid
     in the unit square, which is the assumed mapped
     domain used internally in the FieldData object and
     get the corresponding physical coordinates.
   */

   // Get a grid on the unit square.  The argument sets the size.
   FArrayBox xi;
   getUnitSquareGrid(10, 50, xi);

   const Box& box(xi.box());

   // Get the corresponding physical coordinates.
   FArrayBox physical_coordinates(box,2);
   field_data.getPhysicalCoordinates(xi, physical_coordinates);

   // Write a file containing the physical coordinates that
   // can be plotted with the plot_coords.m script using
   // either MatLab or Octave.  The file name is "coordsX",
   // where X is the block number.
   field_data.writePhysicalCoordinates(physical_coordinates);

   FArrayBox field_unit_vector(box,6);

#ifdef USE_DCT_FIELD  

   // Using the discrete cosine transform (DCT) field representation,
   // whose coefficients are contained in "field_file_name", get the field
   // unit vector (first two components) and the derivatives of its components
   // with respect to the physical coordinates (next 4 components) at
   // this set of physical coordinates.  The numbering of the 6 components
   // is explained at the top of FieldData::convertToMappedDerivatives(),
   // which converts the derivatives in physical space to derivatives
   // in mapped space.

   field_data.getFieldUnitVectorFromDCT(physical_coordinates, field_unit_vector);
   field_data.convertToMappedDerivatives(xi, field_unit_vector);

   FArrayBox magnetic_flux(physical_coordinates.box(),1);
   field_data.getMagneticFluxFromDCT(physical_coordinates, magnetic_flux);

   field_data.writeMagneticFlux(physical_coordinates, magnetic_flux);
#else

   // Using the field data contained in the coordinate mapping
   // file "geometry_file_name" get the field unit vector (first two
   // components) and the derivatives of its components with
   // respect to the mapped coordinates (next 4 components) at
   // this set of physical coordinates.  The numbering of the 6 components
   // is explained at the top of FieldData::getFieldUnitVectorFromMapping().

   field_data.getFieldUnitVectorFromMappingFile(physical_coordinates, field_unit_vector);

#endif

   // Write a file containing the field unit vector that
   // can be plotted with the plot_vectors.m script using
   // either MatLab or Octave.  The file name is "vectorsX",
   // where X is the block number.
   field_data.writeVectors(physical_coordinates, field_unit_vector);

#ifdef CH_MPI
   MPI_Finalize();
#endif

   return 0;
}

void ConstBCFunction::operator()(FArrayBox&           a_state,
                                 const Box&           a_valid,
                                 const ProblemDomain& a_domain,
                                 Real                 a_dx,
                                 bool                 a_homogeneous)
{
  const Box& domainBox = a_domain.domainBox();

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    if (!a_domain.isPeriodic(idir))
    {
      for (SideIterator sit; sit.ok(); ++sit)
      {
        Side::LoHiSide side = sit();

        if (a_valid.sideEnd(side)[idir] == domainBox.sideEnd(side)[idir])
        {
          int  bcType;
          Real bcValue;

          if (side == Side::Lo)
          {
            bcType  = m_loSideType [idir];
            bcValue = m_loSideValue[idir];
          }
          else
          {
            bcType  = m_hiSideType [idir];
            bcValue = m_hiSideValue[idir];
          }

          if (bcType == 0)
          {
            // Neumann BC
            int isign = sign(side);

            Box toRegion = adjCellBox(a_valid, idir, side, 1);
            toRegion &= a_state.box();

            Box fromRegion = toRegion;
            fromRegion.shift(idir, -isign);

            a_state.copy(a_state, fromRegion, 0, toRegion, 0, a_state.nComp());

            if (!a_homogeneous)
            {
              for (BoxIterator bit(toRegion); bit.ok(); ++bit)
              {
                const IntVect& ivTo = bit();
                // IntVect ivClose = ivTo - isign*BASISV(idir);

                for (int icomp = 0; icomp < a_state.nComp(); icomp++)
                {
                  a_state(ivTo, icomp) += Real(isign)*a_dx*bcValue;
                }
              }
            }
          }
          else if (bcType == 1)
          {
            // Dirichlet BC
            int isign = sign(side);

            Box toRegion = adjCellBox(a_valid, idir, side, 1);
            toRegion &= a_state.box();

            for (BoxIterator bit(toRegion); bit.ok(); ++bit)
            {
              const IntVect& ivTo = bit();

              IntVect ivClose = ivTo -   isign*BASISV(idir);
              // IntVect ivFar   = ivTo - 2*isign*BASISV(idir);

              Real inhomogVal = 0.0;

              if (!a_homogeneous)
              {
                inhomogVal = bcValue;
              }

              for (int icomp = 0; icomp < a_state.nComp(); icomp++)
              {
                Real nearVal = a_state(ivClose, icomp);
                // Real farVal  = a_state(ivFar,   icomp);

                Real ghostVal = linearInterp(inhomogVal, nearVal);

                a_state(ivTo, icomp) = ghostVal;
              }
            }
          }
          else
          {
            MayDay::Abort("ConstBCFunction::operator() - unknown BC type");
          }
        } // if ends match
      } // end loop over sides
    } // if not periodic in this direction
  } // end loop over directions
//.........这里部分代码省略.........

void setFunc(const Box&, int m, FArrayBox& F)
{
  F.setVal(procID());
}

void PatchGodunov::PPMNormalPred(FArrayBox&       a_WMinus,
                                 FArrayBox&       a_WPlus,
                                 const Real&      a_dt,
                                 const Real&      a_dx,
                                 const FArrayBox& a_W,
                                 const FArrayBox& a_flat,
                                 const int&       a_dir,
                                 const Box&       a_box)
{
  int numprim = m_gdnvPhysics->numPrimitives();

  Box faceBox = a_box;
  // added by petermc, 22 Sep 2008:
  // for 4th order, need extra faces in all the directions
  if (m_highOrderLimiter) faceBox.grow(1);
  faceBox.surroundingNodes(a_dir);
  FArrayBox WFace(faceBox,numprim);

  // Return WFace on face-centered faceBox.
  m_util.PPMFaceValues(WFace,a_W,numprim,
                       m_useCharLimiting || m_usePrimLimiting,
                       a_dir,faceBox,m_currentTime,m_gdnvPhysics);

  // To save on storage, we use the input values as temporaries for the
  // delta's
  a_WMinus.setVal(0.0);
  a_WPlus .setVal(0.0);

  a_WMinus -= a_W;
  a_WPlus  -= a_W;

  WFace.shiftHalf(a_dir,1);
  a_WMinus += WFace;

  WFace.shift(a_dir,-1);
  a_WPlus  += WFace;

  FArrayBox lambda(a_box, numprim);
  m_gdnvPhysics->charValues(lambda, a_W, a_dir,a_box);

  if (m_useCharLimiting && m_usePrimLimiting)
    {
      MayDay::Error("PatchGodunov::PPMNormalPred: Attempt to limit slopes in primitive AND characteristic coordinates - not implemented");
    }

  // Apply limiter on characteristic or primitive variables. Either
  // way, must end up with characteristic variables to pass to to the
  // normal predictor utility.  Currently, cannot do both.

  // If doing characteristic limiting then transform before limiting
  if (m_useCharLimiting)
    {
      // Transform from primitive to characteristic variables
      m_gdnvPhysics->charAnalysis(a_WMinus,a_W,a_dir,a_box);
      m_gdnvPhysics->charAnalysis(a_WPlus ,a_W,a_dir,a_box);
    }

  if (m_useCharLimiting || m_usePrimLimiting)
    {
      // Do slope limiting
      // m_util.PPMLimiter(a_WMinus,a_WPlus,numprim,a_box);

      // petermc, 4 Sep 2008:  included a_W and a_dir in argument list
      m_util.PPMLimiter(a_WMinus, a_WPlus, a_W, numprim, a_dir, a_box);

      // Do slope flattening
      if (m_useFlattening)
        {
          m_util.applyFlattening(a_WMinus,a_flat,a_box);
          m_util.applyFlattening(a_WPlus ,a_flat,a_box);
        }
    }

  // If not doing characteristic limiting then transform after any limiting
  if (!m_useCharLimiting)
    {
      // Transform from primitive to characteristic variables
      m_gdnvPhysics->charAnalysis(a_WMinus,a_W,a_dir,a_box);
      m_gdnvPhysics->charAnalysis(a_WPlus ,a_W,a_dir,a_box);
    }

  // To the normal prediction in characteristic variables
  m_util.PPMNormalPred(a_WMinus,a_WPlus,lambda,a_dt / a_dx,numprim,a_box);

  // Construct the increments to the primitive variables
  m_gdnvPhysics->charSynthesis(a_WMinus,a_W,a_dir,a_box);
  m_gdnvPhysics->charSynthesis(a_WPlus ,a_W,a_dir,a_box);

  // Apply a physics-dependent post-normal predictor step:
  // For example:
  //   - adjust/bound delta's so constraints on extrapolated primitive
  //     quantities are enforced (density and pressure > 0).
  //   - compute source terms that depend on the spatially varying
  //     coefficients.
  m_gdnvPhysics->postNormalPred(a_WMinus,a_WPlus,a_W,a_dt,a_dx,a_dir,a_box);

  // Compute the state from the increments
  a_WMinus += a_W;
  a_WPlus  += a_W;
}

void PatchGodunov::PLMNormalPred(FArrayBox&       a_WMinus,
                                 FArrayBox&       a_WPlus,
                                 const Real&      a_dt,
                                 const Real&      a_dx,
                                 const FArrayBox& a_W,
                                 const FArrayBox& a_flat,
                                 const int&       a_dir,
                                 const Box&       a_box)
{
  int numprim = m_gdnvPhysics->numPrimitives();

  // This will hold 2nd or 4th order slopes
  FArrayBox dW(a_box,numprim);

  if (m_useFourthOrderSlopes)
    {
      // 2nd order slopes need to be computed over a larger box to accommodate
      // the 4th order slope computation
      Box boxVL = a_box;
      boxVL.grow(a_dir,1);
      boxVL &= m_domain;

      // Compute 2nd order (van Leer) slopes
      FArrayBox dWvL(boxVL, numprim);
      m_util.vanLeerSlopes(dWvL,a_W,numprim,
                           m_useCharLimiting || m_usePrimLimiting,
                           a_dir,boxVL);
      m_gdnvPhysics->getPhysIBC()->setBdrySlopes(dWvL,a_W,a_dir,m_currentTime);

      // Compute 4th order slopes, without limiting.
      m_util.fourthOrderSlopes(dW,a_W,dWvL, numprim, a_dir,a_box);
    }
  else
    {
      // Compute 2nd order (van Leer) slopes
      m_util.vanLeerSlopes(dW,a_W,numprim,
                           m_useCharLimiting || m_usePrimLimiting,
                           a_dir,a_box);
      m_gdnvPhysics->getPhysIBC()->setBdrySlopes(dW,a_W,a_dir,m_currentTime);
    }

  // To save on storage, we use the input values as temporaries for the
  // delta's
  a_WMinus.setVal(0.0);
  a_WPlus .setVal(0.0);

  if (m_useCharLimiting || m_useFourthOrderSlopes)
    {
      // Compute one-sided differences as inputs for limiting.
      m_util.oneSidedDifferences(a_WMinus,a_WPlus,a_W,a_dir,a_box);
    }

  FArrayBox lambda(a_box, numprim);
  m_gdnvPhysics->charValues(lambda, a_W, a_dir,a_box);

  if (m_useCharLimiting && m_usePrimLimiting)
    {
      MayDay::Error("PatchGodunov::PLMNormalPred: Attempt to limit slopes in primitive AND characteristic coordinates - not implemented");
    }

  // Apply limiter on characteristic or primitive variables. Either
  // way, must end up with characteristic variables to pass to to the
  // normal predictor utility.  Currently, cannot do both.

  // If doing characteristic limiting then transform before limiting
  if (m_useCharLimiting)
    {
      // Transform from primitive to characteristic variables
      m_gdnvPhysics->charAnalysis(a_WMinus,a_W,a_dir,a_box);
      m_gdnvPhysics->charAnalysis(a_WPlus ,a_W,a_dir,a_box);
      m_gdnvPhysics->charAnalysis(dW      ,a_W,a_dir,a_box);
    }

  if (m_useCharLimiting || m_usePrimLimiting)
    {
      // Limiting is already done for 2nd order slopes in primitive variables
      // so don't do it again
      if (m_useCharLimiting || m_useFourthOrderSlopes)
        {
          // Do slope limiting
          m_util.slopeLimiter(dW,a_WMinus,a_WPlus,numprim,a_box);
        }

      // Do slope flattening
      if (m_useFlattening)
        {
          m_util.applyFlattening(dW,a_flat,a_box);
        }
    }

  // If not doing characteristic limiting then transform after any limiting
  if (!m_useCharLimiting)
    {
      // Transform from primitive to characteristic variables
      m_gdnvPhysics->charAnalysis(dW,a_W,a_dir,a_box);
    }

  // To the normal prediction in characteristic variables
  m_util.PLMNormalPred(a_WMinus,a_WPlus,dW,lambda,a_dt / a_dx,a_box);

//.........这里部分代码省略.........

// Compute the time-centered values of the primitive variable on the
// interior cell faces of the cell-centered box a_box.
void PatchGodunov::computeWHalf(FluxBox&         a_WHalf,
                                const FArrayBox& a_U,
                                const FArrayBox& a_S,
                                const Real&      a_dt,
                                const Box&       a_box)
{
  CH_assert(isDefined());
  CH_assert(a_box == m_currentBox);

  // Get the number of various variables
  int numPrim = m_gdnvPhysics->numPrimitives();

  // The current box of valid cells
  Box curBox = m_currentBox;

  // Boxes for face-centered state - used for the riemann() and
  // artificialViscosity() calls
  Box faceBox[SpaceDim];

  // Boxes for face-centered fluxes
  Box fluxBox[SpaceDim];

  // Boxes for cell-centered state - used for the updatePrim() calls
  Box ccBox[SpaceDim];

  for (int dir1 = 0; dir1 < SpaceDim; ++dir1)
    {
      // faceBox[dir1] is face-centered in direction "dir1",
      // is valid "curBox" grown by 1 in all directions except "dir1",
      // and stays one cell away, in "dir1", from the domain boundary.
      faceBox[dir1] = curBox; // valid cell-centered box
      faceBox[dir1].grow(1);
      faceBox[dir1] &= m_domain;
      faceBox[dir1].grow(dir1, -1);
      faceBox[dir1].surroundingNodes(dir1);

      // ccBox[dir1] is cell-centered,
      // is valid "curBox" grown by 1 in all directions except "dir1",
      // but only within the domain, "m_domain".
      ccBox[dir1] = curBox;
      ccBox[dir1].grow(1);
      ccBox[dir1].grow(dir1, -1);
      ccBox[dir1] &= m_domain;

      // fluxBox[dir1] is face-centered in direction "dir1",
      // consisting of all those faces of cells of "ccBox[dir1]".
      fluxBox[dir1] = ccBox[dir1];
      fluxBox[dir1].surroundingNodes(dir1);

      // The difference between faceBox[dir1] and fluxBox[dir1] is:
      // if curBox abuts the boundary of the domain in direction "dir1",
      // then fluxBox[dir1] contains faces along that boundary,
      // but faceBox[dir1] does not.
    }

  // cell-centered box:  but restrict it to within domain
  Box WBox = a_U.box();
  WBox &= m_domain;

  // Primitive variables
  FArrayBox W(WBox,numPrim);

  // Calculate the primitive variables W from the conserved variables a_U,
  // on cell-centered WBox.
  m_gdnvPhysics->consToPrim(W, a_U, WBox);

  // slopeBox is the cell-centered box where slopes will be needed:
  // it is one larger than the final update box "curBox" of valid cells,
  // but within domain.
  // On slopeBox we define the FABs flattening, WMinus, and WPlus.
  Box slopeBox = curBox;
  slopeBox.grow(1);
  slopeBox &= m_domain;

  // Compute flattening once for all slopes if needed
  FArrayBox flattening(slopeBox, 1); // cell-centered
  if (m_useFlattening)
    {
      Interval velInt    = m_gdnvPhysics->velocityInterval();
      int pressureIndex  = m_gdnvPhysics->pressureIndex();
      Real smallPressure = m_gdnvPhysics->smallPressure();
      int bulkIndex      = m_gdnvPhysics->bulkModulusIndex();

      m_util.computeFlattening(flattening,
                               W,
                               velInt,
                               pressureIndex,
                               smallPressure,
                               bulkIndex,
                               slopeBox);
    }

  // Intermediate, extrapolated primitive variables
  FArrayBox WMinus[SpaceDim];
  FArrayBox WPlus [SpaceDim];

  // Initial fluxes
  FArrayBox WHalf1[SpaceDim];
//.........这里部分代码省略.........

// Adjust boundary fluxes to account for artificial viscosity
void RampIBC::artViscBC(FArrayBox&       a_F,
                        const FArrayBox& a_U,
                        const FArrayBox& a_divVel,
                        const int&       a_dir,
                        const Real&      a_time)
{
  CH_assert(m_isFortranCommonSet == true);
  CH_assert(m_isDefined == true);

  FArrayBox &divVel = (FArrayBox &)a_divVel;

  Box fluxBox = divVel.box();
  fluxBox.enclosedCells(a_dir);
  fluxBox.grow(a_dir,1);

  Box loBox,hiBox,centerBox,entireBox;
  int hasLo,hasHi;
  loHiCenterFace(loBox,hasLo,hiBox,hasHi,centerBox,entireBox,
                 fluxBox,m_domain,a_dir);

  if (hasLo == 1)
    {
      loBox.shiftHalf(a_dir,1);
      divVel.shiftHalf(a_dir,1);
      a_F.shiftHalf(a_dir,1);

      int loHiSide = -1;

      FORT_RAMPARTVISCF(CHF_FRA(a_F),
                        CHF_CONST_FRA(a_U),
                        CHF_CONST_FRA1(divVel,0),
                        CHF_CONST_INT(loHiSide),
                        CHF_CONST_REAL(m_dx),
                        CHF_CONST_INT(a_dir),
                        CHF_BOX(loBox));

      loBox.shiftHalf(a_dir,-1);
      divVel.shiftHalf(a_dir,-1);
      a_F.shiftHalf(a_dir,-1);
    }

  if (hasHi == 1)
    {
      hiBox.shiftHalf(a_dir,-1);
      divVel.shiftHalf(a_dir,-1);
      a_F.shiftHalf(a_dir,-1);

      int loHiSide = 1;

      FORT_RAMPARTVISCF(CHF_FRA(a_F),
                        CHF_CONST_FRA(a_U),
                        CHF_CONST_FRA1(divVel,0),
                        CHF_CONST_INT(loHiSide),
                        CHF_CONST_REAL(m_dx),
                        CHF_CONST_INT(a_dir),
                        CHF_BOX(hiBox));

      hiBox.shiftHalf(a_dir,1);
      divVel.shiftHalf(a_dir,1);
      a_F.shiftHalf(a_dir,1);
    }
}