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

Functions

void DiffusionSrcForState_N (const Box &bx, const Box &domain, int start_comp, int num_comp, const bool &exp_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 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_z, 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_N()

void DiffusionSrcForState_N ( const Box &  bx,
const Box &  domain,
int  start_comp,
int  num_comp,
const bool &  exp_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 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_z,
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]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
67 {
68  BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N);
69 
70  DiffChoice diffChoice = solverChoice.diffChoice;
71  TurbChoice turbChoice = solverChoice.turbChoice[level];
72 
73  amrex::ignore_unused(use_most);
74 
75  const Real dx_inv = cellSizeInv[0];
76  const Real dy_inv = cellSizeInv[1];
77  const Real dz_inv = cellSizeInv[2];
78 
79  const auto& dom_lo = lbound(domain);
80  const auto& dom_hi = ubound(domain);
81 
82  Real l_inv_theta0 = 1.0 / turbChoice.theta_ref;
83  Real l_abs_g = std::abs(grav_gpu[2]);
84 
85  bool l_use_keqn = ( (turbChoice.les_type == LESType::Deardorff) ||
86  (turbChoice.rans_type == RANSType::kEqn) );
87  bool l_use_mynn = ( (turbChoice.pbl_type == PBLType::MYNN25) || (turbChoice.pbl_type == PBLType::MYNNEDMF) ) ;
88 
89  bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha);
90  bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) ||
91  (turbChoice.les_type == LESType::Deardorff ) ||
92  (turbChoice.rans_type == RANSType::kEqn ) ||
93  (turbChoice.pbl_type == PBLType::MYNN25 ) ||
94  (turbChoice.pbl_type == PBLType::MYNNEDMF ) ||
95  (turbChoice.pbl_type == PBLType::YSU ) );
96 
97  const Box xbx = surroundingNodes(bx,0);
98  const Box ybx = surroundingNodes(bx,1);
99  const Box zbx = surroundingNodes(bx,2);
100 
101  const int end_comp = start_comp + num_comp - 1;
102 
103  // Theta, KE, Scalar
104  Vector<Real> alpha_eff(NPRIMVAR_max, 0.0);
105  if (l_consA) {
106  for (int i = 0; i < NPRIMVAR_max; ++i) {
107  switch (i) {
108  case PrimTheta_comp:
109  alpha_eff[PrimTheta_comp] = diffChoice.alpha_T;
110  break;
111  case PrimScalar_comp:
112  alpha_eff[PrimScalar_comp] = diffChoice.alpha_C;
113  break;
114  case PrimQ1_comp:
115  alpha_eff[PrimQ1_comp] = diffChoice.alpha_C;
116  break;
117  case PrimQ2_comp:
118  alpha_eff[PrimQ2_comp] = diffChoice.alpha_C;
119  break;
120  case PrimQ3_comp:
121  alpha_eff[PrimQ3_comp] = diffChoice.alpha_C;
122  break;
123  case PrimQ4_comp:
124  alpha_eff[PrimQ4_comp] = diffChoice.alpha_C;
125  break;
126  case PrimQ5_comp:
127  alpha_eff[PrimQ5_comp] = diffChoice.alpha_C;
128  break;
129  case PrimQ6_comp:
130  alpha_eff[PrimQ6_comp] = diffChoice.alpha_C;
131  break;
132  default:
133  alpha_eff[i] = 0.0;
134  break;
135  }
136  }
137  } else {
138  for (int i = 0; i < NPRIMVAR_max; ++i) {
139  switch (i) {
140  case PrimTheta_comp:
141  alpha_eff[PrimTheta_comp] = diffChoice.rhoAlpha_T;
142  break;
143  case PrimScalar_comp:
144  alpha_eff[PrimScalar_comp] = diffChoice.rhoAlpha_C;
145  break;
146  case PrimQ1_comp:
147  alpha_eff[PrimQ1_comp] = diffChoice.rhoAlpha_C;
148  break;
149  case PrimQ2_comp:
150  alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C;
151  break;
152  case PrimQ3_comp:
153  alpha_eff[PrimQ3_comp] = diffChoice.rhoAlpha_C;
154  break;
155  case PrimQ4_comp:
156  alpha_eff[PrimQ4_comp] = diffChoice.rhoAlpha_C;
157  break;
158  case PrimQ5_comp:
159  alpha_eff[PrimQ5_comp] = diffChoice.rhoAlpha_C;
160  break;
161  case PrimQ6_comp:
162  alpha_eff[PrimQ6_comp] = diffChoice.rhoAlpha_C;
163  break;
164  default:
165  alpha_eff[i] = 0.0;
166  break;
167  }
168  }
169  }
170 
171  Vector<int> eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h,
172  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h ,
173  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h };
174  Vector<int> eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h,
175  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h ,
176  EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h };
177  Vector<int> eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::Scalar_v,
178  EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v ,
179  EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v };
180 
181  // Device vectors
182  Gpu::AsyncVector<Real> alpha_eff_d;
183  Gpu::AsyncVector<int> eddy_diff_idx_d,eddy_diff_idy_d,eddy_diff_idz_d;
184  alpha_eff_d.resize(alpha_eff.size());
185  eddy_diff_idx_d.resize(eddy_diff_idx.size());
186  eddy_diff_idy_d.resize(eddy_diff_idy.size());
187  eddy_diff_idz_d.resize(eddy_diff_idz.size());
188 
189  Gpu::copy(Gpu::hostToDevice, alpha_eff.begin() , alpha_eff.end() , alpha_eff_d.begin());
190  Gpu::copy(Gpu::hostToDevice, eddy_diff_idx.begin(), eddy_diff_idx.end(), eddy_diff_idx_d.begin());
191  Gpu::copy(Gpu::hostToDevice, eddy_diff_idy.begin(), eddy_diff_idy.end(), eddy_diff_idy_d.begin());
192  Gpu::copy(Gpu::hostToDevice, eddy_diff_idz.begin(), eddy_diff_idz.end(), eddy_diff_idz_d.begin());
193 
194  // Capture pointers for device code
195  Real* d_alpha_eff = alpha_eff_d.data();
196  int* d_eddy_diff_idx = eddy_diff_idx_d.data();
197  int* d_eddy_diff_idy = eddy_diff_idy_d.data();
198  int* d_eddy_diff_idz = eddy_diff_idz_d.data();
199 
200  // Compute fluxes at each face
201  if (l_consA && l_turb) {
202  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
203  {
204  const int qty_index = start_comp + n;
205  const int prim_index = qty_index - 1;
206  // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
207  const int prim_scal_index = prim_index;
208 
209  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
210  BCVars::RhoScalar_bc_comp : qty_index;
211  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
212 
213  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
214  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
215  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
216  ext_dir_on_xlo &= (i == dom_lo.x);
217 
218  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
219  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
220  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
221  ext_dir_on_xlo &= (i == dom_hi.x+1);
222 
223  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
224  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
225  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index])
226  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) );
227 
228  if (ext_dir_on_xlo) {
229  xflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
230  + 3. * cell_prim(i , j, k, prim_index)
231  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
232  } else if (ext_dir_on_xhi) {
233  xflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
234  - 3. * cell_prim(i-1, j, k, prim_index)
235  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
236  } else {
237  xflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) -
238  cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_u(i,j,0);
239  }
240  });
241  ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
242  {
243  const int qty_index = start_comp + n;
244  const int prim_index = qty_index - 1;
245  // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
246  const int prim_scal_index = prim_index;
247 
248  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
249  BCVars::RhoScalar_bc_comp : qty_index;
250  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
251  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
252  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
253  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
254  ext_dir_on_ylo &= (j == dom_lo.y);
255  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
256  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
257  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
258  ext_dir_on_yhi &= (j == dom_hi.y+1);
259 
260  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
261  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
262  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index])
263  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) );
264 
265  if (ext_dir_on_ylo) {
266  yflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
267  + 3. * cell_prim(i, j , k, prim_index)
268  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
269  } else if (ext_dir_on_yhi) {
270  yflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
271  - 3. * cell_prim(i, j-1, k, prim_index)
272  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
273  } else {
274  yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0);
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  const int prim_scal_index = prim_index;
283 
284  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
285  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
286  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index])
287  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) );
288 
289  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
290  BCVars::RhoScalar_bc_comp : qty_index;
291  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
292 
293  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
294  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
295  && k == dom_lo.z);
296  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
297  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
298  && k == dom_hi.z+1);
299  bool most_on_zlo = ( use_most && exp_most &&
300  (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) &&
301  k == dom_lo.z);
302 
303  if (ext_dir_on_zlo) {
304  zflux(i,j,k,qty_index) = -rhoAlpha * ( -(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) ) * dz_inv;
307  } else if (ext_dir_on_zhi) {
308  zflux(i,j,k,qty_index) = -rhoAlpha * ( (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) ) * dz_inv;
311  } else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
312  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
313  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
314  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
315  } else {
316  zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
317  }
318 
319  if (qty_index == RhoTheta_comp) {
320  if (!most_on_zlo) {
321  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
322  }
323  } else if (qty_index == RhoQ1_comp) {
324  if (!most_on_zlo) {
325  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
326  }
327  } else if (qty_index == RhoQ2_comp) {
328  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
329  }
330  });
331  } else if (l_turb) {
332  // with MolecDiffType::Constant or None
333  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
334  {
335  const int qty_index = start_comp + n;
336  const int prim_index = qty_index - 1;
337 
338  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
339  BCVars::RhoScalar_bc_comp : qty_index;
340  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
341 
342  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
343  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
344  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
345  ext_dir_on_xlo &= (i == dom_lo.x);
346 
347  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
348  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
349  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
350  ext_dir_on_xhi &= (i == dom_hi.x+1);
351 
352  Real rhoAlpha = d_alpha_eff[prim_index];
353  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
354  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
355 
356  if (ext_dir_on_xlo) {
357  xflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
358  + 3. * cell_prim(i , j, k, prim_index)
359  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
360  } else if (ext_dir_on_xhi) {
361  xflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
362  - 3. * cell_prim(i-1, j, k, prim_index)
363  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
364  } else {
365  xflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_u(i,j,0);
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  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
374  BCVars::RhoScalar_bc_comp : qty_index;
375  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
376 
377  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
378  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
379  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
380  ext_dir_on_ylo &= (j == dom_lo.y);
381 
382  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
383  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
384  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
385  ext_dir_on_yhi &= (j == dom_hi.y+1);
386 
387  Real rhoAlpha = d_alpha_eff[prim_index];
388  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
389  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
390 
391  if (ext_dir_on_ylo) {
392  yflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
393  + 3. * cell_prim(i, j , k, prim_index)
394  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
395  } else if (ext_dir_on_yhi) {
396  yflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
397  - 3. * cell_prim(i, j-1, k, prim_index)
398  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
399  } else {
400  yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0);
401  }
402  });
403  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
404  {
405  const int qty_index = start_comp + n;
406  const int prim_index = qty_index - 1;
407 
408  Real rhoAlpha = d_alpha_eff[prim_index];
409  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
410  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
411 
412  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
413  BCVars::RhoScalar_bc_comp : qty_index;
414  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
415  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
416  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
417  && k == dom_lo.z);
418  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
419  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
420  && k == dom_hi.z+1);
421  bool most_on_zlo = ( use_most && exp_most &&
422  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) &&
423  k == dom_lo.z);
424 
425  if (ext_dir_on_zlo) {
426  zflux(i,j,k,qty_index) = -rhoAlpha * ( -(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) ) * dz_inv;
429  } else if (ext_dir_on_zhi) {
430  zflux(i,j,k,qty_index) = -rhoAlpha * ( (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) ) * dz_inv;
433  } else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
434  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
435  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
436  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
437  } else {
438  zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
439  }
440 
441  if (qty_index == RhoTheta_comp) {
442  if (!most_on_zlo) {
443  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
444  }
445  } else if (qty_index == RhoQ1_comp) {
446  if (!most_on_zlo) {
447  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
448  }
449  } else if (qty_index == RhoQ2_comp) {
450  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
451  }
452  });
453  } else if(l_consA) {
454  // without an LES/PBL model
455  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
456  {
457  const int qty_index = start_comp + n;
458  const int prim_index = qty_index - 1;
459 
460  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
461  BCVars::RhoScalar_bc_comp : qty_index;
462  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
463 
464  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
465  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
466  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
467  ext_dir_on_xlo &= (i == dom_lo.x);
468 
469  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
470  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
471  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
472  ext_dir_on_xhi &= (i == dom_hi.x+1);
473 
474  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
475  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
476 
477  if (ext_dir_on_xlo) {
478  xflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
479  + 3. * cell_prim(i , j, k, prim_index)
480  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
481  } else if (ext_dir_on_xhi) {
482  xflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
483  - 3. * cell_prim(i-1, j, k, prim_index)
484  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
485  } else {
486  xflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_u(i,j,0);
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  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
495  BCVars::RhoScalar_bc_comp : qty_index;
496  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
497 
498  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
499  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
500  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
501  ext_dir_on_ylo &= (j == dom_lo.y);
502 
503  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
504  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
505  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
506  ext_dir_on_yhi &= (j == dom_hi.y+1);
507 
508  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
509  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
510 
511  if (ext_dir_on_ylo) {
512  yflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
513  + 3. * cell_prim(i, j , k, prim_index)
514  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
515  } else if (ext_dir_on_yhi) {
516  yflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
517  - 3. * cell_prim(i, j-1, k, prim_index)
518  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
519  } else {
520  yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0);
521  }
522  });
523  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
524  {
525  const int qty_index = start_comp + n;
526  const int prim_index = qty_index - 1;
527 
528  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
529  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
530 
531  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
532  BCVars::RhoScalar_bc_comp : qty_index;
533  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
534  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
535  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
536  && k == dom_lo.z);
537  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
538  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
539  && k == dom_hi.z+1);
540  bool most_on_zlo = ( use_most && exp_most &&
541  (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) &&
542  k == dom_lo.z);
543 
544  if (ext_dir_on_zlo) {
545  zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
546  + 3. * cell_prim(i, j, k , prim_index)
547  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
548  } else if (ext_dir_on_zhi) {
549  zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
550  - 3. * cell_prim(i, j, k-1, prim_index)
551  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
552  } else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
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 * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
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  } else {
573  // with MolecDiffType::Constant or None
574  // without an LES/PBL model
575  ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
576  {
577  const int qty_index = start_comp + n;
578  const int prim_index = qty_index - 1;
579 
580  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
581  BCVars::RhoScalar_bc_comp : qty_index;
582  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
583 
584  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
585  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
586  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
587  ext_dir_on_xlo &= (i == dom_lo.x);
588 
589  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
590  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
591  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
592  ext_dir_on_xhi &= (i == dom_hi.x+1);
593 
594  Real rhoAlpha = d_alpha_eff[prim_index];
595 
596  if (ext_dir_on_xlo) {
597  xflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
598  + 3. * cell_prim(i , j, k, prim_index)
599  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
600  } else if (ext_dir_on_xhi) {
601  xflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
602  - 3. * cell_prim(i-1, j, k, prim_index)
603  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
604  } else {
605  xflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_u(i,j,0);
606  }
607  });
608  ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
609  {
610  const int qty_index = start_comp + n;
611  const int prim_index = qty_index - 1;
612 
613  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
614  BCVars::RhoScalar_bc_comp : qty_index;
615  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
616 
617  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
618  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
619  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
620  ext_dir_on_ylo &= (j == dom_lo.y);
621 
622  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
623  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
624  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
625  ext_dir_on_yhi &= (j == dom_hi.y+1);
626 
627  Real rhoAlpha = d_alpha_eff[prim_index];
628 
629  if (ext_dir_on_ylo) {
630  yflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
631  + 3. * cell_prim(i, j , k, prim_index)
632  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
633  } else if (ext_dir_on_yhi) {
634  yflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
635  - 3. * cell_prim(i, j-1, k, prim_index)
636  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
637  } else {
638  yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0);
639  }
640  });
641  ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
642  {
643  const int qty_index = start_comp + n;
644  const int prim_index = qty_index - 1;
645 
646  Real rhoAlpha = d_alpha_eff[prim_index];
647 
648  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
649  BCVars::RhoScalar_bc_comp : qty_index;
650  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
651  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
652  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
653  && k == dom_lo.z);
654  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
655  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
656  && k == dom_hi.z+1);
657  bool most_on_zlo = ( use_most && exp_most &&
658  (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) &&
659  k == dom_lo.z);
660 
661  if (ext_dir_on_zlo) {
662  zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
663  + 3. * cell_prim(i, j, k , prim_index)
664  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
665  } else if (ext_dir_on_zhi) {
666  zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
667  - 3. * cell_prim(i, j, k-1, prim_index)
668  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
669  } else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
670  zflux(i,j,k,qty_index) = hfx_z(i,j,0);
671  } else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
672  zflux(i,j,k,qty_index) = qfx1_z(i,j,0);
673  } else {
674  zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
675  }
676 
677  if (qty_index == RhoTheta_comp) {
678  if (!most_on_zlo) {
679  hfx_z(i,j,k) = zflux(i,j,k,qty_index);
680  }
681  } else if (qty_index == RhoQ1_comp) {
682  if (!most_on_zlo) {
683  qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
684  }
685  } else if (qty_index == RhoQ2_comp) {
686  qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
687  }
688  });
689  }
690 
691  // Use fluxes to compute RHS
692  for (int n(0); n < num_comp; ++n)
693  {
694  int qty_index = start_comp + n;
695  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
696  {
697 
698  cell_rhs(i,j,k,qty_index) -= (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
699  +(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
700  +(zflux(i ,j ,k+1,qty_index) - zflux(i, j, k, qty_index)) * dz_inv; // Diffusive flux in z-dir
701  });
702  }
703 
704  // Using Deardorff (see Sullivan et al 1994)
705  // or k-eqn RANS (see Axell & Liungman 2001)
706  //
707  // Note: At this point, the thermal diffusivity ("Khv" field in ERF), the
708  // subgrid heat flux ("hfx_z" here), and the subgrid dissipation
709  // ("diss" here) have been updated by ComputeTurbulentViscosityLES --
710  // at the beginning of each timestep.
711  // The strain rate magnitude is updated at the beginning of the first
712  // RK stage only, therefore the shear production term also does not
713  // change between RK stages.
714  // The surface heat flux hfx_z(i,j,-1) is updated in MOSTStress at
715  // each RK stage if using the ERF_EXPLICIT_MOST_STRESS path, but that
716  // does not change the buoyancy production term here.
717  if (l_use_keqn && (start_comp <= RhoKE_comp) && (end_comp >= RhoKE_comp)) {
718  int qty_index = RhoKE_comp;
719  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
720  {
721  // Add Buoyancy Source
722  // where the SGS buoyancy flux tau_{theta,i} = -KH * dtheta/dx_i,
723  // such that for dtheta/dz < 0, there is a positive (upward) heat
724  // flux; the TKE buoyancy production is then
725  // B = g/theta_0 * tau_{theta,w}
726  // for a dry atmosphere.
727  // TODO: To account for moisture, the Brunt-Vaisala frequency,
728  // N^2 = g[1/theta * dtheta/dz + ...]
729  // **should** be a function of the water vapor and total water
730  // mixing ratios, depending on whether conditions are saturated or
731  // not (see the WRF model description, Skamarock et al 2019).
732  cell_rhs(i,j,k,qty_index) += l_abs_g * l_inv_theta0 * hfx_z(i,j,k);
733 
734  // TKE shear production
735  // P = -tau_ij * S_ij = 2 * mu_turb * S_ij * S_ij
736  // Note: This assumes that the horizontal and vertical diffusivities
737  // of momentum are equal
738  cell_rhs(i,j,k,qty_index) += 2.0*mu_turb(i,j,k,EddyDiff::Mom_v) * SmnSmn_a(i,j,k);
739 
740  // TKE dissipation
741  cell_rhs(i,j,k,qty_index) -= diss(i,j,k);
742  });
743  }
744 
745  // Using PBL
746  if (l_use_mynn && start_comp <= RhoKE_comp && end_comp >=RhoKE_comp) {
747  int qty_index = RhoKE_comp;
748  auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;
749 
750  const int rhoqv_comp = solverChoice.RhoQv_comp;
751  const int rhoqc_comp = solverChoice.RhoQc_comp;
752  const int rhoqr_comp = solverChoice.RhoQr_comp;
753 
754  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
755  {
756  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
757  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
758  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
759  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) );
760  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
761  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
762 
763  // This computes shear production, buoyancy production, and dissipation terms only.
764  cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
765  mu_turb,cellSizeInv,domain,
766  pbl_mynn_B1_l,tm_arr(i,j,0),
767  rhoqv_comp, rhoqc_comp, rhoqr_comp,
768  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
769  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
770  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
771  use_most);
772  });
773  }
774 }
MolecDiffType::ConstantAlpha
@ ConstantAlpha
DiffusionSrcForState_N
void DiffusionSrcForState_N(const Box &bx, const Box &domain, int start_comp, int num_comp, const bool &exp_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 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_z, 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_N.cpp:40
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
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
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:199
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: