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

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 
166  if (smag2d) {
167  mu_turb(i, j, k, EddyDiff::Mom_v) = 0.0;
168  } else {
169  mu_turb(i, j, k, EddyDiff::Mom_v) = CsDeltaSqr * cell_data(i, j, k, Rho_comp) * std::sqrt(2.0*SmnSmn);
170  }
171  }
172 
173  // =====================================================================
174  // HEAT FLUX CALCULATION
175  // =====================================================================
176  Real dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
177  - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv;
178  hfx_x(i,j,k) = 0.0;
179  hfx_y(i,j,k) = 0.0;
180  hfx_z(i,j,k) = -inv_Pr_t*mu_turb(i,j,k,EddyDiff::Mom_v) * dtheta_dz;
181  });
182  }
183  }
184  // DEARDORFF: Fill Kturb for momentum in horizontal and vertical
185  //***********************************************************************************
186  else if (turbChoice.les_type == LESType::Deardorff)
187  {
188  const Real l_C_k = turbChoice.Ck;
189  const Real l_C_e = turbChoice.Ce;
190  const Real l_C_e_wall = turbChoice.Ce_wall;
191  const Real Ce_lcoeff = amrex::max(0.0, l_C_e - 1.9*l_C_k);
192  const Real l_abs_g = const_grav;
193 
194  const bool use_ref_theta = (turbChoice.theta_ref > 0);
195  const Real l_inv_theta0 = (use_ref_theta) ? 1.0 / turbChoice.theta_ref : 1.0;
196 
197 #ifdef _OPENMP
198 #pragma omp parallel if (Gpu::notInLaunchRegion())
199 #endif
200  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
201  {
202  Box bxcc = mfi.tilebox();
203 
204  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
205  const Array4<Real>& hfx_x = Hfx1.array(mfi);
206  const Array4<Real>& hfx_y = Hfx2.array(mfi);
207  const Array4<Real>& hfx_z = Hfx3.array(mfi);
208  const Array4<Real>& diss = Diss.array(mfi);
209 
210  const Array4<Real const > &cell_data = cons_in.array(mfi);
211 
212  Array4<Real const> mf_u = mapfac[MapFacType::u_x]->const_array(mfi);
213  Array4<Real const> mf_v = mapfac[MapFacType::v_y]->const_array(mfi);
214 
215  Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
216 
217  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
218  {
219  Real dxInv = cellSizeInv[0];
220  Real dyInv = cellSizeInv[1];
221  Real dzInv = cellSizeInv[2];
222  if (use_terrain) {
223  // the terrain grid is only deformed in z for now
224  dzInv /= Compute_h_zeta_AtCellCenter(i,j,k, cellSizeInv, z_nd_arr);
225  }
226  Real Delta;
227  if (isotropic) {
228  Real cellVolMsf = 1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0) * dzInv);
229  Delta = std::cbrt(cellVolMsf);
230  } else {
231  Delta = 1.0 / dzInv;
232  }
233 
234  Real dtheta_dz;
235  if (use_thetav_grad) {
236  dtheta_dz = 0.5 * ( GetThetav(i, j, k+1, cell_data, moisture_indices)
237  - GetThetav(i, j, k-1, cell_data, moisture_indices) )*dzInv;
238  } else if (use_thetal_grad) {
239  dtheta_dz = 0.5 * ( GetThetal(i, j, k+1, cell_data, moisture_indices)
240  - GetThetal(i, j, k-1, cell_data, moisture_indices) )*dzInv;
241  } else {
242  dtheta_dz = 0.5 * ( cell_data(i, j, k+1, RhoTheta_comp) / cell_data(i, j, k+1, Rho_comp)
243  - cell_data(i, j, k-1, RhoTheta_comp) / cell_data(i, j, k-1, Rho_comp) )*dzInv;
244  }
245 
246  // Calculate stratification-dependent mixing length (Deardorff 1980, Eqn. 10a)
247  Real E = amrex::max(cell_data(i,j,k,RhoKE_comp)/cell_data(i,j,k,Rho_comp),Real(0.0));
248  Real stratification = l_abs_g * dtheta_dz * l_inv_theta0;
249  if (!use_ref_theta) {
250  // l_inv_theta0 == 1, divide by actual theta
251  stratification *= cell_data(i,j,k,Rho_comp) /
252  cell_data(i,j,k,RhoTheta_comp);
253  }
254 
255  // Following WRF, the stratification effects are applied to the vertical length scales
256  // in the case of anistropic mixing
257  Real length;
259  if (stratification <= eps) {
260  length = Delta; // cbrt(dx*dy*dz) -or- dz
261  } else {
262  length = 0.76 * std::sqrt(E / amrex::max(stratification,eps));
263  // mixing length should be _reduced_ for stable stratification
264  length = amrex::min(length, Delta);
265  // following WRF, make sure the mixing length isn't too small
266  length = amrex::max(length, 0.001 * Delta);
267  }
268 
269  Real DeltaH = (isotropic) ? length : std::sqrt(1.0 / (dxInv * mf_u(i,j,0) * dyInv * mf_v(i,j,0)));
270 
271  Real Pr_inv_v = (1. + 2.*length/Delta);
272  Real Pr_inv_h = (isotropic) ? Pr_inv_v : inv_Pr_t;
273 
274  // Calculate eddy diffusivities
275  // K = rho * C_k * l * KE^(1/2)
276  mu_turb(i,j,k,EddyDiff::Mom_h) = cell_data(i,j,k,Rho_comp) * l_C_k * DeltaH * std::sqrt(E);
277  mu_turb(i,j,k,EddyDiff::Mom_v) = cell_data(i,j,k,Rho_comp) * l_C_k * length * std::sqrt(E);
278  // KH = (1 + 2*l/delta) * mu_turb
279  mu_turb(i,j,k,EddyDiff::Theta_h) = Pr_inv_h * mu_turb(i,j,k,EddyDiff::Mom_h);
280  mu_turb(i,j,k,EddyDiff::Theta_v) = Pr_inv_v * mu_turb(i,j,k,EddyDiff::Mom_v);
281  // Store lengthscale for TKE source terms
282  mu_turb(i,j,k,EddyDiff::Turb_lengthscale) = length;
283 
284  // Calculate SFS quantities
285  // - dissipation
286  Real Ce;
287  if ((l_C_e_wall > 0) && (k==0)) {
288  Ce = l_C_e_wall;
289  } else {
290  Ce = 1.9*l_C_k + Ce_lcoeff*length / Delta;
291  }
292  diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length;
293 
294  // - heat flux
295  // (Note: If using SurfaceLayer, the value at k=0 will
296  // be overwritten)
297  hfx_x(i,j,k) = 0.0;
298  hfx_y(i,j,k) = 0.0;
299  hfx_z(i,j,k) = -mu_turb(i,j,k,EddyDiff::Theta_v) * dtheta_dz; // (rho*w)' theta' [kg m^-2 s^-1 K]
300  });
301  }
302  }
303 
304  // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0)
305  //***********************************************************************************
306  int ngc(1);
307  // EddyDiff mapping : Theta_h KE_h Scalar_h Q_h
308  Vector<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
309  Gpu::AsyncVector<Real> d_Factors; d_Factors.resize(Factors.size());
310  Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin());
311  Real* fac_ptr = d_Factors.data();
312 
313  const bool use_KE = ( turbChoice.les_type == LESType::Deardorff );
314 
315 #ifdef _OPENMP
316 #pragma omp parallel if (Gpu::notInLaunchRegion())
317 #endif
318  for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
319  {
320  Box bxcc = mfi.tilebox();
321  Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
322  Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1);
323  bxcc.growLo(0,ngc); bxcc.growHi(0,ngc);
324  bxcc.growLo(1,ngc); bxcc.growHi(1,ngc);
325 
326  const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
327 
328  for (auto n = 1; n < (EddyDiff::NumDiffs-1)/2; ++n) {
329  int offset = (EddyDiff::NumDiffs-1)/2;
330  switch (n)
331  {
332  case EddyDiff::KE_h:
333  if (use_KE) {
334  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
335  {
336  int indx = n;
337  int indx_v = indx + offset;
338  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
339  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
340  });
341  }
342  break;
343  default:
344  ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
345  {
346  int indx = n;
347  int indx_v = indx + offset;
348 
349  // NOTE: Theta_h, Theta_v have already been set for Deardorff
350  if (!(indx_v == EddyDiff::Theta_v && use_KE)) {
351  mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1];
352  mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx);
353  }
354  });
355  break;
356  }
357  }
358  }
359 }
@ 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:66
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:39
#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:337
amrex::Real sigma_k
Definition: ERF_TurbStruct.H:353
amrex::Real Sc_t_inv
Definition: ERF_TurbStruct.H:333
StratType strat_type
Definition: ERF_TurbStruct.H:358
amrex::Real Ck
Definition: ERF_TurbStruct.H:342
amrex::Real Cs
Definition: ERF_TurbStruct.H:336
amrex::Real Pr_t_inv
Definition: ERF_TurbStruct.H:329
amrex::Real Ce_wall
Definition: ERF_TurbStruct.H:341
amrex::Real Ce
Definition: ERF_TurbStruct.H:340
bool use_smag_stratification
Definition: ERF_TurbStruct.H:364
amrex::Real theta_ref
Definition: ERF_TurbStruct.H:356
bool mix_isotropic
Definition: ERF_TurbStruct.H:361

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

Referenced by ComputeTurbulentViscosity().

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