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