ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_ComputeTurbulentViscosity.cpp File Reference
#include "ERF_SurfaceLayer.H"
#include "ERF_EddyViscosity.H"
#include "ERF_Diffusion.H"
#include "ERF_PBLModels.H"
#include "ERF_TileNoZ.H"
#include "ERF_TerrainMetrics.H"
#include "ERF_MoistUtils.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 &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, bool use_terrain, 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)
 
void ComputeTurbulentViscosityRANS (Vector< std::unique_ptr< MultiFab >> &, const MultiFab &cons_in, const MultiFab &wdist, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, bool use_terrain, Vector< std::unique_ptr< MultiFab >> &, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< SurfaceLayer > &SurfLayer, const MultiFab *z_0)
 
void ComputeTurbulentViscosity (Real 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 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, bool vert_only)
 

Function Documentation

◆ ComputeTurbulentViscosity()

void ComputeTurbulentViscosity ( Real  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 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,
bool  vert_only 
)

Wrapper to compute turbulent viscosity with LES or PBL.

Parameters
[in]xvelvelocity in x-dir
[in]yvelvelocity in y-dir
[in]Tau_levstrain at this level
[in]cons_incell center conserved quantities
[out]eddyViscosityturbulent viscosity
[in]Hfx1heat flux in x-dir
[in]Hfx2heat flux in y-dir
[in]Hfx3heat flux in z-dir
[in]Dissdissipation of turbulent kinetic energy
[in]geomproblem geometry
[in]mapfacmap factors
[in]turbChoicecontainer with turbulence parameters
[in]mostpointer to Monin-Obukhov class if instantiated
[in]vert_onlyflag for vertical components of eddyViscosity
609 {
610  BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity);
611  //
612  // In LES mode, the turbulent viscosity is isotropic (unless mix_isotropic is set to false), so
613  // the LES model sets both horizontal and vertical viscosities
614  //
615  // In PBL mode, the primary purpose of the PBL model is to control vertical transport, so the PBL model sets the vertical viscosity.
616  // Optionally, the PBL model can be run in conjunction with an LES model that sets the horizontal viscosity
617  // (this isn’t truly LES, but the model form is the same as Smagorinsky).
618  //
619  // ComputeTurbulentViscosityLES populates the LES viscosity for both horizontal and vertical components.
620  // ComputeTurbulentViscosityPBL computes the PBL viscosity just for the vertical component.
621  //
622 
623  TurbChoice turbChoice = solverChoice.turbChoice[level];
624  const Real const_grav = solverChoice.gravity;
625 
626  if (!SurfLayer) {
627  AMREX_ALWAYS_ASSERT(!vert_only);
628  }
629 
630  bool impose_phys_bcs = false;
631 
632  if (turbChoice.les_type != LESType::None) {
633  impose_phys_bcs = true;
635  cons_in, eddyViscosity,
636  Hfx1, Hfx2, Hfx3, Diss,
637  geom, use_terrain_fitted_coords,
638  mapfac, z_phys_nd, turbChoice, const_grav,
639  SurfLayer, solverChoice.moisture_indices);
640  }
641 
642  if (turbChoice.rans_type != RANSType::None) {
643  impose_phys_bcs = true;
645  cons_in, wdist,
646  eddyViscosity,
647  Hfx1, Hfx2, Hfx3, Diss,
648  geom, use_terrain_fitted_coords,
649  mapfac, z_phys_nd, turbChoice, const_grav,
650  SurfLayer, z_0);
651  }
652 
653  if (turbChoice.pbl_type == PBLType::MYJ) {
654  ComputeDiffusivityMYJ(dt, xvel, yvel, cons_in, eddyViscosity,
655  geom, turbChoice, SurfLayer,
656  use_terrain_fitted_coords, use_moisture,
657  level, bc_ptr, vert_only, z_phys_nd,
658  solverChoice.moisture_indices);
659  } else if (turbChoice.pbl_type == PBLType::MYNN25) {
660  ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity,
661  geom, turbChoice, SurfLayer,
662  use_terrain_fitted_coords, use_moisture,
663  level, bc_ptr, vert_only, z_phys_nd,
664  solverChoice.moisture_indices);
665  } else if (turbChoice.pbl_type == PBLType::MYNNEDMF) {
666  ComputeDiffusivityMYNNEDMF(xvel, yvel, cons_in, eddyViscosity,
667  geom, turbChoice, SurfLayer,
668  use_terrain_fitted_coords, use_moisture,
669  level, bc_ptr, vert_only, z_phys_nd,
670  solverChoice.moisture_indices);
671  } else if (turbChoice.pbl_type == PBLType::YSU) {
672  ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity,
673  geom, turbChoice, SurfLayer,
674  use_terrain_fitted_coords, use_moisture,
675  level, bc_ptr, vert_only, z_phys_nd,
676  solverChoice.moisture_indices);
677  } else if (turbChoice.pbl_type == PBLType::MRF) {
678  ComputeDiffusivityMRF(xvel, yvel, cons_in, eddyViscosity,
679  geom, turbChoice, SurfLayer,
680  use_terrain_fitted_coords, use_moisture,
681  level, bc_ptr, vert_only, z_phys_nd,
682  solverChoice.moisture_indices);
683  }
684 
685  //
686  // At all levels we need to fill values outside the physical boundary for the LES coeffs.
687  // In addition, for all cases, if at level > 0, we want to fill fine ghost cell values that
688  // overlie coarse grid cells (and that are not in another fine valid region) with
689  // extrapolated values from the interior, rather than interpolating from the coarser level,
690  // since we may be using a different turbulence model there.
691  //
692  // Note: here "covered" refers to "covered by valid region of another grid at this level"
693  // Note: here "physbnd" refers to "cells outside the domain if not periodic"
694  // Note: here "interior" refers to "valid cells, i.e. inside 'my' grid"
695  //
696  {
697  int is_covered = 0;
698  int is_notcovered = 1;
699  int is_physbnd = 2;
700  int is_interior = 3;
701  iMultiFab cc_mask(eddyViscosity.boxArray(),eddyViscosity.DistributionMap(),1,1);
702  cc_mask.BuildMask(geom.Domain(), geom.periodicity(), is_covered, is_notcovered, is_physbnd, is_interior);
703 
704  Box domain = geom.Domain();
705  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
706  if (geom.isPeriodic(i)) {
707  domain.grow(i,1);
708  }
709  }
710 
711  eddyViscosity.FillBoundary(geom.periodicity());
712 
713  int ncomp = eddyViscosity.nComp();
714 
715 #ifdef _OPENMP
716 #pragma omp parallel if (Gpu::notInLaunchRegion())
717 #endif
718  for (MFIter mfi(eddyViscosity); mfi.isValid(); ++mfi)
719  {
720  Box vbx = mfi.validbox();
721 
722  Box planex_lo = mfi.growntilebox(1); planex_lo.setBig(0, vbx.smallEnd(0)-1);
723  Box planey_lo = mfi.growntilebox(1); planey_lo.setBig(1, vbx.smallEnd(1)-1);
724  Box planez_lo = mfi.growntilebox(1); planez_lo.setBig(2, vbx.smallEnd(2)-1);
725 
726  Box planex_hi = mfi.growntilebox(1); planex_hi.setSmall(0, vbx.bigEnd(0)+1);
727  Box planey_hi = mfi.growntilebox(1); planey_hi.setSmall(1, vbx.bigEnd(1)+1);
728  Box planez_hi = mfi.growntilebox(1); planez_hi.setSmall(2, vbx.bigEnd(2)+1);
729 
730  int i_lo = vbx.smallEnd(0); int i_hi = vbx.bigEnd(0);
731  int j_lo = vbx.smallEnd(1); int j_hi = vbx.bigEnd(1);
732  int k_lo = vbx.smallEnd(2); int k_hi = vbx.bigEnd(2);
733 
734  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
735  const Array4<int>& mask_arr = cc_mask.array(mfi);
736 
737  auto domlo = lbound(domain);
738  auto domhi = ubound(domain);
739 
740  ParallelFor(planex_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
741  {
742  int lj = amrex::min(amrex::max(j, domlo.y), domhi.y);
743  int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
744  if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i_lo,lj,lk) != is_notcovered) ||
745  (mask_arr(i,j,k) == is_physbnd && i < domlo.x && impose_phys_bcs)) {
746  for (int n = 0; n < ncomp; n++) {
747  mu_turb(i,j,k,n) = mu_turb(i_lo,lj,lk,n);
748  }
749  }
750  });
751  ParallelFor(planex_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
752  {
753  int lj = amrex::min(amrex::max(j, domlo.y), domhi.y);
754  int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
755  if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i_hi,lj,lk) != is_notcovered) ||
756  (mask_arr(i,j,k) == is_physbnd && i > domhi.x && impose_phys_bcs)) {
757  for (int n = 0; n < ncomp; n++) {
758  mu_turb(i,j,k,n) = mu_turb(i_hi,lj,lk,n);
759  }
760  }
761  });
762  ParallelFor(planey_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
763  {
764  int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
765  if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j_lo,lk) != is_notcovered) ||
766  (mask_arr(i,j,k) == is_physbnd && j < domlo.y && impose_phys_bcs)) {
767  for (int n = 0; n < ncomp; n++) {
768  mu_turb(i,j,k,n) = mu_turb(i,j_lo,lk,n);
769  }
770  }
771  });
772  ParallelFor(planey_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
773  {
774  int lk = amrex::min(amrex::max(k, domlo.z), domhi.z);
775  if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j_hi,lk) != is_notcovered)||
776  (mask_arr(i,j,k) == is_physbnd && j > domhi.y && impose_phys_bcs)) {
777  for (int n = 0; n < ncomp; n++) {
778  mu_turb(i,j,k,n) = mu_turb(i,j_hi,lk,n);
779  }
780  }
781  });
782  ParallelFor(planez_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
783  {
784  if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j,k_lo) != is_notcovered) ||
785  (mask_arr(i,j,k) == is_physbnd && k < domlo.z && impose_phys_bcs)) {
786  if (mask_arr(i,j,k) == is_physbnd) {
787  }
788  for (int n = 0; n < ncomp; n++) {
789  mu_turb(i,j,k,n) = mu_turb(i,j,k_lo,n);
790  }
791  }
792  });
793  ParallelFor(planez_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
794  {
795  if ((mask_arr(i,j,k) == is_notcovered && mask_arr(i,j,k_hi) != is_notcovered) ||
796  (mask_arr(i,j,k) == is_physbnd && k > domhi.z && impose_phys_bcs)) {
797  for (int n = 0; n < ncomp; n++) {
798  mu_turb(i,j,k,n) = mu_turb(i,j,k_hi,n);
799  }
800  }
801  });
802  } // mfi
803  eddyViscosity.FillBoundary(geom.periodicity());
804  }
805 }
void ComputeDiffusivityMRF(const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
Definition: ERF_ComputeDiffusivityMRF.cpp:12
void ComputeDiffusivityMYJ(Real dt, const MultiFab &xvel, const MultiFab &yvel, MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &, std::unique_ptr< SurfaceLayer > &, bool use_terrain_fitted_coords, bool, int, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
Definition: ERF_ComputeDiffusivityMYJ.cpp:13
void ComputeDiffusivityMYNN25(const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
Definition: ERF_ComputeDiffusivityMYNN25.cpp:13
void ComputeDiffusivityMYNNEDMF(const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
Definition: ERF_ComputeDiffusivityMYNNEDMF.cpp:4177
void ComputeDiffusivityYSU(const MultiFab &xvel, const MultiFab &yvel, const MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool, int level, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
Definition: ERF_ComputeDiffusivityYSU.cpp:11
void ComputeTurbulentViscosityLES(Vector< std::unique_ptr< MultiFab >> &Tau_lev, const MultiFab &cons_in, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, bool use_terrain, 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)
Definition: ERF_ComputeTurbulentViscosity.cpp:27
void ComputeTurbulentViscosity(Real 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 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, bool vert_only)
Definition: ERF_ComputeTurbulentViscosity.cpp:591
void ComputeTurbulentViscosityRANS(Vector< std::unique_ptr< MultiFab >> &, const MultiFab &cons_in, const MultiFab &wdist, MultiFab &eddyViscosity, MultiFab &Hfx1, MultiFab &Hfx2, MultiFab &Hfx3, MultiFab &Diss, const Geometry &geom, bool use_terrain, Vector< std::unique_ptr< MultiFab >> &, const std::unique_ptr< MultiFab > &z_phys_nd, const TurbChoice &turbChoice, const Real const_grav, std::unique_ptr< SurfaceLayer > &SurfLayer, const MultiFab *z_0)
Definition: ERF_ComputeTurbulentViscosity.cpp:370
TurbChoice turbChoice
Definition: ERF_DiffSetup.H:2
amrex::Real Real
Definition: ERF_ShocInterface.H:19
@ xvel
Definition: ERF_IndexDefines.H:141
@ yvel
Definition: ERF_IndexDefines.H:142
amrex::Real gravity
Definition: ERF_DataStruct.H:864
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:826
MoistureComponentIndices moisture_indices
Definition: ERF_DataStruct.H:930
Definition: ERF_TurbStruct.H:41
PBLType pbl_type
Definition: ERF_TurbStruct.H:375
RANSType rans_type
Definition: ERF_TurbStruct.H:372
LESType les_type
Definition: ERF_TurbStruct.H:330

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 &  Hfx1,
MultiFab &  Hfx2,
MultiFab &  Hfx3,
MultiFab &  Diss,
const Geometry &  geom,
bool  use_terrain,
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 
)

Function for computing the turbulent viscosity with LES.

Parameters
[in]Tau_levstrain at this level
[in]cons_incell center conserved quantities
[out]eddyViscosityturbulent viscosity
[in]Hfx1heat flux in x-dir
[in]Hfx2heat flux in y-dir
[in]Hfx3heat flux in z-dir
[in]Dissdissipation of turbulent kinetic energy
[in]geomproblem geometry
[in]mapfacmap factors
[in]turbChoicecontainer with turbulence parameters
36 {
37  const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
38  const Box& domain = geom.Domain();
39 
40  Real inv_Pr_t = turbChoice.Pr_t_inv;
41  Real inv_Sc_t = turbChoice.Sc_t_inv;
42  Real inv_sigma_k = 1.0 / turbChoice.sigma_k;
43 
44  bool use_thetav_grad = (turbChoice.strat_type == StratType::thetav);
45  bool use_thetal_grad = (turbChoice.strat_type == StratType::thetal);
46 
47  bool isotropic = turbChoice.mix_isotropic;
48 
49  // SMAGORINSKY: Fill Kturb for momentum in horizontal and vertical
50  //***********************************************************************************
51  if (turbChoice.les_type == LESType::Smagorinsky)
52  {
53  Real Cs = turbChoice.Cs;
54  bool smag2d = turbChoice.smag2d;
55 
56  // Define variables that need to be captured in the lambda
57  Real l_abs_g = const_grav;
58  const bool use_ref_theta = (turbChoice.theta_ref > 0);
59  bool l_use_moisture = (moisture_indices.qv > 0);
60  int l_rho_qv_comp = moisture_indices.qv;
61  bool l_use_smag_stratification = turbChoice.use_smag_stratification;
62 
63  #ifdef _OPENMP
64  #pragma omp parallel if (Gpu::notInLaunchRegion())
65  #endif
66  for (MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
67  {
68  // NOTE: This gets us the lateral ghost cells for lev>0; which
69  // have been filled from FP Two Levels.
70  Box bxcc = mfi.growntilebox(1) & domain;
71  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
72  const Array4<Real>& hfx_x = Hfx1.array(mfi);
73  const Array4<Real>& hfx_y = Hfx2.array(mfi);
74  const Array4<Real>& hfx_z = Hfx3.array(mfi);
75  const Array4<Real const > &cell_data = cons_in.array(mfi);
76  Array4<Real const> tau11 = Tau_lev[TauType::tau11]->array(mfi);
77  Array4<Real const> tau22 = Tau_lev[TauType::tau22]->array(mfi);
78  Array4<Real const> tau33 = Tau_lev[TauType::tau33]->array(mfi);
79  Array4<Real const> tau12 = Tau_lev[TauType::tau12]->array(mfi);
80  Array4<Real const> tau13 = Tau_lev[TauType::tau13]->array(mfi);
81  Array4<Real const> tau23 = Tau_lev[TauType::tau23]->array(mfi);
82  Array4<Real const> mf_u = mapfac[MapFacType::u_x]->const_array(mfi);
83  Array4<Real const> mf_v = mapfac[MapFacType::v_y]->const_array(mfi);
84  Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
85 
86  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
87  {
88  // =====================================================================
89  // STRAIN RATE CALCULATION
90  // =====================================================================
91  Real SmnSmn;
92  if (smag2d) {
93  SmnSmn = ComputeSmnSmn2D(i,j,k,tau11,tau22,tau12);
94  } else {
95  SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23);
96  }
97 
98  // =====================================================================
99  // GRID SCALE CALCULATION
100  // =====================================================================
101  Real dxInv = cellSizeInv[0];
102  Real dyInv = cellSizeInv[1];
103  Real dzInv = cellSizeInv[2];
104  if (use_terrain) {
105  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
106  }
107 
108  Real Delta;
109  Real DeltaH;
110  if (isotropic) {
111  Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
112  Delta = std::cbrt(cellVolMsf);
113  DeltaH = Delta;
114  } else { // separate horizontal/vertical length scales
115  // Typically, dx and dy >> dz
116  // If using enhanced Smagorinsky, dz is the more meaningful grid scale
117  Delta = 1.0 / dzInv;
118  DeltaH = std::sqrt(1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0)));
119  }
120 
121  // =====================================================================
122  // MIXING LENGTH CALCULATION: STANDARD vs ADJUSTED
123  // =====================================================================
124  Real mixing_length;
125 
126  if (l_use_smag_stratification) {
127  // ENHANCED SMAGORINSKY: Our moist stratification approach
128  Real inv_theta = (use_ref_theta) ? 1.0 / turbChoice.theta_ref
129  : cell_data(i,j,k,Rho_comp) /
130  cell_data(i,j,k,RhoTheta_comp);
131  Real stratification = ComputeStratificationForSmagorinsky(
132  i, j, k, cell_data, dzInv, l_abs_g, inv_theta,
133  l_use_moisture, l_rho_qv_comp, moisture_indices);
134 
135  // Stability-limited mixing length
136  mixing_length = Delta;
137  if (stratification > 1e-10) {
138  Real strain_rate_magnitude = std::sqrt(2.0 * SmnSmn);
139  if (strain_rate_magnitude > 1e-10) {
140  Real velocity_scale = strain_rate_magnitude * Delta;
141  Real buoyancy_length = std::sqrt(velocity_scale * velocity_scale / stratification);
142 
143  mixing_length = amrex::min(Delta, buoyancy_length);
144  mixing_length = amrex::max(mixing_length, 0.001 * Delta);
145  }
146  }
147 
148  } else {
149  // STANDARD SMAGORINSKY: Grid-scale mixing length only
150  mixing_length = Delta; // Always use grid scale
151  }
152 
153  // =====================================================================
154  // EDDY VISCOSITY CALCULATION (same for standard/enhanced)
155  // =====================================================================
156  Real CsDeltaSqr = Cs * Cs * mixing_length * mixing_length;
157 
158  if (isotropic) {
159  mu_turb(i, j, k, EddyDiff::Mom_h) = CsDeltaSqr * cell_data(i, j, k, Rho_comp) * std::sqrt(2.0*SmnSmn);
160  mu_turb(i, j, k, EddyDiff::Mom_v) = mu_turb(i, j, k, EddyDiff::Mom_h);
161 
162  } else { // anisotropic or smag2d
163  Real CsDeltaHSqr = Cs * Cs * DeltaH * DeltaH;
164  mu_turb(i, j, k, EddyDiff::Mom_h) = CsDeltaHSqr * cell_data(i, j, k, Rho_comp) * std::sqrt(2.0*SmnSmn);
165  mu_turb(i, j, k, EddyDiff::Mom_v) = CsDeltaSqr * cell_data(i, j, k, Rho_comp) * std::sqrt(2.0*SmnSmn);
166  }
167 
168  // =====================================================================
169  // HEAT FLUX CALCULATION
170  // =====================================================================
171  Real dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
172  - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv;
173  hfx_x(i,j,k) = 0.0;
174  hfx_y(i,j,k) = 0.0;
175  hfx_z(i,j,k) = -inv_Pr_t*mu_turb(i,j,k,EddyDiff::Mom_v) * dtheta_dz;
176  });
177  }
178  }
179  // DEARDORFF: Fill Kturb for momentum in horizontal and vertical
180  //***********************************************************************************
181  else if (turbChoice.les_type == LESType::Deardorff)
182  {
183  const Real l_C_k = turbChoice.Ck;
184  const Real l_C_e = turbChoice.Ce;
185  const Real l_C_e_wall = turbChoice.Ce_wall;
186  const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k);
187  const Real l_abs_g = const_grav;
188 
189  const bool use_ref_theta = (turbChoice.theta_ref > 0);
190  const Real l_inv_theta0 = (use_ref_theta) ? 1.0 / turbChoice.theta_ref : 1.0;
191 
192 #ifdef _OPENMP
193 #pragma omp parallel if (Gpu::notInLaunchRegion())
194 #endif
195  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
196  {
197  Box bxcc = mfi.tilebox();
198 
199  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
200  const Array4<Real>& hfx_x = Hfx1.array(mfi);
201  const Array4<Real>& hfx_y = Hfx2.array(mfi);
202  const Array4<Real>& hfx_z = Hfx3.array(mfi);
203  const Array4<Real>& diss = Diss.array(mfi);
204 
205  const Array4<Real const > &cell_data = cons_in.array(mfi);
206 
207  Array4<Real const> mf_u = mapfac[MapFacType::u_x]->const_array(mfi);
208  Array4<Real const> mf_v = mapfac[MapFacType::v_y]->const_array(mfi);
209 
210  Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
211 
212  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
213  {
214  Real dxInv = cellSizeInv[0];
215  Real dyInv = cellSizeInv[1];
216  Real dzInv = cellSizeInv[2];
217  if (use_terrain) {
218  // the terrain grid is only deformed in z for now
219  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
220  }
221  Real Delta;
222  if (isotropic) {
223  Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
224  Delta = std::cbrt(cellVolMsf);
225  } else {
226  Delta = 1.0 / dzInv;
227  }
228 
229  Real dtheta_dz;
230  if (use_thetav_grad) {
231  dtheta_dz = 0.5 * ( GetThetav(i, j, k+1, cell_data, moisture_indices)
232  - GetThetav(i, j, k-1, cell_data, moisture_indices) )*dzInv;
233  } else if (use_thetal_grad) {
234  dtheta_dz = 0.5 * ( GetThetal(i, j, k+1, cell_data, moisture_indices)
235  - GetThetal(i, j, k-1, cell_data, moisture_indices) )*dzInv;
236  } else {
237  dtheta_dz = 0.5 * ( cell_data(i, j, k+1, RhoTheta_comp) / cell_data(i, j, k+1, Rho_comp)
238  - cell_data(i, j, k-1, RhoTheta_comp) / cell_data(i, j, k-1, Rho_comp) )*dzInv;
239  }
240 
241  // Calculate stratification-dependent mixing length (Deardorff 1980, Eqn. 10a)
242  Real E = amrex::max(cell_data(i,j,k,RhoKE_comp)/cell_data(i,j,k,Rho_comp),Real(0.0));
243  Real stratification = l_abs_g * dtheta_dz * l_inv_theta0;
244  if (!use_ref_theta) {
245  // l_inv_theta0 == 1, divide by actual theta
246  stratification *= cell_data(i,j,k,Rho_comp) /
247  cell_data(i,j,k,RhoTheta_comp);
248  }
249 
250  // Following WRF, the stratification effects are applied to the vertical length scales
251  // in the case of anistropic mixing
252  Real length;
254  if (stratification <= eps) {
255  length = Delta; // cbrt(dx*dy*dz) -or- dz
256  } else {
257  length = 0.76 * std::sqrt(E / amrex::max(stratification,eps));
258  // mixing length should be _reduced_ for stable stratification
259  length = amrex::min(length, Delta);
260  // following WRF, make sure the mixing length isn't too small
261  length = amrex::max(length, 0.001 * Delta);
262  }
263 
264  Real DeltaH = (isotropic) ? length : std::sqrt(1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0)));
265 
266  Real Pr_inv_v = (1. + 2.*length/Delta);
267  Real Pr_inv_h = (isotropic) ? Pr_inv_v : inv_Pr_t;
268 
269  // Calculate eddy diffusivities
270  // K = rho * C_k * l * KE^(1/2)
271  mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * DeltaH * std::sqrt(E);
272  mu_turb(i,j,k,EddyDiff::Mom_v) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E);
273  // KH = (1 + 2*l/delta) * mu_turb
274  mu_turb(i,j,k,EddyDiff::Theta_h) = Pr_inv_h * mu_turb(i,j,k,EddyDiff::Mom_h);
275  mu_turb(i,j,k,EddyDiff::Theta_v) = Pr_inv_v * mu_turb(i,j,k,EddyDiff::Mom_v);
276  // Store lengthscale for TKE source terms
277  mu_turb(i,j,k,EddyDiff::Turb_lengthscale) = length;
278 
279  // Calculate SFS quantities
280  // - dissipation
281  Real Ce;
282  if ((l_C_e_wall > 0) && (k==0)) {
283  Ce = l_C_e_wall;
284  } else {
285  Ce = 1.9*l_C_k + Ce_lcoeff*length / Delta;
286  }
287  diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length;
288 
289  // - heat flux
290  // (Note: If using SurfaceLayer, the value at k=0 will
291  // be overwritten)
292  hfx_x(i,j,k) = 0.0;
293  hfx_y(i,j,k) = 0.0;
294  hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
295  });
296  }
297  }
298 
299  // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0)
300  //***********************************************************************************
301  int ngc(1);
302  // EddyDiff mapping : Theta_h KE_h Scalar_h Q_h
303  Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
304  Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
305  Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
306  Real* fac_ptr = d_Factors.data();
307 
308  const bool use_KE = ( turbChoice.les_type == LESType::Deardorff );
309 
310 #ifdef _OPENMP
311 #pragma omp parallel if (Gpu::notInLaunchRegion())
312 #endif
313  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
314  {
315  Box bxcc = mfi.tilebox();
316  Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
317  Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1);
318  bxcc.growLo(0,ngc); bxcc.growHi(0,ngc);
319  bxcc.growLo(1,ngc); bxcc.growHi(1,ngc);
320 
321  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
322 
323  for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) {
324  int offset = (EddyDiff::NumDiffs-1)/2;
325  switch (n)
326  {
327  case EddyDiff::KE_h:
328  if (use_KE) {
329  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
330  {
331  int indx = n;
332  int indx_v = indx + offset;
333  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
334  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
335  });
336  }
337  break;
338  default:
339  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
340  {
341  int indx = n;
342  int indx_v = indx + offset;
343 
344  // NOTE: Theta_h, Theta_v have already been set for Deardorff
345  if (!(indx_v == EddyDiff::Theta_v && use_KE)) {
346  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
347  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
348  }
349  });
350  break;
351  }
352  }
353  }
354 }
@ tau12
Definition: ERF_DataStruct.H:30
@ tau23
Definition: ERF_DataStruct.H:30
@ tau33
Definition: ERF_DataStruct.H:30
@ tau22
Definition: ERF_DataStruct.H:30
@ tau11
Definition: ERF_DataStruct.H:30
@ tau13
Definition: ERF_DataStruct.H:30
@ v_y
Definition: ERF_DataStruct.H:23
@ u_x
Definition: ERF_DataStruct.H:22
if(l_use_mynn &&start_comp<=RhoKE_comp &&end_comp >=RhoKE_comp)
Definition: ERF_DiffQKEAdjustment.H:2
Real l_inv_theta0
Definition: ERF_DiffSetup.H:12
Real l_abs_g
Definition: ERF_DiffSetup.H:13
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)
Definition: ERF_EddyViscosity.H:67
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)
Definition: ERF_EddyViscosity.H:40
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetal(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:101
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeStratificationForSmagorinsky(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, amrex::Real dzInv, amrex::Real abs_g, amrex::Real inv_theta0, bool use_moisture, int rho_qv_comp, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:176
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:72
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:47
@ Theta_v
Definition: ERF_IndexDefines.H:176
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:180
@ Mom_h
Definition: ERF_IndexDefines.H:170
@ Mom_v
Definition: ERF_IndexDefines.H:175
@ NumDiffs
Definition: ERF_IndexDefines.H:181
@ Theta_h
Definition: ERF_IndexDefines.H:171
@ KE_h
Definition: ERF_IndexDefines.H:172
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
int qv
Definition: ERF_DataStruct.H:100
bool smag2d
Definition: ERF_TurbStruct.H:342
amrex::Real sigma_k
Definition: ERF_TurbStruct.H:358
amrex::Real Sc_t_inv
Definition: ERF_TurbStruct.H:338
StratType strat_type
Definition: ERF_TurbStruct.H:363
amrex::Real Ck
Definition: ERF_TurbStruct.H:347
amrex::Real Cs
Definition: ERF_TurbStruct.H:341
amrex::Real Pr_t_inv
Definition: ERF_TurbStruct.H:334
amrex::Real Ce_wall
Definition: ERF_TurbStruct.H:346
amrex::Real Ce
Definition: ERF_TurbStruct.H:345
bool use_smag_stratification
Definition: ERF_TurbStruct.H:369
amrex::Real theta_ref
Definition: ERF_TurbStruct.H:361
bool mix_isotropic
Definition: ERF_TurbStruct.H:366

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeTurbulentViscosityRANS()

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

Function for computing the eddy viscosity with RANS.

Parameters
[in]Tau_levstrain at this level
[in]cons_incell center conserved quantities
[out]eddyViscosityturbulent viscosity
[in]Hfx1heat flux in x-dir
[in]Hfx2heat flux in y-dir
[in]Hfx3heat flux in z-dir
[in]Dissdissipation of turbulent kinetic energy
[in]geomproblem geometry
[in]mapfacmap factor
[in]turbChoicecontainer with turbulence parameters
386 {
387  const GpuArray<Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
388  const bool use_SurfLayer = (SurfLayer != nullptr);
389 
390  Real inv_Pr_t = turbChoice.Pr_t_inv;
391  Real inv_Sc_t = turbChoice.Sc_t_inv;
392  Real inv_sigma_k = 1.0 / turbChoice.sigma_k;
393 
394  // One-Equation k model (Axell & Liungman 2001, Environ Fluid Mech)
395  //***********************************************************************************
396  if (turbChoice.rans_type == RANSType::kEqn)
397  {
398  const Real Cmu0 = turbChoice.Cmu0;
399  const Real Cmu0_pow3 = Cmu0 * Cmu0 * Cmu0;
400  const Real inv_Cb_sq = 1.0 / (turbChoice.Cb * turbChoice.Cb);
401  const Real Rt_crit = turbChoice.Rt_crit;
402  const Real Rt_min = turbChoice.Rt_min;
403  const Real l_g_max = turbChoice.l_g_max;
404  const Real abs_g = const_grav;
405 
406  const bool use_ref_theta = (turbChoice.theta_ref > 0);
407  const Real inv_theta0 = (use_ref_theta) ? 1.0 / turbChoice.theta_ref : 1.0;
408 
409 #ifdef _OPENMP
410 #pragma omp parallel if (Gpu::notInLaunchRegion())
411 #endif
412  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
413  {
414  Box bxcc = mfi.tilebox();
415 
416  const Array4<Real const>& d_arr = wdist.const_array(mfi);
417  const Array4<Real const>& z0_arr = (use_SurfLayer) ? z_0->const_array(mfi) : Array4<Real const>{};
418 
419  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
420  const Array4<Real>& hfx_x = Hfx1.array(mfi);
421  const Array4<Real>& hfx_y = Hfx2.array(mfi);
422  const Array4<Real>& hfx_z = Hfx3.array(mfi);
423  const Array4<Real>& diss = Diss.array(mfi);
424 
425  const Array4<Real const>& cell_data = cons_in.array(mfi);
426 
427  const Array4<Real const>& z_nd_arr = z_phys_nd->const_array(mfi);
428 
429  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
430  {
432  Real tke = amrex::max(cell_data(i,j,k,RhoKE_comp)/cell_data(i,j,k,Rho_comp), eps);
433 
434  // Estimate stratification
435  Real dzInv = cellSizeInv[2];
436  if (use_terrain) {
437  // the terrain grid is only deformed in z for now
438  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
439  }
440  Real dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
441  - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv;
442  Real N2 = abs_g * inv_theta0 * dtheta_dz; // Brunt–Väisälä frequency squared
443  if (!use_ref_theta) {
444  // inv_theta0 == 1, divide by actual theta
445  N2 *= cell_data(i,j,k,Rho_comp) /
446  cell_data(i,j,k,RhoTheta_comp);
447  }
448 
449  // Geometric length scale (AL01, Eqn. 22)
450  Real l_g = (z0_arr) ? KAPPA * (d_arr(i, j, k) + z0_arr(i, j, 0))
451  : KAPPA * d_arr(i, j, k);
452 
453  // Enforce a maximum value
454  l_g = l_g_max * l_g / (l_g_max + l_g);
455 
456  // Turbulent Richardson number (AL01, Eqn. 29)
457  // using the old dissipation value
458  Real diss0 = max(diss(i, j, k) / cell_data(i, j, k, Rho_comp),
459  eps);
460  Real Rt = tke*tke * N2 / diss0;
461 
462  // Turbulent length scale
463  Real length;
464  if (std::abs(N2) <= eps) {
465  length = l_g;
466  } else if (N2 > eps) {
467  // Stable (AL01, Eqn. 26)
468  length = std::sqrt(1.0 /
469  (1.0 / (l_g * l_g) + inv_Cb_sq * N2 / tke));
470  } else {
471  // Unstable (AL01, Eqn. 28)
472  length = l_g * std::sqrt(1.0 - Cmu0_pow3*Cmu0_pow3 * inv_Cb_sq * Rt);
473  }
474  mu_turb(i, j, k, EddyDiff::Turb_lengthscale) = length;
475 
476  // Burchard & Petersen smoothing function
477  Rt = (Rt >= Rt_crit) ? Rt
478  : std::max(Rt,
479  Rt - std::pow(Rt - Rt_crit, 2) /
480  (Rt + Rt_min - 2*Rt_crit));
481 
482  // Stability functions
483  // Note: These use the smoothed turbulent Richardson number
484  Real cmu = (Cmu0 + 0.108*Rt)
485  / (1.0 + 0.308*Rt + 0.00837*Rt*Rt); // (AL01, Eqn. 31)
486  Real cmu_prime = Cmu0 / (1 + 0.277*Rt); // (AL01, Eqn. 32)
487 
488  // Calculate eddy diffusivities
489  // K = rho * nu_t = rho * c_mu * tke^(1/2) * length
490  Real nut = cmu * std::sqrt(tke) * length; // eddy viscosity
491  Real nut_prime = cmu_prime / cmu * nut; // eddy diffusivity
492  mu_turb(i, j, k, EddyDiff::Mom_h) = cell_data(i, j, k, Rho_comp) * nut;
493  mu_turb(i, j, k, EddyDiff::Mom_v) = mu_turb(i, j, k, EddyDiff::Mom_h);
494  mu_turb(i, j, k, EddyDiff::Theta_v) = cell_data(i, j, k, Rho_comp) * nut_prime;
495 
496  // Calculate SFS quantities
497  // - dissipation (AL01 Eqn. 19)
498  diss(i, j, k) = cell_data(i, j, k, Rho_comp) * Cmu0_pow3 * std::pow(tke,1.5) / length;
499 
500  // - heat flux
501  // (Note: If using SurfaceLayer, the value at k=0 will
502  // be overwritten)
503  hfx_x(i, j, k) = 0.0;
504  hfx_y(i, j, k) = 0.0;
505  // Note: buoyant production = g/theta0 * hfx == -nut_prime * N^2 (c.f. AL01 Eqn. 15)
506  // = nut_prime * g/theta0 * dtheta/dz
507  // ==> hfx = nut_prime * dtheta/dz
508  // Our convention is such that dtheta/dz < 0 gives a positive
509  // (upward) heat flux.
510  hfx_z(i, j, k) = -mu_turb(i, j, k, EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
511  });
512  }
513  }
514 
515  // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0)
516  //***********************************************************************************
517  int ngc(1);
518  // EddyDiff mapping : Theta_h KE_h Scalar_h Q_h
519  Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
520  Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
521  Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
522  Real* fac_ptr = d_Factors.data();
523 
524  const bool use_KE = ( turbChoice.rans_type == RANSType::kEqn );
525 
526 #ifdef _OPENMP
527 #pragma omp parallel if (Gpu::notInLaunchRegion())
528 #endif
529  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
530  {
531  Box bxcc = mfi.tilebox();
532  Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
533  Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1);
534  bxcc.growLo(0,ngc); bxcc.growHi(0,ngc);
535  bxcc.growLo(1,ngc); bxcc.growHi(1,ngc);
536 
537  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
538 
539  for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) {
540  int offset = (EddyDiff::NumDiffs-1)/2;
541  switch (n)
542  {
543  case EddyDiff::KE_h:
544  if (use_KE) {
545  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
546  {
547  int indx = n;
548  int indx_v = indx + offset;
549  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
550  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
551  });
552  }
553  break;
554  default:
555  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
556  {
557  int indx = n;
558  int indx_v = indx + offset;
559 
560  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
561 
562  // NOTE: Theta_v has already been set for Deardorff
563  if (!(indx_v == EddyDiff::Theta_v && use_KE)) {
564  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
565  }
566  });
567  break;
568  }
569  }
570  }
571 }
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
amrex::Real Rt_min
Definition: ERF_TurbStruct.H:353
amrex::Real Cmu0
Definition: ERF_TurbStruct.H:350
amrex::Real Cb
Definition: ERF_TurbStruct.H:351
amrex::Real Rt_crit
Definition: ERF_TurbStruct.H:352
amrex::Real l_g_max
Definition: ERF_TurbStruct.H:354

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:
Here is the caller graph for this function: