ERF_EddyViscosity.H

Go to the documentation of this file.

/** \file ERF_EddyViscosity.H */
 
#ifndef ERF_EDDY_VISCOSITY_H_
#define ERF_EDDY_VISCOSITY_H_
 
#include "AMReX_BCRec.H"
#include "ERF_SurfaceLayer.H"
#include "ERF_DataStruct.H"
#include "ERF_IndexDefines.H"
#include "ERF_NumericalConstants.H"
#include "ERF_EOS.H"
#include "ERF_EB.H"
 
#ifdef ERF_USE_EAMXX_SHOC
#include "ERF_ShocInterface.H"
#endif
 
/** Compute turbulent viscosity and related flux fields for the configured turbulence model. */
void
ComputeTurbulentViscosity (double dt,
                           const amrex::MultiFab& xvel ,
                           const amrex::MultiFab& yvel ,
                           amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev,
                           amrex::MultiFab& cons_in,
                           const amrex::MultiFab& wdist,
                           amrex::MultiFab& eddyViscosity,
                           amrex::MultiFab& Hfx1,
                           amrex::MultiFab& Hfx2,
                           amrex::MultiFab& Hfx3,
                           amrex::MultiFab& Diss,
                           const amrex::Geometry& geom,
                           amrex::Vector<std::unique_ptr<amrex::MultiFab>>& mapfac,
                           const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                           const std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                           const SolverChoice& solverChoice,
                           std::unique_ptr<SurfaceLayer>& SurfLayer,
                           const amrex::MultiFab* z_0,
                           const bool& use_terrain_fitted_coords,
                           const bool& use_moisture,
                           int level,
                           const amrex::BCRec* bc_ptr,
                           const eb_& ebfact,
                           bool vert_only = false,
                           const amrex::MultiFab* qheating_rates = nullptr,
                           const amrex::MultiFab* terrain_blank = nullptr,
                           const amrex::MultiFab* forest_drag = nullptr);
 
/** Compute LES turbulent viscosity and dissipation.
 *
 * No subgrid heat flux is stored: Hfx3 is z-nodal and the theta diffusion of each RK stage
 * writes its face fluxes there (ERF_ComputeTurbulentViscosity.cpp, ERF_TurbKESources.H).
 */
void
ComputeTurbulentViscosityLES (amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev,
                              const amrex::MultiFab& cons_in,
                              amrex::MultiFab& eddyViscosity,
                              amrex::MultiFab& Diss,
                              const amrex::Geometry& geom,
                              bool use_terrain_fitted_coords,
                              amrex::Vector<std::unique_ptr<amrex::MultiFab>>& mapfac,
                              const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                              const TurbChoice& turbChoice,
                              const amrex::Real const_grav,
                              std::unique_ptr<SurfaceLayer>& SurfLayer,
                              const MoistureComponentIndices& moisture_indices,
                              const amrex::MultiFab* xvel = nullptr,
                              const amrex::MultiFab* yvel = nullptr);
 
/**
 * Compute the RANS eddy viscosity, length scale and dissipation (k-eqn, Axell & Liungman 2001).
 * ib_blank is the solid fraction of the fraction-stress immersed wall law (else nullptr), and
 * flat_wall_diss (erf.if_flat_wall_dissipation) makes its partial wall cells dissipate with the
 * length of flat ground's first cell. Declared here for the unit tests; ComputeTurbulentViscosity
 * is the caller in the code.
 */
void
ComputeTurbulentViscosityRANS (int level,
                               amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev,
                               const amrex::MultiFab& cons_in,
                               const amrex::MultiFab& wdist,
                               amrex::MultiFab& eddyViscosity,
                               amrex::MultiFab& Diss,
                               const amrex::Geometry& geom,
                               bool use_terrain_fitted_coords,
                               const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                               const TurbChoice& turbChoice,
                               const amrex::Real const_grav,
                               std::unique_ptr<SurfaceLayer>& SurfLayer,
                               const amrex::MultiFab* z_0,
                               const CanopyTurbulence& canopy,
                               const amrex::MultiFab* forest_drag,
                               const amrex::MultiFab* ib_blank = nullptr,
                               bool flat_wall_diss = false,
                               amrex::Real ib_small_volfrac = amrex::Real(0.005),
                               amrex::Real ib_z0 = amrex::Real(0.0));
 
/** Compute LES turbulent viscosity and heat fluxes for embedded-boundary cells. */
void
ComputeTurbulentViscosityLES_EB (amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev,
                              const amrex::MultiFab& cons_in,
                              amrex::MultiFab& eddyViscosity,
                              amrex::MultiFab& Hfx1,
                              amrex::MultiFab& Hfx2,
                              amrex::MultiFab& Hfx3,
                              const amrex::Geometry& geom,
                              amrex::Vector<std::unique_ptr<amrex::MultiFab>>& mapfac,
                              const TurbChoice& turbChoice,
                              const amrex::Real const_grav,
                              const SolverChoice& solverChoice,
                              std::unique_ptr<SurfaceLayer>& SurfLayer,
                              const MoistureComponentIndices& moisture_indices,
                              const eb_& ebfact,
                              const amrex::MultiFab* xvel = nullptr,
                              const amrex::MultiFab* yvel = nullptr);
 
/**
 * Compute twice-contracted strain-rate magnitude from stress-like strain components.
 *
 * @param[in] i cell-centered i-index
 * @param[in] j cell-centered j-index
 * @param[in] k cell-centered k-index
 * @param[in] tau11 11 strain component
 * @param[in] tau22 22 strain component
 * @param[in] tau33 33 strain component
 * @param[in] tau12 12 strain component
 * @param[in] tau13 13 strain component
 * @param[in] tau23 23 strain component
 * @return S_mn S_mn at cell centers
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
ComputeSmnSmn (int& i, int& j, int& k,
               const amrex::Array4<amrex::Real const>& tau11,
               const amrex::Array4<amrex::Real const>& tau22,
               const amrex::Array4<amrex::Real const>& tau33,
               const amrex::Array4<amrex::Real const>& tau12,
               const amrex::Array4<amrex::Real const>& tau13,
               const amrex::Array4<amrex::Real const>& tau23)
{
    amrex::Real s11bar = tau11(i,j,k);
    amrex::Real s22bar = tau22(i,j,k);
    amrex::Real s33bar = tau33(i,j,k);
    amrex::Real s12bar = fourth * ( tau12(i  , j  , k  ) + tau12(i  , j+1, k  )
                                  + tau12(i+1, j  , k  ) + tau12(i+1, j+1, k  ) );
    amrex::Real s13bar = fourth * ( tau13(i  , j  , k  ) + tau13(i  , j  , k+1)
                                  + tau13(i+1, j  , k  ) + tau13(i+1, j  , k+1) );
    amrex::Real s23bar = fourth * ( tau23(i  , j  , k  ) + tau23(i  , j  , k+1)
                                  + tau23(i  , j+1, k  ) + tau23(i  , j+1, k+1) );
 
    amrex::Real SmnSmn = s11bar*s11bar + s22bar*s22bar + s33bar*s33bar
                       + two*s12bar*s12bar + two*s13bar*s13bar + two*s23bar*s23bar;
 
    return SmnSmn;
}
 
/** The six strain components ComputeSmnSmn reads, as const arrays of one box. */
struct StrainArrays {
    amrex::Array4<amrex::Real const> tau11, tau22, tau33, tau12, tau13, tau23;
};
 
/** Strain arrays of the box of mfi from the Tau_lev container (tau11..tau23 of TauType). */
inline StrainArrays
strain_arrays (amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Tau_lev, const amrex::MFIter& mfi)
{
    return { Tau_lev[TauType::tau11]->const_array(mfi), Tau_lev[TauType::tau22]->const_array(mfi),
             Tau_lev[TauType::tau33]->const_array(mfi), Tau_lev[TauType::tau12]->const_array(mfi),
             Tau_lev[TauType::tau13]->const_array(mfi), Tau_lev[TauType::tau23]->const_array(mfi) };
}
 
