REMORA_Tagging.cpp

Go to the documentation of this file.

#include <REMORA.H>
#include <REMORA_Derive.H>
 
using namespace amrex;
 
namespace {
/**
 * Copy mf's topmost and bottommost valid z-planes into its z-ghost planes.
 *
 * Several of the fields ErrorEst tags on live in MultiFabs with no ghost cells in z
 * (zvel_new, vec_mskr3d), and the vorticity branch only ever computes valid cells, so
 * mf's k = klo-1 and k = khi+1 planes would otherwise be read uninitialized by the GRAD
 * (adjacent_difference_greater) test, which unconditionally differences k+-1. There is
 * no data below k=0 or above k=N to copy, so impose a zero vertical gradient: GRAD then
 * sees no vertical difference at the surface and bottom cells instead of garbage.
 *
 * @param[inout] mf   single-component MultiFab with at least one ghost cell in z
 */
void
fill_z_ghost_planes (MultiFab& mf)
{
    AMREX_ALWAYS_ASSERT(mf.nComp() == 1 && mf.nGrowVect()[2] >= 1);
 
    // Not tiled: each iteration writes the two z-planes of a whole column, so tiling in
    // z would have several tiles writing the same cell.
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(mf); mfi.isValid(); ++mfi)
    {
        const Box& vbx = mfi.validbox();
        const int klo = vbx.smallEnd(2);
        const int khi = vbx.bigEnd(2);
 
        // Laterally grown as well, so the ghost columns get their z-planes too --
        // GRAD reads i+-1 and j+-1 in those columns as well as k+-1.
        const Box& gbx = mfi.fabbox();
        auto arr = mf.array(mfi);
 
        ParallelFor(makeSlab(gbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
        {
            arr(i,j,klo-1) = arr(i,j,klo);
            arr(i,j,khi+1) = arr(i,j,khi);
        });
    }
}
 
/**
 * Fill mf's lateral ghost cells from the nearest valid cell of the same box.
 *
 * Call this *before* FillBoundary: FillBoundary reads only valid regions, so it
 * overwrites every ghost cell backed by a real same-level or periodic neighbor and
 * leaves only the physical-boundary and coarse/fine ghosts holding the zero-gradient
 * value written here. That keeps the GRAD test's i+-1 / j+-1 reads defined and
 * deterministic everywhere without tagging on a jump we invented.
 *
 * @param[inout] mf   single-component MultiFab
 */
void
fill_lateral_ghosts_zero_grad (MultiFab& mf)
{
    AMREX_ALWAYS_ASSERT(mf.nComp() == 1);
    const IntVect ng = mf.nGrowVect();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& vbx = mfi.validbox();
        const Box  gbx = mfi.growntilebox(IntVect(ng[0],ng[1],0));
        auto arr = mf.array(mfi);
 
        const int ilo = vbx.smallEnd(0); const int ihi = vbx.bigEnd(0);
        const int jlo = vbx.smallEnd(1); const int jhi = vbx.bigEnd(1);
 
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            if (i < ilo || i > ihi || j < jlo || j > jhi) {
                arr(i,j,k) = arr(amrex::min(amrex::max(i,ilo),ihi),
                                 amrex::min(amrex::max(j,jlo),jhi), k);
            }
        });
    }
}
} // namespace
 
/**
 * Apply this criterion to tba, skipping values the land/sea mask makes meaningless.
 * The parameters are documented on the declaration in REMORA_ErrorTag.H.
 *
 * The one thing to know when reading the body: a cell this declines to test is left exactly
 * as it was found, not cleared. That is what lets a static box keep the tags it set over a
 * coast, and it is the whole difference from the derefine criteria this replaced.
 */
void
REMORAErrorTag::operator() (TagBoxArray&    tba,
                            const MultiFab* mf,
                            const MultiFab* mskr3d,
                            char            clearval,
                            char            tagval,
                            Real            time,
                            int             level,
                            const Geometry& geom,
                            const IntVect&  mask_lo,
                            const IntVect&  mask_hi) const
{
    BL_PROFILE("REMORAErrorTag::operator()");
 
    // The tests that read a field cell by cell, and so have something for the mask to guard.
    // RELGRAD and VORT do too, and are absent only because refinement_criteria_setup cannot
    // build them -- so assert rather than let a future one slip silently past the guard.
    const bool masked_test = (m_test == GRAD || m_test == LESS || m_test == GREATER);
 
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mskr3d == nullptr || masked_test ||
                                     (m_test != RELGRAD && m_test != VORT),
                                     "REMORAErrorTag: RELGRAD and VORT have no mask guard; "
                                     "give them one before building them from inputs");
 
    if (mskr3d == nullptr || !masked_test) {
        amrex::AMRErrorTag::operator()(tba, mf, clearval, tagval, time, level, geom);
        return;
    }
 
    AMREX_ALWAYS_ASSERT(mf != nullptr);
 
    // One box index runs the tags, the field and the mask below, so all three have to be
    // distributed the same way over the same boxes.
    AMREX_ALWAYS_ASSERT(mskr3d->boxArray()        == mf->boxArray() &&
                        mskr3d->DistributionMap() == mf->DistributionMap());
 
    AMREX_ALWAYS_ASSERT(mask_lo[0] <= 0 && mask_lo[1] <= 0 &&
                        mask_hi[0] >= 0 && mask_hi[1] >= 0);
 
    // Rejected rather than silently ignored: the mask is column-constant, so a vertical offset
    // could not mean anything, and for the same reason nothing below reads the mask at k+-1 --
    // which is why mskr3d needs no vertical ghost cells.
    AMREX_ALWAYS_ASSERT(mask_lo[2] == 0 && mask_hi[2] == 0);
 
    // Furthest the loop below reaches laterally: the field's own dependence on mskr, plus one
    // for GRAD, which asks the same of the neighbor it differences against.
    const int grad_reach = (m_test == GRAD) ? 1 : 0;
    for (int d = 0; d < 2; ++d) {
        AMREX_ALWAYS_ASSERT(mskr3d->nGrowVect()[d] >=
                            amrex::max(-mask_lo[d], mask_hi[d]) + grad_reach);
    }
 
    if ( (level <  0                 ) || (level >= m_info.m_max_level) ||
         (time  <  m_info.m_min_time ) || (time  >  m_info.m_max_time ) ) {
        return;
    }
 
    auto const& tagma = tba.arrays();
    auto const& datma = mf->const_arrays();
    auto const& mskma = mskr3d->const_arrays();
 
    auto const nvalues = std::ssize(m_value);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nvalues > 0,
                                     "Threshold values not properly set in REMORAErrorTag");
    auto const threshold = m_value[(level < nvalues) ? level : nvalues-1];
    auto const tag_update = m_info.m_derefine ? clearval : tagval;
    auto const volume_weighting = m_info.m_volume_weighting;
    auto const& geomdata = geom.data();
    auto const test = m_test;
    auto const mlo = mask_lo;
    auto const mhi = mask_hi;
 
    ParallelFor(tba, [=] AMREX_GPU_DEVICE (int bi, int i, int j, int k) noexcept
    {
        auto const& msk = mskma[bi];
 
        // Whether the value at index (ii,jj,k) was computed from water alone. A value that read
        // the land side reports on the coast whichever cell index it is filed under, so this
        // asks about every rho-cell the value depends on, not just its own.
        auto value_is_clean = [=] (int ii, int jj) noexcept
        {
            for     (int jo = mlo[1]; jo <= mhi[1]; ++jo) {
                for (int io = mlo[0]; io <= mhi[0]; ++io) {
                    if (msk(ii+io,jj+jo,k) <= Real(0.5)) { return false; }
                }
            }
            return true;
        };
 
        if (!value_is_clean(i,j)) { return; }
 
        auto const& dat = datma[bi];
        bool tag_it = false;
 
        if (test == GRAD)
        {
            // A difference only between two clean values, as AMReX's EB form of this test
            // takes one only across a connected face.
            Real ax = Real(0.0);
            if (value_is_clean(i+1,j)) {
                ax = amrex::max(ax, std::abs(dat(i+1,j,k) - dat(i,j,k)));
            }
            if (value_is_clean(i-1,j)) {
                ax = amrex::max(ax, std::abs(dat(i,j,k) - dat(i-1,j,k)));
            }
 
            Real ay = Real(0.0);
            if (value_is_clean(i,j+1)) {
                ay = amrex::max(ay, std::abs(dat(i,j+1,k) - dat(i,j,k)));
            }
            if (value_is_clean(i,j-1)) {
                ay = amrex::max(ay, std::abs(dat(i,j,k) - dat(i,j-1,k)));
            }
 
            // The mask is constant down a column, so both vertical neighbors of a water cell
            // are water and neither face needs a guard.
            Real az = std::abs(dat(i,j,k+1) - dat(i,j,k));
            az = amrex::max(az, std::abs(dat(i,j,k) - dat(i,j,k-1)));
 
            tag_it = (amrex::max(ax,ay,az) >= threshold);
        }
        else
        {
            const Real vol = volume_weighting
                           ? Geometry::Volume(IntVect{AMREX_D_DECL(i,j,k)}, geomdata)
                           : Real(1.0);
            tag_it = (test == LESS) ? (dat(i,j,k) * vol <= threshold)
                                    : (dat(i,j,k) * vol >= threshold);
        }
 
        if (tag_it) { tagma[bi](i,j,k) = tag_update; }
    });
    Gpu::streamSynchronize();
}
 
/**
 * Function to tag cells for refinement -- this overrides the pure virtual function in AmrCore
 *
 * @param[in]  levc    level of refinement (0 is coarsest level)
 * @param[out] tags    array of tagged cells
 * @param[in]  time    current time
 * @param[in]  ngrow   number of grow cells
*/
void
REMORA::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
{
    const int clearval = TagBox::CLEAR;
    const int   tagval = TagBox::SET;
 
    //
    // This mf must have ghost cells because we may take differences between adjacent values
    //
    std::unique_ptr<MultiFab> mf = std::make_unique<MultiFab>(grids[levc], dmap[levc], 1, 1);
 
    // Any cell-centered tracer may drive refinement, named as it is elsewhere: "temp",
    // "salt", "tracer", "tracer_1", or a biology tracer such as "NO3".
    auto cons_comp_for_field = [this] (const std::string& field) {
        for (int icomp = 0; icomp < ncons; ++icomp) {
            if (cons_names[icomp] == field) { return icomp; }
        }
        return -1;
    };
 
    // Which tracers exist depends on runtime input -- "tracer" only when remora.nscalar > 0,
    // the biology names only with a biology model -- so "use a tracer name" sends a reader
    // hunting for a typo that is not there. Name the ones this run actually has.
    auto valid_field_names = [this] () {
        std::string names;
        for (int icomp = 0; icomp < ncons; ++icomp) { names += cons_names[icomp] + ", "; }
        names += "x_velocity, y_velocity, z_velocity, vorticity, mask";
#ifdef REMORA_USE_PARTICLES
        names += ", <particle>_count";
#endif
        return names;
    };
 
    for (int j=0; j < ref_tags.size(); ++j)
    {
        const int cons_comp = cons_comp_for_field(ref_tags[j].Field());
 
        // Which rho-cells the value this criterion reads depends on, relative to the index it
        // is stored at. Set alongside the fill rather than in a second switch on the field
        // name, so adding a field cannot leave it filled but unguarded. Zero means the value
        // is its own cell and nothing else.
        IntVect mask_lo = IntVect::TheZeroVector();
        IntVect mask_hi = IntVect::TheZeroVector();
 
        if (cons_comp >= 0) {
            FillPatch(levc, time, *cons_new[levc], cons_new, BCVars::cons_bc, BdyVars::t,
                0,true,false);
        }
        // This allows dynamic refinement based on the value of a tracer
        if (cons_comp >= 0)
        {
            // A tracer is masked in place by advance_3d_ml, so its own cell is the whole story.
            MultiFab::Copy(*mf,*cons_new[levc],cons_comp,0,1,1);
        } else if (ref_tags[j].Field() == "x_velocity") {
            FillPatch(levc, time, *xvel_new[levc], xvel_new, xvel_bc(), BdyVars::u,0,true,true);
            MultiFab::Copy(*mf,*xvel_new[levc],0,0,1,1);
            // u is stored at a cell index but lives on that cell's low-x face, and vert_mean_3d
            // multiplies it by msku(i,j) = mskr(i-1,j)*mskr(i,j). So u at a water cell whose
            // i-1 neighbor is land is an exact zero: a mask artifact, not slack water.
            mask_lo = IntVect(-1,0,0);
        } else if (ref_tags[j].Field() == "y_velocity") {
            FillPatch(levc, time, *yvel_new[levc], yvel_new, yvel_bc(), BdyVars::v,0,true,true);
            MultiFab::Copy(*mf,*yvel_new[levc],0,0,1,1);
            // As for u, with mskv(i,j) = mskr(i,j-1)*mskr(i,j).
            mask_lo = IntVect(0,-1,0);
        } else if (ref_tags[j].Field() == "z_velocity") {
            FillPatch(levc, time, *zvel_new[levc], zvel_new, zvel_bc(), BdyVars::null,0,true,true);
            // zvel_new has no ghost cells in z, so we can only ask the copy for lateral ones
            MultiFab::Copy(*mf,*zvel_new[levc],0,0,1,IntVect(1,1,0));
            fill_z_ghost_planes(*mf);
            // zvel_new is identically zero: nothing masks it, and nothing puts a computed
            // value in it -- the vertical velocity the model solves for lives in a scratch
            // array inside advance_3d. So its own cell is as good an answer as any. If it is
            // ever wired up, W is built from Huon and Hvom at i+1 and j+1, making its real
            // dependence the five-point cross.
        } else if (ref_tags[j].Field() == "vorticity") {
            // Fill the ghost cells of the face-based velocities -- including at
            // coarse/fine boundaries, which is what FillPatch's FillPatchTwoLevels
            // interpolates -- since the cell-centered velocities in those ghost cells
            // are read when computing vorticity below
            FillPatch(levc, time, *xvel_new[levc], xvel_new, xvel_bc(), BdyVars::u,0,true,true);
            FillPatch(levc, time, *yvel_new[levc], yvel_new, yvel_bc(), BdyVars::v,0,true,true);
            FillPatch(levc, time, *zvel_new[levc], zvel_new, zvel_bc(), BdyVars::null,0,true,true);
 
            // Vorticity needs no vertical neighbor, and zvel_new has no z ghost cells
            // to average from, so only grow laterally
            MultiFab mf_cc_vel(grids[levc],dmap[levc],3,IntVect(1,1,0));
            average_face_to_cellcenter(mf_cc_vel,0,
                                       Array<const MultiFab*,3>{xvel_new[levc],
                                                                yvel_new[levc],
                                                                zvel_new[levc]},
                                       IntVect(1,1,0));
            // Impose bc's at domain boundaries at all levels
            FillBdyCCVels(levc, mf_cc_vel);
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
            for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
            {
                const Box& bx = mfi.tilebox();
                auto& dfab = (*mf)[mfi];
                auto& sfab = mf_cc_vel[mfi];
                auto pm = vec_pm[levc]->const_array(mfi);
                auto pn = vec_pn[levc]->const_array(mfi);
                auto maskr = vec_mskr[levc]->const_array(mfi);
                derived::remora_dervort(bx, dfab, 0, 1, sfab, pm, pn, maskr, Geom(levc), time, nullptr, levc);
            } // mfi
 
            // remora_dervort only writes valid cells, so fill mf's own ghosts before the
            // tagging criteria difference across them. The zero-gradient pass goes first;
            // FillBoundary then replaces every ghost that has a real same-level or periodic
            // neighbor, leaving only the physical-boundary and coarse/fine ghosts extrapolated.
            fill_lateral_ghosts_zero_grad(*mf);
            mf->FillBoundary(geom[levc].periodicity());
            fill_z_ghost_planes(*mf);
 
            // remora_dervort differences the cell-centered velocities at i+-1 and j+-1 without
            // masking them (see the TODO there), and each is an average of two faces, so this
            // value depends on the whole 3x3 block around it. Narrow once the derive is masked.
            //
            // Untested, unlike the face velocities that DogboneAnalytic_MLdryface covers: the
            // only masked problem available is near-irrotational, so any assertion on it would
            // be pinned to roundoff. Covering it wants a masked case with shear along a coast.
            mask_lo = IntVect(-1,-1,0);
            mask_hi = IntVect( 1, 1,0);
 
        } else if (ref_tags[j].Field() == "mask") {
            // vec_mskr3d has no z ghost cells, so this copy leaves mf's top and bottom
            // ghost planes unwritten while the GRAD test differences in z. The mask is
            // constant down a column, so the zero vertical gradient imposed here is exact.
            MultiFab::Copy(*mf,*vec_mskr3d[levc],0,0,1,IntVect(1,1,0));
            fill_z_ghost_planes(*mf);
        } else if (ref_tags[j].Field() == "h") {
            // Bathymetry is 2D, so spread it down each column of mf. One lateral
            // ghost ring is enough for the GRAD test; vec_h holds NGROW+1. The value
            // is constant down a column, so the zero vertical gradient imposed on
            // the z ghost planes is exact.
            for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
            {
                const Box bx = mfi.growntilebox(IntVect(1,1,0));
                Array4<const Real> const& h   = vec_h[levc]->const_array(mfi);
                Array4<      Real> const& arr = mf->array(mfi);
                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    arr(i,j,k) = h(i,j,0);
                });
            }
            fill_z_ghost_planes(*mf);
#ifdef REMORA_USE_PARTICLES
        } else {
            //
            // This allows dynamic refinement based on the number of particles per cell
            //
            // Note that we must count all the particles in levels both at and above the current,
            //      since otherwise, e.g., if the particles are all at level 1, counting particles at
            //      level 0 will not trigger refinement when regridding so level 1 will disappear,
            //      then come back at the next regridding
            //
            const auto& particles_namelist( particleData.getNames() );
            mf->setVal(zero);
            bool matched_particle_count = false;
            for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++)
            {
                std::string tmp_string(particles_namelist[i]+"_count");
                IntVect rr = IntVect::TheUnitVector();
                if (ref_tags[j].Field() == tmp_string) {
                    matched_particle_count = true;
                    for (int lev = levc; lev <= finest_level; lev++)
                    {
                        MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0);
                        particleData[particles_namelist[i]]->IncrementWithTotal(temp_dat, lev);
 
                        MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0);
 
                        if (lev == levc) {
                            MultiFab::Copy(*mf, temp_dat, 0, 0, 1, 0);
                        } else {
                            for (int d = 0; d < AMREX_SPACEDIM; d++) {
                                rr[d] *= ref_ratio[levc][d];
                            }
                            average_down(temp_dat, temp_dat_crse, 0, 1, rr);
                            MultiFab::Add(*mf, temp_dat_crse, 0, 0, 1, 0);
                        }
                    }
                }
            }
 
            // A box-only indicator carries no field name and tags geometrically, so it
            // never reads mf; only a named field that matched nothing is an error.
            if (!matched_particle_count && !ref_tags[j].Field().empty()) {
                amrex::Abort("Unknown refinement field '" + ref_tags[j].Field() +
                             "'. This run has: " + valid_field_names() + ".");
            }
