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"
 
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");
 
        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 ilo = tbz.smallEnd(0);
            const int jlo = tbz.smallEnd(1);
            const int klo = tbz.smallEnd(2);
            const int khi = tbz.bigEnd(2)-1;
 
            amrex::Real rho_local_sfc;
            if (klo == 0) {
                rho_local_sfc = rho_sfc;
            } else {
                rho_local_sfc = rho_arr(ilo,jlo,klo-1);
            }
 
            amrex::ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                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
 
    } 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
    }
}
 
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);
 
    // We only set them to zero to make sure they are initialized
    amrex::Real T_tr = zero, theta_tr = zero, theta_0 = zero;
 
    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
}
 
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;
       });
    }
}
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);
}
 
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);
}