ERF_ComputeTurbulentViscosity.cpp File Reference

#include <ERF_ABLMost.H>
#include <ERF_EddyViscosity.H>
#include <ERF_Diffusion.H>
#include <ERF_PBLModels.H>
#include <ERF_TileNoZ.H>
#include <ERF_TerrainMetrics.H>

Include dependency graph for ERF_ComputeTurbulentViscosity.cpp:

Functions
void	ComputeTurbulentViscosityLES (const MultiFab &Tau11, const MultiFab &Tau22, const MultiFab &Tau33, const MultiFab &Tau12, const MultiFab &Tau13, const MultiFab &Tau23, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, bool use_terrain, const MultiFab &mapfac_u, const MultiFab &mapfac_v, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< ABLMost > &most, const bool &exp_most)
void	ComputeTurbulentViscosityRANS (const MultiFab &, const MultiFab &, const MultiFab &, const MultiFab &, const MultiFab &, const MultiFab &, const MultiFab &cons_in, const MultiFab &wdist, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, bool use_terrain, const MultiFab &, const MultiFab &, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< ABLMost > &most, const FArrayBox *z_0, const bool &exp_most)
void	ComputeTurbulentViscosity (const MultiFab &xvel, const MultiFab &yvel, const MultiFab &Tau11, const MultiFab &Tau22, const MultiFab &Tau33, const MultiFab &Tau12, const MultiFab &Tau13, const MultiFab &Tau23, const MultiFab &cons_in, const MultiFab &wdist, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, const MultiFab &mapfac_u, const MultiFab &mapfac_v, const std::unique_ptr< MultiFab > &z_phys_nd, const SolverChoice &solverChoice, std::unique_ptr< ABLMost > &most, const amrex::FArrayBox *z_0, const bool &exp_most, const bool &use_terrain_fitted_coords, const bool &use_moisture, int level, const BCRec *bc_ptr, bool vert_only)

Function Documentation

ComputeTurbulentViscosity()

void ComputeTurbulentViscosity	(	const MultiFab &	xvel,
		const MultiFab &	yvel,
		const MultiFab &	Tau11,
		const MultiFab &	Tau22,
		const MultiFab &	Tau33,
		const MultiFab &	Tau12,
		const MultiFab &	Tau13,
		const MultiFab &	Tau23,
		const MultiFab &	cons_in,
		const MultiFab &	wdist,
		MultiFab &	eddyViscosity,
		MultiFab &	Hfx1,
		MultiFab &	Hfx2,
		MultiFab &	Hfx3,
		MultiFab &	Diss,
		const Geometry &	geom,
		const MultiFab &	mapfac_u,
		const MultiFab &	mapfac_v,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const SolverChoice &	solverChoice,
		std::unique_ptr< ABLMost > &	most,
		const amrex::FArrayBox *	z_0,
		const bool &	exp_most,
		const bool &	use_terrain_fitted_coords,
		const bool &	use_moisture,
		int	level,
		const BCRec *	bc_ptr,
		bool	vert_only
	)

Wrapper to compute turbulent viscosity with LES or PBL.

Parameters

[in]	xvel	velocity in x-dir
[in]	yvel	velocity in y-dir
[in]	Tau11	11 strain
[in]	Tau22	22 strain
[in]	Tau33	33 strain
[in]	Tau12	12 strain
[in]	Tau13	13 strain
[in]	Tau23	23 strain
[in]	cons_in	cell center conserved quantities
[out]	eddyViscosity	turbulent viscosity
[in]	Hfx1	heat flux in x-dir
[in]	Hfx2	heat flux in y-dir
[in]	Hfx3	heat flux in z-dir
[in]	Diss	dissipation of turbulent kinetic energy
[in]	geom	problem geometry
[in]	mapfac_u	map factor at x-face
[in]	mapfac_v	map factor at y-face
[in]	turbChoice	container with turbulence parameters
[in]	most	pointer to Monin-Obukhov class if instantiated
[in]	vert_only	flag for vertical components of eddyViscosity

{
    BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity);
    //
    // In LES mode, the turbulent viscosity is isotropic, 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 (most) {
        bool l_use_turb = ( turbChoice.les_type  == LESType::Smagorinsky ||
                            turbChoice.les_type  == LESType::Deardorff   ||
                            turbChoice.rans_type == RANSType::kEqn       ||
                            turbChoice.pbl_type  == PBLType::MYNN25      ||
                            turbChoice.pbl_type  == PBLType::MYNNEDMF    ||
                            turbChoice.pbl_type  == PBLType::YSU );
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_use_turb,
          "A turbulence model must be utilized with MOST boundaries to compute the turbulent viscosity");
    } else {
        AMREX_ALWAYS_ASSERT(!vert_only);
    }
 
    bool impose_phys_bcs = false;
 
    if (turbChoice.les_type != LESType::None) {
        impose_phys_bcs = true;
        ComputeTurbulentViscosityLES(Tau11, Tau22, Tau33,
                                     Tau12, Tau13, Tau23,
                                     cons_in, eddyViscosity,
                                     Hfx1, Hfx2, Hfx3, Diss,
                                     geom, use_terrain_fitted_coords,
                                     mapfac_u, mapfac_v,
                                     z_phys_nd, turbChoice, const_grav,
                                     most, exp_most);
    }
 
    if (turbChoice.rans_type != RANSType::None) {
        impose_phys_bcs = true;
        ComputeTurbulentViscosityRANS(Tau11, Tau22, Tau33,
                                      Tau12, Tau13, Tau23,
                                      cons_in, wdist,
                                      eddyViscosity,
                                      Hfx1, Hfx2, Hfx3, Diss,
                                      geom, use_terrain_fitted_coords,
                                      mapfac_u, mapfac_v,
                                      z_phys_nd, turbChoice, const_grav,
                                      most, z_0, exp_most);
    }
 
    if (turbChoice.pbl_type == PBLType::MYNN25) {
        ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity,
                                 geom, turbChoice, most,
                                 use_terrain_fitted_coords, use_moisture,
                                 level, bc_ptr, vert_only, z_phys_nd,
                                 solverChoice.RhoQv_comp, solverChoice.RhoQc_comp,
                                 solverChoice.RhoQr_comp);
    } else if (turbChoice.pbl_type == PBLType::MYNNEDMF) {
        ComputeDiffusivityMYNNEDMF(xvel, yvel, cons_in, eddyViscosity,
                                   geom, turbChoice, most,
                                   use_terrain_fitted_coords, use_moisture,
                                   level, bc_ptr, vert_only, z_phys_nd,
                                   solverChoice.RhoQv_comp, solverChoice.RhoQc_comp,
                                   solverChoice.RhoQr_comp);
    } else if (turbChoice.pbl_type == PBLType::YSU) {
        ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity,
                              geom, turbChoice, most,
                              use_terrain_fitted_coords, use_moisture,
                              level, bc_ptr, vert_only, z_phys_nd);
    }
 
    //
    // 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,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,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,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,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) == is_physbnd     && k < domlo.z && impose_phys_bcs)) {
                    if (mask_arr(i,j,k) == is_physbnd) {
                    }
                    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) == 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);
                    }
                }
            });
        } // 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	(	const MultiFab &	Tau11,
		const MultiFab &	Tau22,
		const MultiFab &	Tau33,
		const MultiFab &	Tau12,
		const MultiFab &	Tau13,
		const MultiFab &	Tau23,
		const MultiFab &	cons_in,
		MultiFab &	eddyViscosity,
		MultiFab &	Hfx1,
		MultiFab &	Hfx2,
		MultiFab &	Hfx3,
		MultiFab &	Diss,
		const Geometry &	geom,
		bool	use_terrain,
		const MultiFab &	mapfac_u,
		const MultiFab &	mapfac_v,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const TurbChoice &	turbChoice,
		const Real	const_grav,
		std::unique_ptr< ABLMost > &	most,
		const bool &	exp_most
	)

Function for computing the turbulent viscosity with LES.

Parameters

[in]	Tau11	11 strain
[in]	Tau22	22 strain
[in]	Tau33	33 strain
[in]	Tau12	12 strain
[in]	Tau13	13 strain
[in]	Tau23	23 strain
[in]	cons_in	cell center conserved quantities
[out]	eddyViscosity	turbulent viscosity
[in]	Hfx1	heat flux in x-dir
[in]	Hfx2	heat flux in y-dir
[in]	Hfx3	heat flux in z-dir
[in]	Diss	dissipation of turbulent kinetic energy
[in]	geom	problem geometry
[in]	mapfac_u	map factor at x-face
[in]	mapfac_v	map factor at y-face
[in]	turbChoice	container with turbulence parameters

{
    const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
    const Box& domain = geom.Domain();
    const int& klo    = domain.smallEnd(2);
    const bool use_most = (most != nullptr);
 
    Real inv_Pr_t    = turbChoice.Pr_t_inv;
    Real inv_Sc_t    = turbChoice.Sc_t_inv;
    Real inv_sigma_k = 1.0 / turbChoice.sigma_k;
 
    // SMAGORINSKY: Fill Kturb for momentum in horizontal and vertical
    //***********************************************************************************
    if (turbChoice.les_type == LESType::Smagorinsky)
    {
      Real Cs = turbChoice.Cs;
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
      for (MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
      {
          // NOTE: This gets us the lateral ghost cells for lev>0; which
          //       have been filled from FP Two Levels.
          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);
 
          Array4<Real const> tau11 = Tau11.array(mfi);
          Array4<Real const> tau22 = Tau22.array(mfi);
          Array4<Real const> tau33 = Tau33.array(mfi);
          Array4<Real const> tau12 = Tau12.array(mfi);
          Array4<Real const> tau13 = Tau13.array(mfi);
          Array4<Real const> tau23 = Tau23.array(mfi);
 
          Array4<Real const> mf_u = mapfac_u.array(mfi);
          Array4<Real const> mf_v = mapfac_v.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 SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most,exp_most);
              Real dxInv = cellSizeInv[0];
              Real dyInv = cellSizeInv[1];
              Real dzInv = cellSizeInv[2];
              if (use_terrain) {
                  // the terrain grid is only deformed in z for now
                  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
              }
              Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
              Real DeltaMsf   = std::pow(cellVolMsf,1.0/3.0);
              Real CsDeltaSqrMsf = Cs*Cs*DeltaMsf*DeltaMsf;
 
              mu_turb(i, j, k, EddyDiff::Mom_h) = CsDeltaSqrMsf * cell_data(i, j, k, Rho_comp) * std::sqrt(2.0*SmnSmn);
              mu_turb(i, j, k, EddyDiff::Mom_v) = mu_turb(i, j, k, EddyDiff::Mom_h);
 
              // Calculate SFS quantities
              Real dtheta_dz;
              if (use_most && k==klo) {
                  if (exp_most) {
                      dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
                                  - cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp) )*dzInv;
                  } else {
                      dtheta_dz = 0.5 * (-3 * cell_data(i,j,k  ,RhoTheta_comp)
                                            / cell_data(i,j,k  ,Rho_comp)
                                        + 4 * cell_data(i,j,k+1,RhoTheta_comp)
                                            / cell_data(i,j,k+1,Rho_comp)
                                        -     cell_data(i,j,k+2,RhoTheta_comp)
                                            / cell_data(i,j,k+2,Rho_comp) ) * dzInv;
                  }
              } else {
                  dtheta_dz = 0.5 * ( 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;
              }
              // - heat flux
              //   (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will
              //    be overwritten when BCs are applied)
              hfx_x(i,j,k) = 0.0;
              hfx_y(i,j,k) = 0.0;
              hfx_z(i,j,k) = -inv_Pr_t*mu_turb(i,j,k,EddyDiff::Mom_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
          });
      }
    }
    // 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(0.0, l_C_e - 1.9*l_C_k);
        const Real l_abs_g      = const_grav;
        const Real l_inv_theta0 = 1.0 / turbChoice.theta_ref;
 
#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);
            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>& diss    = Diss.array(mfi);
 
            const Array4<Real const > &cell_data = cons_in.array(mfi);
 
            Array4<Real const> mf_u = mapfac_u.array(mfi);
            Array4<Real const> mf_v = mapfac_v.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) {
                    // the terrain grid is only deformed in z for now
                    dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
                }
                Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
                Real DeltaMsf   = std::pow(cellVolMsf,1.0/3.0);
 
                // Calculate stratification-dependent mixing length (Deardorff 1980)
                Real eps       = std::numeric_limits<Real>::epsilon();
                Real dtheta_dz;
                if (use_most && k==klo) {
                    if (exp_most) {
                        dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
                                    - cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp) )*dzInv;
                    } else {
                        dtheta_dz = 0.5 * (-3 * cell_data(i,j,k  ,RhoTheta_comp)
                                              / cell_data(i,j,k  ,Rho_comp)
                                          + 4 * cell_data(i,j,k+1,RhoTheta_comp)
                                              / cell_data(i,j,k+1,Rho_comp)
                                          -     cell_data(i,j,k+2,RhoTheta_comp)
                                              / cell_data(i,j,k+2,Rho_comp) ) * dzInv;
                    }
                } else {
                    dtheta_dz = 0.5 * ( 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;
                }
                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;
                Real length;
                if (stratification <= eps) {
                    length = DeltaMsf;
                } else {
                    length = 0.76 * std::sqrt(E / amrex::max(stratification,eps));
                    // mixing length should be _reduced_ for stable stratification
                    length = amrex::min(length, DeltaMsf);
                    // following WRF, make sure the mixing length isn't too small
                    length = amrex::max(length, 0.001 * DeltaMsf);
                }
 
                // 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 * length * std::sqrt(E);
                mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h);
                // KH = (1 + 2*l/delta) * mu_turb
                mu_turb(i,j,k,EddyDiff::Theta_h) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_h);
                mu_turb(i,j,k,EddyDiff::Theta_v) = mu_turb(i,j,k,EddyDiff::Theta_h);
                // 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 = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf;
                }
                diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length;
 
                // - heat flux
                //   (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will
                //    be overwritten when BCs are applied)
                hfx_x(i,j,k) = 0.0;
                hfx_y(i,j,k) = 0.0;
                hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
            });
        }
    }
 
    // 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);
        bxcc.growLo(0,ngc); bxcc.growHi(0,ngc);
        bxcc.growLo(1,ngc); bxcc.growHi(1,ngc);
 
        const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
 
        for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) {
            int offset = (EddyDiff::NumDiffs-1)/2;
            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;
 
                    // 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,indx);
                    }
                });
                break;
          }
       }
    }
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeTurbulentViscosityRANS()

void ComputeTurbulentViscosityRANS	(	const MultiFab &	,
		const MultiFab &	,
		const MultiFab &	,
		const MultiFab &	,
		const MultiFab &	,
		const MultiFab &	,
		const MultiFab &	cons_in,
		const MultiFab &	wdist,
		MultiFab &	eddyViscosity,
		MultiFab &	Hfx1,
		MultiFab &	Hfx2,
		MultiFab &	Hfx3,
		MultiFab &	Diss,
		const Geometry &	geom,
		bool	use_terrain,
		const MultiFab &	,
		const MultiFab &	,
		const std::unique_ptr< MultiFab > &	z_phys_nd,
		const TurbChoice &	turbChoice,
		const Real	const_grav,
		std::unique_ptr< ABLMost > &	most,
		const FArrayBox *	z_0,
		const bool &	exp_most
	)

Function for computing the eddy viscosity with RANS.

Parameters

[in]	Tau11	11 strain
[in]	Tau22	22 strain
[in]	Tau33	33 strain
[in]	Tau12	12 strain
[in]	Tau13	13 strain
[in]	Tau23	23 strain
[in]	cons_in	cell center conserved quantities
[out]	eddyViscosity	turbulent viscosity
[in]	Hfx1	heat flux in x-dir
[in]	Hfx2	heat flux in y-dir
[in]	Hfx3	heat flux in z-dir
[in]	Diss	dissipation of turbulent kinetic energy
[in]	geom	problem geometry
[in]	mapfac_u	map factor at x-face
[in]	mapfac_v	map factor at y-face
[in]	turbChoice	container with turbulence parameters

{
    const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
    const Box& domain = geom.Domain();
    const int& klo    = domain.smallEnd(2);
    const bool use_most = (most != nullptr);
 
    Real inv_Pr_t    = turbChoice.Pr_t_inv;
    Real inv_Sc_t    = turbChoice.Sc_t_inv;
    Real inv_sigma_k = 1.0 / 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  = 1.0 / (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 abs_g      = const_grav;
        const Real inv_theta0 = 1.0 / turbChoice.theta_ref;
 
#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_most) ? z_0->const_array() : Array4<Real const>{};
 
            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>& 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);
 
            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), eps);
 
                // Estimate stratification
                Real dtheta_dz;
                Real dzInv = cellSizeInv[2];
                if (use_terrain) {
                    // the terrain grid is only deformed in z for now
                    dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
                }
                if (use_most && k==klo) {
                    if (exp_most) {
                        dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
                                    - cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp) )*dzInv;
                    } else {
                        dtheta_dz = 0.5 * (-3 * cell_data(i,j,k  ,RhoTheta_comp)
                                              / cell_data(i,j,k  ,Rho_comp)
                                          + 4 * cell_data(i,j,k+1,RhoTheta_comp)
                                              / cell_data(i,j,k+1,Rho_comp)
                                          -     cell_data(i,j,k+2,RhoTheta_comp)
                                              / cell_data(i,j,k+2,Rho_comp) ) * dzInv;
                    }
                } else {
                    dtheta_dz = 0.5 * ( 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;
                }
                Real N2 = abs_g * inv_theta0 * dtheta_dz; // Brunt–Väisälä frequency squared
 
                // Geometric length scale (AL01, Eqn. 22)
                Real l_g = (z0_arr) ? KAPPA * (d_arr(i, j, k) + z0_arr(i, j, 0))
                                    : KAPPA * d_arr(i, j, k);
 
                // Enforce a maximum value
                l_g = l_g_max * l_g / (l_g_max + l_g);
 
                // Turbulent Richardson number (AL01, Eqn. 29)
                // using the old dissipation value
                Real diss0 = max(diss(i, j, k) / cell_data(i, j, k, Rho_comp),
                                 eps);
                Real Rt = tke*tke * N2 / diss0;
 
                // Turbulent length scale
                Real length;
                if (std::abs(N2) <= eps) {
                    length = l_g;
                } else if (N2 > eps) {
                    // Stable (AL01, Eqn. 26)
                    length = std::sqrt(1.0 /
                            (1.0 / (l_g * l_g) + inv_Cb_sq * N2 / tke));
                } else {
                    // Unstable (AL01, Eqn. 28)
                    length = l_g * std::sqrt(1.0 - Cmu0_pow3*Cmu0_pow3 * inv_Cb_sq * Rt);
                }
                mu_turb(i, j, k, EddyDiff::Turb_lengthscale) = length;
 
                // Burchard & Petersen smoothing function
                Rt = (Rt >= Rt_crit) ? Rt
                                     : std::max(Rt,
                                                Rt - std::pow(Rt - Rt_crit, 2) /
                                                     (Rt + Rt_min - 2*Rt_crit));
 
                // Stability functions
                // Note: These use the smoothed turbulent Richardson number
                Real cmu = (Cmu0 + 0.108*Rt)
                         / (1.0 + 0.308*Rt + 0.00837*Rt*Rt); // (AL01, Eqn. 31)
                Real cmu_prime = Cmu0 / (1 + 0.277*Rt); // (AL01, Eqn. 32)
 
                // Calculate eddy diffusivities
                // K = rho * nu_t = rho * c_mu * tke^(1/2) * length
                Real nut = cmu * std::sqrt(tke) * length; // eddy viscosity
                Real nut_prime = cmu_prime / cmu * nut;   // eddy diffusivity
                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;
 
                // Calculate SFS quantities
                // - dissipation (AL01 Eqn. 19)
                diss(i, j, k) = cell_data(i, j, k, Rho_comp) * Cmu0_pow3 * std::pow(tke,1.5) / length;
 
                // - heat flux
                //   (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will
                //    be overwritten when BCs are applied)
                hfx_x(i, j, k) = 0.0;
                hfx_y(i, j, k) = 0.0;
                // Note: buoyant production = g/theta0 * hfx == -nut_prime * N^2 (c.f. AL01 Eqn. 15)
                //                          = nut_prime * g/theta0 * dtheta/dz
                //                  ==> hfx = nut_prime * dtheta/dz
                //   Our convention is such that dtheta/dz < 0 gives a positive
                //   (upward) heat flux.
                hfx_z(i, j, k) = -mu_turb(i, j, k, EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
            });
        }
    }
 
    // 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.rans_type == RANSType::kEqn );
 
#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);
        bxcc.growLo(0,ngc); bxcc.growHi(0,ngc);
        bxcc.growLo(1,ngc); bxcc.growHi(1,ngc);
 
        const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
 
        for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) {
            int offset = (EddyDiff::NumDiffs-1)/2;
            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;
 
                    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 Deardorff
                    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: