ERF_InitDensityHSE.H

Go to the documentation of this file.

/**
 * Initialize hydrostatically balanced density
 *
 * Calls init_isentropic_hse_terrain() or init_isentropic_hse() for cases with
 * and without terrain, respectively. Hydrostatic equilibrium (HSE) is satisfied discretely.
 * Note that these routines presume that qv==0 when evaluating the EOS. Both
 * density and pressure in HSE are calculated but only the density is used at
 * this point.
*/
 
#include "ERF_HSEUtils.H"
 
/**
 * @brief Initialize density in hydrostatic equilibrium for dry cases.
 *
 * @param[out] rho_hse MultiFab to be filled with HSE density.
 * @param z_phys_nd Unused nodal physical height MultiFab.
 * @param[in] z_phys_cc Cell-centered physical height MultiFab.
 * @param[in] geom Geometry defining the domain.
 * @param[in] stretched_dz_h Vector of stretched grid spacings.
 * @param[in] is_constant_dz Whether the grid spacing is constant.
 * @param[in] is_stretched_dz Whether the grid spacing is stretched.
 */
void
erf_init_dens_hse_dry (amrex::MultiFab& rho_hse,
                       std::unique_ptr<amrex::MultiFab>& /*z_phys_nd*/,
                       std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                       amrex::Geometry const& geom,
                       const amrex::Vector<amrex::Real>& stretched_dz_h,
                       bool is_constant_dz, bool is_stretched_dz) override
{
    const amrex::Real T_sfc    = amrex::Real(300.);
    const amrex::Real rho_sfc  = p_0 / (R_d*T_sfc);
    const amrex::Real Thetabar = T_sfc;
 
    if (!is_constant_dz && !is_stretched_dz) {
 
        const int domlo_z = geom.Domain().smallEnd(2);
        const int domhi_z = geom.Domain().bigEnd(2);
        if (domhi_z > 255) amrex::Abort("1D Arrays are hard-wired to only 256 high");
 
        //
        // A column is integrated upward from the bottom of its box.  A box stacked on another box
        //    of this level (the BoxArray is split in z, e.g. amr.max_grid_size below the number of
        //    cells in z) must continue from the density that box reached, not start afresh from
        //    the ghost cell below it.  The boxes are therefore integrated in bands of equal lowest
        //    index, bottom up, and before each band the cells just below it are filled from the
        //    bands already done.  With a single band nothing is filled.
        //
        const amrex::Vector<int> bands = column_bands(rho_hse.boxArray());
 
        for (const int klo_band : bands)
        {
            if (klo_band != bands[0]) {
                fill_below_band(rho_hse, 0, 1, klo_band, amrex::IntVect(1,1,0), geom);
            }
 
            for ( amrex::MFIter mfi(rho_hse, TileNoZ()); mfi.isValid(); ++mfi )
            {
                amrex::Array4<amrex::Real      > rho_arr  = rho_hse.array(mfi);
                amrex::Array4<amrex::Real const> z_cc_arr = z_phys_cc->const_array(mfi);
 
                // Create a flat box with same horizontal extent but only one cell in vertical
                const amrex::Box& tbz = mfi.nodaltilebox(2);
                amrex::Box b2d = tbz; // Copy constructor
                b2d.grow(0,1); b2d.grow(1,1); // Grow by one in the lateral directions
                b2d.setRange(2,0);
 
                const int klo = tbz.smallEnd(2);
                const int khi = tbz.bigEnd(2)-1;
 
                if (klo != klo_band) { continue; }
 
                amrex::ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
                {
                    //
                    // The density we start each column's integration from.  If this box does
                    //    not reach the bottom of the domain then that is the density in the
                    //    ghost cell below the box: the density the box of this level below it
                    //    reached or, where there is no such box, the one interpolated from the
                    //    coarser level.
                    //
                    // NOTE: we must do this lookup here, inside the kernel, for two reasons:
                    //       rho_arr is device memory, so reading it on the host is invalid in
                    //       a GPU build, and -- more importantly -- each column has its own
                    //       surface density, so we must not use the value at one corner of
                    //       the box for every column.
                    //
                    const amrex::Real rho_local_sfc = (klo == 0) ? rho_sfc : rho_arr(i,j,klo-1);
 
                    amrex::Array1D<amrex::Real,0,255> r;
                    amrex::Array1D<amrex::Real,0,255> p;
                    HSEutils::init_isentropic_hse_terrain(i,j,rho_local_sfc,Thetabar,&(r(0)),&(p(0)),z_cc_arr,klo,khi);
 
                    for (int k = klo; k <= khi; k++) {
                       rho_arr(i,j,k) = r(k);
                    }
 
                    // Impose Neumann conditions at bottom and top of domain boundary
                    if (klo == domlo_z) {
                        rho_arr(i,j,domlo_z-1) = rho_arr(i,j,domlo_z);
                    }
                    if (khi == domhi_z) {
                        rho_arr(i,j,domhi_z+1) = rho_arr(i,j,domhi_z);
                    }
                });
            } // mfi
        } // band
 
    } else {
 
 
        // Note that this integrates over the entire domain height,
        //      but only copies into each box as appropriate.
        // There is no assumption that any one box spans the whole vertical height.
        const int klo = geom.Domain().smallEnd(2);
        const int khi = geom.Domain().bigEnd(2);
 
        // These are at cell centers (unstaggered)
        amrex::Vector<amrex::Real> h_r(khi+2);
        amrex::Vector<amrex::Real> h_p(khi+2);
 
        amrex::Gpu::DeviceVector<amrex::Real> d_r(khi+2);
        amrex::Gpu::DeviceVector<amrex::Real> d_p(khi+2);
 
        if (is_constant_dz) {
            const amrex::Real dz = geom.CellSize()[2];
            HSEutils::init_isentropic_hse_constant_dz(rho_sfc,Thetabar,h_r.data(),h_p.data(),dz,klo,khi);
        } else { // stretched_dz
            HSEutils::init_isentropic_hse_stretched_dz(rho_sfc,Thetabar,h_r.data(),h_p.data(),stretched_dz_h.data(),klo,khi);
        }
 
        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin());
        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin());
 
        amrex::Real* r = d_r.data();
 
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
        for ( amrex::MFIter mfi(rho_hse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
        {
            const amrex::Box& bx = mfi.growntilebox(1);
            const amrex::Array4<amrex::Real> rho_hse_arr = rho_hse[mfi].array();
            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                int kk = std::max(k,0);
                rho_hse_arr(i,j,k) = r[kk];
            });
        } // mfi
    }
}
 
/**
 * @brief Initialize density in hydrostatic equilibrium for moist cases.
 *
 * @param[out] rho_hse MultiFab to be filled with HSE density.
 * @param z_phys_nd Unused nodal physical height MultiFab.
 * @param[in] geom Geometry defining the domain.
 */
void
erf_init_dens_hse_moist (amrex::MultiFab& rho_hse,
                         std::unique_ptr<amrex::MultiFab>& /*z_phys_nd*/,
                         amrex::Geometry const& geom) override
{
    const amrex::Real dz        = geom.CellSize()[2];
    const int khi        = geom.Domain().bigEnd()[2];
 
    // These are at cell centers (unstaggered)
    amrex::Vector<amrex::Real> h_r(khi+2);
    amrex::Vector<amrex::Real> h_p(khi+2);
    amrex::Vector<amrex::Real> h_t(khi+2);
    amrex::Vector<amrex::Real> h_q_v(khi+2);
 
    amrex::ParmParse pp("prob");
    amrex::Real q_t = amrex::Real(0.02);
    pp.query("qt_init",q_t);
 
    amrex::Real eq_pot_temp = amrex::Real(320.0);
    pp.query("eq_pot_temp",eq_pot_temp);
 
    bool use_empirical = false;
    pp.query("use_empirical_psat",use_empirical);
 
    amrex::Real z_tr_1 = -one;
    amrex::Real z_tr_2 = -one;
    pp.query("height", z_tr_1);
    pp.query("z_tr", z_tr_2);
 
    // Used only with T_from_theta_in_moist_init, but kept at physical values:
    // zeros raise an FPE in the speculatively evaluated compute_theta arm.
    amrex::Real T_tr = amrex::Real(300), theta_tr = amrex::Real(300), theta_0 = amrex::Real(300);
 
    bool T_from_theta = false;
    pp.query("T_from_theta_in_moist_init", T_from_theta);
 
    if (T_from_theta) {
        pp.get("theta_tr",theta_tr); pp.get("theta_0",theta_0);
        pp.get("T_tr",T_tr);
    }
 
    HSEutils::init_isentropic_hse_no_terrain( h_t.data(), h_r.data(), h_p.data(),
                                              h_q_v.data(), dz, khi, q_t, eq_pot_temp, use_empirical,
                                              T_from_theta, z_tr_1, z_tr_2,
                                              theta_0, theta_tr, T_tr);
 
    amrex::Gpu::DeviceVector<amrex::Real> d_r(khi+2);
    amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin());
    amrex::Real* r = d_r.data();
 
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
    for ( amrex::MFIter mfi(rho_hse,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const amrex::Box& bx = mfi.growntilebox(1);
        const amrex::Array4<amrex::Real> rho_hse_arr   = rho_hse[mfi].array();
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            int kk = std::max(k,0);
            rho_hse_arr(i,j,k) = r[kk];
        });
    } // mfi
}
 
/**
 * @brief Initialize HSE density as a constant value.
 *
 * @param[out] rho_hse MultiFab to be filled with constant HSE density.
 */
void
erf_init_const_dens_hse (amrex::MultiFab& rho_hse) override
{
    amrex::Real rho_0 = base_parms.rho_0;
    for ( amrex::MFIter mfi(rho_hse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
       amrex::Array4<amrex::Real> rho_hse_arr = rho_hse.array(mfi);
       const amrex::Box& gbx = mfi.growntilebox(1);
       ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
       {
           rho_hse_arr(i,j,k) = rho_0;
       });
    }
}
/**
 * @brief Initialize constant density and potential temperature in HSE.
 *
 * @param[out] rho_hse MultiFab for HSE density.
 * @param[out] p_hse MultiFab for HSE pressure.
 * @param[out] pi_hse MultiFab for HSE internal energy.
 * @param[out] th_hse MultiFab for HSE potential temperature.
 * @param[out] qv_hse MultiFab for HSE water vapor mixing ratio.
 * @param[in] l_rdOcp Constant R_d/Cp.
 */
void
erf_init_const_dens_and_th_hse (amrex::MultiFab& rho_hse, amrex::MultiFab&  p_hse,
                                amrex::MultiFab&  pi_hse, amrex::MultiFab& th_hse,
                                amrex::MultiFab&  qv_hse, amrex::Real l_rdOcp) override
{
    amrex::Real rho_0 = base_parms.rho_0;
    amrex::Real   T_0 = base_parms.T_0;
 
    amrex::ParmParse pp_prob("prob");
    amrex::Real dtheta_dz = zero; pp_prob.query("dtheta_dz",dtheta_dz);
 
    rho_hse.setVal(rho_0);
    th_hse.setVal(T_0);
 
    amrex::Real rt0 = rho_0 * T_0;
 
    amrex::Real p0  = getPgivenRTh(rt0);
    p_hse.setVal(p0);
 
    amrex::Real pi0 = getExnergivenRTh(rt0, l_rdOcp);
 
    pi_hse.setVal(pi0);
 
    // Default to zero background moisture
    qv_hse.setVal(0.0);
}
 
/**
 * @brief Initialize constant density and linearly varying potential temperature in HSE.
 *
 * @param[out] rho_hse MultiFab for HSE density.
 * @param[out] p_hse MultiFab for HSE pressure.
 * @param[out] pi_hse MultiFab for HSE internal energy.
 * @param[out] th_hse MultiFab for HSE potential temperature.
 * @param[out] qv_hse MultiFab for HSE water vapor mixing ratio.
 * @param[in] l_rdOcp Constant R_d/Cp.
 * @param[in] z_phys_cc Cell-centered physical height MultiFab.
 */
void
erf_init_const_dens_and_linear_th_hse (amrex::MultiFab& rho_hse, amrex::MultiFab&  p_hse,
                                       amrex::MultiFab&  pi_hse, amrex::MultiFab& th_hse,
                                       amrex::MultiFab&  qv_hse, amrex::Real l_rdOcp,
                                       std::unique_ptr<amrex::MultiFab>& z_phys_cc) override
{
    amrex::Real rho_0 = base_parms.rho_0;
    amrex::Real   T_0 = base_parms.T_0;
 
    amrex::ParmParse pp_prob("prob");
    amrex::Real dtheta_dz = zero; pp_prob.query("dtheta_dz",dtheta_dz);
 
    rho_hse.setVal(rho_0);
 
    for ( amrex::MFIter mfi(th_hse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
        amrex::Array4<amrex::Real> th_hse_arr  = th_hse.array(mfi);
        amrex::Array4<amrex::Real>  p_hse_arr  = p_hse.array(mfi);
        amrex::Array4<amrex::Real> pi_hse_arr  = pi_hse.array(mfi);
        amrex::Array4<amrex::Real const> z_arr = z_phys_cc->array(mfi);
        const amrex::Box& gbx = mfi.growntilebox(1);
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            th_hse_arr(i,j,k) = T_0 + dtheta_dz * (z_arr(i,j,k)- z_arr(i,j,0));
             p_hse_arr(i,j,k) = getPgivenRTh(rho_0*th_hse_arr(i,j,k));
            pi_hse_arr(i,j,k) = getExnergivenRTh(rho_0*th_hse_arr(i,j,k), l_rdOcp);
        });
    } // mfi
 
    // Default to zero background moisture
    qv_hse.setVal(0.0);
}