REMORA_FillPatch.cpp

Go to the documentation of this file.

#include <REMORA.H>
#include <REMORA_PhysBCFunct.H>
#include <REMORA_IndexDefines.H>
#include <REMORA_TimeInterpolatedData.H>
 
using namespace amrex;
 
PhysBCFunctNoOp null_bc;
 
/**
 * Fill valid and ghost data in the MultiFab "mf"
 * This version fills the MultiFab mf_to_fill in valid regions with
 * the "state data" at the given time;
 * values in mf when it is passed in are *not* used.
 *
 * @param[in   ] lev            level to FillPatch on
 * @param[in   ] time           current time
 * @param[inout] mf_to_fill     MultiFab to fill
 * @param[in   ] mfs            vector over levels of multifabs associated with mf_to_fill
 * @param[in   ] bccomp         index into domain_bcs_type
 * @param[in   ] bdy_var_type   when filling from NetCDF, which boundary data
 * @param[in   ] icomp          component to fill
 * @param[in   ] fill_all       whether to fill all components
 * @param[in   ] fill_set       whether to use FillPatchers
 * @param[in   ] icomp_calc     component to use in RHS boundary calculation
 * @param[in   ] dt_lev         time step at this level
 * @param[in   ] mf_calc        data to use in RHS boundary calculation
 */
void
REMORA::FillPatch (int lev, Real time, MultiFab& mf_to_fill, Vector<MultiFab*> const& mfs,
                  const int  bccomp,
#ifdef REMORA_USE_NETCDF
                  const int bdy_var_type,
#else
                  const int /*bdy_var_type*/,
#endif
                  const int  icomp,
                  const bool fill_all,
                  const bool fill_set,
                  const int  n_not_fill,
#ifdef REMORA_USE_NETCDF
                  const int  icomp_calc,
                  const Real dt_lev,
#else
                  const int /*icomp_calc*/,
                  const Real /*dt_lev*/,
#endif
                  const MultiFab& mf_calc,
                  Vector<MultiFab*> const& mfs_crse_old,
                  Vector<MultiFab*> const& mfs_crse_new)
{
    BL_PROFILE_VAR("REMORA::FillPatch()",REMORA_FillPatch);
    amrex::Interpolater* mapper = nullptr;
 
    Box mf_box(mf_to_fill.boxArray()[0]);
    bool is_2d = mf_box.length(2) == 1;
 
    MultiFab* mask;
 
    //
    // ***************************************************************************
    // The first thing we do is interpolate the momenta on the "valid" faces of
    // the fine grids (where the interface is coarse/fine not fine/fine) -- this
    // will not be over-written below because the FillPatch operators see these as
    // valid faces.
    // ***************************************************************************
    if (lev>0 && fill_set) {
        if (cf_set_width > 0 &&
            mf_box.ixType() == IndexType(IntVect(0,0,0))) {
            FPr_c[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
        } else if (cf_set_width >= 0) {
            if (!is_2d) {
                if (mf_box.ixType() == IndexType(IntVect(1,0,0))) {
                    FPr_u[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                } else if (mf_box.ixType() == IndexType(IntVect(0,1,0))) {
                    FPr_v[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                } else if (mf_box.ixType() == IndexType(IntVect(0,0,1))) {
                    FPr_w[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                }
            } else {
                if (mf_box.ixType() == IndexType(IntVect(1,0,0))) {
                    FPr_ubar[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                } else if (mf_box.ixType() == IndexType(IntVect(0,1,0))) {
                    FPr_vbar[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                }
            }
        }
    }
 
    int ncomp;
    if (fill_all) {
        ncomp = mf_to_fill.nComp();
    } else {
        ncomp = 1;
    }
 
    if (mf_box.ixType() == IndexType(IntVect(0,0,0)))
    {
        mapper = &cell_cons_interp;
        mask   = vec_mskr[lev].get();
    }
    else if (mf_box.ixType() == IndexType(IntVect(1,0,0)))
    {
        mapper = &face_cons_linear_interp;
        mask   = vec_msku[lev].get();
    }
    else if (mf_box.ixType() == IndexType(IntVect(0,1,0)))
    {
        mapper = &face_cons_linear_interp;
        mask   = vec_mskv[lev].get();
    }
    else {
        mapper = &face_cons_linear_interp;
        mask   = vec_mskr[lev].get();
    }
 
    if (lev == 0)
    {
        Vector<MultiFab*> fmf = {mfs[lev], mfs[lev]};
        Vector<Real> ftime    = {t_old[lev], t_new[lev]};
        amrex::FillPatchSingleLevel(mf_to_fill, time, fmf, ftime, icomp, icomp, ncomp,
                                    geom[lev], null_bc, bccomp);
    }
    else
    {
        Vector<MultiFab*> fmf = {mfs[lev], mfs[lev]};
        Vector<Real> ftime    = {t_old[lev], t_new[lev]};
        // A subcycled fine level asks for times inside the parent's step, so the coarse
        // contribution has to be interpolated. Passing one MultiFab twice makes
        // FillPatchTwoLevels return it whatever the time: fine in lockstep, wrong for
        // subcycling.
        Vector<MultiFab*> cmf = {mfs[lev-1], mfs[lev-1]};
        Vector<Real> ctime    = {time, time};
        if (do_substep && int(mfs_crse_old.size()) >= lev && int(mfs_crse_new.size()) >= lev) {
            cmf   = {mfs_crse_old[lev-1], mfs_crse_new[lev-1]};
            ctime = {t_old[lev-1], t_new[lev-1]};
            mfs_crse_old[lev-1]->FillBoundary(geom[lev-1].periodicity());
            mfs_crse_new[lev-1]->FillBoundary(geom[lev-1].periodicity());
        }
 
        mfs[lev-1]->FillBoundary(geom[lev-1].periodicity());
        amrex::FillPatchTwoLevels(mf_to_fill, mf_to_fill.nGrowVect(), IntVect(0,0,0),
                                  time, cmf, ctime, fmf, ftime,
                                  icomp, icomp, ncomp, geom[lev-1], geom[lev],
                                  //null_bc, bccomp, null_bc, bccomp,
                                  refRatio(lev-1),
                                  mapper, domain_bcs_type, bccomp);
    } // lev > 0
 
    // This is currently unconditionally true, but do_bc flag should actually be passed in
    const bool do_bc=true;
    if (do_bc) {
        // ***************************************************************************
        // Physical bc's at domain boundary
        // ***************************************************************************
 
        // Enforce physical boundary conditions
        (*physbcs[lev])(mf_to_fill,*mask,icomp,ncomp,mf_to_fill.nGrowVect(),time,bccomp,n_not_fill,mf_calc,
                *vec_msku[lev].get(), *vec_mskv[lev].get());
 
#ifdef REMORA_USE_NETCDF
        // Fill the data which is stored in the boundary data read from netcdf files
        if ( (solverChoice.boundary_from_netcdf) && (bdy_var_type != BdyVars::null) )
        {
            fill_from_bdyfiles(lev, mf_to_fill,*mask,time,bccomp,bdy_var_type, icomp,icomp_calc,mf_calc,dt_lev);
        }
#endif
        // Fill corners of the domain with periodic data
        if  ( mf_box.ixType() == IndexType(IntVect(0,0,0)) ) {
            mf_to_fill.EnforcePeriodicity(geom[lev].periodicity());
        }
 
        // Also enforce free-slip at top boundary (on xvel or yvel)
        if ( (mf_box.ixType() == IndexType(IntVect(1,0,0))) ||
             (mf_box.ixType() == IndexType(IntVect(0,1,0))) )
        {
            int khi = geom[lev].Domain().bigEnd(2);
            for (MFIter mfi(mf_to_fill); mfi.isValid(); ++mfi)
            {
                Box gbx  = mfi.growntilebox(); // Note this is face-centered since vel is
                gbx.setSmall(2,khi+1);
                if (gbx.ok()) {
                    Array4<Real> vel_arr = mf_to_fill.array(mfi);
                    ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                    {
                        vel_arr(i,j,k) = vel_arr(i,j,khi);
                    });
                }
            }
        }
    }
}
 
/**
 * Fill valid and ghost data in the MultiFab "mf"
 * This version fills the MultiFab mf_to_fill in valid regions with
 * the "state data" at the given time;
 * values in mf when it is passed in are *not* used.
 * Unlike FillPatch, FillPatchNoBC does not apply boundary conditions
 *
 * @param[in   ] lev            level to FillPatch on
 * @param[in   ] time           current time
 * @param[inout] mf_to_fill     MultiFab to fill
 * @param[in   ] mfs            vector over levels of multifabs associated with mf_to_fill
 * @param[in   ] bdy_var_type   when filling from NetCDF, which boundary data
 * @param[in   ] icomp          component to fill
 * @param[in   ] fill_all       whether to fill all components
 * @param[in   ] fill_set       whether to use FillPatchers
 */
void
REMORA::FillPatchNoBC (int lev, Real time, MultiFab& mf_to_fill, Vector<MultiFab*> const& mfs,
#ifdef REMORA_USE_NETCDF
                  const int /*bdy_var_type*/,
#else
                  const int /*bdy_var_type*/,
#endif
                  const int  icomp,
                  const bool fill_all,
                  const bool fill_set,
                  Vector<MultiFab*> const& mfs_crse_old,
                  Vector<MultiFab*> const& mfs_crse_new)
{
    // HACK: Note that this is hacky; should be able to have a single call to FillPatch with a
    // flag for bcs, but for some reason it was acting weird, so we're splitting this out into
    // two functions with repeated code for the time being
    BL_PROFILE_VAR("REMORA::FillPatchNoBC()",REMORA_FillPatch);
    int bccomp;
    amrex::Interpolater* mapper = nullptr;
 
    Box mf_box(mf_to_fill.boxArray()[0]);
    bool is_2d = mf_box.length(2) == 1;
 
    //
    // ***************************************************************************
    // The first thing we do is interpolate the momenta on the "valid" faces of
    // the fine grids (where the interface is coarse/fine not fine/fine) -- this
    // will not be over-written below because the FillPatch operators see these as
    // valid faces.
    // ***************************************************************************
    if (lev>0 && fill_set) {
        if (cf_set_width > 0 &&
            mf_box.ixType() == IndexType(IntVect(0,0,0))) {
            FPr_c[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
        } else if (cf_set_width >= 0) {
            if (!is_2d) {
                if (mf_box.ixType() == IndexType(IntVect(1,0,0))) {
                    FPr_u[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                } else if (mf_box.ixType() == IndexType(IntVect(0,1,0))) {
                    FPr_v[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                } else if (mf_box.ixType() == IndexType(IntVect(0,0,1))) {
                    FPr_w[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                }
            } else {
                if (mf_box.ixType() == IndexType(IntVect(1,0,0))) {
                    FPr_ubar[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                } else if (mf_box.ixType() == IndexType(IntVect(0,1,0))) {
                    FPr_vbar[lev-1].FillSet(mf_to_fill, time, null_bc, domain_bcs_type);
                }
            }
        }
    }
 
    int ncomp;
    if (fill_all) {
        ncomp = mf_to_fill.nComp();
    } else {
        ncomp = 1;
    }
 
    if (mf_box.ixType() == IndexType(IntVect(0,0,0)))
    {
        bccomp = 0;
        mapper = &cell_cons_interp;
    }
    else if (mf_box.ixType() == IndexType(IntVect(1,0,0)))
    {
        bccomp = xvel_bc();
        mapper = &face_cons_linear_interp;
    }
    else if (mf_box.ixType() == IndexType(IntVect(0,1,0)))
    {
        bccomp = yvel_bc();
        mapper = &face_cons_linear_interp;
    }
    else {
        bccomp = zvel_bc();
        mapper = &face_cons_linear_interp;
    }
 
    if (lev == 0)
    {
        Vector<MultiFab*> fmf = {mfs[lev], mfs[lev]};
        Vector<Real> ftime    = {t_old[lev], t_new[lev]};
        amrex::FillPatchSingleLevel(mf_to_fill, time, fmf, ftime, icomp, icomp, ncomp,
                                    geom[lev], null_bc, bccomp);
    }
    else
    {
        Vector<MultiFab*> fmf = {mfs[lev], mfs[lev]};
        Vector<Real> ftime    = {t_old[lev], t_new[lev]};
        // A subcycled fine level asks for times inside the parent's step, so the coarse
        // contribution has to be interpolated. Passing one MultiFab twice makes
        // FillPatchTwoLevels return it whatever the time: fine in lockstep, wrong for
        // subcycling.
        Vector<MultiFab*> cmf = {mfs[lev-1], mfs[lev-1]};
        Vector<Real> ctime    = {time, time};
        if (do_substep && int(mfs_crse_old.size()) >= lev && int(mfs_crse_new.size()) >= lev) {
            cmf   = {mfs_crse_old[lev-1], mfs_crse_new[lev-1]};
            ctime = {t_old[lev-1], t_new[lev-1]};
            mfs_crse_old[lev-1]->FillBoundary(geom[lev-1].periodicity());
            mfs_crse_new[lev-1]->FillBoundary(geom[lev-1].periodicity());
        }
 
        mfs[lev-1]->FillBoundary(geom[lev-1].periodicity());
        amrex::FillPatchTwoLevels(mf_to_fill, mf_to_fill.nGrowVect(), IntVect(0,0,0),
                                  time, cmf, ctime, fmf, ftime,
                                  icomp, icomp, ncomp, geom[lev-1], geom[lev],
                                  //null_bc, bccomp, null_bc, bccomp,
                                  refRatio(lev-1),
                                  mapper, domain_bcs_type, bccomp);
    } // lev > 0
}
 
 
// utility to copy in data from old/new data into a struct that holds data for FillPatching
TimeInterpolatedData
REMORA::GetDataAtTime (int /*lev*/, Real /*time*/)
{
    BL_PROFILE_VAR("GetDataAtTime()",GetDataAtTime);
    TimeInterpolatedData data;
 
// HACK HACK HACK
#if 0
    const Real teps = (t_new[lev] - t_old[lev]) * 1.e-3_rt;
 
    if (time > t_new[lev] - teps && time < t_new[lev] + teps)
    {
        for (int i = 0; i < Vars::NumTypes; ++i) {
            data.add_var(&vars_new[lev][i], data.non_owning);
        }
        data.set_time(t_new[lev]);
    }
    else if (time > t_old[lev] - teps && time < t_old[lev] + teps)
    {
        for (int i = 0; i < Vars::NumTypes; ++i) {
            data.add_var(&vars_old[lev][i], data.non_owning);
        }
        data.set_time(t_old[lev]);
    }
    else if (time > t_old[lev] && time < t_new[lev])
    {
        // do first order interpolation in time between [t_old[lev], t_new[lev]]
        // time interpolation includes the ghost cells
        for (int i = 0; i < Vars::NumTypes; ++i) {
            MultiFab* mf_tmp = new MultiFab(vars_new[lev][i].boxArray(),
                                            vars_new[lev][i].DistributionMap(),
                                            vars_new[lev][i].nComp(), vars_new[lev][i].nGrowVect());
            mf_tmp->setVal(0.0_rt);
 
            const Real dt_fraction = (time - t_old[lev]) / (t_new[lev] - t_old[lev]);
            MultiFab::Saxpy(*mf_tmp, 1.0_rt - dt_fraction, vars_old[lev][i], 0, 0, mf_tmp->nComp(), mf_tmp->nGrowVect());
            MultiFab::Saxpy(*mf_tmp,          dt_fraction, vars_new[lev][i], 0, 0, mf_tmp->nComp(), mf_tmp->nGrowVect());
 
            data.add_var(mf_tmp, data.owning);
        }
        data.set_time(time);
    }
    else
    {
        amrex::Error("Requested data at a time outside the interval [t_old, t_new]");
    }
 
    // We need to make sure to fill these before we compute the viscosity
    for (int i = 0; i < Vars::NumTypes; ++i) {
        data.get_var(Vars::xvel).FillBoundary(geom[lev].periodicity());
        data.get_var(Vars::yvel).FillBoundary(geom[lev].periodicity());
        data.get_var(Vars::zvel).FillBoundary(geom[lev].periodicity());
        data.get_var(Vars::cons).FillBoundary(geom[lev].periodicity());
    }
#endif
    return data;
}
 
 
/**
 * Like FillCoarsePatch but uses piecewise constant interpolater
 *
 * @param[in   ] lev            level to FillPatch on
 * @param[in   ] time           current time
 * @param[inout] mf_to_fill     MultiFab to fill
 * @param[inout] mf_crse        MultiFab of coarse data
 * @param[in   ] bccomp         index into domain_bcs_type
 * @param[in   ] bdy_var_type   when filling from NetCDF, which boundary data
 * @param[in   ] icomp          component to fill
 * @param[in   ] fill_all       whether to fill all components
 * @param[in   ] icomp_calc     component to use in RHS boundary calculation
 * @param[in   ] dt_lev         time step at this level
 * @param[in   ] mf_calc        data to use in RHS boundary calculation
 */
void
REMORA::FillCoarsePatchPC (int lev, Real time, MultiFab* mf_to_fill,
                           MultiFab* mf_crse,
                           const int bccomp,
                           const int bdy_var_type,
                           const int  icomp,
                           const bool fill_all,
                           const int  n_not_fill,
                           const int  icomp_calc,
                           const Real dt_lev,
                           const MultiFab& mf_calc) {
    amrex::Interpolater* mapper = &pc_interp;
    FillCoarsePatchMap(lev, time, mf_to_fill, mf_crse, bccomp, bdy_var_type,icomp,
                    fill_all, n_not_fill, icomp_calc, dt_lev, mf_calc, mapper);
}
 
/**
 * Fill an entire multifab by interpolating from the coarser level --
 * this is used only when a new level of refinement is being created
 * during a run (i.e not at initialization)
 * This will never be used with static refinement.
 *
 * @param[in   ] lev            level to FillPatch on
 * @param[in   ] time           current time
 * @param[inout] mf_to_fill     MultiFab to fill
 * @param[inout] mf_crse        MultiFab of coarse data
 * @param[in   ] bccomp         index into domain_bcs_type
 * @param[in   ] bdy_var_type   when filling from NetCDF, which boundary data
 * @param[in   ] icomp          component to fill
 * @param[in   ] fill_all       whether to fill all components
 * @param[in   ] icomp_calc     component to use in RHS boundary calculation
 * @param[in   ] dt_lev         time step at this level
 * @param[in   ] mf_calc        data to use in RHS boundary calculation
 */
void
REMORA::FillCoarsePatch (int lev, Real time, MultiFab* mf_to_fill,
                         MultiFab* mf_crse,
                         const int bccomp,
                         const int bdy_var_type,
                         const int  icomp,
                         const bool fill_all,
                         const int  n_not_fill,
                         const int  icomp_calc,
                         const Real dt_lev,
                         const MultiFab& mf_calc) {
    amrex::Interpolater* mapper = nullptr;
    FillCoarsePatchMap(lev, time, mf_to_fill, mf_crse, bccomp, bdy_var_type,icomp,
                    fill_all, n_not_fill, icomp_calc, dt_lev, mf_calc, mapper);
}
 
/**
 * Fill an entire multifab by interpolating from the coarser level --
 * this is used only when a new level of refinement is being created
 * during a run (i.e not at initialization)
 * This will never be used with static refinement.
 *
 * @param[in   ] lev            level to FillPatch on
 * @param[in   ] time           current time
 * @param[inout] mf_to_fill     MultiFab to fill
 * @param[inout] mf_crse        MultiFab of coarse data
 * @param[in   ] bccomp         index into domain_bcs_type
 * @param[in   ] bdy_var_type   when filling from NetCDF, which boundary data
 * @param[in   ] icomp          component to fill
 * @param[in   ] fill_all       whether to fill all components
 * @param[in   ] icomp_calc     component to use in RHS boundary calculation
 * @param[in   ] dt_lev         time step at this level
 * @param[in   ] mf_calc        data to use in RHS boundary calculation
 */
void
REMORA::FillCoarsePatchMap (int lev, Real time, MultiFab* mf_to_fill, MultiFab* mf_crse,
                         const int bccomp,
#ifdef REMORA_USE_NETCDF
                         const int bdy_var_type,
#else
                         const int /*bdy_var_type*/,
#endif
                         const int  icomp,
                         const bool fill_all,
                         const int  n_not_fill,
                         const int  icomp_calc,
                         const Real dt_lev,
                         const MultiFab& mf_calc,
                               amrex::Interpolater* mapper)
{
    BL_PROFILE_VAR("FillCoarsePatch()",FillCoarsePatch);
    AMREX_ASSERT(lev > 0);
 
    MultiFab* mask = nullptr;
 
    int ncomp;
    if (fill_all) {
        ncomp = mf_to_fill->nComp();
    } else {
        ncomp = 1;
    }
 
    Box box_mf(mf_to_fill->boxArray()[0]);
 
    if (box_mf.ixType() == IndexType(IntVect(0,0,0)))
    {
        mask = vec_mskr[lev].get();
    }
    else if (box_mf.ixType() == IndexType(IntVect(1,0,0)))
    {
        mask = vec_msku[lev].get();
    }
    else if (box_mf.ixType() == IndexType(IntVect(0,1,0)))
    {
        mask = vec_mskv[lev].get();
    }
    else if (box_mf.ixType() == IndexType(IntVect(0,0,1)))
    {
        mask = vec_mskr[lev].get();
    } else {
          amrex::Abort("Dont recognize this box type in REMORA_FillPatch");
    }
 
    if (mapper == nullptr) {
        if (box_mf.ixType() == IndexType(IntVect(0,0,0)))
        {
            mapper = &cell_cons_interp;
        }
        else if (box_mf.ixType() == IndexType(IntVect(1,0,0)))
        {
            mapper = &face_cons_linear_interp;
        }
        else if (box_mf.ixType() == IndexType(IntVect(0,1,0)))
        {
            mapper = &face_cons_linear_interp;
        }
        else if (box_mf.ixType() == IndexType(IntVect(0,0,1)))
        {
            mapper = &face_cons_linear_interp;
        } else {
              amrex::Abort("Dont recognize this box type in REMORA_FillPatch");
        }
    }
 
    mf_crse->FillBoundary(geom[lev-1].periodicity());
    amrex::InterpFromCoarseLevel(*mf_to_fill, mf_to_fill->nGrowVect(), IntVect(0,0,0),
            *mf_crse, 0, icomp, ncomp, geom[lev-1], geom[lev],
                                 refRatio(lev-1),
                                 mapper, domain_bcs_type, bccomp);
 
    // ***************************************************************************
    // Physical bc's at domain boundary
    // ***************************************************************************
 
    // The radiation conditions compare the field with its previous state (mf_calc),
    // which a level being created does not have. Give them the interpolated field
    // itself: the ghost cells then keep the coarse values and the nudging term,
    // with no time step, is zero. Only the radiation conditions read it, so the
    // copy is made only when one of them is set on this field, and only of the
    // components being filled. It cannot be an alias of the field: the radiation
    // kernels read the along-boundary neighbours of the column they write, so on a
    // GPU a thread could read a cell another thread has already overwritten.
    MultiFab calc_copy;
    const MultiFab* calc = &mf_calc;
    int icomp_calc_used = icomp_calc;
    if (!mf_calc.ok()) {
        bool radiation = false;
        for (int n = 0; n < ncomp; ++n) {
            const BCRec& bcr = domain_bcs_type[bccomp+n];
            for (int dir = 0; dir < 2; ++dir) {
                for (int bc : {bcr.lo(dir), bcr.hi(dir)}) {
                    radiation = radiation || bc == REMORABCType::orlanski_rad
                                          || bc == REMORABCType::orlanski_rad_nudge;
                }
            }
        }
        if (radiation) {
            calc_copy.define(mf_to_fill->boxArray(), mf_to_fill->DistributionMap(),
                             ncomp, mf_to_fill->nGrowVect());
            MultiFab::Copy(calc_copy, *mf_to_fill, icomp, 0, ncomp,
                           mf_to_fill->nGrowVect());
            calc = &calc_copy;
            icomp_calc_used = 0;
        }
    }
 
    // Enforce physical boundary conditions
    (*physbcs[lev])(*mf_to_fill,*mask,icomp,ncomp,mf_to_fill->nGrowVect(),
            time,bccomp,n_not_fill,*calc,
            *vec_msku[lev].get(), *vec_mskv[lev].get());
 
#ifdef REMORA_USE_NETCDF
    // Fill the data which is stored in the boundary data read from netcdf files
    if ( (solverChoice.boundary_from_netcdf) && (bdy_var_type != BdyVars::null) )
    {
        fill_from_bdyfiles(lev, *mf_to_fill,*mask,time,bccomp,bdy_var_type, icomp,icomp_calc_used,*calc,dt_lev);
    }
#else
    amrex::ignore_unused(icomp_calc_used, dt_lev);
#endif
    // Fill corners of the domain with periodic data
    if  ( box_mf.ixType() == IndexType(IntVect(0,0,0)) ) {
        mf_to_fill->EnforcePeriodicity(geom[lev].periodicity());
    } else {
        // The face-centred analogue. InterpFromCoarseLevel fills a periodic ghost from the
        // coarse level, which for a face tangential to the periodic direction does not
        // reproduce the fine value it wraps to: on a new level of Channel_Test, vbar's x
        // ghosts at the seam differed from their images by 2e-4 to 6e-4 where zeta's,
        // corrected above, were exact. The first barotropic substep then sees different vbar
        // at the two copies of the periodic u-face, and that one-shot asymmetry amplifies into
        // a blow-up. Copying the ghosts from this level's own valid data, as every later
        // FillPatch does, removes it.
        mf_to_fill->FillBoundary(geom[lev].periodicity());
    }
 
    // Enforce free-slip at top boundary (on xvel or yvel)
    if ( (box_mf.ixType() == IndexType(IntVect(1,0,0))) ||
         (box_mf.ixType() == IndexType(IntVect(0,1,0))) )
    {
        int khi = geom[lev].Domain().bigEnd(2);
        for (MFIter mfi(*mf_to_fill); mfi.isValid(); ++mfi)
        {
            Box gbx  = mfi.growntilebox(); // Note this is face-centered since vel is
            gbx.setSmall(2,khi+1);
            if (gbx.ok()) {
                Array4<Real> vel_arr = mf_to_fill->array(mfi);
                ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    vel_arr(i,j,k) = vel_arr(i,j,khi);
                });
            }
        }
    }
 
}
 
namespace {
/**
 * Sign for reflecting one cell-centered velocity component into a ghost cell.
 * The normal velocity vanishes at a wall or symmetry plane; only a no-slip wall
 * zeroes the tangential ones too. Anything else copies the adjacent valid cell.
 *
 * @param[in   ] bc        physical BC type for that component on that face
 * @param[in   ] is_normal whether the component is normal to the face
 */
Real
cc_vel_bdy_mult (REMORA_BC bc, bool is_normal)
{
    if (is_normal) {
        return ( (bc == REMORA_BC::slip_wall   ) ||
                 (bc == REMORA_BC::no_slip_wall) ||
                 (bc == REMORA_BC::symmetry    ) ) ? -one : one;
    }
    return (bc == REMORA_BC::no_slip_wall) ? -one : one;
}
} // namespace
 
/**
 * Fill boundary ghost cells for cell-centered velocities. Each component takes
 * the sign implied by its own BC, which with remora.boundary_per_variable can
 * differ between u, v and w on one face.
 *
 * @param[in   ] lev        level to fill at
 * @param[inout] mf_cc_vel  MultiFab of cell-centered velocities
 */
void
REMORA::FillBdyCCVels (int lev, MultiFab& mf_cc_vel)
{
    // Impose bc's at domain boundaries
    Box domain(Geom(lev).Domain());
 
    int ihi = domain.bigEnd(0);
    int jhi = domain.bigEnd(1);
    int khi = domain.bigEnd(2);
 
    // Impose periodicity first
    mf_cc_vel.FillBoundary(geom[lev].periodicity());
 
    // Sign multipliers per face, per velocity component
    GpuArray<GpuArray<Real,AMREX_SPACEDIM>,AMREX_SPACEDIM*2> mult;
    for (OrientationIter oit; oit; ++oit) {
        Orientation ori = oit();
        for (int n = 0; n < AMREX_SPACEDIM; ++n) {
            mult[ori][n] = cc_vel_bdy_mult(phys_bc_type[xvel_bc()+n][ori],
                                           n == ori.coordDir());
        }
    }
 
    // Callers pass different ghost widths (the plotfile ones carry none in z)
    const IntVect ng = mf_cc_vel.nGrowVect();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(mf_cc_vel, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        // Note that we don't fill corners here -- only the cells that share a face
        //      with interior cells -- this is all that is needed to calculate vorticity
        const Box& bx = mfi.tilebox();
        const Array4<Real>& vel_arr = mf_cc_vel.array(mfi);
 
        if (!Geom(lev).isPeriodic(0) && ng[0] > 0) {
            // Low-x side
            if (bx.smallEnd(0) <= domain.smallEnd(0)) {
                auto const m = mult[Orientation(Direction::x,Orientation::low)];
                ParallelFor(makeSlab(bx,0,0), [=] AMREX_GPU_DEVICE(int , int j, int k) noexcept
                {
                    vel_arr(-1,j,k,0) = m[0]*vel_arr(0,j,k,0); // u
                    vel_arr(-1,j,k,1) = m[1]*vel_arr(0,j,k,1); // v
                    vel_arr(-1,j,k,2) = m[2]*vel_arr(0,j,k,2); // w
                });
            }
 
            // High-x side
            if (bx.bigEnd(0) >= domain.bigEnd(0)) {
                auto const m = mult[Orientation(Direction::x,Orientation::high)];
                ParallelFor(makeSlab(bx,0,0), [=] AMREX_GPU_DEVICE(int , int j, int k) noexcept
                {
                    vel_arr(ihi+1,j,k,0) = m[0]*vel_arr(ihi,j,k,0); // u
                    vel_arr(ihi+1,j,k,1) = m[1]*vel_arr(ihi,j,k,1); // v
                    vel_arr(ihi+1,j,k,2) = m[2]*vel_arr(ihi,j,k,2); // w
                });
            }
        } // !periodic
 
        if (!Geom(lev).isPeriodic(1) && ng[1] > 0) {
            // Low-y side
            if (bx.smallEnd(1) <= domain.smallEnd(1)) {
                auto const m = mult[Orientation(Direction::y,Orientation::low)];
                ParallelFor(makeSlab(bx,1,0), [=] AMREX_GPU_DEVICE(int i, int  , int k) noexcept
                {
                    vel_arr(i,-1,k,0) = m[0]*vel_arr(i,0,k,0); // u
                    vel_arr(i,-1,k,1) = m[1]*vel_arr(i,0,k,1); // v
                    vel_arr(i,-1,k,2) = m[2]*vel_arr(i,0,k,2); // w
                });
            }
 
            // High-y side
            if (bx.bigEnd(1) >= domain.bigEnd(1)) {
                auto const m = mult[Orientation(Direction::y,Orientation::high)];
                ParallelFor(makeSlab(bx,1,0), [=] AMREX_GPU_DEVICE(int i, int , int k) noexcept
                {
                    vel_arr(i,jhi+1,k,0) = m[0]*vel_arr(i,jhi,k,0); // u
                    vel_arr(i,jhi+1,k,1) = m[1]*vel_arr(i,jhi,k,1); // v
                    vel_arr(i,jhi+1,k,2) = m[2]*vel_arr(i,jhi,k,2); // w
                });
            }
        } // !periodic
 
        if (!Geom(lev).isPeriodic(2) && ng[2] > 0) {
            // Low-z side
            if (bx.smallEnd(2) <= domain.smallEnd(2)) {
                auto const m = mult[Orientation(Direction::z,Orientation::low)];
                ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
                {
                    vel_arr(i,j,-1,0) = m[0]*vel_arr(i,j,0,0); // u
                    vel_arr(i,j,-1,1) = m[1]*vel_arr(i,j,0,1); // v
                    vel_arr(i,j,-1,2) = m[2]*vel_arr(i,j,0,2); // w
                });
            }
 
            // High-z side
            if (bx.bigEnd(2) >= domain.bigEnd(2)) {
                auto const m = mult[Orientation(Direction::z,Orientation::high)];
                ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
                {
                    vel_arr(i,j,khi+1,0) = m[0]*vel_arr(i,j,khi,0); // u
                    vel_arr(i,j,khi+1,1) = m[1]*vel_arr(i,j,khi,1); // v
                    vel_arr(i,j,khi+1,2) = m[2]*vel_arr(i,j,khi,2); // w
                });
            }
        } // !periodic
    } // MFIter
}