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