ERF_Substep_NS.cpp File Reference

#include <ERF_TI_fast_headers.H>
#include "ERF_Constants.H"

Include dependency graph for ERF_Substep_NS.cpp:

Functions
void	erf_substep_NS (int step, int nrk, int level, int finest_level, Vector< MultiFab > &S_slow_rhs, const Vector< MultiFab > &S_prev, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &qt, const MultiFab &pi_stage, const MultiFab &fast_coeffs, Vector< MultiFab > &S_data, MultiFab &lagged_delta_rt, MultiFab &avg_xmom, MultiFab &avg_ymom, MultiFab &avg_zmom, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const Geometry geom, const Real gravity, amrex::Gpu::DeviceVector< amrex::Real > &stretched_dz_d, const double dtau_d, const Real facinv, Vector< std::unique_ptr< MultiFab >> &mapfac, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine, bool l_use_moisture, bool l_reflux, bool l_real_bc, const amrex::Real *sinesq_stag_d, const Real l_damp_coef, const Real beta_d, const CFFaceFluxes &cf_fluxes)

Function Documentation

erf_substep_NS()

void erf_substep_NS	(	int	step,
		int	nrk,
		int	level,
		int	finest_level,
		Vector< MultiFab > &	S_slow_rhs,
		const Vector< MultiFab > &	S_prev,
		Vector< MultiFab > &	S_stage_data,
		const MultiFab &	S_stage_prim,
		const MultiFab &	qt,
		const MultiFab &	pi_stage,
		const MultiFab &	fast_coeffs,
		Vector< MultiFab > &	S_data,
		MultiFab &	lagged_delta_rt,
		MultiFab &	avg_xmom,
		MultiFab &	avg_ymom,
		MultiFab &	avg_zmom,
		const MultiFab &	cc_src,
		const MultiFab &	xmom_src,
		const MultiFab &	ymom_src,
		const MultiFab &	zmom_src,
		const Geometry	geom,
		const Real	gravity,
		amrex::Gpu::DeviceVector< amrex::Real > &	stretched_dz_d,
		const double	dtau_d,
		const Real	facinv,
		Vector< std::unique_ptr< MultiFab >> &	mapfac,
		YAFluxRegister *	fr_as_crse,
		YAFluxRegister *	fr_as_fine,
		bool	l_use_moisture,
		bool	l_reflux,
		bool	l_real_bc,
		const amrex::Real *	sinesq_stag_d,
		const Real	l_damp_coef,
		const Real	beta_d,
		const CFFaceFluxes &	cf_fluxes
	)

Function for computing the fast RHS with no terrain and variable vertical spacing

Parameters

[in]	step	which fast time step within each Runge-Kutta step
[in]	nrk	which Runge-Kutta step
[in]	level	level of resolution
[in]	finest_level	finest level of resolution
[in]	S_slow_rhs	slow RHS computed in erf_slow_rhs_pre
[in]	S_prev	if step == 0, this is S_old, else the previous fast solution
[in]	S_stage_data	solution at previous RK stage
[in]	S_stage_prim	primitive variables at previous RK stage
[in]	pi_stage	Exner function at previous RK stage
[in]	fast_coeffs	coefficients for the tridiagonal solve used in the fast integrator, and beta_s on the w-faces (component 5)
[out]	S_data	current solution
[in,out]	lagged_delta_rt
[in,out]	avg_xmom	time-averaged x-momentum to be used for updating slow variables
[in,out]	avg_ymom	time-averaged y-momentum to be used for updating slow variables
[in,out]	avg_zmom	time-averaged z-momentum to be used for updating slow variables
[in]	cc_src	source terms for conserved variables
[in]	xmom_src	source terms for x-momentum
[in]	ymom_src	source terms for y-momentum
[in]	zmom_src	source terms for z-momentum
[in]	geom	container for geometric information
[in]	gravity	magnitude of gravity
[in]	stretched_dz_d
[in]	dtau	fast time step
[in]	facinv	inverse factor for time-averaging the momenta
[in]	mapfac	vector of map factors
[in,out]	fr_as_crse	YAFluxRegister at level l at level l / l+1 interface
[in,out]	fr_as_fine	YAFluxRegister at level l at level l-1 / l interface
[in]	l_use_moisture
[in]	l_reflux	should we add fluxes to the FluxRegisters?
[in]	l_real_bc
[in,out]	cf_fluxes	coarse/fine face fluxes: recorded here on a coarse level, imposed on a fine one

{
    //
    // NOTE: for step > 0, S_data and S_prev point to the same MultiFab data!!
    //
 
    BL_PROFILE_REGION("erf_substep_NS()");
 
    Real dtau = static_cast<Real>(dtau_d);
 
    const Box& domain = geom.Domain();
    auto const domlo  = lbound(domain);
    auto const domhi  = ubound(domain);
 
    int ilo = domlo.x;
    int ihi = domhi.x + 1;
    int jlo = domlo.y;
    int jhi = domhi.y + 1;
 
    bool l_rayleigh_impl_for_w = (sinesq_stag_d != nullptr);
 
    const Real* dx = geom.CellSize();
    const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
 
    Real dxi = dxInv[0];
    Real dyi = dxInv[1];
    Real dzi = dxInv[2];
 
    auto dz_ptr = stretched_dz_d.data();
 
    const auto& ba = S_stage_data[IntVars::cons].boxArray();
    const auto& dm = S_stage_data[IntVars::cons].DistributionMap();
 
    MultiFab Delta_rho_theta(        ba                , dm, 1, 1);
    MultiFab Delta_rho_w    (convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,0));
 
    MultiFab     coeff_A_mf(fast_coeffs, make_alias, 0, 1);
    MultiFab inv_coeff_B_mf(fast_coeffs, make_alias, 1, 1);
    MultiFab     coeff_C_mf(fast_coeffs, make_alias, 2, 1);
    MultiFab     coeff_P_mf(fast_coeffs, make_alias, 3, 1);
    MultiFab     coeff_Q_mf(fast_coeffs, make_alias, 4, 1);
 
    // beta_s on the w-faces, from which
    //    beta_1 = myhalf * (one - beta_s) multiplies explicit terms
    //    beta_2 = myhalf * (one + beta_s) multiplies implicit terms
    MultiFab      beta_s_mf(fast_coeffs, make_alias, 5, 1);
 
    // *************************************************************************
    // Set gravity as a vector
    const    Array<Real,AMREX_SPACEDIM> grav{zero, zero, -gravity};
    const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
 
    // This will hold theta extrapolated forward in time
    MultiFab extrap(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(),1,1);
 
    // This will hold the update for (rho) and (rho theta)
    MultiFab temp_rhs(S_stage_data[IntVars::zmom].boxArray(),S_stage_data[IntVars::zmom].DistributionMap(),2,0);
 
    // This will hold the new x- and y-momenta temporarily (so that we don't overwrite values we need when tiling)
    MultiFab temp_cur_xmom(S_stage_data[IntVars::xmom].boxArray(),S_stage_data[IntVars::xmom].DistributionMap(),1,0);
    MultiFab temp_cur_ymom(S_stage_data[IntVars::ymom].boxArray(),S_stage_data[IntVars::ymom].DistributionMap(),1,0);
 
    // We assume that in the first step (nrk == 0) we are only doing one substep.
    AMREX_ALWAYS_ASSERT(nrk > 0 || step == 0);
 
    // *************************************************************************
    // First set up some arrays we'll need
    // *************************************************************************
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Array4<const Real>& prev_cons  = S_prev[IntVars::cons].const_array(mfi);
        const Array4<const Real>& prev_zmom =  S_prev[IntVars::zmom].const_array(mfi);
 
        const Array4<const Real>& stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
        const Array4<const Real>& stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi);
 
        const Array4<Real>& prev_drho_w     = Delta_rho_w.array(mfi);
        const Array4<Real>& prev_drho_theta = Delta_rho_theta.array(mfi);
        const Array4<Real>& lagged_arr      = lagged_delta_rt.array(mfi);
        const Array4<Real>& theta_extrap    = extrap.array(mfi);
        const Array4<const Real>& prim      = S_stage_prim.const_array(mfi);
 
        Box gbx = mfi.growntilebox(1);
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            prev_drho_theta(i,j,k) = prev_cons(i,j,k,RhoTheta_comp) - stage_cons(i,j,k,RhoTheta_comp);
 
            if (step == 0) {
                theta_extrap(i,j,k) = prev_drho_theta(i,j,k);
            } else {
                theta_extrap(i,j,k) = prev_drho_theta(i,j,k) + beta_d *
                  ( prev_drho_theta(i,j,k) - lagged_arr(i,j,k) );
            }
 
            // NOTE: qv is not changing over the fast steps so we use the stage data
            Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero;
            theta_extrap(i,j,k) *= (one + RvoRd*qv);
 
            // We define lagged_delta_rt for our next step as the current delta_rt
            // (after using it above to extrapolate theta for this step)
            lagged_arr(i,j,k) = prev_drho_theta(i,j,k);
        });
 
        // NOTE: We must do this here because for step > 0, prev_zmom and cur_zmom both point to the same data,
        //       so by the time we would use prev_zmom to define zflux, it would have already been over-written.
        //
        // NOTE: grownnodaltilebox, not nodaltilebox(2) grown by hand.  Growing by hand pushes
        //       every tile one cell past its own share of the grid in x and y, so once the grid
        //       is tiled two tiles write the same Delta_rho_w cells -- a data race under OpenMP,
        //       even though both store the same value.  This is identical to the old expression
        //       when there is one tile per grid.
        //
        Box gtbz = mfi.grownnodaltilebox(2,IntVect(1,1,0));
        ParallelFor(gtbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
            prev_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k);
        });
    } // mfi
 
    // *************************************************************************
    // Define updates in the current RK stage
    // *************************************************************************
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        Box tbx = mfi.nodaltilebox(0);
        Box tby = mfi.nodaltilebox(1);
 
        const Array4<Real const>& xmom_src_arr   = xmom_src.const_array(mfi);
        const Array4<Real const>& ymom_src_arr   = ymom_src.const_array(mfi);
 
        const Array4<const Real> & stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
        const Array4<const Real> & stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
        const Array4<const Real> & qt_arr     = qt.const_array(mfi);
 
        const Array4<const Real>& slow_rhs_rho_u = S_slow_rhs[IntVars::xmom].const_array(mfi);
        const Array4<const Real>& slow_rhs_rho_v = S_slow_rhs[IntVars::ymom].const_array(mfi);
 
        const Array4<Real>& temp_cur_xmom_arr  = temp_cur_xmom.array(mfi);
        const Array4<Real>& temp_cur_ymom_arr  = temp_cur_ymom.array(mfi);
 
        const Array4<const Real>& prev_xmom = S_prev[IntVars::xmom].const_array(mfi);
        const Array4<const Real>& prev_ymom = S_prev[IntVars::ymom].const_array(mfi);
 
        // These store the advection momenta which we will use to update the slow variables
        const Array4<      Real>& avg_xmom_arr = avg_xmom.array(mfi);
        const Array4<      Real>& avg_ymom_arr = avg_ymom.array(mfi);
 
        const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
 
        const Array4<Real>& theta_extrap = extrap.array(mfi);
 
        // Map factors
        const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
        const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
        const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
        const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
 
        // *********************************************************************
        // Define updates in the RHS of {x, y, z}-momentum equations
        //
        // NOTE: avg_{x,y}mom is a *flux* that is later handed to the scalar
        //       advection in erf_slow_rhs_post, which applies (mfsq/detJ) to its
        //       divergence.  So what we accumulate here must be in the same
        //       convention as the base value set in AdvectionSrcForRho, namely
        //       ax*rho_u/mf_uy (and ay*rho_v/mf_vx).  Two weightings are needed:
        //
        //       (1) the map factor, matching the density fluxes formed below;
        //
        //       (2) h_zeta = ax = detJ, which for this (laterally homogeneous)
        //           mesh is just dz_ptr[k]/dz.  Note this factor cancels out of
        //           the rho update we do here -- temp_rhs uses a bare dxi and no
        //           1/detJ -- but it does NOT cancel for avg_{x,y}mom, because
        //           the consumer downstream re-applies the 1/detJ that we never
        //           applied.  It is identically 1 for MeshType::ConstantDz.
        // *********************************************************************
        if (nrk == 0 and step == 0) { // prev == stage
            ParallelFor(tbx, tby,
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                Real h_zeta = dz_ptr[k] * dzi;
                Real new_drho_u = dtau * slow_rhs_rho_u(i,j,k) + dtau * xmom_src_arr(i,j,k);;
                avg_xmom_arr(i,j,k) += facinv * new_drho_u * h_zeta / mf_uy(i,j,0);
                temp_cur_xmom_arr(i,j,k) = stage_xmom(i,j,k) + new_drho_u;
            },
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                Real h_zeta = dz_ptr[k] * dzi;
                Real new_drho_v = dtau * slow_rhs_rho_v(i,j,k) + dtau * ymom_src_arr(i,j,k);
                avg_ymom_arr(i,j,k) += facinv * new_drho_v * h_zeta / mf_vx(i,j,0);
                temp_cur_ymom_arr(i,j,k) = stage_ymom(i,j,k) + new_drho_v;
            });
        } else {
            ParallelFor(tbx, tby,
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
                Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? Real(0.) :
                  (theta_extrap(i,j,k) - theta_extrap(i-1,j,k))*dxi;
                gpx *= mf_ux(i,j,0);
 
                Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) : zero;
 
                Real pi_c =  myhalf * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i,j,k,0));
                Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx / (one + q);
 
                Real new_drho_u = prev_xmom(i,j,k) - stage_xmom(i,j,k)
                                + dtau * fast_rhs_rho_u + dtau * slow_rhs_rho_u(i,j,k)
                                + dtau * xmom_src_arr(i,j,k);
 
                Real h_zeta = dz_ptr[k] * dzi;
                avg_xmom_arr(i,j,k) += facinv * new_drho_u * h_zeta / mf_uy(i,j,0);
 
                temp_cur_xmom_arr(i,j,k) = stage_xmom(i,j,k) + new_drho_u;
            },
            [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
                Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? Real(0.) :
                  (theta_extrap(i,j,k) - theta_extrap(i,j-1,k))*dyi;
                gpy *= mf_vy(i,j,0);
 
                Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) : zero;
 
                Real pi_c =  myhalf * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j,k,0));
                Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy / (one + q);
 
                Real new_drho_v = prev_ymom(i,j,k) - stage_ymom(i,j,k)
                                + dtau * fast_rhs_rho_v + dtau * slow_rhs_rho_v(i,j,k)
                                + dtau * ymom_src_arr(i,j,k);
 
                Real h_zeta = dz_ptr[k] * dzi;
                avg_ymom_arr(i,j,k) += facinv * new_drho_v * h_zeta / mf_vx(i,j,0);
 
                temp_cur_ymom_arr(i,j,k) = stage_ymom(i,j,k) + new_drho_v;
            });
        } // nrk > 0 and/or step > 0
 
        // On the coarse/fine faces of a fine level, the coarse level's flux over this substep
        // (CFFaceFluxes) in place of the level's own
        for (int dir = 0; dir < 2; ++dir) {
            if (!cf_fluxes.tgt[dir]) { continue; }
            const Box& fbx = (dir == 0) ? tbx : tby;
            const Array4<const Real>& tgt = cf_fluxes.tgt[dir]->const_array(mfi);
            const Array4<Real>& cur = (dir == 0) ? temp_cur_xmom_arr : temp_cur_ymom_arr;
            const Array4<Real>& avg = (dir == 0) ? avg_xmom_arr      : avg_ymom_arr;
            const Array4<const Real>& mf = (dir == 0) ? mf_uy : mf_vx;
            ParallelFor(fbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                if (tgt(i,j,k,1) > myhalf) {
                    Real new_mom = tgt(i,j,k,0) * mf(i,j,0);
                    Real h_zeta  = dz_ptr[k] * dzi;
                    avg(i,j,k) += facinv * (new_mom - cur(i,j,k)) * h_zeta / mf(i,j,0);
                    cur(i,j,k)  = new_mom;
                }
            });
        }
 
        // On a coarse level, the flux each face applied over this substep
        for (int dir = 0; dir < 2; ++dir) {
            if (!cf_fluxes.rec[dir]) { continue; }
            const Box& fbx = (dir == 0) ? tbx : tby;
            const Array4<Real>& rec = cf_fluxes.rec[dir]->array(mfi);
            const Array4<const Real>& cur = (dir == 0) ? temp_cur_xmom_arr : temp_cur_ymom_arr;
            const Array4<const Real>& mf  = (dir == 0) ? mf_uy : mf_vx;
            ParallelFor(fbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                rec(i,j,k) = cur(i,j,k) / mf(i,j,0);
            });
        }
    } //mfi
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    {
    std::array<FArrayBox,AMREX_SPACEDIM> flux;
    for ( MFIter mfi(S_stage_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
    {
        Box bx  = mfi.tilebox();
        Box tbz = surroundingNodes(bx,2);
 
        Box vbx = mfi.validbox();
        const auto& vbx_hi = ubound(vbx);
 
        const Array4<Real const>& zmom_src_arr   = zmom_src.const_array(mfi);
 
        const Array4<const Real>& stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
        const Array4<const Real>& stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
        const Array4<const Real>& stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi);
        const Array4<const Real> & prim       = S_stage_prim.const_array(mfi);
 
        const Array4<const Real>& prev_drho_theta = Delta_rho_theta.array(mfi);
 
        const Array4<const Real>& prev_cons  = S_prev[IntVars::cons].const_array(mfi);
        const Array4<const Real>& stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
 
        const Array4<const Real>& slow_rhs_cons  = S_slow_rhs[IntVars::cons].const_array(mfi);
        const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi);
 
        const Array4<const Real>& prev_zmom = S_prev[IntVars::zmom].const_array(mfi);
        const Array4<      Real>&  cur_zmom = S_data[IntVars::zmom].array(mfi);
 
        const Array4<Real>& temp_cur_xmom_arr  = temp_cur_xmom.array(mfi);
        const Array4<Real>& temp_cur_ymom_arr  = temp_cur_ymom.array(mfi);
 
        // These store the advection momenta which we will use to update the slow variables
        const Array4<      Real>& avg_zmom_arr = avg_zmom.array(mfi);
 
        // Map factors
        const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
        const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
        const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
        const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
 
        FArrayBox RHS_fab;
        RHS_fab.resize(tbz,1, The_Async_Arena());
 
        FArrayBox soln_fab;
        soln_fab.resize(tbz,1, The_Async_Arena());
 
        auto const& RHS_a  = RHS_fab.array();
        auto const& soln_a = soln_fab.array();
 
        auto const& temp_rhs_arr = temp_rhs.array(mfi);
 
        auto const&     coeffA_a =     coeff_A_mf.array(mfi);
        auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
        auto const&     coeffC_a =     coeff_C_mf.array(mfi);
        auto const&     coeffP_a =     coeff_P_mf.array(mfi);
        auto const&     coeffQ_a =     coeff_Q_mf.array(mfi);
        auto const&     beta_s_a =      beta_s_mf.const_array(mfi);
 
        // *************************************************************************
        // Define flux arrays for use in advection
        // *************************************************************************
        for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
            flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
            flux[dir].setVal<RunOn::Device>(0);
        }
        const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
            flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
 
        // *********************************************************************
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
            Real xflux_lo = (temp_cur_xmom_arr(i  ,j,k) - stage_xmom(i  ,j,k)) / mf_uy(i  ,j,0);
            Real xflux_hi = (temp_cur_xmom_arr(i+1,j,k) - stage_xmom(i+1,j,k)) / mf_uy(i+1,j,0);
            Real yflux_lo = (temp_cur_ymom_arr(i,j  ,k) - stage_ymom(i,j  ,k)) / mf_vx(i,j  ,0);
            Real yflux_hi = (temp_cur_ymom_arr(i,j+1,k) - stage_ymom(i,j+1,k)) / mf_vx(i,j+1,0);
 
            Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
 
            temp_rhs_arr(i,j,k,Rho_comp     ) =  ( xflux_hi - xflux_lo ) * dxi * mfsq
                                               + ( yflux_hi - yflux_lo ) * dyi * mfsq;
            temp_rhs_arr(i,j,k,RhoTheta_comp) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
                                                   xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq +
                                                 ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
                                                   yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) * myhalf;
 
            if (l_reflux) {
                (flx_arr[0])(i,j,k,0) = xflux_lo;
                (flx_arr[0])(i,j,k,1) = (flx_arr[0])(i  ,j,k,0) * myhalf * (prim(i,j,k,0) + prim(i-1,j,k,0));
 
                (flx_arr[1])(i,j,k,0) = yflux_lo;
                (flx_arr[1])(i,j,k,1) = (flx_arr[1])(i,j  ,k,0) * myhalf * (prim(i,j,k,0) + prim(i,j-1,k,0));
 
                if (i == vbx_hi.x) {
                    (flx_arr[0])(i+1,j,k,0) = xflux_hi;
                    (flx_arr[0])(i+1,j,k,1) = (flx_arr[0])(i+1,j,k,0) * myhalf * (prim(i,j,k,0) + prim(i+1,j,k,0));
                }
                if (j == vbx_hi.y) {
                    (flx_arr[1])(i,j+1,k,0) = yflux_hi;
                    (flx_arr[1])(i,j+1,k,1) = (flx_arr[1])(i,j+1,k,0) * myhalf * (prim(i,j,k,0) + prim(i,j+1,k,0));
                }
            }
        });
 
        Box bx_shrunk_in_k = bx;
        int klo = tbz.smallEnd(2);
        int khi = tbz.bigEnd(2);
        bx_shrunk_in_k.setSmall(2,klo+1);
        bx_shrunk_in_k.setBig(2,khi-1);
 
        // Note that the notes use "g" to mean the magnitude of gravity, so it is positive
        // We set grav_gpu[2] to be the vector component which is negative
        // We define halfg to match the notes (which is why we take the absolute value)
        Real halfg = std::abs(myhalf * grav_gpu[2]);
 
        // *********************************************************************
        // fast_loop_on_shrunk
        // *********************************************************************
        //Note we don't act on the bottom or top boundaries of the domain
        ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            Real coeff_P = coeffP_a(i,j,k);
            Real coeff_Q = coeffQ_a(i,j,k);
 
            Real theta_t_lo  = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
            Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k  ,PrimTheta_comp) );
            Real theta_t_hi  = myhalf * ( prim(i,j,k  ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
 
            Real Omega_kp1 = prev_zmom(i,j,k+1) - stage_zmom(i,j,k+1);
            Real Omega_k   = prev_zmom(i,j,k  ) - stage_zmom(i,j,k  );
            Real Omega_km1 = prev_zmom(i,j,k-1) - stage_zmom(i,j,k-1);
 
            // line 2 last two terms (order dtau)
            Real old_drho_k   = prev_cons(i,j,k  ,Rho_comp) - stage_cons(i,j,k  ,Rho_comp);
            Real old_drho_km1 = prev_cons(i,j,k-1,Rho_comp) - stage_cons(i,j,k-1,Rho_comp);
            Real R0_tmp = coeff_P * prev_drho_theta(i,j,k) + coeff_Q * prev_drho_theta(i,j,k-1)
                         - halfg * ( old_drho_k + old_drho_km1 );
 
            // lines 3-5 residuals (order dtau^2) one <-> beta_2
            Real R1_tmp =  halfg * (-slow_rhs_cons(i,j,k  ,Rho_comp) - slow_rhs_cons(i,j,k-1,Rho_comp)
                                    + temp_rhs_arr(i,j,k  ,Rho_comp) +  temp_rhs_arr(i,j,k-1,Rho_comp) )
                + ( coeff_P * (slow_rhs_cons(i,j,k  ,RhoTheta_comp) - temp_rhs_arr(i,j,k  ,RhoTheta_comp)) +
                    coeff_Q * (slow_rhs_cons(i,j,k-1,RhoTheta_comp) - temp_rhs_arr(i,j,k-1,RhoTheta_comp)) );
 
            // Explicit part of the vertical flux on each face is beta_1 * Omega. beta_1 on face k
            // is factored out, and dbeta_1 is how much larger beta_1 is on the neighbouring faces,
            // which is exactly zero (so the roundoff is that of a scalar beta_1) for uniform beta_s
            Real beta_1      = myhalf * (one - beta_s_a(i,j,k));
            Real dbeta_1_kp1 = myhalf * (one - beta_s_a(i,j,k+1)) - beta_1;
            Real dbeta_1_km1 = myhalf * (one - beta_s_a(i,j,k-1)) - beta_1;
 
            // lines 6&7 consolidated (reuse Omega & metrics) (order dtau^2)
            R1_tmp +=  beta_1 * ( ( (Omega_kp1 - Omega_k) / dz_ptr[k] + (Omega_k - Omega_km1) / dz_ptr[k-1] ) * halfg
                                 +(-(Omega_kp1*theta_t_hi  - Omega_k  *theta_t_mid) * coeff_P / dz_ptr[k]
                                   -(Omega_k  *theta_t_mid - Omega_km1*theta_t_lo ) * coeff_Q / dz_ptr[k-1]) )
                     + ( ( dbeta_1_kp1*Omega_kp1 / dz_ptr[k] - dbeta_1_km1*Omega_km1 / dz_ptr[k-1] ) * halfg
                        -  dbeta_1_kp1*Omega_kp1*theta_t_hi * coeff_P / dz_ptr[k]
                        +  dbeta_1_km1*Omega_km1*theta_t_lo * coeff_Q / dz_ptr[k-1] );
 
            // line 1
            Real beta_2 = myhalf * (one + beta_s_a(i,j,k));
            RHS_a(i,j,k) = Omega_k + dtau * (slow_rhs_rho_w(i,j,k) + R0_tmp + dtau * beta_2 * R1_tmp + zmom_src_arr(i,j,k));
 
        }); // bx_shrunk_in_k
 
        Box b2d = tbz;
        b2d.setRange(2,0);
 
        auto const lo = lbound(bx);
        auto const hi = ubound(bx);
 
        ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
        {
            // NOTE: RHS is built on a shrunk box in z, explicitly set
            //       the lo and hi faces here
            // On a box face inside the domain (a c/f face), w = w_old + dtau * slow_rhs.
            // At the domain bottom and top the row is zero: W'' = 0 at the ground and at
            // a wall, and W''_top = W''_top-1 for outflow, the matrix row set in
            // make_fast_coeffs deciding which
            RHS_a (i,j,lo.z) = (lo.z == domlo.z)
                             ? zero
                             : prev_zmom(i,j,lo.z) - stage_zmom(i,j,lo.z)
                             + dtau * slow_rhs_rho_w(i,j,lo.z)
                             + dtau * zmom_src_arr(i,j,lo.z);
 
            RHS_a (i,j,hi.z+1) = (hi.z == domhi.z)
                               ? zero
                               : prev_zmom(i,j,hi.z+1) - stage_zmom(i,j,hi.z+1)
                               + dtau * slow_rhs_rho_w(i,j,hi.z+1)
                               + dtau * zmom_src_arr(i,j,hi.z+1);
        }); // b2d
 
        // (rho w) on the coarse/fine faces of a fine level at the end of this substep, from the coarse level
        if (cf_fluxes.tgt_w) {
            const Array4<const Real>& tgt   = cf_fluxes.tgt[2]->const_array(mfi);
            const Array4<const Real>& tgt_w = cf_fluxes.tgt_w->const_array(mfi);
            ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                if (tgt(i,j,lo.z,1) > myhalf) {
                    RHS_a(i,j,lo.z) = tgt_w(i,j,lo.z) - stage_zmom(i,j,lo.z);
                }
                if (tgt(i,j,hi.z+1,1) > myhalf) {
                    RHS_a(i,j,hi.z+1) = tgt_w(i,j,hi.z+1) - stage_zmom(i,j,hi.z+1);
                }
            });
        }
 
#ifdef AMREX_USE_GPU
        ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
        {
            // Dirichlet condition at k = lo.z
            soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
          cur_zmom(i,j,lo.z) = stage_zmom(i,j,lo.z) + soln_a(i,j,lo.z);
 
          for (int k = lo.z+1; k <= hi.z+1; k++) {
              soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
          }
 
          cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
 
          for (int k = hi.z; k >= lo.z; k--) {
              soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1);
              cur_zmom(i,j,k) = stage_zmom(i,j,k) + soln_a(i,j,k);
          }
        }); // b2d
#else
        for (int j = lo.y; j <= hi.y; ++j) {
            AMREX_PRAGMA_SIMD
            for (int i = lo.x; i <= hi.x; ++i) {
                soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
            }
        }
        for (int k = lo.z+1; k <= hi.z+1; ++k) {
            for (int j = lo.y; j <= hi.y; ++j) {
                AMREX_PRAGMA_SIMD
                for (int i = lo.x; i <= hi.x; ++i) {
                    soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
                }
            }
        }
        for (int j = lo.y; j <= hi.y; ++j) {
            AMREX_PRAGMA_SIMD
            for (int i = lo.x; i <= hi.x; ++i) {
                cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
            }
        }
        for (int k = hi.z; k >= lo.z; --k) {
            for (int j = lo.y; j <= hi.y; ++j) {
                AMREX_PRAGMA_SIMD
                for (int i = lo.x; i <= hi.x; ++i) {
                    soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
                    cur_zmom(i,j,k) = stage_zmom(i,j,k) + soln_a(i,j,k);
                }
            }
        }
#endif
        if (l_rayleigh_impl_for_w) {
            // Every face but the ground, as in erf_substep_T: the top face is no longer
            // pinned with an outflow top, and the faces at the ends of the box are solved for.
            // A coarse/fine face whose (rho w) is imposed from the coarse level's record is not
            // damped: the coarse level has damped it already, and it must match the record
            const bool has_tgt_w = (cf_fluxes.tgt_w != nullptr);
            const Array4<const Real> tgt_flag = (has_tgt_w) ? cf_fluxes.tgt[2]->const_array(mfi) : Array4<const Real>{};
            ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
              const bool imposed = has_tgt_w && (tgt_flag(i,j,k,1) > myhalf);
              if (k > domlo.z && !imposed) {
                  Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
                  cur_zmom(i,j,k) /= (one + damping_coeff);
              }
            });
        }
 
        // **************************************************************************
        // Define updates in the RHS of rho and (rho theta)
        // **************************************************************************
        const Array4<Real>&  prev_drho_w = Delta_rho_w.array(mfi);
        const bool has_tgt_z = (cf_fluxes.tgt[2] != nullptr);
        const bool has_rec_z = (cf_fluxes.rec[2] != nullptr);
        const Array4<const Real> tgt_z = (has_tgt_z) ? cf_fluxes.tgt[2]->const_array(mfi) : Array4<const Real>{};
        const Array4<      Real> rec_z = (has_rec_z) ? cf_fluxes.rec[2]->array(mfi)       : Array4<Real>{};
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            Real zflux_lo = myhalf * (one + beta_s_a(i,j,k  )) * soln_a(i,j,k  )
                          + myhalf * (one - beta_s_a(i,j,k  )) * prev_drho_w(i,j,k  );
            Real zflux_hi = myhalf * (one + beta_s_a(i,j,k+1)) * soln_a(i,j,k+1)
                          + myhalf * (one - beta_s_a(i,j,k+1)) * prev_drho_w(i,j,k+1);
 
            // On the coarse/fine faces of a fine level, the coarse level's flux over this substep
            if (has_tgt_z) {
                Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
                if (tgt_z(i,j,k  ,1) > myhalf) { zflux_lo = tgt_z(i,j,k  ,0) * mfsq - stage_zmom(i,j,k  ); }
                if (tgt_z(i,j,k+1,1) > myhalf) { zflux_hi = tgt_z(i,j,k+1,0) * mfsq - stage_zmom(i,j,k+1); }
            }
 
            // On a coarse level, the flux each face applied over this substep
            if (has_rec_z) {
                Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
                rec_z(i,j,k) = (zflux_lo + stage_zmom(i,j,k)) / mfsq;
                if (k == vbx_hi.z) {
                    rec_z(i,j,k+1) = (zflux_hi + stage_zmom(i,j,k+1)) / mfsq;
                }
            }
 
            avg_zmom_arr(i,j,k)      += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
            if (l_reflux) {
                (flx_arr[2])(i,j,k,0) =        zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
                (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k) + prim(i,j,k-1));
            }
 
            if (k == vbx_hi.z) {
                avg_zmom_arr(i,j,k+1)      += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
                if (l_reflux) {
                    (flx_arr[2])(i,j,k+1,0) =          zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
                    (flx_arr[2])(i,j,k+1,1) = (flx_arr[2])(i,j,k+1,0) * myhalf * (prim(i,j,k) + prim(i,j,k+1));
                }
            }
 
            Real dz_inv = one / dz_ptr[k];
            temp_rhs_arr(i,j,k,Rho_comp     ) += dz_inv * ( zflux_hi - zflux_lo );
            temp_rhs_arr(i,j,k,RhoTheta_comp) += myhalf * dz_inv * ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1))
                                                                   - zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) );
        });
 
        // We only add to the flux registers in the final RK step
        if (l_reflux) {
            int strt_comp_reflux = 0;
            // For now we don't reflux (rho theta) because it seems to create issues at c/f boundaries
            int  num_comp_reflux = 1;
            if (level < finest_level) {
                fr_as_crse->CrseAdd(mfi,
                    {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
                    dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
            }
            if (level > 0) {
                fr_as_fine->FineAdd(mfi,
                    {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
                    dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
            }
 
            // This is necessary here so we don't go on to the next FArrayBox without
            // having finished copying the fluxes into the FluxRegisters (since the fluxes
            // are stored in temporary FArrayBox's)
            Gpu::streamSynchronize();
 
        } // two-way coupling
    } // mfi
    } // OMP
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
 
        const Array4<      Real>&  cur_cons = S_data[IntVars::cons].array(mfi);
        const Array4<const Real>& prev_cons = S_prev[IntVars::cons].const_array(mfi);
        auto const& temp_rhs_arr     = temp_rhs.const_array(mfi);
        auto const& slow_rhs_cons    = S_slow_rhs[IntVars::cons].const_array(mfi);
        const Array4<Real const>& cc_src_arr   = cc_src.const_array(mfi);
 
        if (step == 0) {
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                cur_cons(i,j,k,Rho_comp)      = prev_cons(i,j,k,Rho_comp) +
                                                dtau * (slow_rhs_cons(i,j,k,Rho_comp) - temp_rhs_arr(i,j,k,Rho_comp));
                cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp) +
                                                dtau * (slow_rhs_cons(i,j,k,RhoTheta_comp) - temp_rhs_arr(i,j,k,RhoTheta_comp));
 
                // add in source terms for cell-centered conserved variables
                cur_cons(i,j,k,Rho_comp)      += dtau * cc_src_arr(i,j,k,Rho_comp);
                cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp);
            });
        } else {
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                //
                // We didn't need to set cur_cons = prev_cons above because they point to the same data for step > 0
                //
                cur_cons(i,j,k,Rho_comp)      += dtau * (slow_rhs_cons(i,j,k,Rho_comp) - temp_rhs_arr(i,j,k,Rho_comp));
                cur_cons(i,j,k,RhoTheta_comp) += dtau * (slow_rhs_cons(i,j,k,RhoTheta_comp) - temp_rhs_arr(i,j,k,RhoTheta_comp));
 
                // add in source terms for cell-centered conserved variables
                cur_cons(i,j,k,Rho_comp)      += dtau * cc_src_arr(i,j,k,Rho_comp);
                cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp);
            });
        } // step = 0
 
        const Array4<Real>& cur_xmom = S_data[IntVars::xmom].array(mfi);
        const Array4<Real>& cur_ymom = S_data[IntVars::ymom].array(mfi);
 
        const Array4<Real const>& temp_cur_xmom_arr = temp_cur_xmom.const_array(mfi);
        const Array4<Real const>& temp_cur_ymom_arr = temp_cur_ymom.const_array(mfi);
 
        Box tbx = surroundingNodes(bx,0);
        Box tby = surroundingNodes(bx,1);
 
        ParallelFor(tbx, tby,
        [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            cur_xmom(i,j,k) = temp_cur_xmom_arr(i,j,k);
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            cur_ymom(i,j,k) = temp_cur_ymom_arr(i,j,k);
        });
 
    } // mfi
}

Here is the call graph for this function: