ERF_ComputeDiffusivityMYNNEDMF.cpp File Reference

#include <algorithm>
#include <iostream>
#include <vector>
#include <cmath>
#include <functional>
#include <limits>
#include "ERF_SurfaceLayer.H"
#include "ERF_DirectionSelector.H"
#include "ERF_Diffusion.H"
#include "ERF_Constants.H"
#include "ERF_TurbStruct.H"
#include "ERF_PBLModels.H"

Include dependency graph for ERF_ComputeDiffusivityMYNNEDMF.cpp:

Functions
void	ComputeDiffusivityMYNNEDMF (const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const std::unique_ptr< MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
	Compute eddy viscosity and diffusivity coefficients using the MYNN-EDMF closure. More...

Function Documentation

ComputeDiffusivityMYNNEDMF()

void ComputeDiffusivityMYNNEDMF	(	const MultiFab &	xvel,
		const MultiFab &	yvel,
		const MultiFab &	cons_in,
		MultiFab &	eddyViscosity,
		const Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const std::unique_ptr< MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices
	)

Compute eddy viscosity and diffusivity coefficients using the MYNN-EDMF closure.

Parameters

[in]	xvel	Horizontal x-velocity field.
[in]	yvel	Horizontal y-velocity field.
[in]	cons_in	Field of conservative variables.
[out]	eddyViscosity	MultiFab to be filled with computed eddy diffusivities.
[in]	geom	Grid geometry.
[in]	turbChoice	Turbulence closure options and parameters.
[in]	SurfLayer	Surface layer data for MOST parameters.
[in]	use_terrain_fitted_coords	Flag to use terrain-fitted vertical coordinates.
[in]	use_moisture	Flag to include moisture in buoyancy calculations.
[in]	level	Current level index.
[in]	bc_ptr	Pointer to boundary condition records.
[in]	z_phys_nd	Physical height field at nodes.
[in]	z_phys_cc	Physical height field at cell centers.
[in]	moisture_indices	Indices for mapping moisture variables in the conservative field.

{
    auto mynn     = turbChoice.pbl_mynn;
    const StratType pbl_strat = turbChoice.pbl_strat_type;
    auto level2   = turbChoice.pbl_mynn_level2;
 
    Real Lt_alpha = mynn.Lt_alpha;
 
    // Dirichlet flags to switch derivative stencil
    bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
    bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir) );
    bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
    bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) );
    bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
    bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) );
 
    // Epsilon
    Real eps = std::numeric_limits<Real>::epsilon();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    // NOTE: we must not tile in z here because the body of this loop assumes that each
    //       iterate spans the entire column: it accumulates vertical integrals into a
    //       per-iterate qintegral fab (as in MYNN25)
    for ( MFIter mfi(eddyViscosity,TileNoZ()); mfi.isValid(); ++mfi) {
 
        // NOTE: the valid box, not a grown box, as in MYNN25.  Growing by one in z put
        //       k = -1 and k = nz in the loop, and the vertical-derivative stencil below
        //       reaches k-1 and k+1, so the scheme read cons_in, xvel and yvel two cells
        //       outside the domain -- ghost cells this routine has no guarantee about.
        //       The ghost values it computed were thrown away regardless:
        //       ComputeTurbulentViscosity refills every eddy-viscosity ghost cell after
        //       this routine returns, by FillBoundary and by extrapolation onto the
        //       physical-boundary planes.
        const Box& bx = mfi.tilebox();
        const Array4<Real const>& cell_data = cons_in.array(mfi);
        const Array4<Real      >& K_turb    = eddyViscosity.array(mfi);
        const Array4<Real const>& uvel      = xvel.array(mfi);
        const Array4<Real const>& vvel      = yvel.array(mfi);
 
        // Compute some quantities that are constant in each column: each iterate must
        // hold a whole column for the vertical integrals below to be complete
        const Box& dbx = geom.Domain();
        AMREX_ALWAYS_ASSERT(bx.smallEnd(2) == dbx.smallEnd(2) && bx.bigEnd(2) == dbx.bigEnd(2));
 
        const GeometryData gdata = geom.data();
 
        const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
 
        FArrayBox qintegral(xybx,2,The_Async_Arena());
        FArrayBox qturb(bx,1,The_Async_Arena());
 
        qintegral.setVal<RunOn::Device>(0);
 
        const Array4<Real> qint = qintegral.array();
        const Array4<Real> qvel = qturb.array();
 
        // vertical integrals to compute lengthscale
        if (use_terrain_fitted_coords) {
            const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
            const auto invCellSize = geom.InvCellSizeArray();
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                qvel(i,j,k) = std::sqrt(two * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
                AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > zero, "KE must have a positive value");
 
                const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
                const Real dz   = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr);
                Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz);
                Gpu::Atomic::Add(&qint(i,j,0,1),      qvel(i,j,k)*dz);
            });
        } else {
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                qvel(i,j,k) = std::sqrt(two * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
                AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > zero, "KE must have a positive value");
 
                // Not multiplying by dz: its constant and would fall out when we divide qint0/qint1 anyway
 
                const Real Zval = (k + myhalf)*gdata.CellSize(2);
                Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k));
                Gpu::Atomic::Add(&qint(i,j,0,1),      qvel(i,j,k));
            });
        }
 
        int izmin = geom.Domain().smallEnd(2);
        int izmax = geom.Domain().bigEnd(2);
 
        // Spatially varying MOST
        Real d_kappa   = KAPPA;
        Real d_gravity = CONST_GRAV;
 
        const auto& t_mean_mf = SurfLayer->get_mac_avg(level,5); // theta_v
        const auto& q_mean_mf = SurfLayer->get_mac_avg(level,4); // q_v
        const auto& u_star_mf = SurfLayer->get_u_star(level);
        const auto& t_star_mf = SurfLayer->get_t_star(level);
        const auto& q_star_mf = SurfLayer->get_q_star(level);
 
        const auto& tm_arr     = t_mean_mf->const_array(mfi);
        const auto& qm_arr     = q_mean_mf->const_array(mfi);
        const auto& u_star_arr = u_star_mf->const_array(mfi);
        const auto& t_star_arr = t_star_mf->const_array(mfi);
        const auto& q_star_arr = (use_moisture) ? q_star_mf->const_array(mfi) : Array4<Real>{};
 
        const Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
        const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // Compute some partial derivatives that we will need (second order)
            // U and V derivatives are interpolated to account for staggered grid
            Real dthetadz, dudz, dvdz;
            ComputeVerticalDerivativesPBL(i, j, k,
                                          uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
                                          c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                          u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                          v_ext_dir_on_zlo, v_ext_dir_on_zhi,
                                          dthetadz, dudz, dvdz,
                                          moisture_indices, pbl_strat);
 
            // Spatially varying MOST
            Real theta0 = tm_arr(i,j,0);
            Real qv0    = qm_arr(i,j,0);
            Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
            Real surface_latent_heat{0};
            if (use_moisture) {
                // Compute buoyancy flux (Stull Eqn. 4.4.5d)
                surface_latent_heat = -u_star_arr(i,j,0) * q_star_arr(i,j,0);
                surface_heat_flux *= (one + Real(0.61)*qv0);
                surface_heat_flux += Real(0.61) * theta0 * surface_latent_heat;
            }
 
            Real l_obukhov;
            if (std::abs(surface_heat_flux) > eps && u_star_arr(i,j,0) > eps) {
                l_obukhov = -( theta0 * u_star_arr(i,j,0)*u_star_arr(i,j,0)*u_star_arr(i,j,0) )
                           / ( d_kappa * d_gravity * surface_heat_flux );
            } else {
                l_obukhov = std::numeric_limits<Real>::max();
            }
 
            // Surface-layer length scale (NN09, Eqn. 53)
            AMREX_ASSERT(l_obukhov != 0);
            const Real zval = use_terrain_fitted_coords
                            ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
                            : (k + myhalf)*gdata.CellSize(2);
            const Real zeta = zval/l_obukhov;
            Real l_S;
            if (zeta >= one) {
                l_S = KAPPA*zval/Real(3.7);
            } else if (zeta >= 0) {
                l_S = KAPPA*zval/(1+Real(2.7)*zeta);
            } else {
                l_S = KAPPA*zval*std::pow(one - Real(100.0) * zeta, Real(0.2));
            }
 
            // Replace the resolved gradients in the first cell with the MOST
            // profile gradients; see ApplySurfaceLayerGradientsPBL (ERF #4037)
            if (k == izmin) {
                PBLSurfaceLayerGradient sl;
                sl.u_star  = u_star_arr(i,j,0);
                sl.tstar_v = ComputeVirtualTStarPBL(t_star_arr(i,j,0),
                                                    (use_moisture) ? q_star_arr(i,j,0) : zero,
                                                    theta0, qv0, use_moisture);
                sl.zval    = zval;
                sl.zeta    = zeta;
                ApplySurfaceLayerGradientsPBL(sl, dthetadz, dudz, dvdz);
            }
 
            // ABL-depth length scale (NN09, Eqn. 54)
            Real l_T;
            if (qint(i,j,0,1) > zero) {
                l_T = Lt_alpha*qint(i,j,0,0)/qint(i,j,0,1);
            } else {
                l_T = std::numeric_limits<Real>::max();
            }
 
            // Buoyancy length scale (NN09, Eqn. 55)
            Real l_B;
            if (dthetadz > 0) {
                Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
                if (zeta < 0) {
                    Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T; // velocity scale
                    qc = std::pow(qc,one/three);
                    l_B = (one + Real(5.0)*std::sqrt(qc/(N_brunt_vaisala * l_T))) * qvel(i,j,k)/N_brunt_vaisala;
                } else {
                    l_B = qvel(i,j,k) / N_brunt_vaisala;
                }
            } else {
                l_B = std::numeric_limits<Real>::max();
            }
 
            // Master length scale
            Real Lm;
            if (mynn.config == MYNNConfigType::Chen2021) {
                Lm = std::pow(one/(l_S*l_S) + one/(l_T*l_T) + one/(l_B*l_B), -myhalf);
            } else {
                // NN09, Eqn 52
                Lm = one / (one/l_S + one/l_T + one/l_B);
            }
 
            // Calculate nondimensional production terms
            Real shearProd  = dudz*dudz + dvdz*dvdz;
            Real buoyProd   = -(CONST_GRAV/theta0) * dthetadz;
            Real L2_over_q2 = Lm*Lm/(qvel(i,j,k)*qvel(i,j,k));
            Real GM         = L2_over_q2 * shearProd;
            Real GH         = L2_over_q2 * buoyProd;
 
            // Equilibrium (Level-2) q calculation follows NN09, Appendix 2
            Real Rf  = level2.calc_Rf(GM, GH);
            Real SM2 = level2.calc_SM(Rf);
            Real qe2 = mynn.B1*Lm*Lm*SM2*(one-Rf)*shearProd;
            // Clamp before the sqrt rather than selecting after it: the optimiser
            // evaluates a guarded sqrt of a negative qe2 before the selection, which
            // trips amrex.fpe_trap_invalid. The value is unchanged for every qe2.
            Real qe  = std::sqrt(amrex::max(qe2, amrex::Real(zero)));
 
            // Level 2 limiting (Helfand and Labraga 1988)
            Real alphac  = (qvel(i,j,k) > qe) ? one : qvel(i,j,k) / (qe + eps);
 
            // Level Real(2.5) stability functions
            Real SM, SH, SQ;
            mynn.calc_stability_funcs(SM,SH,SQ,GM,GH,alphac);
 
            // Clip SM, SH, SQ following WRF.  SQ is proportional to the *unclipped*
            // SM (NN09 Eqn. 67 is evaluated before SM is limited), so it needs its own
            // bounds; without them the TKE diffusivity below can go negative and turn
            // the vertical TKE diffusion anti-diffusive.
            SM = amrex::min(amrex::max(SM,mynn.SMmin), mynn.SMmax);
            SH = amrex::min(amrex::max(SH,mynn.SHmin), mynn.SHmax);
            SQ = amrex::min(amrex::max(SQ,mynn.SQmin), mynn.SQmax);
 
            // Finally, compute the eddy viscosity/diffusivities
            const Real rho = cell_data(i,j,k,Rho_comp);
            K_turb(i,j,k,EddyDiff::Mom_v)   = rho * Lm * qvel(i,j,k) * SM;
            K_turb(i,j,k,EddyDiff::Theta_v) = rho * Lm * qvel(i,j,k) * SH;
            K_turb(i,j,k,EddyDiff::KE_v)    = rho * Lm * qvel(i,j,k) * SQ;
 
            // TODO: implement partial-condensation scheme?
            // Currently, implementation matches NN09 without rain (i.e.,
            // the liquid water potential temperature is equal to the
            // potential temperature.
 
            // NN09 gives the total water content flux; this assumes that
            // all the species have the same eddy diffusivity.
            K_turb(i,j,k,EddyDiff::Q_v) = rho * Lm * qvel(i,j,k) * SH;
 
            K_turb(i,j,k,EddyDiff::Turb_lengthscale) = Lm;
        });
    }
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function: