ERF_InitDensityHSE.H File Reference

#include "ERF_HSEUtils.H"

Include dependency graph for ERF_InitDensityHSE.H:

Go to the source code of this file.

Functions
void	erf_init_dens_hse_dry (amrex::MultiFab &rho_hse, std::unique_ptr< amrex::MultiFab > &, 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
	Initialize density in hydrostatic equilibrium for dry cases. More...
void	erf_init_dens_hse_moist (amrex::MultiFab &rho_hse, std::unique_ptr< amrex::MultiFab > &, amrex::Geometry const &geom) override
	Initialize density in hydrostatic equilibrium for moist cases. More...
void	erf_init_const_dens_hse (amrex::MultiFab &rho_hse) override
	Initialize HSE density as a constant value. More...
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
	Initialize constant density and potential temperature in HSE. More...
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
	Initialize constant density and linearly varying potential temperature in HSE. More...

Function Documentation

erf_init_const_dens_and_linear_th_hse()

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

Initialize constant density and linearly varying potential temperature in HSE.

Parameters

[out]	rho_hse	MultiFab for HSE density.
[out]	p_hse	MultiFab for HSE pressure.
[out]	pi_hse	MultiFab for HSE internal energy.
[out]	th_hse	MultiFab for HSE potential temperature.
[out]	qv_hse	MultiFab for HSE water vapor mixing ratio.
[in]	l_rdOcp	Constant R_d/Cp.
[in]	z_phys_cc	Cell-centered physical height MultiFab.

{
    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);
}

Here is the call graph for this function:

erf_init_const_dens_and_th_hse()

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

Initialize constant density and potential temperature in HSE.

Parameters

[out]	rho_hse	MultiFab for HSE density.
[out]	p_hse	MultiFab for HSE pressure.
[out]	pi_hse	MultiFab for HSE internal energy.
[out]	th_hse	MultiFab for HSE potential temperature.
[out]	qv_hse	MultiFab for HSE water vapor mixing ratio.
[in]	l_rdOcp	Constant R_d/Cp.

{
    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);
}

Here is the call graph for this function:

erf_init_const_dens_hse()

void erf_init_const_dens_hse ( amrex::MultiFab &  rho_hse )

override

Initialize HSE density as a constant value.

Parameters

[out]

rho_hse

MultiFab to be filled with constant HSE density.

{
    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;
       });
    }
}

Here is the call graph for this function:

erf_init_dens_hse_dry()

void erf_init_dens_hse_dry ( amrex::MultiFab &  rho_hse, std::unique_ptr< amrex::MultiFab > &  , 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

Initialize density in hydrostatic equilibrium for dry cases.

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.

Parameters

[out]	rho_hse	MultiFab to be filled with HSE density.
	z_phys_nd	Unused nodal physical height MultiFab.
[in]	z_phys_cc	Cell-centered physical height MultiFab.
[in]	geom	Geometry defining the domain.
[in]	stretched_dz_h	Vector of stretched grid spacings.
[in]	is_constant_dz	Whether the grid spacing is constant.
[in]	is_stretched_dz	Whether the grid spacing is stretched.

{
    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
    }
}

Here is the call graph for this function:

erf_init_dens_hse_moist()

void erf_init_dens_hse_moist ( amrex::MultiFab &  rho_hse, std::unique_ptr< amrex::MultiFab > &  , amrex::Geometry const &  geom  )

override

Initialize density in hydrostatic equilibrium for moist cases.

Parameters

[out]	rho_hse	MultiFab to be filled with HSE density.
	z_phys_nd	Unused nodal physical height MultiFab.
[in]	geom	Geometry defining the domain.

{
    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
}

Here is the call graph for this function: