C++ BoxArray类代码示例

MultiFab*
Nyx::build_fine_mask()
{
    BL_ASSERT(level > 0); // because we are building a mask for the coarser level

    if (fine_mask != 0) return fine_mask;

    BoxArray baf = parent->boxArray(level);
    baf.coarsen(crse_ratio);

    const BoxArray& bac = parent->boxArray(level-1);
    fine_mask = new MultiFab(bac,parent->DistributionMap(level-1), 1,0);
    fine_mask->setVal(1.0);

#ifdef _OPENMP
#pragma omp parallel
#endif
    for (MFIter mfi(*fine_mask); mfi.isValid(); ++mfi)
    {
        FArrayBox& fab = (*fine_mask)[mfi];

        std::vector< std::pair<int,Box> > isects = baf.intersections(fab.box());

        for (int ii = 0; ii < isects.size(); ii++)
        {
            fab.setVal(0.0,isects[ii].second,0);
        }
    }

    return fine_mask;
}

//
// What's the slowest way I can think of to compute all the norms??
//
Real
MFNorm (const MultiFab& mfab, 
        const int       exponent,
        const int       srcComp,
        const int       numComp,
        const int       numGrow)
{
    BL_ASSERT (numGrow <= mfab.nGrow());
    BoxArray boxes = mfab.boxArray();
    boxes.grow(numGrow);
    //
    // Get a copy of the multifab
    //
    MultiFab mftmp(mfab.boxArray(), numComp, 0);
    MultiFab::Copy(mftmp,mfab,srcComp,0,numComp,numGrow);
    //
    // Calculate the Norms
    //
    Real myNorm = 0;
    if ( exponent == 0 )
    {
        for ( MFIter mftmpmfi(mftmp); mftmpmfi.isValid(); ++mftmpmfi)
        {
            mftmp[mftmpmfi].abs(boxes[mftmpmfi.index()], 0, numComp);
            myNorm = std::max(myNorm, mftmp[mftmpmfi].norm(0, 0, numComp));
        }
	ParallelDescriptor::ReduceRealMax(myNorm);

    } else if ( exponent == 1 )
    {
        for ( MFIter mftmpmfi(mftmp); mftmpmfi.isValid(); ++mftmpmfi)
        {
            mftmp[mftmpmfi].abs(boxes[mftmpmfi.index()], 0, numComp);

            myNorm += mftmp[mftmpmfi].norm(1, 0, numComp);
        }
	ParallelDescriptor::ReduceRealSum(myNorm);

    } else if ( exponent == 2 )
    {
        for ( MFIter mftmpmfi(mftmp); mftmpmfi.isValid(); ++mftmpmfi)
        {
            mftmp[mftmpmfi].abs(boxes[mftmpmfi.index()], 0, numComp);

            myNorm += pow(mftmp[mftmpmfi].norm(2, 0, numComp), 2);
        }
	ParallelDescriptor::ReduceRealSum(myNorm);
        myNorm = sqrt( myNorm );

    } else {

        BoxLib::Error("Invalid exponent to norm function");
    }
    
    return myNorm;
}

BoxList::BoxList (const BoxArray &ba)
    :
    lbox(),
    btype()
{
    if (ba.size() > 0)
        btype = ba[0].ixType();
    for (int i = 0, N = ba.size(); i < N; ++i)
        push_back(ba[i]);
}

int
iMultiFab::norm0 (int comp, const BoxArray& ba, bool local) const
{
    int nm0 = -std::numeric_limits<int>::max();

#ifdef _OPENMP
#pragma omp parallel
#endif
    {
	std::vector< std::pair<int,Box> > isects;
	int priv_nm0 = -std::numeric_limits<int>::max();

	for (MFIter mfi(*this); mfi.isValid(); ++mfi)
	{
	    ba.intersections(mfi.validbox(),isects);
	    
	    for (int i = 0, N = isects.size(); i < N; i++)
	    {
		priv_nm0 = std::max(priv_nm0, get(mfi).norm(isects[i].second, 0, comp, 1));
	    }
	}
#ifdef _OPENMP
#pragma omp critical (imultifab_norm0_ba)
#endif
	{
	    nm0 = std::max(nm0, priv_nm0);
	}
    }
 
    if (!local)
	ParallelDescriptor::ReduceIntMax(nm0);
 
    return nm0;
}

MultiFab_C_to_F::MultiFab_C_to_F (const Geometry& geom,
				  const DistributionMapping& dmap,
				  const BoxArray& ba)
	  
{
    BL_ASSERT(count == 0);
    count++;

    int nb = ba.size();
    int dm = BL_SPACEDIM;

    std::vector<int> lo(nb*dm);
    std::vector<int> hi(nb*dm);

    for ( int i = 0; i < nb; ++i ) {
	const Box& bx = BoxLib::enclosedCells(ba[i]);
        for ( int j = 0; j < dm; ++j ) {
	    lo[j + i*dm] = bx.smallEnd(j);
	    hi[j + i*dm] = bx.bigEnd(j);
	}
    }

    const Box& domain = geom.Domain();

    int pm[dm];
    for ( int i = 0; i < dm; ++i ) {
	pm[i] = geom.isPeriodic(i)? 1 : 0;
    }

    const Array<int>& pmap = dmap.ProcessorMap();

    build_layout_from_c(nb, dm, &lo[0], &hi[0], 
			domain.loVect(), domain.hiVect(), 
			pm, pmap.dataPtr());
}

MCMultiGrid::MCMultiGrid (MCLinOp &_lp)
    :
    initialsolution(0),
    Lp(_lp)
{
    Initialize();

    maxiter = def_maxiter;
    numiter = def_numiter;
    nu_0 = def_nu_0;
    nu_1 = def_nu_1;
    nu_2 = def_nu_2;
    nu_f = def_nu_f;
    usecg = def_usecg;
    verbose = def_verbose;
    rtol_b = def_rtol_b;
    atol_b = def_atol_b;
    nu_b = def_nu_b;
    numLevelsMAX = def_numLevelsMAX;
    numlevels = numLevels();
    numcomps = _lp.numberComponents();
    if ( ParallelDescriptor::IOProcessor() && (verbose > 2) )
    {
	BoxArray tmp = Lp.boxArray();
	std::cout << "MCMultiGrid: numlevels = " << numlevels 
		  << ": ngrid = " << tmp.size() << ", npts = [";
	for ( int i = 0; i < numlevels; ++i ) 
        {
	    if ( i > 0 ) tmp.coarsen(2);
	    std::cout << tmp.d_numPts() << " ";
        }
	std::cout << "]" << '\n';

	std::cout << "MCMultiGrid: " << numlevels
	     << " multigrid levels created for this solve" << '\n';
    }

    if ( ParallelDescriptor::IOProcessor() && (verbose > 4) )
    {
	std::cout << "Grids: " << '\n';
	BoxArray tmp = Lp.boxArray();
	for (int i = 0; i < numlevels; ++i)
	{
            Orientation face(0, Orientation::low);
            const DistributionMapping& map = Lp.bndryData().bndryValues(face).DistributionMap();
	    if (i > 0)
		tmp.coarsen(2);
	    std::cout << " Level: " << i << '\n';
	    for (int k = 0; k < tmp.size(); k++)
	    {
		const Box& b = tmp[k];
		std::cout << "  [" << k << "]: " << b << "   ";
		for (int j = 0; j < BL_SPACEDIM; j++)
		    std::cout << b.length(j) << ' ';
                std::cout << ":: " << map[k] << '\n';
	    }
	}
    }
}

void
BndryRegister::setBoxes (const BoxArray& _grids)
{
    BL_ASSERT(grids.size() == 0);
    BL_ASSERT(_grids.size() > 0);
    BL_ASSERT(_grids[0].cellCentered());

    grids.define(_grids);
    //
    // Check that bndry regions are not allocated.
    //
    for (int k = 0; k < 2*BL_SPACEDIM; k++)
        BL_ASSERT(bndry[k].size() == 0);
}

bool
ParticleBase::FineToCrse (const ParticleBase&                p,
                          int                                flev,
                          const ParGDBBase*                  gdb,
                          const Array<IntVect>&              fcells,
                          const BoxArray&                    fvalid,
                          const BoxArray&                    compfvalid_grown,
                          Array<IntVect>&                    ccells,
                          Array<Real>&                       cfracs,
                          Array<int>&                        which,
                          Array<int>&                        cgrid,
                          Array<IntVect>&                    pshifts,
                          std::vector< std::pair<int,Box> >& isects)
{
    BL_ASSERT(gdb != 0);
    BL_ASSERT(flev > 0);
    //
    // We're in AssignDensity(). We want to know whether or not updating
    // with a particle we'll cross a fine->crse boundary.  Note that crossing
    // a periodic boundary, where the periodic shift lies in our valid region,
    // is not considered a Fine->Crse crossing.
    //
    const int M = fcells.size();

    which.resize(M);
    cgrid.resize(M);
    ccells.resize(M);
    cfracs.resize(M);

    for (int i = 0; i < M; i++)
    {
        cgrid[i] = -1;
        which[i] =  0;
    }

    const Box& ibx = BoxLib::grow(gdb->ParticleBoxArray(flev)[p.m_grid],-1);

    BL_ASSERT(ibx.ok());

    if (ibx.contains(p.m_cell))
        //
        // We're strictly contained in our valid box.
        // We can't cross a fine->crse boundary.
        //
        return false;

    if (!compfvalid_grown.contains(p.m_cell))
        //
        // We're strictly contained in our "valid" region. Note that the valid
        // region contains any periodically shifted ghost cells that intersect
        // valid region.
        //
        return false;
    //
    // Otherwise ...
    //
    const Geometry& cgm = gdb->Geom(flev-1);
    const IntVect&  rr  = gdb->refRatio(flev-1);
    const BoxArray& cba = gdb->ParticleBoxArray(flev-1);

    ParticleBase::CIC_Cells_Fracs(p, cgm.ProbLo(), cgm.CellSize(), cfracs, ccells);

    bool result = false;

    for (int i = 0; i < M; i++)
    {
        IntVect ccell_refined = ccells[i]*rr;
        //
        // We've got to protect against the case when we're at the low
        // end of the domain because coarsening & refining don't work right
        // when indices go negative.
        //
        for (int dm = 0; dm < BL_SPACEDIM; dm++)
            ccell_refined[dm] = std::max(ccell_refined[dm], -1);

        if (!fvalid.contains(ccell_refined))
        {
            result   = true;
            which[i] = 1;

            Box cbx(ccells[i],ccells[i]);
    
            if (!cgm.Domain().contains(ccells[i]))
            {
                //
                // We must be at a periodic boundary.
                // Find valid box into which we can be periodically shifted.
                //
                BL_ASSERT(cgm.isAnyPeriodic());

                cgm.periodicShift(cbx, cgm.Domain(), pshifts);

                BL_ASSERT(pshifts.size() == 1);

                cbx -= pshifts[0];

                ccells[i] -= pshifts[0];
                BL_ASSERT(cbx.ok());
                BL_ASSERT(cgm.Domain().contains(cbx));
            }
//.........这里部分代码省略.........

void
BndryRegister::defineDoit (Orientation _face,
                           IndexType   _typ,
                           int         _in_rad,
                           int         _out_rad,
                           int         _extent_rad,
                           BoxArray&   fsBA)
{
    BL_PROFILE("BndryRegister::defineDoit()");

    BL_ASSERT(grids.size() > 0);

    const int coord_dir = _face.coordDir();
    const int lo_side   = _face.isLow();
    //
    // Build the BoxArray on which to define the FabSet on this face.
    //
    const int N = grids.size();

    fsBA.resize(N);

#ifdef _OPENMP
#pragma omp parallel for
#endif
    for (int idx = 0; idx < N; ++idx)
    {
        Box b;
        //
        // First construct proper box for direction normal to face.
        //
        if (_out_rad > 0)
        {
            if (_typ.ixType(coord_dir) == IndexType::CELL)
                b = BoxLib::adjCell(grids[idx], _face, _out_rad);
            else
                b = BoxLib::bdryNode(grids[idx], _face, _out_rad);

            if (_in_rad > 0)
                b.grow(_face.flip(), _in_rad);
        }
        else
        {
            if (_in_rad > 0)
            {
                if (_typ.ixType(coord_dir) == IndexType::CELL)
                    b = BoxLib::adjCell(grids[idx], _face, _in_rad);
                else
                    b = BoxLib::bdryNode(grids[idx], _face, _in_rad);

                b.shift(coord_dir, lo_side ? _in_rad : -_in_rad);
            }
            else
                BoxLib::Error("BndryRegister::define(): strange values for in_rad, out_rad");
        }
        //
        // Now alter box in all other index directions.
        //
        for (int dir = 0; dir < BL_SPACEDIM; dir++)
        {
            if (dir == coord_dir)
                continue;
            if (_typ.ixType(dir) == IndexType::NODE)
                b.surroundingNodes(dir);
            if (_extent_rad > 0)
                b.grow(dir,_extent_rad);
        }

        BL_ASSERT(b.ok());

        fsBA.set(idx,b);
    }

    BL_ASSERT(fsBA.ok());
}

//.........这里部分代码省略.........
		mw.dataPtr(),ARLIM(mw.loVect()),ARLIM(mw.hiVect()),
		ms.dataPtr(),ARLIM(ms.loVect()),ARLIM(ms.hiVect()),
		mt.dataPtr(),ARLIM(mt.loVect()),ARLIM(mt.hiVect()),
		mb.dataPtr(),ARLIM(mb.loVect()),ARLIM(mb.hiVect()),
		denfab.dataPtr(), 
		ARLIM(denfab.loVect()), ARLIM(denfab.hiVect()),
		exttdptr, ARLIM(fslo), ARLIM(fshi),
		tdfab.dataPtr(),ARLIM(tdfab.loVect()),ARLIM(tdfab.hiVect()),
		inout.box(gn).loVect(), inout.box(gn).hiVect(),
		&nc, h[level]);
#endif
	}
    }

#if 0
  // This "probably" works, but is not strictly needed just because of the way Bill
  // coded up the tangential derivative stuff.  It's handy code though, so I want to
  // keep it around/

  // Clean up corners:
  // The problem here is that APPLYBC fills only grow cells normal to the boundary.
  // As a result, any corner cell on the boundary (either coarse-fine or fine-fine)
  // is not filled.  For coarse-fine, the operator adjusts itself, sliding away from
  // the box edge to avoid referencing that corner point.  On the physical boundary
  // though, the corner point is needed.  Particularly if a fine-fine boundary intersects
  // the physical boundary, since we want the stencil to be independent of the box
  // blocking.  FillBoundary operations wont fix the problem because the "good"
  // data we need is living in the grow region of adjacent fabs.  So, here we play
  // the usual games to treat the newly filled grow cells as "valid" data.

  // Note that we only need to do something where the grids touch the physical boundary.

  const Geometry& geomlev = geomarray[level];
  const BoxArray& grids = inout.boxArray();
  const Box& domain = geomlev.Domain();
  int nGrow = 1;
  int src_comp = 0;
  int num_comp = BL_SPACEDIM;


  // Lets do a quick check to see if we need to do anything at all here
  BoxArray BIGba = BoxArray(grids).grow(nGrow);

  if (! (domain.contains(BIGba.minimalBox())) ) {

    BoxArray boundary_pieces;
    Array<int> proc_idxs;
    Array<Array<int> > old_to_new(grids.size());
    const DistributionMapping& dmap=inout.DistributionMap();

    for (int d=0; d<BL_SPACEDIM; ++d) {
      if (! (geomlev.isPeriodic(d)) ) {

        BoxArray gba = BoxArray(grids).grow(d,nGrow);
        for (int i=0; i<gba.size(); ++i) {
          BoxArray new_pieces = BoxLib::boxComplement(gba[i],domain);
          int size_new = new_pieces.size();
          if (size_new>0) {
            int size_old = boundary_pieces.size();
            boundary_pieces.resize(size_old+size_new);
            proc_idxs.resize(boundary_pieces.size());
            for (int j=0; j<size_new; ++j) {
              boundary_pieces.set(size_old+j,new_pieces[j]);
              proc_idxs[size_old+j] = dmap[i];
              old_to_new[i].push_back(size_old+j);
            }

Real
Adv::advance (Real time,
              Real dt,
              int  iteration,
              int  ncycle)
{
    for (int k = 0; k < NUM_STATE_TYPE; k++) {
        state[k].allocOldData();
        state[k].swapTimeLevels(dt);
    }

    MultiFab& S_new = get_new_data(State_Type);

    const Real prev_time = state[State_Type].prevTime();
    const Real cur_time = state[State_Type].curTime();
    const Real ctr_time = 0.5*(prev_time + cur_time);

    const Real* dx = geom.CellSize();
    const Real* prob_lo = geom.ProbLo();

    //
    // Get pointers to Flux registers, or set pointer to zero if not there.
    //
    FluxRegister *fine    = 0;
    FluxRegister *current = 0;
    
    int finest_level = parent->finestLevel();

    if (do_reflux && level < finest_level) {
	fine = &getFluxReg(level+1);
	fine->setVal(0.0);
    }

    if (do_reflux && level > 0) {
	current = &getFluxReg(level);
    }

    MultiFab fluxes[BL_SPACEDIM];

    if (do_reflux)
    {
	for (int j = 0; j < BL_SPACEDIM; j++)
	{
	    BoxArray ba = S_new.boxArray();
	    ba.surroundingNodes(j);
	    fluxes[j].define(ba, NUM_STATE, 0, Fab_allocate);
	}
    }

    // State with ghost cells
    MultiFab Sborder(grids, NUM_STATE, NUM_GROW);
    FillPatch(*this, Sborder, NUM_GROW, time, State_Type, 0, NUM_STATE);

#ifdef _OPENMP
#pragma omp parallel
#endif
    {
	FArrayBox flux[BL_SPACEDIM], uface[BL_SPACEDIM];

	for (MFIter mfi(S_new, true); mfi.isValid(); ++mfi)
	{
	    const Box& bx = mfi.tilebox();

	    const FArrayBox& statein = Sborder[mfi];
	    FArrayBox& stateout      =   S_new[mfi];

	    // Allocate fabs for fluxes and Godunov velocities.
	    for (int i = 0; i < BL_SPACEDIM ; i++) {
		const Box& bxtmp = BoxLib::surroundingNodes(bx,i);
		flux[i].resize(bxtmp,NUM_STATE);
		uface[i].resize(BoxLib::grow(bxtmp,1),1);
	    }

	    get_face_velocity(level, ctr_time,
			      D_DECL(BL_TO_FORTRAN(uface[0]),
				     BL_TO_FORTRAN(uface[1]),
				     BL_TO_FORTRAN(uface[2])),
			      dx, prob_lo);

            advect(time, bx.loVect(), bx.hiVect(),
		   BL_TO_FORTRAN_3D(statein), 
		   BL_TO_FORTRAN_3D(stateout),
		   D_DECL(BL_TO_FORTRAN_3D(uface[0]),
			  BL_TO_FORTRAN_3D(uface[1]),
			  BL_TO_FORTRAN_3D(uface[2])),
		   D_DECL(BL_TO_FORTRAN_3D(flux[0]), 
			  BL_TO_FORTRAN_3D(flux[1]), 
			  BL_TO_FORTRAN_3D(flux[2])), 
		   dx, dt);

	    if (do_reflux) {
		for (int i = 0; i < BL_SPACEDIM ; i++)
		    fluxes[i][mfi].copy(flux[i],mfi.nodaltilebox(i));	  
	    }
	}
    }

    if (do_reflux) {
	if (current) {
	    for (int i = 0; i < BL_SPACEDIM ; i++)
//.........这里部分代码省略.........

bool
BoxList::contains (const BoxArray&  ba) const
{
    BoxArray tba(*this);
    return ba.contains(tba);
}

static
BoxArray
GetBndryCells (const BoxArray& ba,
               int             ngrow,
               const Geometry& geom)
{
    //
    // First get list of all ghost cells.
    //
    BoxList gcells, bcells;

    for (int i = 0; i < ba.size(); ++i)
	gcells.join(BoxLib::boxDiff(BoxLib::grow(ba[i],ngrow),ba[i]));
    //
    // Now strip out intersections with original BoxArray.
    //
    for (BoxList::const_iterator it = gcells.begin(); it != gcells.end(); ++it)
    {
        std::vector< std::pair<int,Box> > isects = ba.intersections(*it);

        if (isects.empty())
            bcells.push_back(*it);
        else
        {
            //
            // Collect all the intersection pieces.
            //
            BoxList pieces;
            for (int i = 0; i < isects.size(); i++)
                pieces.push_back(isects[i].second);
            BoxList leftover = BoxLib::complementIn(*it,pieces);
            bcells.catenate(leftover);
        }
    }
    //
    // Now strip out overlaps.
    //
    gcells.clear();
    gcells = BoxLib::removeOverlap(bcells);
    bcells.clear();

    if (geom.isAnyPeriodic())
    {
        Array<IntVect> pshifts(27);

        const Box& domain = geom.Domain();

        for (BoxList::const_iterator it = gcells.begin(); it != gcells.end(); ++it)
        {
            if (!domain.contains(*it))
            {
                //
                // Add in periodic ghost cells shifted to valid region.
                //
                geom.periodicShift(domain, *it, pshifts);

                for (int i = 0; i < pshifts.size(); i++)
                {
                    const Box& shftbox = *it + pshifts[i];

                    const Box& ovlp = domain & shftbox;
                    BoxList bl = BoxLib::complementIn(ovlp,BoxList(ba));
                    bcells.catenate(bl);
                }
            }
        }

        gcells.catenate(bcells);
    }

    return BoxArray(gcells);
}

int
main (int   argc,
      char* argv[])
{
    BoxLib::Initialize(argc,argv);

    if (argc < 2)
        print_usage(argc,argv);

    ParmParse pp;

    if (pp.contains("help"))
        print_usage(argc,argv);

    bool verbose = false;
    pp.query("verbose",verbose);

    if (verbose>1)
        AmrData::SetVerbose(true);

    std::string infile; pp.get("infile",infile);
    std::string outfile_DEF;

    std::string outType = "tec";
#ifdef USE_TEC_BIN_IO
    bool doBin = true;
    pp.query("doBin",doBin);
    outfile_DEF = infile+(doBin ? ".plt" : ".dat" );
#else
    bool doBin=false;
    outfile_DEF = infile+".dat";
#endif

    // 
    bool connect_cc = true; pp.query("connect_cc",connect_cc);

    std::string outfile(outfile_DEF); pp.query("outfile",outfile);
    DataServices::SetBatchMode();
    Amrvis::FileType fileType(Amrvis::NEWPLT);
    DataServices dataServices(infile, fileType);

    if (!dataServices.AmrDataOk())
        DataServices::Dispatch(DataServices::ExitRequest, NULL);

    AmrData& amrData = dataServices.AmrDataRef();

    const Array<std::string>& names = amrData.PlotVarNames();

    Array<int> comps;
    if (int nc = pp.countval("comps"))
    {
        comps.resize(nc);
        pp.getarr("comps",comps,0,nc);
    }
    else
    {
        int sComp = 0;
        pp.query("sComp",sComp);
        int nComp = amrData.NComp();
        pp.query("nComp",nComp);
        BL_ASSERT(sComp+nComp <= amrData.NComp());
        comps.resize(nComp);
        for (int i=0; i<nComp; ++i)
            comps[i] = sComp + i;
    }

    Box subbox;
    if (int nx=pp.countval("box"))
    {
        Array<int> barr;
        pp.getarr("box",barr,0,nx);
        int d=BL_SPACEDIM;
        BL_ASSERT(barr.size()==2*d);
        subbox=Box(IntVect(D_DECL(barr[0],barr[1],barr[2])),
                   IntVect(D_DECL(barr[d],barr[d+1],barr[d+2]))) & amrData.ProbDomain()[0];

    } else {

        subbox = amrData.ProbDomain()[0];
    }

    int finestLevel = amrData.FinestLevel(); pp.query("finestLevel",finestLevel);
    int Nlev = finestLevel + 1;
    Array<BoxArray> gridArray(Nlev);
    Array<Box> subboxArray(Nlev);

    int nGrowPer = 0; pp.query("nGrowPer",nGrowPer);
    PArray<Geometry> geom(Nlev);

    Array<Real> LO(BL_SPACEDIM,0);
    Array<Real> HI(BL_SPACEDIM,1);
    RealBox rb(LO.dataPtr(),HI.dataPtr());
    int coordSys = 0;
    Array<int> isPer(BL_SPACEDIM,0);

    for (int lev=0; lev<Nlev; ++lev)
    {
        subboxArray[lev]
            = (lev==0 ? subbox
               : Box(subboxArray[lev-1]).refine(amrData.RefRatio()[lev-1]));
//.........这里部分代码省略.........

void main_main ()
{
    // What time is it now?  We'll use this to compute total run time.
    Real strt_time = ParallelDescriptor::second();

    // AMREX_SPACEDIM: number of dimensions
    int n_cell, max_grid_size;
    Vector<int> bc_lo(AMREX_SPACEDIM,0);
    Vector<int> bc_hi(AMREX_SPACEDIM,0);

    // inputs parameters
    {
        // ParmParse is way of reading inputs from the inputs file
        ParmParse pp;

        // We need to get n_cell from the inputs file - this is the
        // number of cells on each side of a square (or cubic) domain.
        pp.get("n_cell", n_cell);

        // The domain is broken into boxes of size max_grid_size
        max_grid_size = 32;
        pp.query("max_grid_size", max_grid_size);

    }

    Vector<int> is_periodic(AMREX_SPACEDIM,0);
    for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
      is_periodic[idim] = 1;
    }

    // make BoxArray and Geometry
    BoxArray ba;
    Geometry geom;
    {
        IntVect dom_lo(AMREX_D_DECL(       0,        0,        0));
        IntVect dom_hi(AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1));
        Box domain(dom_lo, dom_hi);

        // Initialize the boxarray "ba" from the single box "bx"
        ba.define(domain);

        // Break up boxarray "ba" into chunks no larger than
        // "max_grid_size" along a direction
        ba.maxSize(max_grid_size);

        // This defines the physical box, [0, 1] in each direction.
        RealBox real_box({AMREX_D_DECL(0.0, 0.0, 0.0)},
                         {AMREX_D_DECL(1.0, 1.0, 1.0)});

        // This defines a Geometry object
        geom.define(domain, &real_box,
                    CoordSys::cartesian, is_periodic.data());
    }

    // Nghost = number of ghost cells for each array
    int Nghost = 0;

    // do the runtime parameter initializations and microphysics inits
    if (ParallelDescriptor::IOProcessor()) {
      std::cout << "reading extern runtime parameters ..." << std::endl;
    }

    ParmParse ppa("amr");

    std::string probin_file = "probin";

    ppa.query("probin_file", probin_file);

    const int probin_file_length = probin_file.length();
    Vector<int> probin_file_name(probin_file_length);

    for (int i = 0; i < probin_file_length; i++)
      probin_file_name[i] = probin_file[i];

    init_unit_test(probin_file_name.dataPtr(), &probin_file_length);

    // Ncomp = number of components for each array
    int Ncomp = -1;
    init_variables();
    get_ncomp(&Ncomp);

    int name_len = -1;
    get_name_len(&name_len);

    // get the variable names
    Vector<std::string> varnames;

    for (int i=0; i<Ncomp; i++) {
      char* cstring[name_len+1];
      get_var_name(cstring, &i);
      std::string name(*cstring);
      varnames.push_back(name);
    }

    // time = starting time in the simulation
    Real time = 0.0;

    // How Boxes are distrubuted among MPI processes
    DistributionMapping dm(ba);

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