#else
        } else if (!ref_tags[j].Field().empty()) {
            // mf is uninitialized until something writes it, so an unrecognized named
            // field would otherwise tag on garbage. A box-only indicator has no field
            // name, tags geometrically, and never reads mf.
            amrex::Abort("Unknown refinement field '" + ref_tags[j].Field() +
                         "'. This run has: " + valid_field_names() + ".");
#endif
        }
 
        // Two criteria want the unguarded test. One keyed on the mask is asking where the
        // coast is: guarding it would leave no water-water face across which the mask varies,
        // so it would never tag. And with no mask every cell is water, so the guard could not
        // fire -- skipping it keeps an unmasked run on exactly AMReX's own code path.
        const bool unguarded = (ref_tags[j].Field() == "mask") ||
                               (solverChoice.mask_type == MaskType::none);
        const MultiFab* mskr3d_for_tag = unguarded ? nullptr : vec_mskr3d[levc].get();
 
        // time counts from remora.start_time, but the window is on the model clock. Checked
        // here, after the field is filled, so a misnamed field still aborts outside it.
        if (!ref_tags[j].ActiveAt(model_time(time))) {
            continue;
        }
 
        ref_tags[j](tags,mf.get(),mskr3d_for_tag,clearval,tagval,time,levc,geom[levc],
                    mask_lo,mask_hi);
    }
 
    // Promote any tagged cell to a full local z-column.
    for (MFIter mfi(tags, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.validbox();
        auto const& tag = tags.array(mfi);
 
        const int klo = bx.smallEnd(2);
        const int khi = bx.bigEnd(2);
 
        amrex::ParallelFor(makeSlab(bx, 2, 0),
        [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
        {
            bool refine_col = false;
            for (int k = klo; k <= khi; ++k) {
                refine_col = refine_col || (tag(i,j,k) != TagBox::CLEAR);
            }
 
            if (refine_col) {
                for (int k = klo; k <= khi; ++k) {
                    tag(i,j,k) = TagBox::SET;
                }
            }
        });
    }
}
 
/**
 * Function to define the refinement criteria based on user input
*/
void
REMORA::refinement_criteria_setup ()
{
    if (max_level > 0)
    {
        ParmParse pp(pp_prefix);
        Vector<std::string> refinement_indicators;
        pp.queryarr("refinement_indicators",refinement_indicators,0,pp.countval("refinement_indicators"));
        for (int i=0; i<refinement_indicators.size(); ++i)
        {
            std::string ref_prefix = pp_prefix + "." + refinement_indicators[i];
 
            ParmParse ppr(ref_prefix);
            RealBox realbox;
            int lev_for_box;
 
            int num_real_lo      = ppr.countval("in_box_lo");
            int num_indx_lo      = ppr.countval("in_box_lo_indices");
            int num_indx_lo_crse = ppr.countval("in_box_lo_indices_crse");
 
            int num_real_hi      = ppr.countval("in_box_hi");
            int num_indx_hi      = ppr.countval("in_box_hi_indices");
            int num_indx_hi_crse = ppr.countval("in_box_hi_indices_crse");
 
            AMREX_ALWAYS_ASSERT( (num_real_lo      == num_real_hi)      && (num_real_lo      == 0 || num_real_lo      >= 2) );
            AMREX_ALWAYS_ASSERT( (num_indx_lo      == num_indx_hi)      && (num_indx_lo      == 0 || num_indx_lo      >= 2) );
            AMREX_ALWAYS_ASSERT( (num_indx_lo_crse == num_indx_hi_crse) && (num_indx_lo_crse == 0 || num_indx_lo_crse >= 2) );
 
            // Problem low and high (in real not index space) are the same at all levels
            if ( !((num_real_lo >= AMREX_SPACEDIM-1 && num_indx_lo == 0 && num_indx_lo_crse == 0) ||
                   (num_indx_lo >= AMREX_SPACEDIM-1 && num_real_lo == 0 && num_indx_lo_crse == 0) ||
                   (num_indx_lo ==              0   && num_real_lo == 0 && num_indx_lo_crse == 0) ||
                   (num_indx_lo_crse >= AMREX_SPACEDIM-1 && num_real_lo == 0 && num_indx_lo == 0)
                ) )
            {
                amrex::Abort("Must only specify box for refinement using real OR index space with fine/coarse grid indices");
            }
 
            if (num_real_lo > 0) {
                std::vector<Real> box_lo(3), box_hi(3);
                ppr.get("max_level",lev_for_box);
                if (lev_for_box > 0 && lev_for_box <= max_level)
                {
                    ppr.getarr("in_box_lo",box_lo,0,2);
                    ppr.getarr("in_box_hi",box_hi,0,2);
                    box_lo[2] = geom[0].ProbLo(2);
                    box_hi[2] = geom[0].ProbHi(2);
                    realbox = RealBox(&(box_lo[0]),&(box_hi[0]));
 
                    amrex::Print() << "Reading " << realbox << " at level " << lev_for_box << std::endl;
                    num_boxes_at_level[lev_for_box] += 1;
 
                    const auto* dx  = geom[lev_for_box].CellSize();
                    const Real* plo = geom[lev_for_box].ProbLo();
 
                    int ilo = static_cast<int>((box_lo[0] - plo[0])/dx[0]);
                    int jlo = static_cast<int>((box_lo[1] - plo[1])/dx[1]);
                    int klo = static_cast<int>((box_lo[2] - plo[2])/dx[2]);
                    int ihi = static_cast<int>((box_hi[0] - plo[0])/dx[0]-1);
                    int jhi = static_cast<int>((box_hi[1] - plo[1])/dx[1]-1);
                    int khi = static_cast<int>((box_hi[2] - plo[2])/dx[2]-1);
 
                    Box bx_old(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi));
 
                    int mod_ilo = ilo%ref_ratio[lev_for_box-1][0];
                    int mod_jlo = jlo%ref_ratio[lev_for_box-1][1];
 
                    int mod_ihi = (ihi+1)%ref_ratio[lev_for_box-1][0];
                    int mod_jhi = (jhi+1)%ref_ratio[lev_for_box-1][1];
 
                    if (mod_ilo != 0) {
                        ilo -= mod_ilo;
                    }
                    if (mod_jlo != 0) {
                        jlo -= mod_jlo;
                    }
                    if (mod_ihi != 0) {
                        ihi += ref_ratio[lev_for_box-1][0] - mod_ihi;
                    }
                    if (mod_jhi != 0) {
                        jhi += ref_ratio[lev_for_box-1][1] - mod_jhi;
                    }
                    Box bx(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi));
                    if (mod_ilo !=0 || mod_jlo !=0 || mod_ihi != 0 || mod_jhi != 0) {
                        amrex::Print() << "Fine box on level " << lev_for_box << " adjusted from " << bx_old << " to " << bx << " to make it valid for refinement." << std::endl;
                    }
                    boxes_at_level[lev_for_box].push_back(bx);
                    amrex::Print() << "Saving in 'boxes at level' as " << bx << std::endl;
                } // lev
 
            } else if (num_indx_lo > 0) {
 
                std::vector<int> box_lo(3), box_hi(3);
                ppr.get("max_level",lev_for_box);
                if (lev_for_box > 0 && lev_for_box <= max_level)
                {
                    if (n_error_buf[0] != IntVect::TheZeroVector()) {
                        amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
                    }
 
                    ppr.getarr("in_box_lo_indices",box_lo,0,num_indx_lo);
                    ppr.getarr("in_box_hi_indices",box_hi,0,num_indx_hi);
 
                    if (num_indx_lo < AMREX_SPACEDIM) {
                        box_lo[2] = geom[lev_for_box].Domain().smallEnd(2);
                        box_hi[2] = geom[lev_for_box].Domain().bigEnd(2);
                    }
 
                    Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2]));
                    const Box& domain = geom[lev_for_box].Domain();
 
                    if (!domain.contains(bx)) {
                        amrex::Print() << "\n";
                        amrex::Print() << "Box specified       is " << bx << std::endl;
                        amrex::Print() << "But domain at level is " << domain << std::endl;
                        amrex::Error("Specified box doesn't fit in the domain");
                    }
 
                    const auto* dx  = geom[lev_for_box].CellSize();
                    const Real* plo = geom[lev_for_box].ProbLo();
                    realbox = RealBox(plo[0]+ box_lo[0]   *dx[0], plo[1]+ box_lo[1]   *dx[1], plo[2]+ box_lo[2]   *dx[2],
                                      plo[0]+(box_hi[0]+1)*dx[0], plo[1]+(box_hi[1]+1)*dx[1], plo[2]+(box_hi[2]+1)*dx[2]);
 
                    Print() << "Reading " << bx << " at level " << lev_for_box << std::endl;
                    num_boxes_at_level[lev_for_box] += 1;
 
                    if(box_lo[0]%ref_ratio[lev_for_box-1][0] != 0){
                        amrex::Print()<< "Requested ilo in x-direction : " << box_lo[0] << std::endl;
                        amrex::Print() << "ilo = " << box_lo[0] << " is not divisible by ref_ratio in x direction = " <<
                                          ref_ratio[lev_for_box-1][0] << std::endl;
                        amrex::Error("Adjust in_box_lo_indices in x-direction to be divisible by ref_ratio and try again");
                    }
                    if((box_hi[0]+1)%ref_ratio[lev_for_box-1][0] != 0){
                        amrex::Print()<< "Requested ihi in x-direction : " << box_hi[0] << std::endl;
                        amrex::Print() << "ihi+1 = " << box_hi[0]+1 << " is not divisible by ref_ratio in x direction = " <<
                                          ref_ratio[lev_for_box-1][0] << std::endl;
                        amrex::Error("Adjust in_box_hi_indices in x-direction to be divisible by ref_ratio and try again");
                    }
                     if(box_lo[1]%ref_ratio[lev_for_box-1][1] != 0){
                        amrex::Print()<< "Requested jlo in y-direction : " << box_lo[1] << std::endl;
                        amrex::Print() << "jlo = " << box_lo[1] << " is not divisible by ref_ratio in y direction = " <<
                                          ref_ratio[lev_for_box-1][1] << std::endl;
                        amrex::Error("Adjust in_box_lo_indices in y-direction to be divisible by ref_ratio and try again");
                    }
                    if((box_hi[1]+1)%ref_ratio[lev_for_box-1][1] != 0){
                        amrex::Print()<< "Requested jhi in y-direction : " << box_hi[1] << std::endl;
                        amrex::Print() << "jhi+1 = " << box_hi[1]+1 << " is not divisible by ref_ratio in y direction = " <<
                                          ref_ratio[lev_for_box-1][1] << std::endl;
                        amrex::Error("Adjust in_box_hi_indices in y-direction to be divisible by ref_ratio and try again");
                    }
                    if(box_lo[2]%ref_ratio[lev_for_box-1][2] != 0){
                        amrex::Print()<< "Requested klo in z-direction : " << box_lo[2] << std::endl;
                        amrex::Print() << "klo = " << box_lo[2] << " is not   divisible by ref_ratio in z direction = " <<
                                          ref_ratio[lev_for_box-1][2] << std::endl;
                        amrex::Error("Adjust in_box_lo_indices in z-direction to be divisible by ref_ratio and try again");
                    }
                    if((box_hi[2]+1)%ref_ratio[lev_for_box-1][2] != 0){
                        amrex::Print()<< "Requested khi in z-direction : " << box_hi[2] << std::endl;
                        amrex::Print() << "khi+1 = " << box_hi[2]+1 << " is not divisible by ref_ratio in z direction = " <<
                                          ref_ratio[lev_for_box-1][2] << std::endl;
                        amrex::Error("Adjust in_box_hi_indices in z-direction to be divisible by ref_ratio and try again");
                    }
 
                    boxes_at_level[lev_for_box].push_back(bx);
                    Print() << "Saving in 'boxes at level' as " << bx << std::endl;
                } // lev
 
            } else if (num_indx_lo_crse > 0) {
 
                std::vector<int> box_lo(3), box_hi(3);
                ppr.get("max_level",lev_for_box);
                if (lev_for_box > 0 && lev_for_box <= max_level)
                {
                    if (n_error_buf[0] != IntVect::TheZeroVector()) {
                        amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
                    }
 
                    ppr.getarr("in_box_lo_indices_crse",box_lo,0,num_indx_lo_crse);
                    ppr.getarr("in_box_hi_indices_crse",box_hi,0,num_indx_hi_crse);
 
                    if (num_indx_lo_crse < AMREX_SPACEDIM) {
                        box_lo[2] = geom[lev_for_box-1].Domain().smallEnd(2);
                        box_hi[2] = geom[lev_for_box-1].Domain().bigEnd(2);
                    }
 
                    Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2]));
 
                    if (!geom[lev_for_box-1].Domain().contains(bx)) {
                        amrex::Print() << "\n";
                        amrex::Print() << "(Coarse) Box specified       is " << bx << std::endl;
                        amrex::Print() << "But (coarse) domain at level is " << geom[lev_for_box-1].Domain() << std::endl;
                        amrex::Error("Specified box doesn't fit in the domain");
                    }
 
                    bx.refine(ref_ratio[lev_for_box-1]);
 
                    const auto* dx  = geom[lev_for_box-1].CellSize();
 
                    const Real* plo = geom[lev_for_box].ProbLo();
                    realbox = RealBox(plo[0]+ box_lo[0]   *dx[0], plo[1]+ box_lo[1]   *dx[1], plo[2]+ box_lo[2]   *dx[2],
                                      plo[0]+(box_hi[0]+1)*dx[0], plo[1]+(box_hi[1]+1)*dx[1], plo[2]+(box_hi[2]+1)*dx[2]);
 
                    Print() << "Reading " << bx << " at level " << lev_for_box << std::endl;
                    num_boxes_at_level[lev_for_box] += 1;
 
                    boxes_at_level[lev_for_box].push_back(bx);
                    Print() << "Saving in 'boxes at level' as " << bx << std::endl;
                } // lev
            }
 
            AMRErrorTagInfo info;
 
            if (realbox.ok()) {
                info.SetRealBox(realbox);
            }
            // The window is on the model clock. It stays out of info, whose times are Real
            // and are compared against elapsed time; ErrorEst applies it instead.
            double ref_min_time = std::numeric_limits<double>::lowest();
            double ref_max_time = std::numeric_limits<double>::max();
            ppr.query("start_time",ref_min_time);
            ppr.query("end_time",ref_max_time);
            if (ppr.countval("max_level") > 0) {
                int ref_max_level; ppr.get("max_level",ref_max_level);
                info.SetMaxLevel(ref_max_level);
            }
 
            if (ppr.countval("value_greater")) {
                int num_val = ppr.countval("value_greater");
                Vector<Real> value(num_val);
                ppr.getarr("value_greater",value,0,num_val);
                std::string field; ppr.get("field_name",field);
                ref_tags.push_back(REMORAErrorTag(value,AMRErrorTag::GREATER,field,info));
            }
            else if (ppr.countval("value_less")) {
                int num_val = ppr.countval("value_less");
                Vector<Real> value(num_val);
                ppr.getarr("value_less",value,0,num_val);
                std::string field; ppr.get("field_name",field);
                ref_tags.push_back(REMORAErrorTag(value,AMRErrorTag::LESS,field,info));
            }
            else if (ppr.countval("adjacent_difference_greater")) {
                int num_val = ppr.countval("adjacent_difference_greater");
                Vector<Real> value(num_val);
                ppr.getarr("adjacent_difference_greater",value,0,num_val);
                std::string field; ppr.get("field_name",field);
                ref_tags.push_back(REMORAErrorTag(value,AMRErrorTag::GRAD,field,info));
            }
            else if (realbox.ok())
            {
                ref_tags.push_back(REMORAErrorTag(info));
            } else {
                Abort(std::string("Unrecognized refinement indicator for " + refinement_indicators[i]).c_str());
            }
            ref_tags.back().SetModelTimeWindow(ref_min_time, ref_max_time);
        } // loop over criteria
    } // if max_level > 0
}