ShocMoments Class Reference

#include <ERF_ShocMoments.H>

Static Public Member Functions
static void	calc_var_or_covar (const ShocColumnData &col, amrex::Real tunefac, const amrex::FArrayBox &isotropy_zi, const amrex::FArrayBox &tkh_zi, const amrex::FArrayBox &invar1, const amrex::FArrayBox &invar2, amrex::FArrayBox &outvar)
static void	calc_vertflux (const ShocColumnData &col, const amrex::FArrayBox &tkh_zi, const amrex::FArrayBox &invar, amrex::FArrayBox &vertflux)
static void	diagnose_second_moments (ShocColumnData &col, const ShocRuntimeOptions &opts)
static void	clip_third_moments (const ShocColumnData &col, const amrex::FArrayBox &w_sec_zi, amrex::FArrayBox &w3)
static void	diagnose_third_moments (ShocColumnData &col, const ShocRuntimeOptions &opts)
static void	diagnose_moments (ShocColumnData &col, const ShocRuntimeOptions &opts)

Member Function Documentation

calc_var_or_covar()

void ShocMoments::calc_var_or_covar ( const ShocColumnData &  col, amrex::Real  tunefac, const amrex::FArrayBox &  isotropy_zi, const amrex::FArrayBox &  tkh_zi, const amrex::FArrayBox &  invar1, const amrex::FArrayBox &  invar2, amrex::FArrayBox &  outvar  )

static

{
    auto out = outvar.array();
    const auto iso = isotropy_zi.const_array();
    const auto tkh = tkh_zi.const_array();
    const auto v1 = invar1.const_array();
    const auto v2 = invar2.const_array();
    const auto zt = col.zt.const_array();
    const auto zi = col.zi.const_array();
    const auto layout = col.layout;
    const Box col_box(IntVect(0,0,0), IntVect(layout.ncell - 1, 0, 0));
    ParallelFor(col_box, [=] AMREX_GPU_DEVICE (int ic, int, int) noexcept
    {
        for (int k = 1; k < layout.nlev; ++k) {
            const Real dz_zi = cell_spacing_at_iface(zt, zi, layout, ic, k);
            const Real grid_dz2 = 1.0_rt / (dz_zi * dz_zi);
            out(ic,k,0) = tunefac * (iso(ic,k,0) * tkh(ic,k,0)) * grid_dz2 *
                          (v1(ic,k-1,0) - v1(ic,k,0)) * (v2(ic,k-1,0) - v2(ic,k,0));
        }
    });
}

Here is the call graph for this function:

calc_vertflux()

void ShocMoments::calc_vertflux ( const ShocColumnData &  col, const amrex::FArrayBox &  tkh_zi, const amrex::FArrayBox &  invar, amrex::FArrayBox &  vertflux  )

static

{
    auto out = vertflux.array();
    const auto tkh = tkh_zi.const_array();
    const auto v = invar.const_array();
    const auto zt = col.zt.const_array();
    const auto zi = col.zi.const_array();
    const auto layout = col.layout;
    const Box col_box(IntVect(0,0,0), IntVect(layout.ncell - 1, 0, 0));
    ParallelFor(col_box, [=] AMREX_GPU_DEVICE (int ic, int, int) noexcept
    {
        for (int k = 1; k < layout.nlev; ++k) {
            const Real dz_zi = cell_spacing_at_iface(zt, zi, layout, ic, k);
            // E3SM uses top-down vertical indexing. ERF columns are bottom-up,
            // so the interface gradient must use upper-minus-lower to preserve
            // the same physical downgradient flux sign.
            out(ic,k,0) = -(tkh(ic,k,0) / dz_zi) * (v(ic,k,0) - v(ic,k-1,0));
        }
    });
}

Here is the call graph for this function:

clip_third_moments()

void ShocMoments::clip_third_moments ( const ShocColumnData &  col, const amrex::FArrayBox &  w_sec_zi, amrex::FArrayBox &  w3  )

static

{
    auto w3 = w3_fab.array();
    const auto wsec = w_sec_zi.const_array();
 
    const auto layout = col.layout;
    const Box col_box(IntVect(0,0,0), IntVect(layout.ncell - 1, 0, 0));
    ParallelFor(col_box, [=] AMREX_GPU_DEVICE (int ic, int, int) noexcept
    {
        for (int k = 0; k <= layout.nlev; ++k) {
            const Real wsec3 = wsec(ic,k,0) * wsec(ic,k,0) * wsec(ic,k,0);
            const Real clip_cond = shoc_w3clip() * std::sqrt(amrex::max(0.0_rt, 2.0_rt * wsec3));
            // E3SM/EAMxx replaces an out-of-range w3 with a fixed positive
            // constant (0.02). That discards the sign of the third moment,
            // so a strongly negatively skewed column is handed back a
            // positively skewed w3, and it can also move |w3| further away
            // from the bound it is supposed to enforce when clip_cond is
            // small. Limit the magnitude to clip_cond instead and keep the
            // sign, which is continuous in w3 and respects the bound.
            if (amrex::Math::abs(w3(ic,k,0)) > clip_cond) {
                w3(ic,k,0) = std::copysign(clip_cond, w3(ic,k,0));
            }
        }
    });
}

Here is the call graph for this function:

diagnose_moments()

void ShocMoments::diagnose_moments ( ShocColumnData &  col, const ShocRuntimeOptions &  opts  )

static

{
    diagnose_second_moments(col, opts);
    diagnose_third_moments(col, opts);
}

Referenced by ShocDiagnostics::diagnose_post_implicit().

Here is the caller graph for this function:

diagnose_second_moments()

void ShocMoments::diagnose_second_moments ( ShocColumnData &  col, const ShocRuntimeOptions &  opts  )

static

{
    FArrayBox isotropy_zi, tkh_zi, tk_zi;
    const Box iface_box(IntVect(0,0,0), IntVect(col.layout.ncell - 1, col.layout.nlev, 0));
    isotropy_zi.resize(iface_box, 1, The_Async_Arena());
    tkh_zi.resize(iface_box, 1, The_Async_Arena());
    tk_zi.resize(iface_box, 1, The_Async_Arena());
 
    interpolate_cc_to_iface(col, col.isotropy, isotropy_zi, 0.0);
    interpolate_cc_to_iface(col, col.tkh, tkh_zi, 0.0);
    interpolate_cc_to_iface(col, col.tk, tk_zi, 0.0);
 
    auto thl_sec = col.thl_sec.array();
    auto qw_sec = col.qw_sec.array();
    auto qwthl_sec = col.qwthl_sec.array();
    auto w_sec = col.w_sec.array();
 
    const auto tke = col.tke.const_array();
    const auto layout = col.layout;
    const Box col_box(IntVect(0,0,0), IntVect(layout.ncell - 1, 0, 0));
    ParallelFor(col_box, [=] AMREX_GPU_DEVICE (int ic, int, int) noexcept
    {
        for (int k = 0; k < layout.nlev; ++k) {
            w_sec(ic,k,0) = opts.shoc_1p5tke ? 0.0_rt : opts.w2tune * (2.0_rt / 3.0_rt) * tke(ic,k,0);
        }
    });
 
    if (opts.shoc_1p5tke) {
        ParallelFor(col_box, [=] AMREX_GPU_DEVICE (int ic, int, int) noexcept
        {
            for (int k = 0; k <= layout.nlev; ++k) {
                thl_sec(ic,k,0) = 0.0_rt;
                qw_sec(ic,k,0) = 0.0_rt;
                qwthl_sec(ic,k,0) = 0.0_rt;
            }
        });
    } else {
        calc_var_or_covar(col, opts.thl2tune, isotropy_zi, tkh_zi, col.thetal, col.thetal, col.thl_sec);
        calc_var_or_covar(col, opts.qw2tune, isotropy_zi, tkh_zi, col.qw, col.qw, col.qw_sec);
        calc_var_or_covar(col, opts.qwthl2tune, isotropy_zi, tkh_zi, col.thetal, col.qw, col.qwthl_sec);
    }
 
    calc_vertflux(col, tkh_zi, col.thetal, col.wthl_sec);
    calc_vertflux(col, tkh_zi, col.qw, col.wqw_sec);
    calc_vertflux(col, tkh_zi, col.tke, col.wtke_sec);
    calc_vertflux(col, tk_zi, col.u, col.uw_sec);
    calc_vertflux(col, tk_zi, col.v, col.vw_sec);
 
    apply_second_moment_boundary_conditions(col);
    apply_top_taper_to_second_moments(col, opts);
}

Here is the call graph for this function:

diagnose_third_moments()

void ShocMoments::diagnose_third_moments ( ShocColumnData &  col, const ShocRuntimeOptions &  opts  )

static

{
    FArrayBox isotropy_zi, brunt_zi, w_sec_zi, thetal_zi;
    const Box iface_box(IntVect(0,0,0), IntVect(col.layout.ncell - 1, col.layout.nlev, 0));
    isotropy_zi.resize(iface_box, 1, The_Async_Arena());
    brunt_zi.resize(iface_box, 1, The_Async_Arena());
    w_sec_zi.resize(iface_box, 1, The_Async_Arena());
    thetal_zi.resize(iface_box, 1, The_Async_Arena());
 
    interpolate_cc_to_iface(col, col.isotropy, isotropy_zi, 0.0);
    interpolate_cc_to_iface(col, col.brunt, brunt_zi, shoc_large_neg());
    interpolate_cc_to_iface(col, col.w_sec, w_sec_zi, (2.0_rt / 3.0_rt) * shoc_min_tke());
    interpolate_cc_to_iface(col, col.thetal, thetal_zi, 0.0);
 
    auto w3 = col.w3.array();
    const auto w_sec = col.w_sec.const_array();
    const auto thl_sec = col.thl_sec.const_array();
    const auto wthl_sec = col.wthl_sec.const_array();
    const auto tke = col.tke.const_array();
    const auto dz = col.dz.const_array();
    const auto isotropy_i = isotropy_zi.const_array();
    const auto brunt_i = brunt_zi.const_array();
    const auto w_sec_i = w_sec_zi.const_array();
    const auto thetal_i = thetal_zi.const_array();
    const auto zt = col.zt.const_array();
    const auto zi = col.zi.const_array();
    const auto layout = col.layout;
    const Box col_box(IntVect(0,0,0), IntVect(layout.ncell - 1, 0, 0));
    ParallelFor(col_box, [=] AMREX_GPU_DEVICE (int ic, int, int) noexcept
    {
        w3(ic,0,0) = 0.0_rt;
        w3(ic,layout.nlev,0) = 0.0_rt;
 
        for (int k = 1; k < layout.nlev; ++k) {
            if (opts.shoc_1p5tke) {
                w3(ic,k,0) = 0.0_rt;
                continue;
            }
 
            const Real dz_zt_k = amrex::max(dz(ic,k,0), eps());
            const Real dz_zt_km1 = amrex::max(dz(ic,k-1,0), eps());
            const Real dz_zi = cell_spacing_at_iface(zt, zi, layout, ic, k);
            const Real thedz = 1.0_rt / dz_zi;
            const Real thedz2 = 1.0_rt / (dz_zt_k + dz_zt_km1);
            const Real iso = isotropy_i(ic,k,0);
            const Real isosq = iso * iso;
            const Real buoy_sgs2 = isosq * brunt_i(ic,k,0);
            const Real bet2 = CONST_GRAV / amrex::max(thetal_i(ic,k,0), eps());
 
            // E3SM's top-down indexing forms above-minus-below centered
            // differences here. In ERF's bottom-up ordering, the upper
            // neighbor is k+1 and the lower neighbor is k-1, while the cell
            // just above interface k is k and the one below is k-1.
            const Real thl_sec_diff = thl_sec(ic,amrex::min(k+1, layout.nlev),0) - thl_sec(ic,k-1,0);
            const Real wthl_sec_diff = wthl_sec(ic,amrex::min(k+1, layout.nlev),0) - wthl_sec(ic,k-1,0);
            const Real wsec_diff = w_sec(ic,k,0) - w_sec(ic,k-1,0);
            const Real tke_diff = tke(ic,k,0) - tke(ic,k-1,0);
 
            const Real c = opts.c_diag_3rd_mom;
            const Real a0 = (0.52_rt / (c * c)) / amrex::max(c - 2.0_rt, eps());
            const Real a1 = 0.87_rt / (c * c);
            const Real a2 = 0.5_rt / c;
            const Real a3 = 0.6_rt / (c * amrex::max(c - 2.0_rt, eps()));
            const Real a4 = 2.4_rt / (3.0_rt * c + 5.0_rt);
            const Real a5 = 0.6_rt / (c * (3.0_rt + 5.0_rt * c));
            const Real bet2_sq = bet2 * bet2;
            const Real bet2_cu = bet2_sq * bet2;
            const Real iso_cu = isosq * iso;
 
            const Real f0 = thedz2 * bet2_cu * isosq * isosq *
                            wthl_sec(ic,k,0) * thl_sec_diff;
            const Real f1 = thedz2 * bet2_sq * iso_cu *
                            (wthl_sec(ic,k,0) * wthl_sec_diff +
                            0.5_rt * w_sec_i(ic,k,0) * thl_sec_diff);
            const Real f2 = thedz * bet2 * isosq * wthl_sec(ic,k,0) * wsec_diff +
                            2.0_rt * thedz2 * bet2 * isosq * w_sec_i(ic,k,0) * wthl_sec_diff;
            const Real f3 = thedz2 * bet2 * isosq * w_sec_i(ic,k,0) * wthl_sec_diff +
                            thedz * bet2 * isosq * (wthl_sec(ic,k,0) * tke_diff);
            const Real f4 = thedz * iso * w_sec_i(ic,k,0) * (wsec_diff + tke_diff);
            const Real f5 = thedz * iso * w_sec_i(ic,k,0) * wsec_diff;
 
            const Real denom0 = signed_denominator(1.0_rt - (a1 + a3) * buoy_sgs2);
            const Real omega0 = a4 / signed_denominator(1.0_rt - a5 * buoy_sgs2);
            const Real omega1 = omega0 / (2.0_rt * c);
            const Real omega2 = omega1 * f3 + 1.25_rt * omega0 * f4;
            const Real x0 = (a2 * buoy_sgs2 * (1.0_rt - a3 * buoy_sgs2)) / denom0;
            const Real y0 = (2.0_rt * a2 * buoy_sgs2 * x0) / signed_denominator(1.0_rt - a3 * buoy_sgs2);
            const Real x1 = (a0 * f0 + a1 * f1 + a2 * (1.0_rt - a3 * buoy_sgs2) * f2) / denom0;
            const Real y1 = (2.0_rt * a2 * (buoy_sgs2 * x1 + (a0 / amrex::max(a1, eps())) * f0 + f1)) /
                            signed_denominator(1.0_rt - a3 * buoy_sgs2);
            const Real aa0 = omega0 * x0 + omega1 * y0;
            const Real aa1 = omega0 * x1 + omega1 * y1 + omega2;
 
            const Real denom = signed_denominator(c - 1.2_rt * x0 + aa0);
            Real w3_val = (aa1 - 1.2_rt * x1 - 1.5_rt * f5) / denom;
            w3(ic,k,0) = w3_val;
        }
    });
 
    clip_third_moments(col, w_sec_zi, col.w3);
    apply_top_taper_to_third_moments(col, opts);
}

Here is the call graph for this function:

---

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

• Source/PBL/Shoc/ERF_ShocMoments.H
• Source/PBL/Shoc/ERF_ShocMoments.cpp