ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_DiffusionSrcForState_T.cpp File Reference
#include <ERF_Diffusion.H>
#include <ERF_EddyViscosity.H>
#include <ERF_TerrainMetrics.H>
#include <ERF_PBLModels.H>
Include dependency graph for ERF_DiffusionSrcForState_T.cpp:

Functions

void DiffusionSrcForState_T (const Box &bx, const Box &domain, int start_comp, int num_comp, const bool &exp_most, const bool &rot_most, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &cell_data, const Array4< const Real > &cell_prim, const Array4< Real > &cell_rhs, const Array4< Real > &xflux, const Array4< Real > &yflux, const Array4< Real > &zflux, const Array4< const Real > &z_nd, const Array4< const Real > &ax, const Array4< const Real > &ay, const Array4< const Real > &az, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &SmnSmn_a, const Array4< const Real > &mf_m, const Array4< const Real > &mf_u, const Array4< const Real > &mf_v, Array4< Real > &hfx_x, Array4< Real > &hfx_y, Array4< Real > &hfx_z, Array4< Real > &qfx1_x, Array4< Real > &qfx1_y, Array4< Real > &qfx1_z, Array4< Real > &qfx2_z, Array4< Real > &diss, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const int level, const Array4< const Real > &tm_arr, const GpuArray< Real, AMREX_SPACEDIM > grav_gpu, const BCRec *bc_ptr, const bool use_most)
 

Function Documentation

◆ DiffusionSrcForState_T()

void DiffusionSrcForState_T ( const Box &  bx,
const Box &  domain,
int  start_comp,
int  num_comp,
const bool &  exp_most,
const bool &  rot_most,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  cell_prim,
const Array4< Real > &  cell_rhs,
const Array4< Real > &  xflux,
const Array4< Real > &  yflux,
const Array4< Real > &  zflux,
const Array4< const Real > &  z_nd,
const Array4< const Real > &  ax,
const Array4< const Real > &  ay,
const Array4< const Real > &  az,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
const Array4< const Real > &  SmnSmn_a,
const Array4< const Real > &  mf_m,
const Array4< const Real > &  mf_u,
const Array4< const Real > &  mf_v,
Array4< Real > &  hfx_x,
Array4< Real > &  hfx_y,
Array4< Real > &  hfx_z,
Array4< Real > &  qfx1_x,
Array4< Real > &  qfx1_y,
Array4< Real > &  qfx1_z,
Array4< Real > &  qfx2_z,
Array4< Real > &  diss,
const Array4< const Real > &  mu_turb,
const SolverChoice &  solverChoice,
const int  level,
const Array4< const Real > &  tm_arr,
const GpuArray< Real, AMREX_SPACEDIM >  grav_gpu,
const BCRec *  bc_ptr,
const bool  use_most 
)

Function for computing the scalar RHS for diffusion operator without terrain.

Parameters
[in]bxcell center box to loop over
[in]domainbox of the whole domain
[in]start_compstarting component index
[in]num_compnumber of components
[in]uvelocity in x-dir
[in]vvelocity in y-dir
[in]cell_dataconserved cell center vars
[in]cell_primprimitive cell center vars
[out]cell_rhsRHS for cell center vars
[in]xfluxflux in x-dir
[in]yfluxflux in y-dir
[in]zfluxflux in z-dir
[in]z_ndphysical z height
[in]detJJacobian determinant
[in]cellSizeInvinverse cell size array
[in]SmnSmn_astrain rate magnitude
[in]mf_mmap factor at cell center
[in]mf_umap factor at x-face
[in]mf_vmap factor at y-face
[in,out]hfx_zheat flux in z-dir
[in,out]qfx1_zheat flux in z-dir
[out]qfx2_zheat flux in z-dir
[in]dissdissipation of TKE
[in]mu_turbturbulent viscosity
[in]diffChoicecontainer of diffusion parameters
[in]turbChoicecontainer of turbulence parameters
[in]tm_arrtheta mean array
[in]grav_gpugravity vector
[in]bc_ptrcontainer with boundary conditions
[in]use_mostwhether we have turned on MOST BCs
80 {
81  BL_PROFILE_VAR("DiffusionSrcForState_T()",DiffusionSrcForState_T);
82 
83  DiffChoice diffChoice = solverChoice.diffChoice;
84  TurbChoice turbChoice = solverChoice.turbChoice[level];
85 
86  amrex::ignore_unused(use_most);
87 
88  const Real dx_inv = cellSizeInv[0];
89  const Real dy_inv = cellSizeInv[1];
90  const Real dz_inv = cellSizeInv[2];
91 
92  const auto& dom_hi = ubound(domain);
93 
94  Real l_inv_theta0 = 1.0 / turbChoice.theta_ref;
95  Real l_abs_g = std::abs(grav_gpu[2]);
96 
97  bool l_use_keqn = ( (turbChoice.les_type == LESType::Deardorff) ||
98  (turbChoice.rans_type == RANSType::kEqn) );
99  bool l_use_mynn = ( (turbChoice.pbl_type == PBLType::MYNN25) || (turbChoice.pbl_type == PBLType::MYNNEDMF) ) ;
100 
101  bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha);
102  bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) ||
103  (turbChoice.les_type == LESType::Deardorff ) ||
104  (turbChoice.rans_type == RANSType::kEqn ) ||
105  (turbChoice.pbl_type == PBLType::MYNN25 ) ||
106  (turbChoice.pbl_type == PBLType::MYNNEDMF ) ||
107  (turbChoice.pbl_type == PBLType::YSU ) );
108 
109  const Box xbx = surroundingNodes(bx,0);
110  const Box ybx = surroundingNodes(bx,1);
111  const Box zbx = surroundingNodes(bx,2);
112  Box zbx3 = zbx;
113 
114  const int end_comp = start_comp + num_comp - 1;
115 
116  // Theta, KE, Scalar
117  Vector<Real> alpha_eff(NPRIMVAR_max, 0.0);
118  if (l_consA) {
119  for (int i = 0; i < NPRIMVAR_max; ++i) {
120  switch (i) {
121  case PrimTheta_comp:
122  alpha_eff[PrimTheta_comp] = diffChoice.alpha_T;
123  break;
124  case PrimScalar_comp:
125  alpha_eff[PrimScalar_comp] = diffChoice.alpha_C;
126  break;
127  case PrimQ1_comp:
128  alpha_eff[PrimQ1_comp] = diffChoice.alpha_C;
129  break;
130  case PrimQ2_comp:
131  alpha_eff[PrimQ2_comp] = diffChoice.alpha_C;
132  break;
133  case PrimQ3_comp:
134  alpha_eff[PrimQ3_comp] = diffChoice.alpha_C;
135  break;
136  case PrimQ4_comp:
137  alpha_eff[PrimQ4_comp] = diffChoice.alpha_C;
138  break;
139  case PrimQ5_comp:
140  alpha_eff[PrimQ5_comp] = diffChoice.alpha_C;
141  break;
142  case PrimQ6_comp:
143  alpha_eff[PrimQ6_comp] = diffChoice.alpha_C;
144  break;
145  default:
146  alpha_eff[i] = 0.0;
147  break;
148  }
149  }
150  } else {
151  for (int i = 0; i < NPRIMVAR_max; ++i) {
152  switch (i) {
153  case PrimTheta_comp:
154  alpha_eff[PrimTheta_comp] = diffChoice.rhoAlpha_T;
155  break;
156  case PrimScalar_comp:
157  alpha_eff[PrimScalar_comp] = diffChoice.rhoAlpha_C;
158  break;
159  case PrimQ1_comp:
160  alpha_eff[PrimQ1_comp] = diffChoice.rhoAlpha_C;
161  break;
162  case PrimQ2_comp:
163  alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C;
164  break;
165  case PrimQ3_comp:
166  alpha_eff[PrimQ3_comp] = diffChoice.rhoAlpha_C;
167  break;
168  case PrimQ4_comp:
169  alpha_eff[PrimQ4_comp] = diffChoice.rhoAlpha_C;
170  break;
171  case PrimQ5_comp:
172  alpha_eff[PrimQ5_comp] = diffChoice.rhoAlpha_C;
173  break;
174  case PrimQ6_comp:
175  alpha_eff[PrimQ6_comp] = diffChoice.rhoAlpha_C;
176  break;
177  default:
178  alpha_eff[i] = 0.0;
179  break;
180  }
181  }
182  }
183 
184  Vector<int> eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h,
185  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h ,
186  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h };
187  Vector<int> eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h,
188  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h ,
189  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h };
190  Vector<int> eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::Scalar_v,
191  EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v ,
192  EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v };
193 
194  // Device vectors
195  Gpu::AsyncVector<Real> alpha_eff_d;
196  Gpu::AsyncVector<int> eddy_diff_idx_d,eddy_diff_idy_d,eddy_diff_idz_d;
197  alpha_eff_d.resize(alpha_eff.size());
198  eddy_diff_idx_d.resize(eddy_diff_idx.size());
199  eddy_diff_idy_d.resize(eddy_diff_idy.size());
200  eddy_diff_idz_d.resize(eddy_diff_idz.size());
201 
202  Gpu::copy(Gpu::hostToDevice, alpha_eff.begin() , alpha_eff.end() , alpha_eff_d.begin());
203  Gpu::copy(Gpu::hostToDevice, eddy_diff_idx.begin(), eddy_diff_idx.end(), eddy_diff_idx_d.begin());
204  Gpu::copy(Gpu::hostToDevice, eddy_diff_idy.begin(), eddy_diff_idy.end(), eddy_diff_idy_d.begin());
205  Gpu::copy(Gpu::hostToDevice, eddy_diff_idz.begin(), eddy_diff_idz.end(), eddy_diff_idz_d.begin());
206 
207  // Capture pointers for device code
208  Real* d_alpha_eff = alpha_eff_d.data();
209  int* d_eddy_diff_idx = eddy_diff_idx_d.data();
210  int* d_eddy_diff_idy = eddy_diff_idy_d.data();
211  int* d_eddy_diff_idz = eddy_diff_idz_d.data();
212 
213  // Constant alpha & Turb model
214  if (l_consA && l_turb) {
215  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
216  {
217  const int qty_index = start_comp + n;
218  const int prim_index = qty_index - 1;
219  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
220 
221  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
222  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
223  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index])
224  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) );
225 
226  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
227  Real met_h_zeta = ax(i,j,k);
228 
229  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
230  BCVars::RhoScalar_bc_comp : qty_index;
231  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
232 
233  bool most_on_zlo = ( use_most && rot_most &&
234  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
235 
236  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index)
237  - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) );
238  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
239 
240  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
241  xflux(i,j,k,qty_index) = hfx_x(i,j,0);
242  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
243  xflux(i,j,k,qty_index) = qfx1_x(i,j,0);
244  } else {
245  xflux(i,j,k,qty_index) = -rhoAlpha * mf_u(i,j,0) * ( GradCx - (met_h_xi/met_h_zeta)*GradCz );
246  }
247  });
248  ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
249  {
250  const int qty_index = start_comp + n;
251  const int prim_index = qty_index - 1;
252  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
253 
254  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
255  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
256  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index])
257  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) );
258 
259  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
260  Real met_h_zeta = ay(i,j,k);
261 
262  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
263  BCVars::RhoScalar_bc_comp : qty_index;
264  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
265  bool most_on_zlo = ( use_most && rot_most &&
266  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
267 
268  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index)
269  - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) );
270  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
271 
272  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
273  yflux(i,j,k,qty_index) = hfx_y(i,j,0);
274  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
275  yflux(i,j,k,qty_index) = qfx1_y(i,j,0);
276  } else {
277  yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz );
278  }
279  });
280  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
281  {
282  const int qty_index = start_comp + n;
283  const int prim_index = qty_index - 1;
284  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
285 
286  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
287  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
288  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index])
289  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) );
290 
291  Real met_h_zeta = az(i,j,k);
292 
293  Real GradCz;
294  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
295  BCVars::RhoScalar_bc_comp : qty_index;
296  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
297  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
298  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
299  && k == 0);
300  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) ||
301  (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim) )
302  && k == dom_hi.z+1);
303  bool most_on_zlo = ( use_most && exp_most &&
304  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
305 
306  if (ext_dir_on_zlo) {
307  GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
308  + 3. * cell_prim(i, j, k , prim_index)
309  - (1./3.) * cell_prim(i, j, k+1, prim_index) );
310  } else if (ext_dir_on_zhi) {
311  GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index)
312  - 3. * cell_prim(i, j, k-1, prim_index)
313  + (1./3.) * cell_prim(i, j, k-2, prim_index) );
314  } else {
315  GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
316  }
317 
318  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
319  // set the exact value from MOST, don't need finite diff
320  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
321  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
322  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
323  } else {
324  zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
325  }
326 
327  if (qty_index == RhoTheta_comp) {
328  if (!most_on_zlo) {
329  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
330  }
331  } else if (qty_index == RhoQ1_comp) {
332  if (!most_on_zlo) {
333  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
334  }
335  } else if (qty_index == RhoQ2_comp) {
336  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
337  }
338  });
339  // Constant rho*alpha & Turb model
340  } else if (l_turb) {
341  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
342  {
343  const int qty_index = start_comp + n;
344  const int prim_index = qty_index - 1;
345 
346  Real rhoAlpha = d_alpha_eff[prim_index];
347  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
348  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
349 
350  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
351  Real met_h_zeta = ax(i,j,k);
352 
353  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
354  BCVars::RhoScalar_bc_comp : qty_index;
355  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
356  bool most_on_zlo = ( use_most && rot_most &&
357  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
358 
359  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index)
360  - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) );
361  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
362 
363  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
364  xflux(i,j,k,qty_index) = hfx_x(i,j,0);
365  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
366  xflux(i,j,k,qty_index) = qfx1_x(i,j,0);
367  } else {
368  xflux(i,j,k,qty_index) = -rhoAlpha * mf_u(i,j,0) * ( GradCx - (met_h_xi/met_h_zeta)*GradCz );
369  }
370  });
371  ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
372  {
373  const int qty_index = start_comp + n;
374  const int prim_index = qty_index - 1;
375 
376  Real rhoAlpha = d_alpha_eff[prim_index];
377  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
378  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
379 
380  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
381  Real met_h_zeta = ay(i,j,k);
382 
383  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
384  BCVars::RhoScalar_bc_comp : qty_index;
385  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
386  bool most_on_zlo = ( use_most && rot_most &&
387  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
388 
389  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index)
390  - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) );
391  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
392 
393  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
394  yflux(i,j,k,qty_index) = hfx_y(i,j,0);
395  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
396  yflux(i,j,k,qty_index) = qfx1_y(i,j,0);
397  } else {
398  yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz );
399  }
400  });
401  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
402  {
403  const int qty_index = start_comp + n;
404  const int prim_index = qty_index - 1;
405 
406  Real rhoAlpha = d_alpha_eff[prim_index];
407 
408  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
409  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
410 
411  Real met_h_zeta = az(i,j,k);
412 
413  Real GradCz;
414  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
415  BCVars::RhoScalar_bc_comp : qty_index;
416  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
417  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
418  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
419  && k == 0);
420  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) ||
421  (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim))
422  && k == dom_hi.z+1);
423  bool most_on_zlo = ( use_most && exp_most &&
424  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
425 
426  if (ext_dir_on_zlo) {
427  GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
428  + 3. * cell_prim(i, j, k , prim_index)
429  - (1./3.) * cell_prim(i, j, k+1, prim_index) );
430  } else if (ext_dir_on_zhi) {
431  GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index)
432  - 3. * cell_prim(i, j, k-1, prim_index)
433  + (1./3.) * cell_prim(i, j, k-2, prim_index) );
434  } else {
435  GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
436  }
437 
438  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
439  // set the exact value from MOST, don't need finite diff
440  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
441  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
442  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
443  } else {
444  zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
445  }
446 
447  if (qty_index == RhoTheta_comp) {
448  if (!most_on_zlo) {
449  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
450  }
451  } else if (qty_index == RhoQ1_comp) {
452  if (!most_on_zlo) {
453  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
454  }
455  } else if (qty_index == RhoQ2_comp) {
456  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
457  }
458  });
459  // Constant alpha & no LES/PBL model
460  } else if(l_consA) {
461  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
462  {
463  const int qty_index = start_comp + n;
464  const int prim_index = qty_index - 1;
465 
466  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
467  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
468 
469  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
470  Real met_h_zeta = ax(i,j,k);
471 
472  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
473  BCVars::RhoScalar_bc_comp : qty_index;
474  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
475  bool most_on_zlo = ( use_most && rot_most &&
476  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
477 
478  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index)
479  - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) );
480  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
481 
482  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
483  xflux(i,j,k,qty_index) = hfx_x(i,j,0);
484  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
485  xflux(i,j,k,qty_index) = qfx1_x(i,j,0);
486  } else {
487  xflux(i,j,k,qty_index) = -rhoAlpha * mf_u(i,j,0) * ( GradCx - (met_h_xi/met_h_zeta)*GradCz );
488  }
489  });
490  ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
491  {
492  const int qty_index = start_comp + n;
493  const int prim_index = qty_index - 1;
494 
495  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
496  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
497 
498  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
499  Real met_h_zeta = ay(i,j,k);
500 
501  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
502  BCVars::RhoScalar_bc_comp : qty_index;
503  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
504  bool most_on_zlo = ( use_most && rot_most &&
505  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
506 
507  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index)
508  - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) );
509  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
510 
511  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
512  yflux(i,j,k,qty_index) = hfx_y(i,j,0);
513  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
514  yflux(i,j,k,qty_index) = qfx1_y(i,j,0);
515  } else {
516  yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz );
517  }
518  });
519  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
520  {
521  const int qty_index = start_comp + n;
522  const int prim_index = qty_index - 1;
523 
524  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
525  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
526 
527  Real met_h_zeta = az(i,j,k);
528 
529  Real GradCz;
530  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
531  BCVars::RhoScalar_bc_comp : qty_index;
532  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
533  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
534  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
535  && k == 0);
536  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) ||
537  (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim))
538  && k == dom_hi.z+1);
539  bool most_on_zlo = ( use_most && exp_most &&
540  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
541 
542  if (ext_dir_on_zlo) {
543  GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
544  + 3. * cell_prim(i, j, k , prim_index)
545  - (1./3.) * cell_prim(i, j, k+1, prim_index) );
546  } else if (ext_dir_on_zhi) {
547  GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index)
548  - 3. * cell_prim(i, j, k-1, prim_index)
549  + (1./3.) * cell_prim(i, j, k-2, prim_index) );
550  } else {
551  GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
552  }
553 
554  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
555  // set the exact value from MOST, don't need finite diff
556  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
557  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
558  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
559  } else {
560  zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
561  }
562 
563  if (qty_index == RhoTheta_comp) {
564  if (!most_on_zlo) {
565  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
566  }
567  } else if (qty_index == RhoQ1_comp) {
568  if (!most_on_zlo) {
569  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
570  }
571  } else if (qty_index == RhoQ2_comp) {
572  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
573  }
574  });
575  // Constant rho*alpha & no LES/PBL model
576  } else {
577  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
578  {
579  const int qty_index = start_comp + n;
580  const int prim_index = qty_index - 1;
581 
582  Real rhoAlpha = d_alpha_eff[prim_index];
583 
584  Real met_h_xi,met_h_zeta;
585  met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
586  met_h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd);
587 
588  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
589  BCVars::RhoScalar_bc_comp : qty_index;
590  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
591  bool most_on_zlo = ( use_most && rot_most &&
592  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
593 
594  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index)
595  - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) );
596  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
597 
598  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
599  xflux(i,j,k,qty_index) = hfx_x(i,j,0);
600  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
601  xflux(i,j,k,qty_index) = qfx1_x(i,j,0);
602  } else {
603  xflux(i,j,k,qty_index) = -rhoAlpha * mf_u(i,j,0) * ( GradCx - (met_h_xi/met_h_zeta)*GradCz );
604  }
605  });
606  ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
607  {
608  const int qty_index = start_comp + n;
609  const int prim_index = qty_index - 1;
610 
611  Real rhoAlpha = d_alpha_eff[prim_index];
612 
613  Real met_h_eta,met_h_zeta;
614  met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
615  met_h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd);
616 
617  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
618  BCVars::RhoScalar_bc_comp : qty_index;
619  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
620  bool most_on_zlo = ( use_most && rot_most &&
621  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
622 
623  Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index)
624  - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) );
625  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
626 
627  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
628  yflux(i,j,k,qty_index) = hfx_y(i,j,0);
629  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
630  yflux(i,j,k,qty_index) = qfx1_y(i,j,0);
631  } else {
632  yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz );
633  }
634  });
635  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
636  {
637  const int qty_index = start_comp + n;
638  const int prim_index = qty_index - 1;
639 
640  Real rhoAlpha = d_alpha_eff[prim_index];
641 
642  Real met_h_zeta;
643  met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd);
644 
645  Real GradCz;
646  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
647  BCVars::RhoScalar_bc_comp : qty_index;
648  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
649  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
650  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
651  && k == 0);
652  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) ||
653  (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim))
654  && k == dom_hi.z+1);
655  bool most_on_zlo = ( use_most && exp_most &&
656  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0);
657 
658  if (ext_dir_on_zlo) {
659  GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
660  + 3. * cell_prim(i, j, k , prim_index)
661  - (1./3.) * cell_prim(i, j, k+1, prim_index) );
662  } else if (ext_dir_on_zhi) {
663  GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index)
664  - 3. * cell_prim(i, j, k-1, prim_index)
665  + (1./3.) * cell_prim(i, j, k-2, prim_index) );
666  } else {
667  GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
668  }
669 
670  if (most_on_zlo && (qty_index == RhoTheta_comp)) {
671  // set the exact value from MOST, don't need finite diff
672  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
673  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
674  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
675  } else {
676  zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
677  }
678 
679  if (qty_index == RhoTheta_comp) {
680  if (!most_on_zlo) {
681  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
682  }
683  } else if (qty_index == RhoQ1_comp) {
684  if (!most_on_zlo) {
685  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
686  }
687  } else if (qty_index == RhoQ2_comp) {
688  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
689  }
690  });
691  }
692 
693  // Linear combinations for z-flux with terrain
694  //-----------------------------------------------------------------------------------
695  // Extrapolate top and bottom cells
696  {
697  Box planexy = zbx; planexy.setBig(2, planexy.smallEnd(2) );
698  int k_lo = zbx.smallEnd(2); int k_hi = zbx.bigEnd(2);
699  zbx3.growLo(2,-1); zbx3.growHi(2,-1);
700  ParallelFor(planexy, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int , int n) noexcept
701  {
702  const int qty_index = start_comp + n;
703  Real met_h_xi,met_h_eta;
704 
705  { // Bottom face
706  met_h_xi = Compute_h_xi_AtKface (i,j,k_lo,cellSizeInv,z_nd);
707  met_h_eta = Compute_h_eta_AtKface(i,j,k_lo,cellSizeInv,z_nd);
708 
709  Real xfluxlo = 0.5 * ( xflux(i , j , k_lo , qty_index) + xflux(i+1, j , k_lo , qty_index) );
710  Real xfluxhi = 0.5 * ( xflux(i , j , k_lo+1, qty_index) + xflux(i+1, j , k_lo+1, qty_index) );
711  Real xfluxbar = 1.5*xfluxlo - 0.5*xfluxhi;
712 
713  Real yfluxlo = 0.5 * ( yflux(i , j , k_lo , qty_index) + yflux(i , j+1, k_lo , qty_index) );
714  Real yfluxhi = 0.5 * ( yflux(i , j , k_lo+1, qty_index) + yflux(i , j+1, k_lo+1, qty_index) );
715  Real yfluxbar = 1.5*yfluxlo - 0.5*yfluxhi;
716 
717  zflux(i,j,k_lo,qty_index) -= met_h_xi*xfluxbar + met_h_eta*yfluxbar;
718  }
719 
720  { // Top face
721  met_h_xi = Compute_h_xi_AtKface (i,j,k_hi,cellSizeInv,z_nd);
722  met_h_eta = Compute_h_eta_AtKface(i,j,k_hi,cellSizeInv,z_nd);
723 
724  Real xfluxlo = 0.5 * ( xflux(i , j , k_hi-2, qty_index) + xflux(i+1, j , k_hi-2, qty_index) );
725  Real xfluxhi = 0.5 * ( xflux(i , j , k_hi-1, qty_index) + xflux(i+1, j , k_hi-1, qty_index) );
726  Real xfluxbar = 1.5*xfluxhi - 0.5*xfluxlo;
727 
728  Real yfluxlo = 0.5 * ( yflux(i , j , k_hi-2, qty_index) + yflux(i , j+1, k_hi-2, qty_index) );
729  Real yfluxhi = 0.5 * ( yflux(i , j , k_hi-1, qty_index) + yflux(i , j+1, k_hi-1, qty_index) );
730  Real yfluxbar = 1.5*yfluxhi - 0.5*yfluxlo;
731 
732  zflux(i,j,k_hi,qty_index) -= met_h_xi*xfluxbar + met_h_eta*yfluxbar;
733  }
734  });
735  }
736  // Average interior cells
737  ParallelFor(zbx3, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
738  {
739  const int qty_index = start_comp + n;
740 
741  Real met_h_xi,met_h_eta;
742  met_h_xi = Compute_h_xi_AtKface (i,j,k,cellSizeInv,z_nd);
743  met_h_eta = Compute_h_eta_AtKface(i,j,k,cellSizeInv,z_nd);
744 
745  Real xfluxbar = 0.25 * ( xflux(i , j , k , qty_index) + xflux(i+1, j , k , qty_index)
746  + xflux(i , j , k-1, qty_index) + xflux(i+1, j , k-1, qty_index) );
747  Real yfluxbar = 0.25 * ( yflux(i , j , k , qty_index) + yflux(i , j+1, k , qty_index)
748  + yflux(i , j , k-1, qty_index) + yflux(i , j+1, k-1, qty_index) );
749 
750  zflux(i,j,k,qty_index) -= met_h_xi*xfluxbar + met_h_eta*yfluxbar;
751  });
752  // Multiply h_zeta by x/y-fluxes
753  ParallelFor(xbx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
754  {
755  const int qty_index = start_comp + n;
756  Real met_h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd);
757  xflux(i,j,k,qty_index) *= met_h_zeta;
758  });
759  ParallelFor(ybx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
760  {
761  const int qty_index = start_comp + n;
762  Real met_h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd);
763  yflux(i,j,k,qty_index) *= met_h_zeta;
764  });
765 
766 
767  // Use fluxes to compute RHS
768  //-----------------------------------------------------------------------------------
769  for (int n(0); n < num_comp; ++n)
770  {
771  int qty_index = start_comp + n;
772  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
773  {
774  Real stateContrib = (xflux(i+1,j ,k ,qty_index) - xflux(i, j, k, qty_index)) * dx_inv * mf_m(i,j,0) // Diffusive flux in x-dir
775  +(yflux(i ,j+1,k ,qty_index) - yflux(i, j, k, qty_index)) * dy_inv * mf_m(i,j,0) // Diffusive flux in y-dir
776  +(zflux(i ,j ,k+1,qty_index) - zflux(i, j, k, qty_index)) * dz_inv; // Diffusive flux in z-dir
777 
778  stateContrib /= detJ(i,j,k);
779 
780  cell_rhs(i,j,k,qty_index) -= stateContrib;
781  });
782  }
783 
784  // Using Deardorff (see Sullivan et al 1994)
785  // or k-eqn RANS (see Axell & Liungman 2001)
786  //
787  // Note: At this point, the thermal diffusivity ("Khv" field in ERF), the
788  // subgrid heat flux ("hfx_z" here), and the subgrid dissipation
789  // ("diss" here) have been updated by ComputeTurbulentViscosityLES --
790  // at the beginning of each timestep.
791  // The strain rate magnitude is updated at the beginning of the first
792  // RK stage only, therefore the shear production term also does not
793  // change between RK stages.
794  // The surface heat flux hfx_z(i,j,-1) is updated in MOSTStress at
795  // each RK stage if using the ERF_EXPLICIT_MOST_STRESS path, but that
796  // does not change the buoyancy production term here.
797  if (l_use_keqn && (start_comp <= RhoKE_comp) && (end_comp >= RhoKE_comp)) {
798  int qty_index = RhoKE_comp;
799  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
800  {
801  // Add Buoyancy Source
802  // where the SGS buoyancy flux tau_{theta,i} = -KH * dtheta/dx_i,
803  // such that for dtheta/dz < 0, there is a positive (upward) heat
804  // flux; the TKE buoyancy production is then
805  // B = g/theta_0 * tau_{theta,w}
806  // for a dry atmosphere.
807  // TODO: To account for moisture, the Brunt-Vaisala frequency,
808  // N^2 = g[1/theta * dtheta/dz + ...]
809  // **should** be a function of the water vapor and total water
810  // mixing ratios, depending on whether conditions are saturated or
811  // not (see the WRF model description, Skamarock et al 2019).
812  cell_rhs(i,j,k,qty_index) += l_abs_g * l_inv_theta0 * hfx_z(i,j,k);
813 
814  // TKE shear production
815  // P = -tau_ij * S_ij = 2 * mu_turb * S_ij * S_ij
816  // Note: This assumes that the horizontal and vertical diffusivities
817  // of momentum are equal
818  cell_rhs(i,j,k,qty_index) += 2.0*mu_turb(i,j,k,EddyDiff::Mom_v) * SmnSmn_a(i,j,k);
819 
820  // TKE dissipation
821  cell_rhs(i,j,k,qty_index) -= diss(i,j,k);
822  });
823  }
824 
825  // Using PBL
826  if (l_use_mynn && start_comp <= RhoKE_comp && end_comp >=RhoKE_comp) {
827  int qty_index = RhoKE_comp;
828  auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;
829 
830  const int rhoqv_comp = solverChoice.RhoQv_comp;
831  const int rhoqc_comp = solverChoice.RhoQc_comp;
832  const int rhoqr_comp = solverChoice.RhoQr_comp;
833 
834  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
835  {
836  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
837  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
838  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
839  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) );
840  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
841  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
842 
843  // This computes shear production, buoyancy production, and dissipation terms only.
844  cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
845  mu_turb,cellSizeInv,domain,
846  pbl_mynn_B1_l,tm_arr(i,j,0),
847  rhoqv_comp, rhoqc_comp, rhoqr_comp,
848  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
849  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
850  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
851  use_most);
852  });
853  }
854 }
MolecDiffType::ConstantAlpha
@ ConstantAlpha
DiffusionSrcForState_T
void DiffusionSrcForState_T(const Box &bx, const Box &domain, int start_comp, int num_comp, const bool &exp_most, const bool &rot_most, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &cell_data, const Array4< const Real > &cell_prim, const Array4< Real > &cell_rhs, const Array4< Real > &xflux, const Array4< Real > &yflux, const Array4< Real > &zflux, const Array4< const Real > &z_nd, const Array4< const Real > &ax, const Array4< const Real > &ay, const Array4< const Real > &az, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &SmnSmn_a, const Array4< const Real > &mf_m, const Array4< const Real > &mf_u, const Array4< const Real > &mf_v, Array4< Real > &hfx_x, Array4< Real > &hfx_y, Array4< Real > &hfx_z, Array4< Real > &qfx1_x, Array4< Real > &qfx1_y, Array4< Real > &qfx1_z, Array4< Real > &qfx2_z, Array4< Real > &diss, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const int level, const Array4< const Real > &tm_arr, const GpuArray< Real, AMREX_SPACEDIM > grav_gpu, const BCRec *bc_ptr, const bool use_most)
Definition: ERF_DiffusionSrcForState_T.cpp:43
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
PrimQ1_comp
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:53
NPRIMVAR_max
#define NPRIMVAR_max
Definition: ERF_IndexDefines.H:33
PrimQ2_comp
#define PrimQ2_comp
Definition: ERF_IndexDefines.H:54
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
PrimQ4_comp
#define PrimQ4_comp
Definition: ERF_IndexDefines.H:56
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
PrimQ6_comp
#define PrimQ6_comp
Definition: ERF_IndexDefines.H:58
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
PrimQ5_comp
#define PrimQ5_comp
Definition: ERF_IndexDefines.H:57
PrimTheta_comp
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:50
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
PrimScalar_comp
#define PrimScalar_comp
Definition: ERF_IndexDefines.H:52
PrimQ3_comp
#define PrimQ3_comp
Definition: ERF_IndexDefines.H:55
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
ComputeQKESourceTerms
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeQKESourceTerms(int i, int j, int k, const amrex::Array4< const amrex::Real > &uvel, const amrex::Array4< const amrex::Real > &vvel, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &K_turb, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Box &domain, amrex::Real pbl_mynn_B1_l, const amrex::Real theta_mean, const int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp, bool c_ext_dir_on_zlo, bool c_ext_dir_on_zhi, bool u_ext_dir_on_zlo, bool u_ext_dir_on_zhi, bool v_ext_dir_on_zlo, bool v_ext_dir_on_zhi, const bool use_most=false, const amrex::Real met_h_zeta=1.0)
Definition: ERF_PBLModels.H:197
Compute_h_xi_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtIface(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:108
Compute_h_zeta_AtKface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtKface(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:180
Compute_h_zeta_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(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:94
Compute_h_xi_AtKface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtKface(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:195
Compute_h_zeta_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(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:137
Compute_h_eta_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(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:165
Compute_h_eta_AtKface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtKface(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:209
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:80
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:88
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:76
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:87
ERFBCType::foextrap
@ foextrap
Definition: ERF_IndexDefines.H:190
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:191
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:193
EddyDiff::Theta_v
@ Theta_v
Definition: ERF_IndexDefines.H:168
EddyDiff::Scalar_v
@ Scalar_v
Definition: ERF_IndexDefines.H:170
EddyDiff::Q_v
@ Q_v
Definition: ERF_IndexDefines.H:171
EddyDiff::Q_h
@ Q_h
Definition: ERF_IndexDefines.H:166
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:167
EddyDiff::Scalar_h
@ Scalar_h
Definition: ERF_IndexDefines.H:165
EddyDiff::KE_v
@ KE_v
Definition: ERF_IndexDefines.H:169
EddyDiff::Theta_h
@ Theta_h
Definition: ERF_IndexDefines.H:163
EddyDiff::KE_h
@ KE_h
Definition: ERF_IndexDefines.H:164
DiffChoice
Definition: ERF_DiffStruct.H:19
DiffChoice::rhoAlpha_C
amrex::Real rhoAlpha_C
Definition: ERF_DiffStruct.H:92
DiffChoice::rhoAlpha_T
amrex::Real rhoAlpha_T
Definition: ERF_DiffStruct.H:91
DiffChoice::alpha_T
amrex::Real alpha_T
Definition: ERF_DiffStruct.H:84
DiffChoice::alpha_C
amrex::Real alpha_C
Definition: ERF_DiffStruct.H:85
DiffChoice::molec_diff_type
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:81
MYNNLevel25::B1
amrex::Real B1
Definition: ERF_MYNNStruct.H:43
SolverChoice::RhoQr_comp
int RhoQr_comp
Definition: ERF_DataStruct.H:783
SolverChoice::diffChoice
DiffChoice diffChoice
Definition: ERF_DataStruct.H:674
SolverChoice::RhoQc_comp
int RhoQc_comp
Definition: ERF_DataStruct.H:777
SolverChoice::RhoQv_comp
int RhoQv_comp
Definition: ERF_DataStruct.H:776
SolverChoice::turbChoice
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:676
TurbChoice
Definition: ERF_TurbStruct.H:31
TurbChoice::pbl_mynn
MYNNLevel25 pbl_mynn
Definition: ERF_TurbStruct.H:242
TurbChoice::pbl_type
PBLType pbl_type
Definition: ERF_TurbStruct.H:240
TurbChoice::rans_type
RANSType rans_type
Definition: ERF_TurbStruct.H:237
TurbChoice::les_type
LESType les_type
Definition: ERF_TurbStruct.H:204
TurbChoice::theta_ref
amrex::Real theta_ref
Definition: ERF_TurbStruct.H:234
Here is the call graph for this function: