SHOCInterface Class Reference

#include <ERF_ShocInterface.H>

Collaboration diagram for SHOCInterface:

Classes
struct	Buffer
struct	SHOCPostprocess
struct	SHOCPreprocess

Public Member Functions
	SHOCInterface (const int &lev, SolverChoice &sc)
void	set_grids (int &level, const amrex::BoxArray &ba, amrex::Geometry &geom, amrex::MultiFab *cons, amrex::MultiFab *xvel, amrex::MultiFab *yvel, amrex::MultiFab *zvel, amrex::Real *w_subsid, amrex::MultiFab *tau13, amrex::MultiFab *tau23, amrex::MultiFab *hfx3, amrex::MultiFab *qfx3, amrex::MultiFab *eddyDiffs, amrex::MultiFab *z_phys)
void	initialize_impl ()
void	run_impl (const double dt)
void	finalize_impl (const double dt)
void	alloc_buffers ()
void	dealloc_buffers ()
void	mf_to_kokkos_buffers ()
void	kokkos_buffers_to_mf (const double dt)
std::string	name () const
void	set_eddy_diffs ()
void	set_diff_stresses ()
void	add_fast_tend (amrex::Vector< amrex::MultiFab > &S_rhs)
void	add_slow_tend (const amrex::MFIter &mfi, const amrex::Box &tbx, const amrex::Array4< amrex::Real > &cc_rhs_arr)

Protected Member Functions
void	check_flux_state_consistency (const double dt)
void	apply_turbulent_mountain_stress ()
void	set_computed_group_impl ()
size_t	requested_buffer_size_in_bytes () const
void	init_buffers ()
Int	wsm_num_slots () const

Protected Attributes
amrex::Vector< int >	m_col_offsets
Int	m_num_cols = 0
Int	m_num_layers = 0
Int	m_npbl
Int	m_nadv
Int	m_num_tracers = 3
Int	m_num_vel_comp = 2
Int	hdtime
int	m_lev
int	m_step
amrex::Geometry	m_geom
amrex::BoxArray	m_ba
amrex::MultiFab *	m_cons = nullptr
amrex::MultiFab *	m_xvel = nullptr
amrex::MultiFab *	m_yvel = nullptr
amrex::MultiFab *	m_zvel = nullptr
amrex::Real *	m_w_subsid = nullptr
amrex::MultiFab *	m_tau13 = nullptr
amrex::MultiFab *	m_tau23 = nullptr
amrex::MultiFab *	m_hfx3 = nullptr
amrex::MultiFab *	m_qfx3 = nullptr
amrex::MultiFab *	m_mu = nullptr
amrex::MultiFab *	m_z_phys = nullptr
amrex::MultiFab	c_tend
amrex::MultiFab	u_tend
amrex::MultiFab	v_tend
bool	m_first_step = true
Buffer	m_buffer
SHF::SHOCInput	input
SHF::SHOCInputOutput	input_output
SHF::SHOCOutput	output
SHF::SHOCHistoryOutput	history_output
SHF::SHOCRuntime	runtime_options
SHOCPreprocess	shoc_preprocess
SHOCPostprocess	shoc_postprocess
ekat::WorkspaceManager< Spack, KT::Device >	workspace_mgr
view_1d	tot_buff_view
bool	apply_tms = false
bool	check_flux_state = false
bool	extra_shoc_diags = false
bool	column_conservation_check = false
view_2d	omega
view_1d	surf_sens_flux
sview_2d	surf_mom_flux
view_1d	surf_evap
view_2d	T_mid
view_2d	qv
view_1d	surf_drag_coeff_tms
view_2d	p_mid
view_2d	p_int
view_2d	pseudo_dens
view_1d	phis
view_3d	horiz_wind
view_2d	sgs_buoy_flux
view_2d	tk
view_2d	cldfrac_liq
view_2d	tke
view_2d	qc
view_1d	pblh
view_2d	inv_qc_relvar
view_2d	tkh
view_2d	w_sec
view_2d	cldfrac_liq_prev
view_1d	ustar
view_1d	obklen
view_2d	brunt
view_2d	shoc_mix
view_2d	isotropy
view_2d	shoc_cond
view_2d	shoc_evap
view_2d	wthl_sec
view_2d	thl_sec
view_2d	wqw_sec
view_2d	qw_sec
view_2d	uw_sec
view_2d	vw_sec
view_2d	w3
view_3d_strided	qtracers
view_1d	vapor_flux
view_1d	water_flux
view_1d	ice_flux
view_1d	heat_flux

Static Protected Attributes
static constexpr Int	shoc_extra_trac_slots = 3

Private Types
using	SHF = scream::shoc::Functions< Real, KokkosDefaultDevice >
using	PF = scream::PhysicsFunctions< KokkosDefaultDevice >
using	C = scream::physics::Constants< Real >
using	KT = ekat::KokkosTypes< KokkosDefaultDevice >
using	SC = scream::shoc::Constants< Real >
using	Spack = typename SHF::Spack
using	IntSmallPack = typename SHF::IntSmallPack
using	Smask = typename SHF::Smask
using	view_1d_int = typename KT::template view_1d< Int >
using	view_1d = typename SHF::view_1d< Real >
using	view_1d_const = typename SHF::view_1d< const Real >
using	view_2d = typename SHF::view_2d< SHF::Spack >
using	view_2d_const = typename SHF::view_2d< const Spack >
using	sview_2d = typename KT::template view_2d< Real >
using	sview_2d_const = typename KT::template view_2d< const Real >
using	view_3d = typename SHF::view_3d< Spack >
using	view_3d_const = typename SHF::view_3d< const Spack >
using	view_3d_strided = typename SHF::view_3d_strided< Spack >
using	WSM = ekat::WorkspaceManager< Spack, KT::Device >
template<typename ScalarT >
using	uview_1d = ekat::Unmanaged< typename KT::template view_1d< ScalarT > >
template<typename ScalarT >
using	uview_2d = ekat::Unmanaged< typename KT::template view_2d< ScalarT > >

Member Typedef Documentation

C

using SHOCInterface::C = scream::physics::Constants<Real>

private

IntSmallPack

using SHOCInterface::IntSmallPack = typename SHF::IntSmallPack

private

KT

using SHOCInterface::KT = ekat::KokkosTypes<KokkosDefaultDevice>

private

PF

using SHOCInterface::PF = scream::PhysicsFunctions<KokkosDefaultDevice>

private

SC

using SHOCInterface::SC = scream::shoc::Constants<Real>

private

SHF

using SHOCInterface::SHF = scream::shoc::Functions<Real, KokkosDefaultDevice>

private

Smask

using SHOCInterface::Smask = typename SHF::Smask

private

Spack

using SHOCInterface::Spack = typename SHF::Spack

private

sview_2d

using SHOCInterface::sview_2d = typename KT::template view_2d<Real>

private

sview_2d_const

using SHOCInterface::sview_2d_const = typename KT::template view_2d<const Real>

private

uview_1d

template<typename ScalarT >

using SHOCInterface::uview_1d = ekat::Unmanaged<typename KT::template view_1d<ScalarT> >

private

uview_2d

template<typename ScalarT >

using SHOCInterface::uview_2d = ekat::Unmanaged<typename KT::template view_2d<ScalarT> >

private

view_1d

using SHOCInterface::view_1d = typename SHF::view_1d<Real>

private

view_1d_const

using SHOCInterface::view_1d_const = typename SHF::view_1d<const Real>

private

view_1d_int

using SHOCInterface::view_1d_int = typename KT::template view_1d<Int>

private

view_2d

using SHOCInterface::view_2d = typename SHF::view_2d<SHF::Spack>

private

view_2d_const

using SHOCInterface::view_2d_const = typename SHF::view_2d<const Spack>

private

view_3d

using SHOCInterface::view_3d = typename SHF::view_3d<Spack>

private

view_3d_const

using SHOCInterface::view_3d_const = typename SHF::view_3d<const Spack>

private

view_3d_strided

using SHOCInterface::view_3d_strided = typename SHF::view_3d_strided<Spack>

private

WSM

using SHOCInterface::WSM = ekat::WorkspaceManager<Spack, KT::Device>

private

Constructor & Destructor Documentation

SHOCInterface()

SHOCInterface::SHOCInterface	(	const int &	lev,
		SolverChoice &	sc
	)

{
    //
    // Defaults set from E3SM/components/eamxx/cime_config/namelist_defaults_eamxx.xml
    //
    // Turn off SGS variability in SHOC, effectively reducing it to a Real(1.5) TKE closure?
    bool def_shoc_1p5tke = false;
    // Minimum value of stability correction
    Real def_lambda_low = Real(0.001);
    // Maximum value of stability correction
    Real def_lambda_high = Real(0.04);
    // Slope of change from lambda_low to lambda_high
    Real def_lambda_slope = Real(2.65);
    // stability threshold for which to apply more stability correction
    Real def_lambda_thresh = Real(0.02);
    // Temperature variance tuning factor
    Real def_thl2tune = one;
    // Moisture variance tuning factor
    Real def_qw2tune = one;
    // Temperature moisture covariance
    Real def_qwthl2tune = one;
    // Vertical velocity variance
    Real def_w2tune = one;
    // Length scale factor
    Real def_length_fac = Real(0.5);
    // Third moment vertical velocity damping factor
    Real def_c_diag_3rd_mom = Real(7.0);
    // Eddy diffusivity coefficient for heat
    Real def_coeff_kh = Real(0.1);
    // Eddy diffusivity coefficient for momentum
    Real def_coeff_km = Real(0.1);
 
    runtime_options.lambda_low    = def_lambda_low;
    runtime_options.lambda_high   = def_lambda_high;
    runtime_options.lambda_slope  = def_lambda_slope;
    runtime_options.lambda_thresh = def_lambda_thresh;
 
    runtime_options.thl2tune   = def_thl2tune;
    runtime_options.qwthl2tune = def_qwthl2tune;
    runtime_options.qw2tune    = def_qw2tune;
    runtime_options.w2tune     = def_w2tune;
 
    runtime_options.length_fac     = def_length_fac;
    runtime_options.c_diag_3rd_mom = def_c_diag_3rd_mom;
    runtime_options.Ckh            = def_coeff_kh;
    runtime_options.Ckm            = def_coeff_km;
    runtime_options.shoc_1p5tke    = def_shoc_1p5tke;
    runtime_options.extra_diags    = extra_shoc_diags;
 
    // Construct parser object for following reads
    ParmParse pp("erf.shoc");
 
    // Parse runtime inputs at start up
    pp.queryAdd("lambda_low"      , runtime_options.lambda_low    );
    pp.queryAdd("lambda_high"     , runtime_options.lambda_high   );
    pp.queryAdd("lambda_slope"    , runtime_options.lambda_slope  );
    pp.queryAdd("lambda_thresh"   , runtime_options.lambda_thresh );
    pp.queryAdd("thl2tune"        , runtime_options.thl2tune      );
    pp.queryAdd("qw2tune"         , runtime_options.qw2tune       );
    pp.queryAdd("qwthl2tune"      , runtime_options.qwthl2tune    );
    pp.queryAdd("w2tune"          , runtime_options.w2tune        );
    pp.queryAdd("length_fac"      , runtime_options.length_fac    );
    pp.queryAdd("c_diag_3rd_mom"  , runtime_options.c_diag_3rd_mom);
    pp.queryAdd("coeff_kh"        , runtime_options.Ckh           );
    pp.queryAdd("coeff_km"        , runtime_options.Ckm           );
    pp.queryAdd("shoc_1p5tke"     , runtime_options.shoc_1p5tke   );
    pp.queryAdd("extra_shoc_diags", runtime_options.extra_diags   );
 
    // Set to default but allow us to change it through the inputs file
    pp.queryAdd("apply_tms", apply_tms);
    pp.queryAdd("check_flux_state", check_flux_state);
    pp.queryAdd("extra_shoc_diags", extra_shoc_diags);
    pp.queryAdd("column_conservation_check", column_conservation_check);
}

Here is the call graph for this function:

Member Function Documentation

add_fast_tend()

void SHOCInterface::add_fast_tend

(

amrex::Vector< amrex::MultiFab > &

S_rhs

)

{
    for (MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
        const auto& vbx_cc  = mfi.validbox();
        const auto& vbx_x   = convert(vbx_cc,IntVect(1,0,0));
        const auto& vbx_y   = convert(vbx_cc,IntVect(0,1,0));
 
        const Array4<const Real>& c_arr = m_cons->const_array(mfi);
 
        const Array4<Real>& cc_rhs_arr = S_rhs[IntVars::cons].array(mfi);
        const Array4<Real>& ru_rhs_arr = S_rhs[IntVars::xmom].array(mfi);
        const Array4<Real>& rv_rhs_arr = S_rhs[IntVars::ymom].array(mfi);
 
        const Array4<const Real>& c_tend_arr = c_tend.const_array(mfi);
        const Array4<const Real>& u_tend_arr = u_tend.const_array(mfi);
        const Array4<const Real>& v_tend_arr = v_tend.const_array(mfi);
 
        ParallelFor(vbx_cc, vbx_x, vbx_y,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            cc_rhs_arr(i,j,k,RhoTheta_comp) += c_arr(i,j,k,Rho_comp) * c_tend_arr(i,j,k,RhoTheta_comp);
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            Real rho_on_u_face = Real(0.5) * (c_arr(i-1,j,k,Rho_comp)
                                              + c_arr(i,j,k,Rho_comp));
            ru_rhs_arr(i,j,k) += rho_on_u_face * u_tend_arr(i,j,k);
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            Real rho_on_v_face = Real(0.5) * (c_arr(i,j-1,k,Rho_comp)
                                              + c_arr(i,j,k,Rho_comp));
            rv_rhs_arr(i,j,k) += rho_on_v_face * v_tend_arr(i,j,k);
        });
    }
}

Here is the call graph for this function:

add_slow_tend()

void SHOCInterface::add_slow_tend	(	const amrex::MFIter &	mfi,
		const amrex::Box &	tbx,
		const amrex::Array4< amrex::Real > &	cc_rhs_arr
	)

{
    bool moist = (m_cons->nComp() > RhoQ1_comp);
 
    const Array4<const Real>& c_arr = m_cons->const_array(mfi);
 
    const Array4<const Real>& c_tend_arr = c_tend.const_array(mfi);
 
    ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
    {
        cc_rhs_arr(i,j,k,RhoKE_comp) += c_arr(i,j,k,Rho_comp) * c_tend_arr(i,j,k,RhoKE_comp);
        if (moist) {
            cc_rhs_arr(i,j,k,RhoQ1_comp) += c_arr(i,j,k,Rho_comp) * c_tend_arr(i,j,k,RhoQ1_comp);
            cc_rhs_arr(i,j,k,RhoQ2_comp) += c_arr(i,j,k,Rho_comp) * c_tend_arr(i,j,k,RhoQ2_comp);
        }
    });
}

Here is the call graph for this function:

alloc_buffers()

void SHOCInterface::alloc_buffers

(

)

{
    // Interface data structures
    //=======================================================
    omega            = view_2d("Omega"             , m_num_cols, m_num_layers  );
    surf_sens_flux   = view_1d("Sfc sens flux"     , m_num_cols);
    surf_mom_flux    = sview_2d("Sfc mom flux"     , m_num_cols, m_num_vel_comp);
    surf_evap        = view_1d("Sfc evap"          , m_num_cols);
    T_mid            = view_2d("T_mid"             , m_num_cols, m_num_layers  );
    qv               = view_2d("Qv"                , m_num_cols, m_num_layers  );
    surf_drag_coeff_tms = view_1d("surf_drag_coeff", m_num_cols);
 
    // Input data structures
    //=======================================================
    p_mid            = view_2d("P_mid"             , m_num_cols, m_num_layers  );
    p_int            = view_2d("P_int"             , m_num_cols, m_num_layers+1);
    pseudo_dens      = view_2d("Pseudo density"    , m_num_cols, m_num_layers  );
    phis             = view_1d("Phis"              , m_num_cols);
 
    // Input/Output data structures
    //=======================================================
    horiz_wind       = view_3d("horiz_wind"        , m_num_cols, m_num_vel_comp, m_num_layers);
    cldfrac_liq      = view_2d("Cld_frac_liq"      , m_num_cols, m_num_layers  );
    tke              = view_2d("Tke"               , m_num_cols, m_num_layers  );
    qc               = view_2d("Qc"                , m_num_cols, m_num_layers  );
 
    // Output data structures
    //=======================================================
    pblh             = view_1d("pbl_height"        , m_num_cols);
    inv_qc_relvar    = view_2d("inv_qc_relvar"     , m_num_cols, m_num_layers);
    tkh              = view_2d("eddy_diff_heat"    , m_num_cols, m_num_layers);
    w_sec            = view_2d("w_sec"             , m_num_cols, m_num_layers);
    cldfrac_liq_prev = view_2d("cld_frac_liq_prev" , m_num_cols, m_num_layers  );
    ustar            = view_1d("ustar"             , m_num_cols);
    obklen           = view_1d("obklen"            , m_num_cols);
 
    // Extra diagnostic data structures
    //=======================================================
    if (extra_shoc_diags) {
        brunt            = view_2d("brunt"    , m_num_cols, m_num_layers);
        shoc_mix         = view_2d("shoc_mix" , m_num_cols, m_num_layers);
        isotropy         = view_2d("isotropy" , m_num_cols, m_num_layers);
        shoc_cond        = view_2d("shoc_cond", m_num_cols, m_num_layers);
        shoc_evap        = view_2d("shoc_evap", m_num_cols, m_num_layers);
 
        wthl_sec         = view_2d("wthl_sec" , m_num_cols, m_num_layers+1);
        thl_sec          = view_2d("thl_sec"  , m_num_cols, m_num_layers+1);
        wqw_sec          = view_2d("wqw_sec"  , m_num_cols, m_num_layers+1);
        qw_sec           = view_2d("qw_sec"   , m_num_cols, m_num_layers+1);
        uw_sec           = view_2d("uw_sec"   , m_num_cols, m_num_layers+1);
        vw_sec           = view_2d("vw_sec"   , m_num_cols, m_num_layers+1);
        w3               = view_2d("w3"       , m_num_cols, m_num_layers+1);
    }
 
    // Tracer data structures
    //=======================================================
    // NOTE: Use layoutright format
    Kokkos::LayoutStride layout(m_num_cols   , m_num_tracers*m_num_layers, // stride for dim0
                                m_num_tracers, m_num_layers,               // stride for dim1
                                m_num_layers , 1                           // stride for dim2
                                );
    qtracers             = view_3d_strided("Qtracers"  , layout);
 
    // Boundary flux data structures
    //=======================================================
    if (column_conservation_check) {
        vapor_flux       = view_1d("vapor_flux", m_num_cols);
        water_flux       = view_1d("water_flux", m_num_cols);
        ice_flux         = view_1d("ice_flux"  , m_num_cols);
        heat_flux        = view_1d("heat_flux" , m_num_cols);
    }
}

apply_turbulent_mountain_stress()

void SHOCInterface::apply_turbulent_mountain_stress ( )

protected

{
    auto rrho_i   = m_buffer.rrho_i;
    auto upwp_sfc = m_buffer.upwp_sfc;
    auto vpwp_sfc = m_buffer.vpwp_sfc;
    auto surf_drag_coeff_tms_d = surf_drag_coeff_tms;
    auto horiz_wind_d = horiz_wind;
 
    const int nlev_v  = (m_num_layers-1)/Spack::n;
    const int nlev_p  = (m_num_layers-1)%Spack::n;
    const int nlevi_v = m_num_layers/Spack::n;
    const int nlevi_p = m_num_layers%Spack::n;
 
    Kokkos::parallel_for("apply_tms", KT::RangePolicy(0, m_num_cols), KOKKOS_LAMBDA (const int i)
    {
        upwp_sfc(i) -= surf_drag_coeff_tms_d(i)*horiz_wind_d(i,0,nlev_v)[nlev_p]/rrho_i(i,nlevi_v)[nlevi_p];
        vpwp_sfc(i) -= surf_drag_coeff_tms_d(i)*horiz_wind_d(i,1,nlev_v)[nlev_p]/rrho_i(i,nlevi_v)[nlevi_p];
    });
}

check_flux_state_consistency()

void SHOCInterface::check_flux_state_consistency ( const double  dt )

protected

{
    using PC  = scream::physics::Constants<Real>;
    using RU  = ekat::ReductionUtils<KT::ExeSpace>;
    using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
 
    const Real gravit = PC::gravit;
    const Real qmin   = 1e-12; // minimum permitted constituent concentration (kg/kg)
 
    const auto& pseudo_density = pseudo_dens;
    auto qv_d = qv;
    auto surf_evap_d = surf_evap;
 
    const auto nlevs           = m_num_layers;
    const auto nlev_packs      = ekat::npack<Spack>(nlevs);
    const auto last_pack_idx   = (nlevs-1)/Spack::n;
    const auto last_pack_entry = (nlevs-1)%Spack::n;
    const auto policy          = TPF::get_default_team_policy(m_num_cols, nlev_packs);
    Kokkos::parallel_for("check_flux_state_consistency",
                         policy,
                         KOKKOS_LAMBDA (const KT::MemberType& team)
    {
        const auto i = team.league_rank();
 
        const auto& pseudo_density_i = ekat::subview(pseudo_density, i);
        const auto& qv_i             = ekat::subview(qv_d, i);
 
        // reciprocal of pseudo_density at the bottom layer
        const auto rpdel = one/pseudo_density_i(last_pack_idx)[last_pack_entry];
 
        // Check if the negative surface latent heat flux can exhaust
        // the moisture in the lowest model level. If so, apply fixer.
        const auto condition = surf_evap_d(i) - (qmin - qv_i(last_pack_idx)[last_pack_entry])/(dt*gravit*rpdel);
        if (condition < 0) {
            const auto cc = std::fabs(surf_evap_d(i)*dt*gravit);
 
            auto tracer_mass = [&](const int k)
            {
                return qv_i(k)*pseudo_density_i(k);
            };
            Real mm = RU::view_reduction(team, 0, nlevs, tracer_mass);
 
            EKAT_KERNEL_ASSERT_MSG(mm >= cc, "Error! Total mass of column vapor should be greater than mass of surf_evap.\n");
 
            Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_packs), [&](const int& k)
            {
                const auto adjust = cc*qv_i(k)*pseudo_density_i(k)/mm;
                qv_i(k) = (qv_i(k)*pseudo_density_i(k) - adjust)/pseudo_density_i(k);
            });
 
            surf_evap_d(i) = 0;
        }
    });
}

dealloc_buffers()

void SHOCInterface::dealloc_buffers

(

)

{
    // Contiguous memory buffer view
    //=======================================================
    tot_buff_view    = view_1d();
 
    // Interface data structures
    //=======================================================
    omega            = view_2d();
    surf_sens_flux   = view_1d();
    surf_mom_flux    = sview_2d();
    surf_evap        = view_1d();
    T_mid            = view_2d();
    qv               = view_2d();
    surf_drag_coeff_tms = view_1d();
 
    // Input data structures
    //=======================================================
    p_mid            = view_2d();
    p_int            = view_2d();
    pseudo_dens      = view_2d();
    phis             = view_1d();
 
    // Input/Output data structures
    //=======================================================
    horiz_wind       = view_3d();
    cldfrac_liq      = view_2d();
    tke              = view_2d();
    qc               = view_2d();
 
    // Output data structures
    //=======================================================
    pblh             = view_1d();
    inv_qc_relvar    = view_2d();
    tkh              = view_2d();
    w_sec            = view_2d();
    cldfrac_liq_prev = view_2d();
    ustar            = view_1d();
    obklen           = view_1d();
 
    // Extra diagnostic data structures
    //=======================================================
    if (extra_shoc_diags) {
        brunt            = view_2d();
        shoc_mix         = view_2d();
        isotropy         = view_2d();
        shoc_cond        = view_2d();
        shoc_evap        = view_2d();
 
        wthl_sec         = view_2d();
        thl_sec          = view_2d();
        wqw_sec          = view_2d();
        qw_sec           = view_2d();
        uw_sec           = view_2d();
        vw_sec           = view_2d();
        w3               = view_2d();
    }
 
    // Tracer data structures
    //=======================================================
    qtracers             = view_3d_strided();
 
    // Boundary flux data structures
    //=======================================================
    if (column_conservation_check) {
        vapor_flux       = view_1d();
        water_flux       = view_1d();
        ice_flux         = view_1d();
        heat_flux        = view_1d();
    }
}

finalize_impl()

void SHOCInterface::finalize_impl

(

const double

dt

)

{
    // Do nothing (per SHOCMacrophysics::finalize_impl())
 
    // Fill the AMReX MFs from Kokkos Views
    kokkos_buffers_to_mf(dt);
 
    // Deallocate the buffer arrays
    dealloc_buffers();
}

init_buffers()

void SHOCInterface::init_buffers ( )

protected

{
 
    // Buffer of contiguous memory
    auto buffer_size = requested_buffer_size_in_bytes();
    tot_buff_view = view_1d("contiguous shoc_buffer",buffer_size);
    Real* mem = reinterpret_cast<Real*>(tot_buff_view.data());
 
    // 1d scalar views
    using scalar_view_t = decltype(m_buffer.wpthlp_sfc);
    scalar_view_t* _1d_scalar_view_ptrs[Buffer::num_1d_scalar_ncol] =
        {&m_buffer.wpthlp_sfc, &m_buffer.wprtp_sfc, &m_buffer.upwp_sfc, &m_buffer.vpwp_sfc
#ifdef SCREAM_SHOC_SMALL_KERNELS
         , &m_buffer.se_b, &m_buffer.ke_b, &m_buffer.wv_b, &m_buffer.wl_b
         , &m_buffer.se_a, &m_buffer.ke_a, &m_buffer.wv_a, &m_buffer.wl_a
         , &m_buffer.kbfs, &m_buffer.ustar2, &m_buffer.wstar
#endif
        };
    for (int i = 0; i < Buffer::num_1d_scalar_ncol; ++i) {
        *_1d_scalar_view_ptrs[i] = scalar_view_t(mem, m_num_cols);
        mem += _1d_scalar_view_ptrs[i]->size();
    }
 
    Spack* s_mem = reinterpret_cast<Spack*>(mem);
 
    // 2d packed views
    const int nlev_packs       = ekat::npack<Spack>(m_num_layers);
    const int nlevi_packs      = ekat::npack<Spack>(m_num_layers+1);
    const int num_tracer_packs = ekat::npack<Spack>(m_num_tracers);
 
    m_buffer.pref_mid = decltype(m_buffer.pref_mid)(s_mem, nlev_packs);
    s_mem += m_buffer.pref_mid.size();
 
    using spack_2d_view_t = decltype(m_buffer.z_mid);
    spack_2d_view_t* _2d_spack_mid_view_ptrs[Buffer::num_2d_vector_mid] =
        {&m_buffer.z_mid, &m_buffer.rrho, &m_buffer.thv, &m_buffer.dz, &m_buffer.zt_grid, &m_buffer.wm_zt, &m_buffer.unused,
         &m_buffer.inv_exner, &m_buffer.thlm, &m_buffer.qw, &m_buffer.dse, &m_buffer.tke_copy, &m_buffer.qc_copy,
         &m_buffer.shoc_ql2, &m_buffer.shoc_mix, &m_buffer.isotropy, &m_buffer.w_sec, &m_buffer.wqls_sec, &m_buffer.brunt
#ifdef SCREAM_SHOC_SMALL_KERNELS
         , &m_buffer.rho_zt, &m_buffer.shoc_qv, &m_buffer.tabs, &m_buffer.dz_zt
#endif
        };
 
    spack_2d_view_t* _2d_spack_int_view_ptrs[Buffer::num_2d_vector_int] =
        {&m_buffer.z_int, &m_buffer.rrho_i, &m_buffer.zi_grid, &m_buffer.thl_sec, &m_buffer.qw_sec,
         &m_buffer.qwthl_sec, &m_buffer.wthl_sec, &m_buffer.wqw_sec, &m_buffer.wtke_sec, &m_buffer.uw_sec,
         &m_buffer.vw_sec, &m_buffer.w3
#ifdef SCREAM_SHOC_SMALL_KERNELS
         , &m_buffer.dz_zi
#endif
        };
 
    for (int i = 0; i < Buffer::num_2d_vector_mid; ++i) {
        *_2d_spack_mid_view_ptrs[i] = spack_2d_view_t(s_mem, m_num_cols, nlev_packs);
        s_mem += _2d_spack_mid_view_ptrs[i]->size();
    }
 
    for (int i = 0; i < Buffer::num_2d_vector_int; ++i) {
        *_2d_spack_int_view_ptrs[i] = spack_2d_view_t(s_mem, m_num_cols, nlevi_packs);
        s_mem += _2d_spack_int_view_ptrs[i]->size();
    }
    m_buffer.wtracer_sfc = decltype(m_buffer.wtracer_sfc)(s_mem, m_num_cols, num_tracer_packs);
    s_mem += m_buffer.wtracer_sfc.size();
 
    // WSM data
    m_buffer.wsm_data = s_mem;
}

initialize_impl()

void SHOCInterface::initialize_impl

(

)

{
    // Alias local variables from temporary buffer
    auto z_mid       = m_buffer.z_mid;
    auto z_int       = m_buffer.z_int;
    auto wpthlp_sfc  = m_buffer.wpthlp_sfc;
    auto wprtp_sfc   = m_buffer.wprtp_sfc;
    auto upwp_sfc    = m_buffer.upwp_sfc;
    auto vpwp_sfc    = m_buffer.vpwp_sfc;
    auto rrho        = m_buffer.rrho;
    auto rrho_i      = m_buffer.rrho_i;
    auto thv         = m_buffer.thv;
    auto dz          = m_buffer.dz;
    auto zt_grid     = m_buffer.zt_grid;
    auto zi_grid     = m_buffer.zi_grid;
    auto wtracer_sfc = m_buffer.wtracer_sfc;
    auto wm_zt       = m_buffer.wm_zt;
    auto inv_exner   = m_buffer.inv_exner;
    auto thlm        = m_buffer.thlm;
    auto qw          = m_buffer.qw;
    auto dse         = m_buffer.dse;
    auto tke_copy    = m_buffer.tke_copy;
    auto qc_copy     = m_buffer.qc_copy;
    auto shoc_ql2    = m_buffer.shoc_ql2;
 
    // For now, set z_int(i,nlevs) = z_surf = 0
    const Real z_surf = Real(0.);
 
    // Set preprocess variables
    shoc_preprocess.set_variables(m_num_cols, m_num_layers, z_surf,
                                  T_mid, p_mid, p_int, pseudo_dens, omega, phis, surf_sens_flux, surf_evap,
                                  surf_mom_flux, qtracers, qv, qc, qc_copy, tke, tke_copy, z_mid, z_int,
                                  dse, rrho, rrho_i, thv, dz, zt_grid, zi_grid, wpthlp_sfc, wprtp_sfc, upwp_sfc,
                                  vpwp_sfc, wtracer_sfc, wm_zt, inv_exner, thlm, qw, cldfrac_liq, cldfrac_liq_prev);
 
    // Input Variables:
    input.zt_grid     = shoc_preprocess.zt_grid;
    input.zi_grid     = shoc_preprocess.zi_grid;
    input.pres        = p_mid;
    input.presi       = p_int;
    input.pdel        = pseudo_dens;
    input.thv         = shoc_preprocess.thv;
    input.w_field     = shoc_preprocess.wm_zt;
    input.wthl_sfc    = shoc_preprocess.wpthlp_sfc;
    input.wqw_sfc     = shoc_preprocess.wprtp_sfc;
    input.uw_sfc      = shoc_preprocess.upwp_sfc;
    input.vw_sfc      = shoc_preprocess.vpwp_sfc;
    input.wtracer_sfc = shoc_preprocess.wtracer_sfc;
    input.inv_exner   = shoc_preprocess.inv_exner;
    input.phis        = phis;
 
    // Input/Output Variables
    input_output.host_dse     = shoc_preprocess.shoc_s;
    input_output.tke          = shoc_preprocess.tke_copy;
    input_output.thetal       = shoc_preprocess.thlm;
    input_output.qw           = shoc_preprocess.qw;
    input_output.horiz_wind   = horiz_wind;
    input_output.wthv_sec     = sgs_buoy_flux;
    input_output.qtracers     = shoc_preprocess.qtracers;
    input_output.tk           = tk;
    input_output.shoc_cldfrac = cldfrac_liq;
    input_output.shoc_ql      = qc_copy;
 
    // Output Variables
    output.pblh     = pblh;
    output.shoc_ql2 = shoc_ql2;
    output.tkh      = tkh;
    output.ustar    = ustar;
    output.obklen   = obklen;
 
    // Output (diagnostic)
    history_output.shoc_mix  = m_buffer.shoc_mix;
    history_output.isotropy  = m_buffer.isotropy;
    if (extra_shoc_diags) {
        history_output.shoc_cond = shoc_cond;
        history_output.shoc_evap = shoc_evap;
    } else {
        history_output.shoc_cond = m_buffer.unused;
        history_output.shoc_evap = m_buffer.unused;
    }
    history_output.w_sec     = w_sec;
    history_output.thl_sec   = m_buffer.thl_sec;
    history_output.qw_sec    = m_buffer.qw_sec;
    history_output.qwthl_sec = m_buffer.qwthl_sec;
    history_output.wthl_sec  = m_buffer.wthl_sec;
    history_output.wqw_sec   = m_buffer.wqw_sec;
    history_output.wtke_sec  = m_buffer.wtke_sec;
    history_output.uw_sec    = m_buffer.uw_sec;
    history_output.vw_sec    = m_buffer.vw_sec;
    history_output.w3        = m_buffer.w3;
    history_output.wqls_sec  = m_buffer.wqls_sec;
    history_output.brunt     = m_buffer.brunt;
 
#ifdef SCREAM_SHOC_SMALL_KERNELS
    temporaries.se_b = m_buffer.se_b;
    temporaries.ke_b = m_buffer.ke_b;
    temporaries.wv_b = m_buffer.wv_b;
    temporaries.wl_b = m_buffer.wl_b;
    temporaries.se_a = m_buffer.se_a;
    temporaries.ke_a = m_buffer.ke_a;
    temporaries.wv_a = m_buffer.wv_a;
    temporaries.wl_a = m_buffer.wl_a;
    temporaries.kbfs = m_buffer.kbfs;
    temporaries.ustar2 = m_buffer.ustar2;
    temporaries.wstar = m_buffer.wstar;
 
    temporaries.rho_zt = m_buffer.rho_zt;
    temporaries.shoc_qv = m_buffer.shoc_qv;
    temporaries.tabs = m_buffer.tabs;
    temporaries.dz_zt = m_buffer.dz_zt;
    temporaries.dz_zi = m_buffer.dz_zi;
#endif
 
    // Set postprocess variables
    shoc_postprocess.set_variables(m_num_cols, m_num_layers,
                                   rrho, qv, qw, qc, qc_copy, tke,
                                   tke_copy, qtracers, shoc_ql2,
                                   cldfrac_liq, inv_qc_relvar,
                                   T_mid, dse, z_mid, phis);
 
    // Setup WSM for internal local variables
    using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
    const auto nlev_packs  = ekat::npack<Spack>(m_num_layers);
    const auto nlevi_packs = ekat::npack<Spack>(m_num_layers+1);
    const auto default_policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
    const Int  num_slots      = wsm_num_slots();
 
    // The WSM data sits at the tail of the contiguous buffer, so the slots we
    // hand it here must not exceed what requested_buffer_size_in_bytes() set aside.
    const size_t wsm_offset = static_cast<size_t>(reinterpret_cast<Real*>(m_buffer.wsm_data)
                                                  - tot_buff_view.data());
    const size_t wsm_reals  = WSM::get_total_bytes_needed(nlevi_packs, num_slots,
                                                          default_policy)/sizeof(Real);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(wsm_offset + wsm_reals <= tot_buff_view.size(),
                                     "SHOC WSM slots exceed the reserved buffer; "
                                     "wsm_num_slots() and requested_buffer_size_in_bytes() disagree");
 
    workspace_mgr.setup(m_buffer.wsm_data, nlevi_packs, num_slots, default_policy);
 
    // NOTE: Vertical indices were permuted, so top and bottom are correct
    // Maximum number of levels in pbl from surface
    const int ntop_shoc = 0;
    const int nbot_shoc = m_num_layers;
    auto p_mid_d = p_mid;
    view_1d pref_mid("pref_mid", m_num_layers);
    Spack* s_mem = reinterpret_cast<Spack*>(pref_mid.data());
    SHF::view_1d<Spack> pref_mid_um(s_mem, m_num_layers);
    const auto policy     =  TPF::get_default_team_policy(m_num_cols, nlev_packs);
    Kokkos::parallel_for("pref_mid",
                         policy,
                         KOKKOS_LAMBDA (const KT::MemberType& team)
    {
        const auto i = team.league_rank();
        if (i==0) {
            const auto& pmid_i = ekat::subview(p_mid_d, i);
            Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_packs), [&](const int& k)
            {
                pref_mid_um(k) = pmid_i(k);
            });
        }
    });
    Kokkos::fence();
    m_npbl = SHF::shoc_init(nbot_shoc,ntop_shoc,pref_mid_um);
 
    // Cell length for input dx and dy
    view_1d cell_length_x("cell_length_x", m_num_cols);
    view_1d cell_length_y("cell_length_y", m_num_cols);
    Kokkos::deep_copy(cell_length_x, m_geom.CellSize(0));
    Kokkos::deep_copy(cell_length_y, m_geom.CellSize(1));
    input.dx = cell_length_x;
    input.dy = cell_length_y;
}

Here is the call graph for this function:

kokkos_buffers_to_mf()

void SHOCInterface::kokkos_buffers_to_mf

(

const double

dt

)

{
    //
    // Expose for device capture
    //
 
    // Buffer data structures
    //=======================================================
    auto thlm_d = m_buffer.thlm;
 
    // Interface data structures
    //=======================================================
    auto T_mid_d = T_mid;
    auto qv_d = qv;
 
    // Input/Output data structures
    //=======================================================
    auto horiz_wind_d = horiz_wind;
    auto tke_d = tke;
    auto qc_d  = qc;
 
    bool moist = (m_cons->nComp() > RhoQ1_comp);
    for (MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
        // NOTE: No ghost cells when going back to MFs
        const auto& vbx_cc = mfi.validbox();
        const auto& vbx_x  = convert(vbx_cc,IntVect(1,0,0));
        const auto& vbx_y  = convert(vbx_cc,IntVect(0,1,0));
 
        // NOTE: Grown box only for mapping
        const auto& gbx  = mfi.tilebox(IntVect(0,0,0),IntVect(1,1,0));
        const int nx     = gbx.length(0);
        const int imin   = gbx.smallEnd(0);
        const int jmin   = gbx.smallEnd(1);
        const int kmax   = gbx.bigEnd(2);
        const int offset = m_col_offsets[mfi.index()];
 
        const Array4<const Real>& cons_arr = m_cons->const_array(mfi);
        const Array4<const Real>& u_arr    = m_xvel->const_array(mfi);
        const Array4<const Real>& v_arr    = m_yvel->const_array(mfi);
 
        const Array4<Real>& c_tend_arr     = c_tend.array(mfi);
        const Array4<Real>& u_tend_arr     = u_tend.array(mfi);
        const Array4<Real>& v_tend_arr     = v_tend.array(mfi);
 
        // Interpolate SHOC cell-centered wind increments to faces. Reconstructing the
        // face winds first would apply a horizontal filter even for a zero increment.
        ParallelFor(vbx_cc, vbx_x, vbx_y,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // NOTE: k gets permuted with ilay
            // map [i,j,k] 0-based to [icol, ilay] 0-based
            const int icol   = (j-jmin)*nx + (i-imin) + offset;
            const int ilay   = kmax - k;
 
            // Density at CC
            Real r  = cons_arr(i,j,k,Rho_comp);
 
            // Theta at CC (eamxx_common_physics_functions_impl.hpp L123)
            Real Th = thlm_d(icol,ilay)[0] / ( one - (one / T_mid_d(icol,ilay)[0])
                                             * (C::LatVap/C::Cpair) * qc_d(icol,ilay)[0] );
 
            // Populate the tendencies
            c_tend_arr(i,j,k,RhoTheta_comp)  = ( Th                  - cons_arr(i,j,k,RhoTheta_comp)/r ) / dt;
            c_tend_arr(i,j,k,RhoKE_comp)     = ( tke_d(icol,ilay)[0] - cons_arr(i,j,k,RhoKE_comp   )/r ) / dt;
            if (moist) {
                c_tend_arr(i,j,k,RhoQ1_comp) = ( qv_d(icol,ilay)[0] - cons_arr(i,j,k,RhoQ1_comp)/r ) / dt;
                c_tend_arr(i,j,k,RhoQ2_comp) = ( qc_d(icol,ilay)[0] - cons_arr(i,j,k,RhoQ2_comp)/r ) / dt;
            }
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // NOTE: k gets permuted with ilay
            // map [i,j,k] 0-based to [icol, ilay] 0-based
            const int icol   = (j-jmin)*nx + (i-imin) + offset;
            const int ilay   = kmax - k;
 
            int icolim = (j-jmin)*nx + (i-1-imin) + offset;
            Real du_cc = horiz_wind_d(icol,0,ilay)[0]
                       - Real(0.5) * (u_arr(i,j,k) + u_arr(i+1,j,k));
            Real du_cc_im = horiz_wind_d(icolim,0,ilay)[0]
                          - Real(0.5) * (u_arr(i-1,j,k) + u_arr(i,j,k));
            u_tend_arr(i,j,k) = Real(0.5) * (du_cc_im + du_cc) / dt;
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // NOTE: k gets permuted with ilay
            // map [i,j,k] 0-based to [icol, ilay] 0-based
            const int icol   = (j-jmin)*nx + (i-imin) + offset;
            const int ilay   = kmax - k;
 
            int icoljm = (j-1-jmin)*nx + (i-imin) + offset;
            Real dv_cc = horiz_wind_d(icol,1,ilay)[0]
                       - Real(0.5) * (v_arr(i,j,k) + v_arr(i,j+1,k));
            Real dv_cc_jm = horiz_wind_d(icoljm,1,ilay)[0]
                          - Real(0.5) * (v_arr(i,j-1,k) + v_arr(i,j,k));
            v_tend_arr(i,j,k) = Real(0.5) * (dv_cc_jm + dv_cc) / dt;
        });
    }
}

Here is the call graph for this function:

mf_to_kokkos_buffers()

void SHOCInterface::mf_to_kokkos_buffers

(

)

{
    // FillBoundary for internal ghost cells for u/v averaging
    m_tau13->FillBoundary(m_geom.periodicity());
    m_tau23->FillBoundary(m_geom.periodicity());
    m_hfx3->FillBoundary(m_geom.periodicity());
    m_qfx3->FillBoundary(m_geom.periodicity());
 
    //
    // Expose for device capture
    //
 
    // Interface data structures
    //=======================================================
    auto omega_d = omega;
    auto surf_sens_flux_d = surf_sens_flux;
    auto surf_mom_flux_d = surf_mom_flux;
    auto surf_evap_d = surf_evap;
    auto T_mid_d = T_mid;
    auto qv_d = qv;
    auto surf_drag_coeff_tms_d = surf_drag_coeff_tms;
 
    // Input data structures
    //=======================================================
    auto p_mid_d = p_mid;
    auto p_int_d = p_int;
    auto pseudo_dens_d = pseudo_dens;
    auto phis_d = phis;
 
    // Input/Output data structures
    //=======================================================
    auto horiz_wind_d = horiz_wind;
    auto cldfrac_liq_d = cldfrac_liq;
    auto tke_d = tke;
    auto qc_d = qc;
 
    // Enforce the correct grid heights and density
    //=======================================================
    auto dz_d = m_buffer.dz;
 
    // Subsidence pointer to device vector data
    Real* w_sub = m_w_subsid;
 
    int  nlay  = m_num_layers;
    Real dz    = m_geom.CellSize(2);
    bool moist = (m_cons->nComp() > RhoQ1_comp);
    auto ProbLoArr = m_geom.ProbLoArray();
 
    auto domain = m_geom.Domain();
    int ilo     = domain.smallEnd(0);
    int ihi     = domain.bigEnd(0);
    int jlo     = domain.smallEnd(1);
    int jhi     = domain.bigEnd(1);
 
    for (MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
        // NOTE: Grown box to get ghost cells in views
        const auto& gbx  = mfi.tilebox(IntVect(0,0,0),IntVect(1,1,0));
        const int nx     = gbx.length(0);
        const int imin   = gbx.smallEnd(0);
        const int jmin   = gbx.smallEnd(1);
        const int kmax   = gbx.bigEnd(2);
        const int offset = m_col_offsets[mfi.index()];
 
        const Array4<const Real>& cons_arr = m_cons->const_array(mfi);
 
        const Array4<const Real>& u_arr    = m_xvel->const_array(mfi);
        const Array4<const Real>& v_arr    = m_yvel->const_array(mfi);
        const Array4<const Real>& w_arr    = m_zvel->const_array(mfi);
 
        const Array4<const Real>& t13_arr   = m_tau13->const_array(mfi);
        const Array4<const Real>& t23_arr   = m_tau23->const_array(mfi);
        const Array4<const Real>& hfx3_arr  = m_hfx3->const_array(mfi);
        const Array4<const Real>& qfx3_arr  = m_qfx3->const_array(mfi);
 
        const Array4<const Real>& z_arr    = (m_z_phys) ? m_z_phys->const_array(mfi) :
                                                          Array4<const Real>{};
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // NOTE: k gets permuted with ilay
            // map [i,j,k] 0-based to [icol, ilay] 0-based
            const int icol   = (j-jmin)*nx + (i-imin) + offset;
            const int ilay   = kmax - k;
            const int ilayi  = kmax + 1 - k;
 
            // EOS input (at CC)
            Real r  = cons_arr(i,j,k,Rho_comp);
            Real rt = cons_arr(i,j,k,RhoTheta_comp);
            Real qv = (moist) ? cons_arr(i,j,k,RhoQ1_comp)/r : Real(0.);
            Real qc = (moist) ? cons_arr(i,j,k,RhoQ2_comp)/r : Real(0.);
 
            // EOS avg to z-face
            Real r_lo  = cons_arr(i,j,k-1,Rho_comp);
            Real rt_lo = cons_arr(i,j,k-1,RhoTheta_comp);
            Real qv_lo = (moist) ? cons_arr(i,j,k-1,RhoQ1_comp)/r_lo : Real(0.);
            Real rt_avg = Real(0.5) * (rt + rt_lo);
            Real qv_avg = Real(0.5) * (qv + qv_lo);
 
            // Delta z
            Real delz = (z_arr) ? Real(0.25) * ( (z_arr(i  ,j  ,k+1) - z_arr(i  ,j  ,k))
                                               + (z_arr(i+1,j  ,k+1) - z_arr(i+1,j  ,k))
                                               + (z_arr(i  ,j+1,k+1) - z_arr(i  ,j+1,k))
                                               + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz;
 
            // W at cc (cannot be 0?; inspection of shoc code...)
            Real w_cc = Real(0.5) * (w_arr(i,j,k) + w_arr(i,j,k+1));
            w_cc += (w_sub) ? Real(0.5) * (w_sub[k] + w_sub[k+1]) : Real(0.);
            Real w_limited = std::copysign(std::max(std::fabs(w_cc),Real(1.0e-6)),w_cc);
 
            // Input/Output data structures
            //=======================================================
            horiz_wind_d(icol,0,ilay)  = Real(0.5) * (u_arr(i,j,k) + u_arr(i+1,j  ,k));
            horiz_wind_d(icol,1,ilay)  = Real(0.5) * (v_arr(i,j,k) + v_arr(i  ,j+1,k));
            cldfrac_liq_d(icol,ilay)   = (qc>Real(0.)) ? Real(1.) : Real(0.);
            tke_d(icol,ilay)           = std::max(cons_arr(i,j,k,RhoKE_comp)/r, Real(0.));
            qc_d(icol,ilay)            = qc;
 
            // Interface data structures
            //=======================================================
            // eamxx_common_physics_functions_impl.hpp: calculate_vertical_velocity
            omega_d(icol,ilay)         = -w_limited * r * CONST_GRAV;
            if (k==0) {
                int ii  = std::min(std::max(i,ilo),ihi);
                int jj  = std::min(std::max(j,jlo),jhi);
 
                surf_mom_flux_d(icol,0)  = Real(0.5) * (t13_arr(ii,jj,k) + t13_arr(ii+1,jj  ,k));
                surf_mom_flux_d(icol,1)  = Real(0.5) * (t23_arr(ii,jj,k) + t23_arr(ii  ,jj+1,k));
                // No unit conversion to W/m^2 (ERF_ShocInterface.H L224)
                surf_sens_flux_d(icol)   = hfx3_arr(ii,jj,k);
                surf_evap_d(icol)        = (moist) ? qfx3_arr(ii,jj,k) : Real(0.);
                // EAMxx TMS expects rho * Cd * |U| [kg/(m^2 s)]. Back it out
                // from the conservative surface-stress and wind magnitudes.
                Real wsp = std::sqrt( horiz_wind_d(icol,0,ilay)[0]*horiz_wind_d(icol,0,ilay)[0]
                               + horiz_wind_d(icol,1,ilay)[0]*horiz_wind_d(icol,1,ilay)[0] );
                Real stress_mag = std::sqrt( surf_mom_flux_d(icol,0)*surf_mom_flux_d(icol,0)
                                           + surf_mom_flux_d(icol,1)*surf_mom_flux_d(icol,1) );
                surf_drag_coeff_tms_d(icol) = (wsp > Real(1.0e-8)) ? stress_mag / wsp : Real(0.);
            }
            T_mid_d(icol,ilay)          = getTgivenRandRTh(r, rt, qv);
            qv_d(icol,ilay)             = qv;
 
            // Input data structures
            //=======================================================
            p_mid_d(icol,ilay)       = getPgivenRTh(rt, qv);
            p_int_d(icol,ilayi)      = getPgivenRTh(rt_avg, qv_avg);
            // eamxx_common_physics_functions_impl.hpp: calculate_density
            pseudo_dens_d(icol,ilay) = r * CONST_GRAV * delz;
            // Enforce the grid spacing
            dz_d(icol,ilay) = delz;
            // Surface geopotential
            if (k==0) {
                Real z = (z_arr) ? Real(0.125) * ( (z_arr(i  ,j  ,k+1) + z_arr(i  ,j  ,k))
                                                 + (z_arr(i+1,j  ,k+1) + z_arr(i+1,j  ,k))
                                                 + (z_arr(i  ,j+1,k+1) + z_arr(i  ,j+1,k))
                                                 + (z_arr(i+1,j+1,k+1) + z_arr(i+1,j+1,k)) ) : ProbLoArr[2];
                phis_d(icol) = CONST_GRAV * z;
            }
 
            if (ilay==0) {
                Real r_hi  = cons_arr(i,j,k+1,Rho_comp);
                Real rt_hi = cons_arr(i,j,k+1,RhoTheta_comp);
                Real qv_hi = (moist) ? std::max(cons_arr(i,j,k+1,RhoQ1_comp)/r_hi,Real(0.)) : Real(0.);
                rt_avg = Real(0.5) * (rt + rt_hi);
                qv_avg = Real(0.5) * (qv + qv_hi);
                p_int_d(icol,0) = getPgivenRTh(rt_avg, qv_avg);
            }
        });
    }
}

Here is the call graph for this function:

name()

std::string SHOCInterface::name ( ) const

inline

{ return "shoc"; }

requested_buffer_size_in_bytes()

size_t SHOCInterface::requested_buffer_size_in_bytes ( ) const

protected

{
    using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
 
    const int nlev_packs       = ekat::npack<Spack>(m_num_layers);
    const int nlevi_packs      = ekat::npack<Spack>(m_num_layers+1);
    const int num_tracer_packs = ekat::npack<Spack>(m_num_tracers);
 
    // Number of Reals needed by local views in the interface
    const size_t interface_request = Buffer::num_1d_scalar_ncol*m_num_cols*sizeof(Real) +
                                     Buffer::num_1d_scalar_nlev*nlev_packs*sizeof(Spack) +
                                     Buffer::num_2d_vector_mid*m_num_cols*nlev_packs*sizeof(Spack) +
                                     Buffer::num_2d_vector_int*m_num_cols*nlevi_packs*sizeof(Spack) +
                                     Buffer::num_2d_vector_tr*m_num_cols*num_tracer_packs*sizeof(Spack);
 
    // Number of Reals needed by the WorkspaceManager passed to shoc_main
    const auto policy        = TPF::get_default_team_policy(m_num_cols, nlev_packs);
    const size_t wsm_request = WSM::get_total_bytes_needed(nlevi_packs, wsm_num_slots(), policy);
 
    return ( (interface_request + wsm_request)/sizeof(Real) );
}

run_impl()

void SHOCInterface::run_impl

(

const double

dt

)

{
    using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
 
    EKAT_REQUIRE_MSG (dt<=300,
                      "Error! SHOC is intended to run with a timestep no longer than 5 minutes.\n"
                      "       Please, reduce timestep (perhaps increasing subcycling iterations).\n");
 
    const auto nlev_packs     = ekat::npack<Spack>(m_num_layers);
    const auto scan_policy    = TPF::get_thread_range_parallel_scan_team_policy(m_num_cols, nlev_packs);
    const auto default_policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
 
    // Preprocessing of SHOC inputs. Kernel contains a parallel_scan,
    // so a special TeamPolicy is required.
    Kokkos::parallel_for("shoc_preprocess",
                         scan_policy,
                         shoc_preprocess);
    Kokkos::fence();
 
    auto wtracer_sfc = shoc_preprocess.wtracer_sfc;
    Kokkos::deep_copy(wtracer_sfc, 0);
 
    if (apply_tms) {
        apply_turbulent_mountain_stress();
    }
 
    if (check_flux_state) {
        check_flux_state_consistency(dt);
    }
 
    // For now set the host timestep to the shoc timestep. This forces
    // number of SHOC timesteps (nadv) to be one
    // TODO: input parameter?
    hdtime = dt;
    m_nadv = std::max(static_cast<int>(round(hdtime/dt)),1);
 
    // Reset internal WSM variables.
    workspace_mgr.reset_internals();
 
    // Run shoc main
    SHF::shoc_main(m_num_cols, m_num_layers, m_num_layers+1, m_npbl, m_nadv, m_num_tracers, dt,
                   workspace_mgr,runtime_options,input,input_output,output,history_output
#ifdef SCREAM_SHOC_SMALL_KERNELS
                   , temporaries
#endif
                   );
 
    // Postprocessing of SHOC outputs
    Kokkos::parallel_for("shoc_postprocess",
                         default_policy,
                         shoc_postprocess);
    Kokkos::fence();
 
    // Extra SHOC output diagnostics
    if (runtime_options.extra_diags) {
        Kokkos::deep_copy(shoc_mix,history_output.shoc_mix);
        Kokkos::deep_copy(brunt,history_output.brunt);
        Kokkos::deep_copy(w3,history_output.w3);
        Kokkos::deep_copy(isotropy,history_output.isotropy);
        Kokkos::deep_copy(wthl_sec,history_output.wthl_sec);
        Kokkos::deep_copy(wqw_sec,history_output.wqw_sec);
        Kokkos::deep_copy(uw_sec,history_output.uw_sec);
        Kokkos::deep_copy(vw_sec,history_output.vw_sec);
        Kokkos::deep_copy(qw_sec,history_output.qw_sec);
        Kokkos::deep_copy(thl_sec,history_output.thl_sec);
    }
}

set_computed_group_impl()

void SHOCInterface::set_computed_group_impl ( )

protected

set_diff_stresses()

void SHOCInterface::set_diff_stresses

(

)

{
    for (MFIter mfi(*m_hfx3); mfi.isValid(); ++mfi) {
        const auto& vbx_cc  = mfi.validbox();
        const auto& vbx_xz  = convert(vbx_cc,IntVect(1,0,1));
        const auto& vbx_yz  = convert(vbx_cc,IntVect(0,1,1));
 
        const Array4<Real>& hfx_arr = m_hfx3->array(mfi);
        const Array4<Real>& qfx_arr = m_qfx3->array(mfi);
 
        const Array4<Real>& t13_arr = m_tau13->array(mfi);
        const Array4<Real>& t23_arr = m_tau23->array(mfi);
 
        ParallelFor(vbx_cc, vbx_xz, vbx_yz,
                    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            hfx_arr(i,j,k) = Real(0.);
            qfx_arr(i,j,k) = Real(0.);
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            t13_arr(i,j,k) = Real(0.);
        },
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            t23_arr(i,j,k) = Real(0.);
        });
    }
}

Here is the call graph for this function:

set_eddy_diffs()

void SHOCInterface::set_eddy_diffs

(

)

{
    //
    // Expose for device capture
    //
 
    // Input/Output data structures
    //=======================================================
    auto tk_d = tk;
 
    // NOTE: Loop on grown box to fill ghost cells but limit
    //       to valid box where views are defined.
    for (MFIter mfi(*m_mu); mfi.isValid(); ++mfi) {
        const auto& gbx_cc = mfi.growntilebox();
        const auto& vbx_cc = mfi.validbox();
 
        // NOTE: Grown box only for mapping
        const auto& gbx  = mfi.tilebox(IntVect(0,0,0),IntVect(1,1,0));
        const int nx     = gbx.length(0);
        const int imin   = gbx.smallEnd(0);
        const int jmin   = gbx.smallEnd(1);
        const int kmax   = gbx.bigEnd(2);
        const int offset = m_col_offsets[mfi.index()];
 
        // Limiting to validbox
        const int iminv  = vbx_cc.smallEnd(0);
        const int imaxv  = vbx_cc.bigEnd(0);
        const int jminv  = vbx_cc.smallEnd(1);
        const int jmaxv  = vbx_cc.bigEnd(1);
        const int kminv  = vbx_cc.smallEnd(2);
        const int kmaxv  = vbx_cc.bigEnd(2);
 
        const Array4<Real>&       mu_arr   = m_mu->array(mfi);
        const Array4<const Real>& cons_arr = m_cons->const_array(mfi);
 
        ParallelFor(gbx_cc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // Limiting
            int ii = std::min(std::max(i,iminv),imaxv);
            int jj = std::min(std::max(j,jminv),jmaxv);
            int kk = std::min(std::max(k,kminv),kmaxv);
 
            // NOTE: k gets permuted with ilay
            // map [i,j,k] 0-based to [icol, ilay] 0-based
            const int icol   = (jj-jmin)*nx + (ii-imin) + offset;
            const int ilay   = kmax - kk;
 
            // NOTE: SHOC provides kinematic diffusivity, while eddyDiffs stores
            //       density-weighted diffusivity. Set mom_v for tau_33; all
            //       other vertical components are 0.
            mu_arr(i,j,k,EddyDiff::Mom_v)   = cons_arr(ii,jj,kk,Rho_comp) * tk_d(icol,ilay)[0];
            mu_arr(i,j,k,EddyDiff::Theta_v) = Real(0.);
            mu_arr(i,j,k,EddyDiff::KE_v)    = Real(0.);
            mu_arr(i,j,k,EddyDiff::Q_v)     = Real(0.);
        });
    }
 
    // Correct the internal ghost cells that have foextrap
    m_mu->FillBoundary(m_geom.periodicity());
}

Here is the call graph for this function:

set_grids()

void SHOCInterface::set_grids	(	int &	level,
		const amrex::BoxArray &	ba,
		amrex::Geometry &	geom,
		amrex::MultiFab *	cons,
		amrex::MultiFab *	xvel,
		amrex::MultiFab *	yvel,
		amrex::MultiFab *	zvel,
		amrex::Real *	w_subsid,
		amrex::MultiFab *	tau13,
		amrex::MultiFab *	tau23,
		amrex::MultiFab *	hfx3,
		amrex::MultiFab *	qfx3,
		amrex::MultiFab *	eddyDiffs,
		amrex::MultiFab *	z_phys
	)

{
    // Set data members that may change
    m_lev     = level;
    m_geom    = geom;
    m_cons    = cons;
    m_xvel    = xvel;
    m_yvel    = yvel;
    m_zvel    = zvel;
    m_w_subsid = w_subsid;
    m_tau13   = tau13;
    m_tau23   = tau23;
    m_hfx3    = hfx3;
    m_qfx3    = qfx3;
    m_mu      = eddyDiffs;
    m_z_phys  = z_phys;
 
    // Ensure the boxes span klo -> khi
    int klo = geom.Domain().smallEnd(2);
    int khi = geom.Domain().bigEnd(2);
 
    // Reset vector of offsets for columnar data
    m_num_layers = geom.Domain().length(2);
 
    int num_cols = 0;
    m_col_offsets.clear();
    m_col_offsets.resize(int(ba.size()));
    for (MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
        // NOTE: Get lateral ghost cells for CC <--> FC
        const auto& gbx = mfi.tilebox(IntVect(0,0,0),IntVect(1,1,0));
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE((klo == gbx.smallEnd(2)) &&
                                         (khi == gbx.bigEnd(2)),
                                         "Vertical decomposition with shoc is not allowed.");
        int nx = gbx.length(0);
        int ny = gbx.length(1);
        m_col_offsets[mfi.index()] = num_cols;
        num_cols += nx * ny;
    }
 
    // Resize the Kokkos variables that persist in memory
    if (num_cols != m_num_cols) {
        sgs_buoy_flux = view_2d();
        tk            = view_2d();
        sgs_buoy_flux = view_2d("sgs_buoy_flux", num_cols, m_num_layers);
        tk            = view_2d("eddy_diff_mom", num_cols, m_num_layers);
    }
    m_num_cols = num_cols;
 
    // Allocate the tendency MultiFabs
    c_tend.define(m_cons->boxArray(), m_cons->DistributionMap(), m_cons->nComp(), 0);
    u_tend.define(m_xvel->boxArray(), m_xvel->DistributionMap(), m_xvel->nComp(), 0);
    v_tend.define(m_yvel->boxArray(), m_yvel->DistributionMap(), m_yvel->nComp(), 0);
 
    // Allocate the buffer arrays in ERF
    alloc_buffers();
 
    // Allocate the m_buffer struct
    init_buffers();
 
    // Fill the KOKKOS Views from AMReX MFs
    mf_to_kokkos_buffers();
}

Here is the call graph for this function:

wsm_num_slots()

Int SHOCInterface::wsm_num_slots ( ) const

inlineprotected

    {
        const Int n_wind_slots = ekat::npack<Spack>(m_num_vel_comp)*Spack::n;
        const Int n_trac_slots = ekat::npack<Spack>(m_num_tracers +
                                                    shoc_extra_trac_slots)*Spack::n;
        return ( 14 + (n_wind_slots + n_trac_slots) );
    }

Member Data Documentation

apply_tms

bool SHOCInterface::apply_tms = false

protected

brunt

view_2d SHOCInterface::brunt

protected

c_tend

amrex::MultiFab SHOCInterface::c_tend

protected

check_flux_state

bool SHOCInterface::check_flux_state = false

protected

cldfrac_liq

view_2d SHOCInterface::cldfrac_liq

protected

cldfrac_liq_prev

view_2d SHOCInterface::cldfrac_liq_prev

protected

column_conservation_check

bool SHOCInterface::column_conservation_check = false

protected

extra_shoc_diags

bool SHOCInterface::extra_shoc_diags = false

protected

hdtime

Int SHOCInterface::hdtime

protected

heat_flux

view_1d SHOCInterface::heat_flux

protected

history_output

SHF::SHOCHistoryOutput SHOCInterface::history_output

protected

horiz_wind

view_3d SHOCInterface::horiz_wind

protected

ice_flux

view_1d SHOCInterface::ice_flux

protected

input

SHF::SHOCInput SHOCInterface::input

protected

input_output

SHF::SHOCInputOutput SHOCInterface::input_output

protected

inv_qc_relvar

view_2d SHOCInterface::inv_qc_relvar

protected

isotropy

view_2d SHOCInterface::isotropy

protected

m_ba

amrex::BoxArray SHOCInterface::m_ba

protected

m_buffer

Buffer SHOCInterface::m_buffer

protected

m_col_offsets

amrex::Vector<int> SHOCInterface::m_col_offsets

protected

m_cons

amrex::MultiFab* SHOCInterface::m_cons = nullptr

protected

m_first_step

bool SHOCInterface::m_first_step = true

protected

m_geom

amrex::Geometry SHOCInterface::m_geom

protected

m_hfx3

amrex::MultiFab* SHOCInterface::m_hfx3 = nullptr

protected

m_lev

int SHOCInterface::m_lev

protected

m_mu

amrex::MultiFab* SHOCInterface::m_mu = nullptr

protected

m_nadv

Int SHOCInterface::m_nadv

protected

m_npbl

Int SHOCInterface::m_npbl

protected

m_num_cols

Int SHOCInterface::m_num_cols = 0

protected

m_num_layers

Int SHOCInterface::m_num_layers = 0

protected

m_num_tracers

Int SHOCInterface::m_num_tracers = 3

protected

Referenced by wsm_num_slots().

m_num_vel_comp

Int SHOCInterface::m_num_vel_comp = 2

protected

Referenced by wsm_num_slots().

m_qfx3

amrex::MultiFab* SHOCInterface::m_qfx3 = nullptr

protected

m_step

int SHOCInterface::m_step

protected

m_tau13

amrex::MultiFab* SHOCInterface::m_tau13 = nullptr

protected

m_tau23

amrex::MultiFab* SHOCInterface::m_tau23 = nullptr

protected

m_w_subsid

amrex::Real* SHOCInterface::m_w_subsid = nullptr

protected

m_xvel

amrex::MultiFab* SHOCInterface::m_xvel = nullptr

protected

m_yvel

amrex::MultiFab* SHOCInterface::m_yvel = nullptr

protected

m_z_phys

amrex::MultiFab* SHOCInterface::m_z_phys = nullptr

protected

m_zvel

amrex::MultiFab* SHOCInterface::m_zvel = nullptr

protected

obklen

view_1d SHOCInterface::obklen

protected

omega

view_2d SHOCInterface::omega

protected

output

SHF::SHOCOutput SHOCInterface::output

protected

p_int

view_2d SHOCInterface::p_int

protected

p_mid

view_2d SHOCInterface::p_mid

protected

pblh

view_1d SHOCInterface::pblh

protected

phis

view_1d SHOCInterface::phis

protected

pseudo_dens

view_2d SHOCInterface::pseudo_dens

protected

qc

view_2d SHOCInterface::qc

protected

qtracers

view_3d_strided SHOCInterface::qtracers

protected

qv

view_2d SHOCInterface::qv

protected

qw_sec

view_2d SHOCInterface::qw_sec

protected

runtime_options

SHF::SHOCRuntime SHOCInterface::runtime_options

protected

sgs_buoy_flux

view_2d SHOCInterface::sgs_buoy_flux

protected

shoc_cond

view_2d SHOCInterface::shoc_cond

protected

shoc_evap

view_2d SHOCInterface::shoc_evap

protected

shoc_extra_trac_slots

constexpr Int SHOCInterface::shoc_extra_trac_slots = 3

staticconstexprprotected

Referenced by wsm_num_slots().

shoc_mix

view_2d SHOCInterface::shoc_mix

protected

shoc_postprocess

SHOCPostprocess SHOCInterface::shoc_postprocess

protected

shoc_preprocess

SHOCPreprocess SHOCInterface::shoc_preprocess

protected

surf_drag_coeff_tms

view_1d SHOCInterface::surf_drag_coeff_tms

protected

surf_evap

view_1d SHOCInterface::surf_evap

protected

surf_mom_flux

sview_2d SHOCInterface::surf_mom_flux

protected

surf_sens_flux

view_1d SHOCInterface::surf_sens_flux

protected

T_mid

view_2d SHOCInterface::T_mid

protected

thl_sec

view_2d SHOCInterface::thl_sec

protected

tk

view_2d SHOCInterface::tk

protected

tke

view_2d SHOCInterface::tke

protected

tkh

view_2d SHOCInterface::tkh

protected

tot_buff_view

view_1d SHOCInterface::tot_buff_view

protected

u_tend

amrex::MultiFab SHOCInterface::u_tend

protected

ustar

view_1d SHOCInterface::ustar

protected

uw_sec

view_2d SHOCInterface::uw_sec

protected

v_tend

amrex::MultiFab SHOCInterface::v_tend

protected

vapor_flux

view_1d SHOCInterface::vapor_flux

protected

vw_sec

view_2d SHOCInterface::vw_sec

protected

w3

view_2d SHOCInterface::w3

protected

w_sec

view_2d SHOCInterface::w_sec

protected

water_flux

view_1d SHOCInterface::water_flux

protected

workspace_mgr

ekat::WorkspaceManager<Spack, KT::Device> SHOCInterface::workspace_mgr

protected

wqw_sec

view_2d SHOCInterface::wqw_sec

protected

wthl_sec

view_2d SHOCInterface::wthl_sec

protected

---

The documentation for this class was generated from the following files:

• Source/PhysicsInterfaces/Shoc/ERF_ShocInterface.H
• Source/PhysicsInterfaces/Shoc/ERF_ShocInterface.cpp