/**
 * Compute twice-contracted strain-rate magnitude for embedded-boundary cells.
 *
 * @param[in] i cell-centered i-index
 * @param[in] j cell-centered j-index
 * @param[in] k cell-centered k-index
 * @param[in] tau11 11 strain component
 * @param[in] tau22 22 strain component
 * @param[in] tau33 33 strain component
 * @param[in] tau12 12 strain component
 * @param[in] tau13 13 strain component
 * @param[in] tau23 23 strain component
 * @param[in] c_cflag cell-centered EB flags
 * @param[in] u_cflag x-face EB flags
 * @param[in] v_cflag y-face EB flags
 * @param[in] w_cflag z-face EB flags
 * @return S_mn S_mn at uncovered cell centers
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
ComputeSmnSmn_EB (int& i, int& j, int& k,
               const amrex::Array4<amrex::Real const>& tau11,
               const amrex::Array4<amrex::Real const>& tau22,
               const amrex::Array4<amrex::Real const>& tau33,
               const amrex::Array4<amrex::Real const>& tau12,
               const amrex::Array4<amrex::Real const>& tau13,
               const amrex::Array4<amrex::Real const>& tau23,
               const amrex::Array4<const amrex::EBCellFlag>& c_cflag,
               const amrex::Array4<const amrex::EBCellFlag>& u_cflag,
               const amrex::Array4<const amrex::EBCellFlag>& v_cflag,
               const amrex::Array4<const amrex::EBCellFlag>& w_cflag)
{
    amrex::Real SmnSmn = zero;
    if (!c_cflag(i,j,k).isCovered()) {
        amrex::Real s11bar = tau11(i,j,k);
        amrex::Real s22bar = tau22(i,j,k);
        amrex::Real s33bar = tau33(i,j,k);
        amrex::Real s12bar = zero;
        amrex::Real s13bar = zero;
        amrex::Real s23bar = zero;
 
        amrex::Real count_s12 = zero;
        if (!u_cflag(i,j,k).isCovered() || !v_cflag(i,j,k).isCovered()) {
            s12bar += tau12(i,j,k);
            count_s12 += one;
        }
        if (!u_cflag(i+1,j,k).isCovered() || !v_cflag(i,j,k).isCovered()) {
            s12bar += tau12(i+1,j,k);
            count_s12 += one;
        }
        if (!u_cflag(i,j,k).isCovered() || !v_cflag(i,j+1,k).isCovered()) {
            s12bar += tau12(i,j+1,k);
            count_s12 += one;
        }
        if (!u_cflag(i+1,j,k).isCovered() || !v_cflag(i,j+1,k).isCovered()) {
            s12bar += tau12(i+1,j+1,k);
            count_s12 += one;
        }
        if (count_s12 > zero) {
            s12bar /= count_s12;
        }
 
        amrex::Real count_s13 = zero;
        if (!u_cflag(i,j,k).isCovered() || !w_cflag(i,j,k).isCovered()) {
            s13bar += tau13(i, j, k);
            count_s13 += one;
        }
        if (!u_cflag(i,j,k).isCovered() || !w_cflag(i,j,k+1).isCovered()) {
            s13bar += tau13(i, j, k+1);
            count_s13 += one;
        }
        if (!u_cflag(i+1,j,k).isCovered() || !w_cflag(i,j,k).isCovered()) {
            s13bar += tau13(i+1, j, k);
            count_s13 += one;
        }
        if (!u_cflag(i+1,j,k).isCovered() || !w_cflag(i,j,k+1).isCovered()) {
            s13bar += tau13(i+1, j, k+1);
            count_s13 += one;
        }
        if (count_s13 > zero) {
            s13bar /= count_s13;
        }
 
        amrex::Real count_s23 = zero;
        if (!v_cflag(i,j,k).isCovered() || !w_cflag(i,j,k).isCovered()) {
            s23bar += tau23(i, j, k);
            count_s23 += one;
        }
        if (!v_cflag(i,j,k).isCovered() || !w_cflag(i,j,k+1).isCovered()) {
            s23bar += tau23(i, j, k+1);
            count_s23 += one;
        }
        if (!v_cflag(i,j+1,k).isCovered() || !w_cflag(i,j,k).isCovered()) {
            s23bar += tau23(i, j+1, k);
            count_s23 += one;
        }
        if (!v_cflag(i,j+1,k).isCovered() || !w_cflag(i,j,k+1).isCovered()) {
            s23bar += tau23(i, j+1, k+1);
            count_s23 += one;
        }
        if (count_s23 > zero) {
            s23bar /= count_s23;
        }
 
        SmnSmn = s11bar*s11bar + s22bar*s22bar + s33bar*s33bar
               + two*s12bar*s12bar + two*s13bar*s13bar + two*s23bar*s23bar;
    }
 
    return SmnSmn;
}
 
/**
 * Compute horizontal two-dimensional strain-rate magnitude.
 *
 * @param[in] i cell-centered i-index
 * @param[in] j cell-centered j-index
 * @param[in] k cell-centered k-index
 * @param[in] tau11 11 strain component
 * @param[in] tau22 22 strain component
 * @param[in] tau12 12 strain component
 * @return horizontal S_mn S_mn contribution
 */
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
ComputeSmnSmn2D (int& i, int& j, int& k,
                 const amrex::Array4<amrex::Real const>& tau11,
                 const amrex::Array4<amrex::Real const>& tau22,
                 const amrex::Array4<amrex::Real const>& tau12)
{
    amrex::Real sdiff = tau11(i,j,k) - tau22(i,j,k);
    amrex::Real s12bar = fourth * ( tau12(i  , j  , k  ) + tau12(i  , j+1, k  )
                                  + tau12(i+1, j  , k  ) + tau12(i+1, j+1, k  ) );
    return myhalf * sdiff*sdiff + two*s12bar*s12bar;
}
#endif