ERF_ComputeTurbulentViscosity.cpp File Reference

#include "ERF_ImmersedWallCell.H"
#include "ERF_SurfaceLayer.H"
#include "ERF_Constants.H"
#include "ERF_EddyViscosity.H"
#include "ERF_RANSClosure.H"
#include "ERF_CanopyTurbulence.H"
#include "ERF_Diffusion.H"
#include "ERF_PBLModels.H"
#include "ERF_TileNoZ.H"
#include "ERF_TerrainMetrics.H"
#include "ERF_TerrainDiffusionLimits.H"
#include "ERF_MoistUtils.H"
#include "ERF_RichardsonNumber.H"

Include dependency graph for ERF_ComputeTurbulentViscosity.cpp:

Functions
void	ComputeTurbulentViscosityLES (Vector< std::unique_ptr< MultiFab >> &Tau_lev, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Diss, const Geometry &geom, bool use_terrain_fitted_coords, Vector< std::unique_ptr< MultiFab >> &mapfac, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< SurfaceLayer > &, const MoistureComponentIndices &moisture_indices, const MultiFab *xvel, const MultiFab *yvel)
void	ComputeTurbulentViscosityLES_EB (Vector< std::unique_ptr< MultiFab >> &Tau_lev, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, const Geometry &geom, Vector< std::unique_ptr< MultiFab >> &mapfac, const TurbChoice &turbChoice, const Real const_grav, [[maybe_unused]] const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &, const MoistureComponentIndices &moisture_indices, const eb_ &ebfact, const MultiFab *xvel, const MultiFab *yvel)
void	ComputeTurbulentViscosityRANS (int level, Vector< std::unique_ptr< MultiFab >> &Tau_lev, const MultiFab &cons_in, const MultiFab &wdist, MultiFab &eddyViscosity, MultiFab &Diss, const Geometry &geom, bool use_terrain_fitted_coords, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< SurfaceLayer > &SurfLayer, const MultiFab *z_0, const CanopyTurbulence &canopy, const MultiFab *forest_drag, const MultiFab *ib_blank, bool flat_wall_diss, Real ib_small_volfrac, Real ib_z0)
void	ComputeTurbulentViscosity (double dt, const MultiFab &xvel, const MultiFab &yvel, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab &cons_in, const MultiFab &wdist, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, Vector< std::unique_ptr< MultiFab >> &mapfac, const std::unique_ptr< MultiFab > &z_phys_nd, const std::unique_ptr< MultiFab > &z_phys_cc, const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, const MultiFab *z_0, const bool &use_terrain_fitted_coords, const bool &use_moisture, int level, const BCRec *bc_ptr, const eb_ &ebfact, bool vert_only, const MultiFab *qheating_rates, const MultiFab *terrain_blank, const MultiFab *forest_drag)

Function Documentation

ComputeTurbulentViscosity()

void ComputeTurbulentViscosity	(	double	dt,
		const MultiFab &	xvel,
		const MultiFab &	yvel,
		Vector< std::unique_ptr< MultiFab >> &	Tau_lev,
		MultiFab &	cons_in,
		const MultiFab &	wdist,
		MultiFab &	eddyViscosity,
		MultiFab &	Hfx1,
		MultiFab &	Hfx2,
		MultiFab &	Hfx3,
		MultiFab &	Diss,
		const Geometry &	geom,
		Vector< std::unique_ptr< MultiFab >> &	mapfac,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const std::unique_ptr< MultiFab > &	z_phys_cc,
		const SolverChoice &	solverChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		const MultiFab *	z_0,
		const bool &	use_terrain_fitted_coords,
		const bool &	use_moisture,
		int	level,
		const BCRec *	bc_ptr,
		const eb_ &	ebfact,
		bool	vert_only,
		const MultiFab *	qheating_rates,
		const MultiFab *	terrain_blank,
		const MultiFab *	forest_drag
	)

Wrapper to compute turbulent viscosity with LES or PBL.

Parameters

[in]	dt	time step
[in]	xvel	velocity in x-dir
[in]	yvel	velocity in y-dir
[in]	Tau_lev	strain at this level
[in]	cons_in	cell center conserved quantities
[in]	wdist	wall distance
[out]	eddyViscosity	turbulent viscosity
[out]	Hfx1	heat flux in x-dir
[out]	Hfx2	heat flux in y-dir
[out]	Hfx3	heat flux in z-dir
[out]	Diss	dissipation of turbulent kinetic energy
[in]	geom	problem geometry
[in]	mapfac	map factors
[in]	z_phys_nd	nodal physical z coordinates
[in]	z_phys_cc	cell-centered physical z coordinates
[in]	solverChoice	container with solver, diffusion, and turbulence parameters
[in]	SurfLayer	optional surface-layer model
[in]	z_0	roughness length
[in]	use_terrain_fitted_coords	flag for terrain-fitted coordinates
[in]	use_moisture	flag for moisture physics
[in]	level	AMR level
[in]	bc_ptr	boundary condition records
[in]	ebfact	EB factories for cell- and face-centered variables
[in]	vert_only	flag for vertical components of eddyViscosity
[in]	qheating_rates	radiation heating rates (SW, LW components)
[in]	terrain_blank	immersed forcing mask, or null
[in]	forest_drag	forest drag field C_d LAD, or null without a forest

{
    BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity);
    //
    // In LES mode, the turbulent viscosity is isotropic (unless mix_isotropic is set to false), so
    // the LES model sets both horizontal and vertical viscosities
    //
    // In PBL mode, the primary purpose of the PBL model is to control vertical transport, so the PBL model sets the vertical viscosity.
    // Optionally, the PBL model can be run in conjunction with an LES model that sets the horizontal viscosity
    // (this isn’t truly LES, but the model form is the same as Smagorinsky).
    //
    // ComputeTurbulentViscosityLES populates the LES viscosity for both horizontal and vertical components.
    // ComputeTurbulentViscosityPBL computes the PBL viscosity just for the vertical component.
    //
 
    TurbChoice turbChoice = solverChoice.turbChoice[level];
    const Real const_grav = solverChoice.gravity;
 
    if (!SurfLayer) {
        AMREX_ALWAYS_ASSERT(!vert_only);
    }
 
    bool impose_phys_bcs = true;
 
    if (turbChoice.les_type != LESType::None) {
        if (solverChoice.terrain_type == TerrainType::EB) {
            ComputeTurbulentViscosityLES_EB(Tau_lev,
                                        cons_in, eddyViscosity,
                                        Hfx1, Hfx2, Hfx3,
                                        geom,
                                        mapfac, turbChoice, const_grav,
                                        solverChoice,
                                        SurfLayer, solverChoice.moisture_indices,
                                        ebfact, &xvel, &yvel);
        } else {
            ComputeTurbulentViscosityLES(Tau_lev,
                                        cons_in, eddyViscosity,
                                        Diss,
                                        geom, use_terrain_fitted_coords,
                                        mapfac, z_phys_nd, turbChoice, const_grav,
                                        SurfLayer, solverChoice.moisture_indices,
                                        &xvel, &yvel);
        }
    }
 
    if (turbChoice.rans_type != RANSType::None) {
        ComputeTurbulentViscosityRANS(level, Tau_lev, cons_in, wdist,
                                      eddyViscosity,
                                      Diss,
                                      geom, use_terrain_fitted_coords,
                                      z_phys_nd, turbChoice, const_grav,
                                      SurfLayer, z_0,
                                      solverChoice.canopyTurb, forest_drag,
                                      (solverChoice.if_fraction_stress) ? terrain_blank : nullptr,
                                      solverChoice.if_flat_wall_dissipation,
                                      solverChoice.if_small_volfrac,
                                      solverChoice.if_z0);
    }
 
    if (turbChoice.pbl_type == PBLType::MYJ) {
        ComputeDiffusivityMYJ(dt, xvel, yvel, cons_in, eddyViscosity,
                              geom, turbChoice, SurfLayer,
                              use_terrain_fitted_coords, use_moisture,
                              level, bc_ptr, vert_only, z_phys_nd, z_phys_cc,
                              solverChoice.moisture_indices);
    } else if (turbChoice.pbl_type == PBLType::MYNN25) {
        ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity,
                                 geom, turbChoice, SurfLayer,
                                 use_terrain_fitted_coords, use_moisture,
                                 level, bc_ptr, vert_only, z_phys_nd, z_phys_cc,
                                 solverChoice.moisture_indices);
    } else if (turbChoice.pbl_type == PBLType::MYNNEDMF) {
        ComputeDiffusivityMYNNEDMF(xvel, yvel, cons_in, eddyViscosity,
                                   geom, turbChoice, SurfLayer,
                                   use_terrain_fitted_coords, use_moisture,
                                   level, bc_ptr, vert_only, z_phys_nd, z_phys_cc,
                                   solverChoice.moisture_indices);
    } else if (turbChoice.pbl_type == PBLType::YSU) {
        ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity,
                              geom, turbChoice, SurfLayer,
                              use_terrain_fitted_coords, use_moisture,
                              level, bc_ptr, vert_only, z_phys_nd, z_phys_cc,
                              solverChoice.moisture_indices);
    } else if (turbChoice.pbl_type == PBLType::MRF) {
        ComputeDiffusivityMRF(xvel, yvel, cons_in, eddyViscosity,
                              geom, turbChoice, SurfLayer,
                              use_terrain_fitted_coords, use_moisture,
                              level, bc_ptr, vert_only, z_phys_nd, z_phys_cc,
                              solverChoice.moisture_indices, terrain_blank);
    } else if (turbChoice.pbl_type == PBLType::YSUNew) {
        ComputeDiffusivityYSUNew(xvel, yvel, cons_in, eddyViscosity,
                                 geom, turbChoice, SurfLayer,
                                 use_terrain_fitted_coords, use_moisture,
                                 level, bc_ptr, vert_only, z_phys_nd, z_phys_cc,
                                 solverChoice.moisture_indices,
                                 qheating_rates, terrain_blank);
    } else if (turbChoice.uses_shoc_family()) {
        // NOTE: Nothing to do here. The SHOC class handles setting the vertical
        //       components of eddyDiffs in slow RHS pre.
    }
 
    //
    // At all levels we need to fill values outside the physical boundary for the LES coeffs.
    // In addition, for all cases, if at level > 0, we want to fill fine ghost cell values that
    // overlie coarse grid cells (and that are not in another fine valid region) with
    // extrapolated values from the interior, rather than interpolating from the coarser level,
    // since we may be using a different turbulence model there.
    //
    // Note: here "covered" refers to "covered by valid region of another grid at this level"
    // Note: here "physbnd" refers to "cells outside the domain if not periodic"
    // Note: here "interior" refers to "valid cells, i.e. inside 'my' grid"
    //
    {
        int is_covered    = 0;
        int is_notcovered = 1;
        int is_physbnd    = 2;
        int is_interior   = 3;
        iMultiFab cc_mask(eddyViscosity.boxArray(),eddyViscosity.DistributionMap(),1,1);
        cc_mask.BuildMask(geom.Domain(), geom.periodicity(), is_covered, is_notcovered, is_physbnd, is_interior);
 
        Box domain = geom.Domain();
        for (int i = 0; i < AMREX_SPACEDIM; ++i) {
            if (geom.isPeriodic(i)) {
                domain.grow(i,1);
            }
        }
 
        eddyViscosity.FillBoundary(geom.periodicity());
 
        int ncomp  = eddyViscosity.nComp();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for (MFIter mfi(eddyViscosity); mfi.isValid(); ++mfi)
        {
            Box vbx = mfi.validbox();
 
            Box planex_lo = mfi.growntilebox(1); planex_lo.setBig(0, vbx.smallEnd(0)-1);
            Box planey_lo = mfi.growntilebox(1); planey_lo.setBig(1, vbx.smallEnd(1)-1);
            Box planez_lo = mfi.growntilebox(1); planez_lo.setBig(2, vbx.smallEnd(2)-1);
 
            Box planex_hi = mfi.growntilebox(1); planex_hi.setSmall(0, vbx.bigEnd(0)+1);
            Box planey_hi = mfi.growntilebox(1); planey_hi.setSmall(1, vbx.bigEnd(1)+1);
            Box planez_hi = mfi.growntilebox(1); planez_hi.setSmall(2, vbx.bigEnd(2)+1);
 
            int i_lo   = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
            int j_lo   = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1);
            int k_lo   = vbx.smallEnd(2); int k_hi = vbx.bigEnd(2);
 
            const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
            const Array4<int>& mask_arr = cc_mask.array(mfi);
 
            auto domlo = lbound(domain);
            auto domhi = ubound(domain);
 
            ParallelFor(planex_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                int lj = amrex::min(amrex::max(j, domlo.y), domhi.y);
                int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
                if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i_lo,lj,lk) != is_notcovered) ||
                    (mask_arr(i,j,k) == is_physbnd    && i < domlo.x && impose_phys_bcs)) {
                    for (int n = 0; n < ncomp; n++) {
                        mu_turb(i,j,k,n) = mu_turb(i_lo,lj,lk,n);
                    }
                }
            });
            ParallelFor(planex_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                int lj = amrex::min(amrex::max(j, domlo.y), domhi.y);
                int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
                if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i_hi,lj,lk) != is_notcovered) ||
                    (mask_arr(i,j,k) == is_physbnd    && i > domhi.x && impose_phys_bcs)) {
                    for (int n = 0; n < ncomp; n++) {
                        mu_turb(i,j,k,n) = mu_turb(i_hi,lj,lk,n);
                    }
                }
            });
            ParallelFor(planey_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
                if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j_lo,lk) != is_notcovered) ||
                    (mask_arr(i,j,k) == is_physbnd    && j < domlo.y && impose_phys_bcs)) {
                    for (int n = 0; n < ncomp; n++) {
                        mu_turb(i,j,k,n) = mu_turb(i,j_lo,lk,n);
                    }
                }
            });
            ParallelFor(planey_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
                if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j_hi,lk) != is_notcovered)||
                    (mask_arr(i,j,k) == is_physbnd    && j > domhi.y && impose_phys_bcs)) {
                    for (int n = 0; n < ncomp; n++) {
                        mu_turb(i,j,k,n) = mu_turb(i,j_hi,lk,n);
                    }
                }
            });
            ParallelFor(planez_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j,k_lo) != is_notcovered) ||
                    (mask_arr(i,j,k) == is_physbnd    && k < domlo.z && impose_phys_bcs)) {
                    for (int n = 0; n < ncomp; n++) {
                        mu_turb(i,j,k,n) = mu_turb(i,j,k_lo,n);
                    }
                }
            });
            ParallelFor(planez_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j,k_hi) != is_notcovered) ||
                    (mask_arr(i,j,k) == is_physbnd    && k > domhi.z && impose_phys_bcs)) {
                    for (int n = 0; n < ncomp; n++) {
                        mu_turb(i,j,k,n) = mu_turb(i,j,k_hi,n);
                    }
                }
            });
            //
            // Edge and corner ghost cells at a C/F boundary are skipped above when every face
            // neighbour toward the valid box is itself a C/F ghost cell, e.g. (i_lo-1,j_lo-1).
            // Fill these from a cell that is valid at this level (in this grid or another),
            // stepping toward the valid box in two directions (xy, xz, yz) before all three.
            // The sources are never written in this block, and taking covered cells as well as
            // our own keeps the result independent of how the level is divided into grids.
            //
            ParallelFor(mfi.growntilebox(1), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (mask_arr(i,j,k) != is_notcovered) { return; }
 
                int si = amrex::min(amrex::max(i, i_lo), i_hi);
                int sj = amrex::min(amrex::max(j, j_lo), j_hi);
                int sk = amrex::min(amrex::max(k, k_lo), k_hi);
                bool gx = (si != i);
                bool gy = (sj != j);
                bool gz = (sk != k);
 
                // Already filled by one of the plane kernels above
                if (gx && mask_arr(si,j,k) != is_notcovered) { return; }
                if (gy && mask_arr(i,sj,k) != is_notcovered) { return; }
                if (gz && mask_arr(i,j,sk) != is_notcovered) { return; }
 
                int ii = si; int jj = sj; int kk = sk;
                if        (gx && gy && mask_arr(si,sj,k ) != is_notcovered) {
                    kk = k;
                } else if (gx && gz && mask_arr(si,j ,sk) != is_notcovered) {
                    jj = j;
                } else if (gy && gz && mask_arr(i ,sj,sk) != is_notcovered) {
                    ii = i;
                }
                for (int n = 0; n < ncomp; n++) {
                    mu_turb(i,j,k,n) = mu_turb(ii,jj,kk,n);
                }
            });
        } // mfi
 
        eddyViscosity.FillBoundary(geom.periodicity());
    }
}

Referenced by ERF::advance_dycore().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeTurbulentViscosityLES()

void ComputeTurbulentViscosityLES	(	Vector< std::unique_ptr< MultiFab >> &	Tau_lev,
		const MultiFab &	cons_in,
		MultiFab &	eddyViscosity,
		MultiFab &	Diss,
		const Geometry &	geom,
		bool	use_terrain_fitted_coords,
		Vector< std::unique_ptr< MultiFab >> &	mapfac,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const TurbChoice &	turbChoice,
		const Real	const_grav,
		std::unique_ptr< SurfaceLayer > &	,
		const MoistureComponentIndices &	moisture_indices,
		const MultiFab *	xvel,
		const MultiFab *	yvel
	)

Function for computing the turbulent viscosity with LES.

Parameters

[in]	Tau_lev	strain at this level
[in]	cons_in	cell center conserved quantities
[out]	eddyViscosity	turbulent viscosity
[out]	Diss	dissipation of turbulent kinetic energy
[in]	geom	problem geometry
[in]	use_terrain_fitted_coords	flag for terrain-fitted coordinates
[in]	mapfac	map factors
[in]	z_phys_nd	nodal physical z coordinates
[in]	turbChoice	container with turbulence parameters
[in]	const_grav	gravitational acceleration
[in]	moisture_indices	moisture component indices
[in]	xvel	x-direction velocity (for moist Ri correction)
[in]	yvel	y-direction velocity (for moist Ri correction)

{
    const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
    const Box& domain = geom.Domain();
 
    Real inv_Pr_t    = turbChoice.Pr_t_inv;
    Real inv_Sc_t    = turbChoice.Sc_t_inv;
    Real inv_sigma_k = one / turbChoice.sigma_k;
 
    bool use_thetav_grad = (turbChoice.strat_type == StratType::thetav);
    bool use_thetal_grad = (turbChoice.strat_type == StratType::thetal);
 
    bool isotropic = turbChoice.mix_isotropic;
 
    // SMAGORINSKY: Fill Kturb for momentum in horizontal and vertical
    //***********************************************************************************
    if (turbChoice.les_type == LESType::Smagorinsky)
    {
        Real Cs = turbChoice.Cs;
        bool smag2d = turbChoice.smag2d;
 
        // Opt-in limits of WRF's smag2d_km on the horizontal viscosity (see
        // Smag2DLimitedViscosity).  The slope limiter needs the terrain slope, so it acts
        // on terrain-fitted meshes only; SolverChoice aborts if it is asked for elsewhere.
        bool l_smag2d_slope = smag2d && turbChoice.smag2d_slope_limiter && use_terrain_fitted_coords;
        Real l_smag2d_cap   = (smag2d) ? turbChoice.smag2d_kh_cap : zero;
        bool l_smag2d_limit = l_smag2d_slope || (l_smag2d_cap > zero);
 
        // Define variables required inside device lambdas (scalars only)
        Real l_abs_g = const_grav;
        Real l_Ri_crit = turbChoice.Ri_crit;
        bool l_use_Ri_corr = turbChoice.use_Ri_correction;
        bool l_has_xvel = (xvel != nullptr);
        bool l_has_yvel = (yvel != nullptr);
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for (MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            Box bxcc  = mfi.growntilebox(1) & domain;
            // NOTE: the closures deliberately store no subgrid heat flux.  hfx_z (Hfx3)
            //       is z-nodal, (0,0,1), while the closure's -K dtheta/dz is cell-centred,
            //       and the theta diffusion of every RK stage overwrites all z-faces (the
            //       surface layer the bottom face) before anything reads them, so the value
            //       written here was both misplaced and unused.  The TKE buoyancy source
            //       reads the face fluxes of that diffusion (ERF_TurbKESources.H).
            const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
            const Array4<Real const >& cell_data = cons_in.array(mfi);
            const StrainArrays tau = strain_arrays(Tau_lev, mfi);
            const Array4<Real const> tau11 = tau.tau11, tau22 = tau.tau22, tau33 = tau.tau33;
            const Array4<Real const> tau12 = tau.tau12, tau13 = tau.tau13, tau23 = tau.tau23;
            Array4<Real const> mf_u = mapfac[MapFacType::u_x]->const_array(mfi);
            Array4<Real const> mf_v = mapfac[MapFacType::v_y]->const_array(mfi);
            Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
 
            Array4<Real const> u_arr = (l_has_xvel) ? xvel->const_array(mfi) : Array4<Real const>{};
            Array4<Real const> v_arr = (l_has_yvel) ? yvel->const_array(mfi) : Array4<Real const>{};
 
            ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                // =====================================================================
                // 1. STRAIN RATE MAGNITUDE CALCULATION
                // =====================================================================
                Real SmnSmn;
                if (smag2d) {
                    SmnSmn = ComputeSmnSmn2D(i,j,k,tau11,tau22,tau12);
                } else {
                    SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23);
                }
                Real strain_rate_magnitude = std::sqrt(two * SmnSmn);
 
                // =====================================================================
                // 2. GRID SCALE CALCULATION (filter width Δ)
                // =====================================================================
                Real dxInv = cellSizeInv[0];
                Real dyInv = cellSizeInv[1];
                Real dzInv = cellSizeInv[2];
                if (use_terrain_fitted_coords) {
                    dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
                }
 
                Real Delta;
                Real DeltaH;
                if (isotropic) {
                    Real cellVolMsf = one / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
                    Delta = std::cbrt(cellVolMsf);
                    DeltaH = Delta;
                } else {
                    Delta = one / dzInv;
                    DeltaH = std::sqrt(one / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0)));
                }
 
                Real rho = cell_data(i, j, k, Rho_comp);
                Real CsDeltaSqr_h = Cs * Cs * DeltaH * DeltaH;
                Real CsDeltaSqr_v = Cs * Cs * Delta * Delta;
 
                Real nu_turb_base_h = CsDeltaSqr_h * strain_rate_magnitude;
                Real nu_turb_base_v = CsDeltaSqr_v * strain_rate_magnitude;
 
                // WRF smag2d_km limits on K_h only, applied before the scalar diffusivities
                // (Theta_h, Scalar_h, Q_h) are derived from Mom_h below; Mom_v is unchanged.
                if (l_smag2d_limit) {
                    Real alpha = one;
                    if (l_smag2d_slope) {
                        alpha = TerrainSlopeFactor(ComputeTerrainCellDrops(i,j,k,z_nd_arr));
                    }
                    nu_turb_base_h = Smag2DLimitedViscosity(nu_turb_base_h, strain_rate_magnitude,
                                                            DeltaH, alpha, l_smag2d_cap, l_smag2d_slope);
                }
 
                Real stability_factor = one;
 
                if (l_use_Ri_corr && l_has_xvel && l_has_yvel) {
                    Real N2 = ComputeN2(i, j, k, dzInv, l_abs_g, cell_data, moisture_indices);
                    Real S2_vert = ComputeVerticalShear2(i, j, k, dzInv, u_arr, v_arr);
                    Real Ri = ComputeRichardson(N2, S2_vert);
                    stability_factor = StabilityFunction(Ri, l_Ri_crit);
                }
 
                if (isotropic) {
                    mu_turb(i, j, k, EddyDiff::Mom_h) = rho * nu_turb_base_h * stability_factor;
                    mu_turb(i, j, k, EddyDiff::Mom_v) = rho * nu_turb_base_v * stability_factor;
                } else {
                    mu_turb(i, j, k, EddyDiff::Mom_h) = rho * nu_turb_base_h;
                    mu_turb(i, j, k, EddyDiff::Mom_v) = rho * nu_turb_base_v * stability_factor;
                }
            });
        }
    }
    // DEARDORFF: Fill Kturb for momentum in horizontal and vertical
    //***********************************************************************************
    else if (turbChoice.les_type == LESType::Deardorff)
    {
        const Real l_C_k        = turbChoice.Ck;
        const Real l_C_e        = turbChoice.Ce;
        const Real l_C_e_wall   = turbChoice.Ce_wall;
        const Real Ce_lcoeff    = amrex::max(amrex::Real(zero), l_C_e - Real(1.9)*l_C_k);
        const Real l_abs_g      = const_grav;
 
        // Clamped divisor: the select is if-converted, so 1/theta_ref runs even
        //   when theta_ref = 0 and would trip fpe_trap_zero (see ERF_SetupDiff.H)
        const bool use_ref_theta = (turbChoice.theta_ref > 0);
        const Real inv_theta_ref = one / amrex::max(turbChoice.theta_ref, std::numeric_limits<Real>::min());
        const Real l_inv_theta0  = (use_ref_theta) ? inv_theta_ref : one;
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            Box bxcc  = mfi.tilebox();
 
            // NOTE: the closures deliberately store no subgrid heat flux.  hfx_z (Hfx3)
            //       is z-nodal, (0,0,1), while the closure's -K dtheta/dz is cell-centred,
            //       and the theta diffusion of every RK stage overwrites all z-faces (the
            //       surface layer the bottom face) before anything reads them, so the value
            //       written here was both misplaced and unused.  The TKE buoyancy source
            //       reads the face fluxes of that diffusion (ERF_TurbKESources.H).
            const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
            const Array4<Real>& diss    = Diss.array(mfi);
 
            const Array4<Real const > &cell_data = cons_in.array(mfi);
 
            Array4<Real const> mf_u = mapfac[MapFacType::u_x]->const_array(mfi);
            Array4<Real const> mf_v = mapfac[MapFacType::v_y]->const_array(mfi);
 
            Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
 
            ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                Real dxInv = cellSizeInv[0];
                Real dyInv = cellSizeInv[1];
                Real dzInv = cellSizeInv[2];
                if (use_terrain_fitted_coords) {
                    // the terrain grid is only deformed in z for now
                    dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
                }
                Real Delta;
                if (isotropic) {
                    Real cellVolMsf = one / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
                    Delta = std::cbrt(cellVolMsf);
                } else {
                    Delta = one / dzInv;
                }
 
                Real dtheta_dz;
                if (use_thetav_grad) {
                    dtheta_dz = myhalf * ( GetThetav(i, j, k+1, cell_data, moisture_indices)
                                          -GetThetav(i, j, k-1, cell_data, moisture_indices) )*dzInv;
                } else if (use_thetal_grad) {
                    dtheta_dz = myhalf * ( GetThetal(i, j, k+1, cell_data, moisture_indices)
                                          -GetThetal(i, j, k-1, cell_data, moisture_indices) )*dzInv;
                } else {
                    dtheta_dz = myhalf * ( cell_data(i, j, k+1, RhoTheta_comp) / cell_data(i, j, k+1, Rho_comp)
                                          -cell_data(i, j, k-1, RhoTheta_comp) / cell_data(i, j, k-1, Rho_comp) )*dzInv;
                }
 
                // Calculate stratification-dependent mixing length (Deardorff 1980, Eqn. 10a)
                Real E              = amrex::max(cell_data(i,j,k,RhoKE_comp)/cell_data(i,j,k,Rho_comp),Real(0.0));
                Real stratification = l_abs_g * dtheta_dz * l_inv_theta0;
                if (!use_ref_theta) {
                    // l_inv_theta0 == 1, divide by actual theta
                    stratification *= cell_data(i,j,k,Rho_comp) /
                                      cell_data(i,j,k,RhoTheta_comp);
                }
 
                // Following WRF, the stratification effects are applied to the vertical length scales
                // in the case of anisotropic mixing
                Real length;
                Real eps       = std::numeric_limits<Real>::epsilon();
                if (stratification <= eps) {
                    length = Delta;  // cbrt(dx*dy*dz) -or- dz
                } else {
                    length = Real(0.76) * std::sqrt(E / amrex::max(stratification,eps));
                    // mixing length should be _reduced_ for stable stratification
                    length = amrex::min(length, Delta);
                    // following WRF, make sure the mixing length isn't too small
                    length = amrex::max(length, Real(0.001) * Delta);
                }
 
                Real DeltaH = (isotropic) ? length : std::sqrt(one / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0)));
 
                Real Pr_inv_v = (one + two*length/Delta);
                Real Pr_inv_h  = (isotropic) ? Pr_inv_v : inv_Pr_t;
 
                // Calculate eddy diffusivities
                // K = rho * C_k * l * KE^(1/2)
                mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * DeltaH * std::sqrt(E);
                mu_turb(i,j,k,EddyDiff::Mom_v) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E);
                // KH = (1 + 2*l/delta) * mu_turb
                mu_turb(i,j,k,EddyDiff::Theta_h) = Pr_inv_h * mu_turb(i,j,k,EddyDiff::Mom_h);
                mu_turb(i,j,k,EddyDiff::Theta_v) = Pr_inv_v * mu_turb(i,j,k,EddyDiff::Mom_v);
                // Store lengthscale for TKE source terms
                mu_turb(i,j,k,EddyDiff::Turb_lengthscale) = length;
 
                // Calculate SFS quantities
                // - dissipation
                Real Ce;
                if ((l_C_e_wall > 0) && (k==0)) {
                    Ce = l_C_e_wall;
                } else {
                    Ce = Real(1.9)*l_C_k + Ce_lcoeff*length / Delta;
                }
                diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,Real(1.5)) / length;
            });
        }
    }
 
    // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0)
    //***********************************************************************************
    int ngc(1);
    // EddyDiff mapping :   Theta_h     KE_h       Scalar_h    Q_h
    Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
    Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
    Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
    Real* fac_ptr = d_Factors.data();
 
    const bool use_KE = ( turbChoice.les_type == LESType::Deardorff );
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        Box bxcc   = mfi.tilebox();
        Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
        Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1);
        //
        // NOTE: growntilebox, not the tilebox grown by hand.  Growing by hand pushes every tile
        //       ngc cells past its own share of the grid, so once the grid is tiled two tiles
        //       store to the same mu_turb cells -- a data race under OpenMP, even though both
        //       store the same value.  The union of the tiles, and the values stored, are
        //       unchanged, and this is identical to the old expression with one tile per grid.
        //
        bxcc = mfi.growntilebox(IntVect(ngc,ngc,0));
 
        const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
 
        for (auto n = 1; n < EddyDiff::Mom_v; ++n) {
            int offset = EddyDiff::Mom_v;
            switch (n)
            {
             case EddyDiff::KE_h:
                if (use_KE) {
                   ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                   {
                       int indx   = n;
                       int indx_v = indx + offset;
                       mu_turb(i,j,k,indx)   = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
                       mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,EddyDiff::Mom_v) * fac_ptr[indx-1];
                   });
                }
                break;
            default:
                ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                {
                    int indx   = n;
                    int indx_v = indx + offset;
 
                    // NOTE: Theta_h, Theta_v have already been set for Deardorff
                    if (!(indx_v == EddyDiff::Theta_v && use_KE)) {
                        mu_turb(i,j,k,indx)   = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
                        mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,EddyDiff::Mom_v) * fac_ptr[indx-1];
                    }
                });
                break;
          }
       }
    }
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeTurbulentViscosityLES_EB()

void ComputeTurbulentViscosityLES_EB	(	Vector< std::unique_ptr< MultiFab >> &	Tau_lev,
		const MultiFab &	cons_in,
		MultiFab &	eddyViscosity,
		MultiFab &	Hfx1,
		MultiFab &	Hfx2,
		MultiFab &	Hfx3,
		const Geometry &	geom,
		Vector< std::unique_ptr< MultiFab >> &	mapfac,
		const TurbChoice &	turbChoice,
		const Real	const_grav,
		[[maybe_unused] ] const SolverChoice &	solverChoice,
		std::unique_ptr< SurfaceLayer > &	,
		const MoistureComponentIndices &	moisture_indices,
		const eb_ &	ebfact,
		const MultiFab *	xvel,
		const MultiFab *	yvel
	)

{
    const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
    Real dxInv = cellSizeInv[0];
    Real dyInv = cellSizeInv[1];
    Real dzInv = cellSizeInv[2];
 
    const Box& domain = geom.Domain();
 
    Real inv_Pr_t    = turbChoice.Pr_t_inv;
    Real inv_Sc_t    = turbChoice.Sc_t_inv;
    Real inv_sigma_k = one / turbChoice.sigma_k;
 
    bool isotropic = turbChoice.mix_isotropic;
 
    AMREX_ASSERT(turbChoice.les_type == LESType::Smagorinsky);
 
    // SMAGORINSKY: Fill Kturb for momentum in horizontal and vertical
    //***********************************************************************************
    if (turbChoice.les_type == LESType::Smagorinsky)
    {
        Real Cs = turbChoice.Cs;
 
        // Define variables required inside device lambdas (scalars only)
        Real l_abs_g = const_grav;
        Real l_Ri_crit = turbChoice.Ri_crit;
        bool l_use_Ri_corr = turbChoice.use_Ri_correction;
        bool l_has_xvel = (xvel != nullptr);
        bool l_has_yvel = (yvel != nullptr);
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for (MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            Box bxcc  = mfi.growntilebox(1) & domain;
            const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
            const Array4<Real>& hfx_x   = Hfx1.array(mfi);
            const Array4<Real>& hfx_y   = Hfx2.array(mfi);
            const Array4<Real>& hfx_z   = Hfx3.array(mfi);
            const Array4<Real const >& cell_data = cons_in.array(mfi);
            const StrainArrays tau = strain_arrays(Tau_lev, mfi);
            const Array4<Real const> tau11 = tau.tau11, tau22 = tau.tau22, tau33 = tau.tau33;
            const Array4<Real const> tau12 = tau.tau12, tau13 = tau.tau13, tau23 = tau.tau23;
            Array4<Real const> mf_u = mapfac[MapFacType::u_x]->const_array(mfi);
            Array4<Real const> mf_v = mapfac[MapFacType::v_y]->const_array(mfi);
 
            Array4<Real const> u_arr = (l_has_xvel) ? xvel->const_array(mfi) : Array4<Real const>{};
            Array4<Real const> v_arr = (l_has_yvel) ? yvel->const_array(mfi) : Array4<Real const>{};
 
            Array4<const EBCellFlag> c_cflag = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
            Array4<const EBCellFlag> u_cflag = (ebfact.get_u_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
            Array4<const EBCellFlag> v_cflag = (ebfact.get_v_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
            Array4<const EBCellFlag> w_cflag = (ebfact.get_w_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
            Array4<const Real      > u_vfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
            Array4<const Real      > v_vfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
 
            ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (!c_cflag(i,j,k).isCovered()) {
                    // =====================================================================
                    // 1. STRAIN RATE MAGNITUDE CALCULATION
                    // =====================================================================
                    Real SmnSmn = zero;
                    if (c_cflag(i,j,k).isRegular()) {
                        SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23);
                    } else if (c_cflag(i,j,k).isSingleValued()) {
                        SmnSmn = ComputeSmnSmn_EB(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,c_cflag,u_cflag,v_cflag,w_cflag);
                    }
                    Real strain_rate_magnitude = std::sqrt(two * SmnSmn);
 
                    // =====================================================================
                    // 2. GRID SCALE CALCULATION (filter width Δ)
                    // =====================================================================
                    Real Delta;
                    Real DeltaH;
                    if (isotropic) {
                        Real cellVolMsf = one / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
                        Delta = std::cbrt(cellVolMsf);
                        DeltaH = Delta;
                    } else {
                        Delta = one / dzInv;
                        DeltaH = std::sqrt(one / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0)));
                    }
 
                    Real rho = cell_data(i, j, k, Rho_comp);
                    Real CsDeltaSqr_h = Cs * Cs * DeltaH * DeltaH;
                    Real CsDeltaSqr_v = Cs * Cs * Delta * Delta;
 
                    Real nu_turb_base_h = CsDeltaSqr_h * strain_rate_magnitude;
                    Real nu_turb_base_v = CsDeltaSqr_v * strain_rate_magnitude;
 
                    Real stability_factor = one;
 
                    if (l_use_Ri_corr && l_has_xvel && l_has_yvel) {
                        Real N2 = zero;
                        Real S2_vert = zero;
                        // Covered and multi-valued cells fall through both branches below
                        // with no shear computed, so start from "not valid" and let the
                        // branch that does the work say otherwise.
                        bool shear_is_valid = false;
                        if (c_cflag(i,j,k).isRegular()) {
                            N2 = ComputeN2(i, j, k, dzInv, l_abs_g, cell_data, moisture_indices);
                            S2_vert = ComputeVerticalShear2(i, j, k, dzInv, u_arr, v_arr);
                            shear_is_valid = true;
                        } else if (c_cflag(i,j,k).isSingleValued()) {
                            N2 = ComputeN2_EB(i, j, k, c_cflag, dzInv, l_abs_g, cell_data, moisture_indices);
                            S2_vert = ComputeVerticalShear2_EB(i, j, k, c_cflag, u_vfrac, v_vfrac, dzInv, u_arr, v_arr, shear_is_valid);
                        }
                        // When the EB geometry leaves no usable shear stencil we keep
                        // stability_factor = 1, i.e. no stratification damping: without a
                        // trustworthy gradient we would rather leave the Smagorinsky
                        // viscosity alone than suppress it on the strength of a Ri we do
                        // not believe. This is a deliberate choice -- damping to 0 instead
                        // would be just as defensible -- so change it knowingly.
                        if (shear_is_valid) {
                            Real Ri = ComputeRichardson(N2, S2_vert);
                            stability_factor = StabilityFunction(Ri, l_Ri_crit);
                        }
                    }
 
                    if (isotropic) {
                        mu_turb(i, j, k, EddyDiff::Mom_h) = rho * nu_turb_base_h * stability_factor;
                        mu_turb(i, j, k, EddyDiff::Mom_v) = rho * nu_turb_base_v * stability_factor;
                    } else {
                        mu_turb(i, j, k, EddyDiff::Mom_h) = rho * nu_turb_base_h;
                        mu_turb(i, j, k, EddyDiff::Mom_v) = rho * nu_turb_base_v * stability_factor;
                    }
 
                    Real dtheta_dz = zero;
 
                    if (c_cflag(i,j,k).isSingleValued() && c_cflag(i,j,k+1).isCovered()) {
                        amrex::Real theta    = cell_data(i, j, k,   RhoTheta_comp) / rho;
                        amrex::Real theta_km1 = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
                        amrex::Real theta_km2 = cell_data(i,j,k-2,RhoTheta_comp)/cell_data(i,j,k-2,Rho_comp);
                        dtheta_dz = (three*theta - Real(4.)*theta_km1 + theta_km2) * myhalf * dzInv;
                    } else if (c_cflag(i,j,k).isSingleValued() && c_cflag(i,j,k-1).isCovered()) {
                        amrex::Real theta    = cell_data(i, j, k,   RhoTheta_comp) / rho;
                        amrex::Real theta_kp1 = cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp);
                        amrex::Real theta_kp2 = cell_data(i,j,k+2,RhoTheta_comp)/cell_data(i,j,k+2,Rho_comp);
                        dtheta_dz = (-theta_kp2 + Real(4.)*theta_kp1 - three*theta) * myhalf * dzInv;
                    } else if (!c_cflag(i,j,k).isCovered()) {
                        amrex::Real theta_km1 = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
                        amrex::Real theta_kp1 = cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp);
                        dtheta_dz = myhalf * (theta_kp1 - theta_km1) * dzInv;
                    }
 
                    // The EB diffusion does not write hfx_z, so on EB terrain this is
                    // the only writer of the SGS heat flux output (the surface layer
                    // overwrites the bottom).  Elsewhere the theta diffusion writes the
                    // face fluxes and the closures do not store a heat flux.  This write
                    // keeps the z-nodal/cell-centred mismatch noted above; moving the EB
                    // heat-flux output to the EB face fluxes belongs with that diffusion.
                    hfx_x(i,j,k) = zero;
                    hfx_y(i,j,k) = zero;
                    hfx_z(i,j,k) = -inv_Pr_t * mu_turb(i,j,k,EddyDiff::Mom_v) * dtheta_dz;
 
                } else {
 
                    hfx_x(i,j,k) = zero;
                    hfx_y(i,j,k) = zero;
                    hfx_z(i,j,k) = zero;
                }// End of if not covered
            });
        }
    }
 
    // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0)
    //***********************************************************************************
    int ngc(1);
    // EddyDiff mapping :   Theta_h     KE_h       Scalar_h    Q_h
    Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
    Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
    Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
    Real* fac_ptr = d_Factors.data();
 
    const bool use_KE = ( turbChoice.les_type == LESType::Deardorff );
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        Box bxcc   = mfi.tilebox();
        Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
        Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1);
        //
        // NOTE: growntilebox, not the tilebox grown by hand.  Growing by hand pushes every tile
        //       ngc cells past its own share of the grid, so once the grid is tiled two tiles
        //       store to the same mu_turb cells -- a data race under OpenMP, even though both
        //       store the same value.  The union of the tiles, and the values stored, are
        //       unchanged, and this is identical to the old expression with one tile per grid.
        //
        bxcc = mfi.growntilebox(IntVect(ngc,ngc,0));
 
        const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
 
        for (auto n = 1; n < EddyDiff::Mom_v; ++n) {
            int offset = EddyDiff::Mom_v;
            switch (n)
            {
             case EddyDiff::KE_h:
                if (use_KE) {
                   ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                   {
                       int indx   = n;
                       int indx_v = indx + offset;
                       mu_turb(i,j,k,indx)   = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
                       mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,EddyDiff::Mom_v) * fac_ptr[indx-1];
                   });
                }
                break;
            default:
                ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                {
                    int indx   = n;
                    int indx_v = indx + offset;
 
                    // NOTE: Theta_h, Theta_v have already been set for Deardorff
                    if (!(indx_v == EddyDiff::Theta_v && use_KE)) {
                        mu_turb(i,j,k,indx)   = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
                        mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,EddyDiff::Mom_v) * fac_ptr[indx-1];
                    }
                });
                break;
          }
       }
    }
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeTurbulentViscosityRANS()

void ComputeTurbulentViscosityRANS	(	int	level,
		Vector< std::unique_ptr< MultiFab >> &	Tau_lev,
		const MultiFab &	cons_in,
		const MultiFab &	wdist,
		MultiFab &	eddyViscosity,
		MultiFab &	Diss,
		const Geometry &	geom,
		bool	use_terrain_fitted_coords,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const TurbChoice &	turbChoice,
		const Real	const_grav,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		const MultiFab *	z_0,
		const CanopyTurbulence &	canopy,
		const MultiFab *	forest_drag,
		const MultiFab *	ib_blank,
		bool	flat_wall_diss,
		Real	ib_small_volfrac,
		Real	ib_z0
	)

Function for computing the eddy viscosity with RANS.

Parameters

[in]	Tau_lev	strain at this level (the realizability cap reads S_ij)
[in]	cons_in	cell center conserved quantities
[in]	wdist	wall distance
[out]	eddyViscosity	turbulent viscosity
[out]	Diss	dissipation of turbulent kinetic energy
[in]	geom	problem geometry
[in]	use_terrain_fitted_coords	flag for terrain-fitted coordinates
[in]	z_phys_nd	nodal physical z coordinates
[in]	turbChoice	container with turbulence parameters
[in]	const_grav	gravitational acceleration
[in]	SurfLayer	optional surface-layer model
[in]	z_0	roughness length
[in]	canopy	canopy terms of the k-eqn model (erf.forest_canopy_*)
[in]	forest_drag	forest drag field C_d LAD, or null without a forest
[in]	ib_blank	solid fraction for the fraction-stress immersed wall law, or null without it
[in]	flat_wall_diss	partial wall cells dissipate with flat ground's first-cell length (erf.if_flat_wall_dissipation)
[in]	ib_small_volfrac	solid fractions below this count as fluid (eb2.small_volfrac)
[in]	ib_z0	roughness length of the immersed wall law (erf.if_z0), used without a surface layer

{
    const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
    const bool use_SurfLayer = (SurfLayer != nullptr);
 
    // Optional cap of the geometric length from the diagnosed PBL height
    const bool lscale_from_pblh = turbChoice.rans_lscale_from_pblh;
    if (lscale_from_pblh) {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(use_SurfLayer && SurfLayer->computes_pblh(),
            "erf.rans_lscale_from_pblh needs zlo.type = surface_layer and erf.most.pblh_calc = MYNN25");
    }
    const MultiFab* pblh_mf = (lscale_from_pblh) ? SurfLayer->get_pblh(level) : nullptr;
 
    Real inv_Pr_t    = turbChoice.Pr_t_inv;
    Real inv_Sc_t    = turbChoice.Sc_t_inv;
    Real inv_sigma_k = one / turbChoice.sigma_k;
 
    // One-Equation k model (Axell & Liungman 2001, Environ Fluid Mech)
    //***********************************************************************************
    if (turbChoice.rans_type == RANSType::kEqn)
    {
        const Real Cmu0       = turbChoice.Cmu0;
        const Real Cmu0_pow3  = Cmu0 * Cmu0 * Cmu0;
        const Real inv_Cb_sq  = one / (turbChoice.Cb * turbChoice.Cb);
        const Real Rt_crit    = turbChoice.Rt_crit;
        const Real Rt_min     = turbChoice.Rt_min;
        const Real l_g_max    = turbChoice.l_g_max;
        const Real l_min      = turbChoice.rans_lscale_min;
        const Real abs_g      = const_grav;
        // floor on k: erf.tke_floor if set, otherwise machine epsilon
        const Real tke_floor  = amrex::max(turbChoice.tke_floor, std::numeric_limits<Real>::epsilon());
 
        // Canopy cap of the turbulent length scale, L <= alpha / (C_d LAD)
        const bool use_canopy_length = canopy.enabled && (canopy.length_alpha > zero);
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!use_canopy_length || forest_drag,
            "erf.forest_canopy_tke needs the forest drag field");
        const Real canopy_alpha = canopy.length_alpha;
        // Realizability cap nu_t <= C_r k / S (AL01::realizable_nut); off at 0
        const Real Cr = turbChoice.rans_realizability;
        const bool use_cap = (Cr > zero);
 
        // Clamped divisor: the select is if-converted, so 1/theta_ref runs even
        //   when theta_ref = 0 and would trip fpe_trap_zero (see ERF_SetupDiff.H)
        const bool use_ref_theta = (turbChoice.theta_ref > 0);
        const Real inv_theta_ref = one / amrex::max(turbChoice.theta_ref, std::numeric_limits<Real>::min());
        const Real inv_theta0    = (use_ref_theta) ? inv_theta_ref : one;
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
        for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            Box bxcc = mfi.tilebox();
 
            const Array4<Real const>& d_arr  = wdist.const_array(mfi);
            const Array4<Real const>& z0_arr = (use_SurfLayer) ? z_0->const_array(mfi) : Array4<Real const>{};
            const Array4<Real const>& pblh_arr = (pblh_mf) ? pblh_mf->const_array(mfi) : Array4<Real const>{};
            const Array4<Real const>& f_drag_arr = (use_canopy_length) ? forest_drag->const_array(mfi)
                                                                        : Array4<Real const>{};
 
            // NOTE: the closures deliberately store no subgrid heat flux.  hfx_z (Hfx3)
            //       is z-nodal, (0,0,1), while the closure's -K dtheta/dz is cell-centred,
            //       and the theta diffusion of every RK stage overwrites all z-faces (the
            //       surface layer the bottom face) before anything reads them, so the value
            //       written here was both misplaced and unused.  The TKE buoyancy source
            //       reads the face fluxes of that diffusion (ERF_TurbKESources.H).
            const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
            const Array4<Real>& diss    = Diss.array(mfi);
 
            const Array4<Real const>& cell_data = cons_in.array(mfi);
 
            const Array4<Real const>& z_nd_arr = z_phys_nd->const_array(mfi);
            // fraction-stress immersed wall law: solid fraction, to find the wall cell
            const Array4<Real const>& ib_arr = (ib_blank) ? ib_blank->const_array(mfi) : Array4<Real const>{};
            const int ib_klo = geom.Domain().smallEnd(2);
            const int ib_khi = geom.Domain().bigEnd(2);
            const Real ib_small = ib_small_volfrac;
 
            // Strain components for the cap only; the arrays stay empty otherwise. The k-eqn
            // closure cannot run with EB terrain (SolverChoice stops the run), so the plain
            // ComputeSmnSmn is the right form here.
            const StrainArrays tau = (use_cap) ? strain_arrays(Tau_lev, mfi) : StrainArrays{};
 
            ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                Real eps = std::numeric_limits<Real>::epsilon();
                Real tke = amrex::max(cell_data(i,j,k,RhoKE_comp)/cell_data(i,j,k,Rho_comp), tke_floor);
 
                // Estimate stratification
                Real dzInv = cellSizeInv[2];
                if (use_terrain_fitted_coords) {
                    // the terrain grid is only deformed in z for now
                    dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
                }
                // In the wall cell of the fraction-stress immersed wall law the cell below is
                // solid: read the wall cell itself there, as the first-order extrapolated ghost
                // of the zlo face does on flat ground, not the solid's temperature
                int km1 = k-1;
                if (ib_arr) {
                    const Real b_m = if_rans_beta(ib_arr, i, j, k-1, ib_klo, ib_khi, ib_small);
                    const Real b_0 = if_rans_beta(ib_arr, i, j, k,   ib_klo, ib_khi, ib_small);
                    const Real b_p = if_rans_beta(ib_arr, i, j, k+1, ib_klo, ib_khi, ib_small);
                    if (ib_wall::is_wall_cell(b_m, b_0, b_p, ib_small)) { km1 = k; }
                }
                Real dtheta_dz = myhalf * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
                                          - cell_data(i,j,km1,RhoTheta_comp)/cell_data(i,j,km1,Rho_comp) )*dzInv;
                Real N2 = abs_g * inv_theta0 * dtheta_dz; // Brunt–Väisälä frequency squared
                if (!use_ref_theta) {
                    // inv_theta0 == 1, divide by actual theta
                    N2 *= cell_data(i,j,k,Rho_comp) /
                          cell_data(i,j,k,RhoTheta_comp);
                }
 
                // Geometric length scale (AL01, Eqn. 22)
                // the surface layer's roughness, or the immersed wall law's (fraction_stress) when
                // there is no surface layer at the bottom
                Real l_g = (z0_arr) ? KAPPA * (d_arr(i, j, k) + z0_arr(i, j, 0))
                         : (ib_arr) ? KAPPA * (d_arr(i, j, k) + ib_z0)
                                    : KAPPA * d_arr(i, j, k);
 
                // Enforce a maximum value: fixed, or kappa * 0.1 * zi clamped to
                // [rans_lscale_min, max_geom_lscale] when the PBL height is diagnosed
                Real l_cap = l_g_max;
                if (pblh_arr) {
                    l_cap = amrex::min(l_g_max, amrex::max(l_min, KAPPA * Real(0.1) * pblh_arr(i,j,0)));
                }
                l_g = AL01::geom_length(l_g, l_cap);
 
                // Turbulent length scale (neutral / stable Eq. 26 / unstable Eq. 28,
                // the latter bounded through the smoothed Rt)
                Real length = AL01::turb_length(l_g, N2, tke, Cmu0_pow3, inv_Cb_sq, Rt_crit, Rt_min, eps);
                // In the canopy the eddies are limited by the foliage: cap L before the
                // dissipation, Rt and the stability functions are evaluated from it
                if (f_drag_arr) {
                    length = canopy_turbulence::limited_length(length, f_drag_arr(i,j,k), canopy_alpha);
                }
                mu_turb(i, j, k, EddyDiff::Turb_lengthscale) = length;
 
                // Dissipation rate (AL01, Eqn. 19)
                diss(i, j, k) = AL01::dissipation(cell_data(i, j, k, Rho_comp), Cmu0_pow3, tke, length);
 
                // erf.if_flat_wall_dissipation (kynema-sgf KLAxell.flat_wall_dissipation): a partial
                // wall cell of the fraction-stress law has its length at the fluid centroid,
                // (1 - beta) dz / 2 from the surface, and a nearly solid one then dissipates TKE
                // several times faster than flat ground's first cell, faster than the wall value
                // restores it. Keep the centroid length for the viscosity, but dissipate with the
                // length the same formula gives at dz / 2, the distance of flat ground's first cell
                if (flat_wall_diss && ib_arr) {
                    const Real bk = ib_arr(i,j,k);
                    if (bk >= ib_small && bk < one - ib_small) {   // a partial cell as ib_wall::is_wall_cell sees it
                        Real l_flat = (z0_arr) ? KAPPA * (myhalf / cellSizeInv[2] + z0_arr(i, j, 0))
                                               : KAPPA * (myhalf / cellSizeInv[2] + ib_z0);
                        l_flat = AL01::geom_length(l_flat, l_cap);
                        const Real len_flat = AL01::turb_length(l_flat, N2, tke, Cmu0_pow3, inv_Cb_sq, Rt_crit, Rt_min, eps);
                        diss(i, j, k) = AL01::dissipation(cell_data(i, j, k, Rho_comp), Cmu0_pow3, tke, len_flat);
                    }
                }
 
                // Turbulent Richardson number (AL01, Eqn. 29 combined with Eqn. 19),
                // smoothed below Rt_crit (Burchard & Petersen)
                Real Rt = AL01::smooth_Rt(AL01::richardson(length, N2, tke, Cmu0_pow3), Rt_crit, Rt_min);
 
                // Stability functions (AL01, Eqns. 31 and 32), using the smoothed Rt
                Real cmu       = AL01::cmu(Rt, Cmu0);
                Real cmu_prime = AL01::cmu_prime(Rt, Cmu0);
                mu_turb(i, j, k, EddyDiff::RANS_Rt)        = Rt;
                mu_turb(i, j, k, EddyDiff::RANS_cmu)       = cmu;
                mu_turb(i, j, k, EddyDiff::RANS_cmu_prime) = cmu_prime;
 
                // Calculate eddy diffusivities
                // K = rho * nu_t = rho * c_mu * tke^(1/2) * length
                Real nut = cmu * std::sqrt(tke) * length; // eddy viscosity
 
                // Realizability cap nu_t <= C_r k / S; cap_ratio = nu_t / (C_r k / S)
                Real cap_ratio = zero;
                if (use_cap) {
                    Real S = std::sqrt(two * ComputeSmnSmn(i,j,k,tau.tau11,tau.tau22,tau.tau33,tau.tau12,tau.tau13,tau.tau23));
                    nut = AL01::realizable_nut(nut, tke, S, Cr, cap_ratio);
                }
                mu_turb(i, j, k, EddyDiff::RANS_cap_ratio) = cap_ratio;
 
                Real nut_prime = cmu_prime / cmu * nut;   // eddy diffusivity, capped with nu_t
                mu_turb(i, j, k, EddyDiff::Mom_h) = cell_data(i, j, k, Rho_comp) * nut;
                mu_turb(i, j, k, EddyDiff::Mom_v) = mu_turb(i, j, k, EddyDiff::Mom_h);
                mu_turb(i, j, k, EddyDiff::Theta_v) = cell_data(i, j, k, Rho_comp) * nut_prime;
            });
        }
    }
 
    // Fill the remaining eddy-diffusivity components from Mom_h. Ghost cells
    // are filled afterwards by the FillBoundary and physical-boundary
    // extrapolation in ComputeTurbulentViscosity, so only the valid cells are
    // set here.
    //***********************************************************************************
    // EddyDiff mapping :   Theta_h     KE_h       Scalar_h    Q_h
    Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
    Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
    Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
    Real* fac_ptr = d_Factors.data();
 
    const bool use_KE = ( turbChoice.rans_type == RANSType::kEqn );
    // With erf.rans_consistent_diffusivities every heat, scalar and moisture
    // diffusivity is rho * cmu' * sqrt(k) * L, i.e. the Theta_v value set by
    // the closure; without it the horizontal heat and all scalar diffusivities
    // are the eddy viscosity over Pr_t / Sc_t (the historical behaviour).
    const bool consistent = use_KE && turbChoice.rans_consistent_diffusivities;
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        Box bxcc = mfi.tilebox();
 
        const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
 
        for (auto n = 1; n < EddyDiff::Mom_v; ++n) {
            int offset = EddyDiff::Mom_v;
            switch (n)
            {
             case EddyDiff::KE_h:
                if (use_KE) {
                   ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                   {
                       int indx   = n;
                       int indx_v = indx + offset;
                       mu_turb(i,j,k,indx)   = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
                       mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
                   });
                }
                break;
            default:
                ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                {
                    int indx   = n;
                    int indx_v = indx + offset;
 
                    if (consistent) {
                        mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Theta_v);
                        if (indx_v != EddyDiff::Theta_v) {
                            mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,EddyDiff::Theta_v);
                        }
                    } else {
                        mu_turb(i,j,k,indx)   = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
 
                        // NOTE: Theta_v has already been set for the closure
                        if (!(indx_v == EddyDiff::Theta_v && use_KE)) {
                            mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
                        }
                    }
                });
                break;
          }
       }
    }
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function: