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 "ERF_SetupDiff.H"
#include "ERF_AddTKESources.H"
#include "ERF_AddQKESources.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 &rotate, 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 > &z_cc, const Array4< const Real > &ax, const Array4< const Real > &ay, const Array4< const Real > &, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &SmnSmn_a, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_ux, const Array4< const Real > &mf_vx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vy, 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_SurfLayer, const Vector< std::unique_ptr< SurfaceLayer >> &SurfLayer, const Real implicit_fac)
 

Function Documentation

◆ DiffusionSrcForState_T()

void DiffusionSrcForState_T ( const Box &  bx,
const Box &  domain,
int  start_comp,
int  num_comp,
const bool &  rotate,
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 > &  z_cc,
const Array4< const Real > &  ax,
const Array4< const Real > &  ay,
const Array4< const Real > &  ,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
const Array4< const Real > &  SmnSmn_a,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_ux,
const Array4< const Real > &  mf_vx,
const Array4< const Real > &  mf_my,
const Array4< const Real > &  mf_uy,
const Array4< const Real > &  mf_vy,
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_SurfLayer,
const Vector< std::unique_ptr< SurfaceLayer >> &  SurfLayer,
const Real  implicit_fac 
)

Function for computing the scalar RHS for diffusion operator with terrain-fitted coordinates.

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]rotateflag to rotate terrain-aligned fluxes
[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]z_cccell-centered physical z height
[in]axarea fraction of x-faces
[in]ayarea fraction of y-faces
[in]detJJacobian determinant
[in]cellSizeInvinverse cell size array
[in]SmnSmn_astrain rate magnitude
[in]mf_mxx map factor at cell centers
[in]mf_uxx map factor at x-faces
[in]mf_vxx map factor at y-faces
[in]mf_myy map factor at cell centers
[in]mf_uyy map factor at x-faces
[in]mf_vyy map factor at y-faces
[in,out]hfx_xheat flux in x-dir
[in,out]hfx_yheat flux in y-dir
[in,out]hfx_zheat flux in z-dir
[in,out]qfx1_xheat flux in x-dir
[in,out]qfx1_yheat flux in y-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]solverChoicecontainer of solver and diffusion parameters
[in]levelAMR level
[in]tm_arrtheta mean array
[in]grav_gpugravity vector
[in]bc_ptrcontainer with boundary conditions
[in]use_SurfLayerwhether we have turned on subgrid diffusion
[in]implicit_fac– factor of implicitness for vertical differences only
97 {
98  BL_PROFILE_VAR("DiffusionSrcForState_T()",DiffusionSrcForState_T);
99 
100  const Real explicit_fac = one - implicit_fac;
101 
102 #include "ERF_SetupDiff.H"
103  Real l_abs_g = std::abs(grav_gpu[2]);
104 
105  const Real dz_inv = cellSizeInv[2];
106 
107  // We need to grow these boxes in the vertical direction when tiling so that we can access xflux and yflux
108  // to modify zflux
109  Box xbx_g1(xbx); Box ybx_g1(ybx);
110  if (xbx_g1.smallEnd(2) != dom_lo.z) xbx_g1.growLo(2,1);
111  if (ybx_g1.smallEnd(2) != dom_lo.z) ybx_g1.growLo(2,1);
112  if (xbx_g1.bigEnd(2) != dom_hi.z) xbx_g1.growHi(2,1);
113  if (ybx_g1.bigEnd(2) != dom_hi.z) ybx_g1.growHi(2,1);
114 
115  for (int n(0); n<num_comp; ++n) {
116  const int qty_index = start_comp + n;
117 
118  // Constant alpha & Turb model
119  if (l_consA && l_turb) {
120  ParallelFor(xbx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
121  {
122  const int prim_index = qty_index - 1;
123  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
124 
125  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
126  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
127  rhoAlpha += myhalf * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index])
128  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) );
129 
130  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
131 
132  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
133  BCVars::RhoScalar_bc_comp : qty_index;
134  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
135 
136  bool SurfLayer_on_xlo = ( SurfLayer_xlo && i == dom_lo.x);
137  bool SurfLayer_on_xhi = ( SurfLayer_xhi && i == dom_hi.x + 1);
138  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
139 
140  Real idz_hi = one / (z_cc(i ,j,k+1) - z_cc(i ,j,k-1));
141  Real idz_lo = one / (z_cc(i-1,j,k+1) - z_cc(i-1,j,k-1));
142  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i-1, j, k+1, prim_index)*idz_lo
143  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i-1, j, k-1, prim_index)*idz_lo );
144  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
145 
146  if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoTheta_comp)) {
147  xflux(i,j,k) = hfx_x(i,j,k);
148  } else if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoQ1_comp)) {
149  xflux(i,j,k) = qfx1_x(i,j,k);
150  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
151  xflux(i,j,k) = hfx_x(i,j,0);
152  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
153  xflux(i,j,k) = qfx1_x(i,j,0);
154  } else {
155  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * ( GradCx - met_h_xi*GradCz );
156  }
157 
158  });
159  ParallelFor(ybx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
160  {
161  const int prim_index = qty_index - 1;
162  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
163 
164  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
165  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
166  rhoAlpha += myhalf * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index])
167  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) );
168 
169  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
170 
171  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
172  BCVars::RhoScalar_bc_comp : qty_index;
173  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
174  bool SurfLayer_on_ylo = ( SurfLayer_ylo && j == dom_lo.y);
175  bool SurfLayer_on_yhi = ( SurfLayer_yhi && j == dom_hi.y + 1);
176  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
177 
178  Real idz_hi = one / (z_cc(i,j ,k+1) - z_cc(i,j ,k-1));
179  Real idz_lo = one / (z_cc(i,j-1,k+1) - z_cc(i,j-1,k-1));
180  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i, j-1, k+1, prim_index)*idz_lo
181  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i, j-1, k-1, prim_index)*idz_lo );
182  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
183 
184  if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoTheta_comp)) {
185  yflux(i,j,k) = hfx_y(i,j,k);
186  } else if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoQ1_comp)) {
187  yflux(i,j,k) = qfx1_y(i,j,k);
188  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
189  yflux(i,j,k) = hfx_y(i,j,0);
190  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
191  yflux(i,j,k) = qfx1_y(i,j,0);
192  } else {
193  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * ( GradCy - met_h_eta*GradCz );
194  }
195 
196  });
197  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
198  {
199  const int prim_index = qty_index - 1;
200  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
201 
202  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
203  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
204  rhoAlpha += myhalf * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index])
205  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) );
206 
207  Real GradCz;
208  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
209  BCVars::RhoScalar_bc_comp : qty_index;
210  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
211  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
212  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
213  && k == dom_lo.z);
214  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
215  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim) )
216  && k == dom_hi.z+1);
217  bool SurfLayer_on_zlo = ( SurfLayer_zlo && k == dom_lo.z);
218  bool SurfLayer_on_zhi = ( SurfLayer_zhi && k == dom_hi.z + 1);
219 
220  if (ext_dir_on_zlo) {
221  // Third order stencil with variable dz
222  Real zm = Compute_Z_AtWFace(i,j,k+1,z_nd);
223  Real dz0 = zm - Compute_Z_AtWFace(i,j,k,z_nd);
224  Real dz1 = Compute_Z_AtWFace(i,j,k+2,z_nd) - zm;
225  Real idz0 = one / dz0;
226  Real f = (dz1 / dz0) + two;
227  Real f2 = f*f;
228  Real c3 = two / (f - f2);
229  Real c2 = -f2*c3;
230  Real c1 = -(one-f2)*c3;
231  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
232  + c2 * cell_prim(i, j, k , prim_index)
233  + c3 * cell_prim(i, j, k+1, prim_index) );
234  } else if (ext_dir_on_zhi) {
235  // Third order stencil with variable dz
236  Real zm = Compute_Z_AtWFace(i,j,k-1,z_nd);
237  Real dz0 = Compute_Z_AtWFace(i,j,k,z_nd) - zm;
238  Real dz1 = zm - Compute_Z_AtWFace(i,j,k-2,z_nd);
239  Real idz0 = one / dz0;
240  Real f = (dz1 / dz0) + two;
241  Real f2 = f*f;
242  Real c3 = two / (f - f2);
243  Real c2 = -f2*c3;
244  Real c1 = -(one-f2)*c3;
245  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
246  + c2 * cell_prim(i, j, k-1, prim_index)
247  + c3 * cell_prim(i, j, k-2, prim_index) ) );
248  } else {
249  Real met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd);
250  GradCz = (dz_inv/met_h_zeta) * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
251  }
252 
253  if (SurfLayer_on_zlo || SurfLayer_on_zhi) {
254  if (qty_index == RhoTheta_comp) {
255  zflux(i,j,k) = hfx_z(i,j,k);
256  } else if (qty_index == RhoQ1_comp) {
257  zflux(i,j,k) = qfx1_z(i,j,k);
258  } else {
259  zflux(i,j,k) = zero;
260  }
261  } else {
262  zflux(i,j,k) = -rhoAlpha * GradCz;
263  }
264 
265  if (qty_index == RhoTheta_comp) {
266  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
267  hfx_z(i,j,k) = zflux(i,j,k);
268  }
269  } else if (qty_index == RhoQ1_comp) {
270  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
271  qfx1_z(i,j,k) = zflux(i,j,k);
272  }
273  } else if (qty_index == RhoQ2_comp) {
274  qfx2_z(i,j,k) = zflux(i,j,k);
275  }
276  });
277  // Constant rho*alpha & Turb model
278  } else if (l_turb) {
279  ParallelFor(xbx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
280  {
281  const int prim_index = qty_index - 1;
282 
283  Real rhoAlpha = d_alpha_eff[prim_index];
284  rhoAlpha += myhalf * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
285  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
286 
287  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
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  bool SurfLayer_on_xlo = ( SurfLayer_xlo && i == dom_lo.x);
293  bool SurfLayer_on_xhi = ( SurfLayer_xhi && i == dom_hi.x + 1);
294  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
295 
296  Real idz_hi = one / (z_cc(i ,j,k+1) - z_cc(i ,j,k-1));
297  Real idz_lo = one / (z_cc(i-1,j,k+1) - z_cc(i-1,j,k-1));
298  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i-1, j, k+1, prim_index)*idz_lo
299  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i-1, j, k-1, prim_index)*idz_lo );
300  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
301 
302  if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoTheta_comp)) {
303  xflux(i,j,k) = hfx_x(i,j,k);
304  } else if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoQ1_comp)) {
305  xflux(i,j,k) = qfx1_x(i,j,k);
306  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
307  xflux(i,j,k) = hfx_x(i,j,0);
308  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
309  xflux(i,j,k) = qfx1_x(i,j,0);
310  } else {
311  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * ( GradCx - met_h_xi*GradCz );
312  }
313 
314  });
315  ParallelFor(ybx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
316  {
317  const int prim_index = qty_index - 1;
318 
319  Real rhoAlpha = d_alpha_eff[prim_index];
320  rhoAlpha += myhalf * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
321  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
322 
323  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
324 
325  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
326  BCVars::RhoScalar_bc_comp : qty_index;
327  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
328  bool SurfLayer_on_ylo = ( SurfLayer_ylo && j == dom_lo.y);
329  bool SurfLayer_on_yhi = ( SurfLayer_yhi && j == dom_hi.y + 1);
330  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
331 
332  Real idz_hi = one / (z_cc(i,j ,k+1) - z_cc(i,j ,k-1));
333  Real idz_lo = one / (z_cc(i,j-1,k+1) - z_cc(i,j-1,k-1));
334  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i, j-1, k+1, prim_index)*idz_lo
335  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i, j-1, k-1, prim_index)*idz_lo );
336  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
337 
338  if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoTheta_comp)) {
339  yflux(i,j,k) = hfx_y(i,j,k);
340  } else if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoQ1_comp)) {
341  yflux(i,j,k) = qfx1_y(i,j,k);
342  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
343  yflux(i,j,k) = hfx_y(i,j,0);
344  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
345  yflux(i,j,k) = qfx1_y(i,j,0);
346  } else {
347  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * ( GradCy - met_h_eta*GradCz );
348  }
349 
350  });
351  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
352  {
353  const int prim_index = qty_index - 1;
354 
355  Real rhoAlpha = d_alpha_eff[prim_index];
356  rhoAlpha += myhalf * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
357  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
358 
359  Real GradCz;
360  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
361  BCVars::RhoScalar_bc_comp : qty_index;
362  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
363  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
364  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
365  && k == dom_lo.z);
366  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
367  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
368  && k == dom_hi.z+1);
369  bool SurfLayer_on_zlo = ( SurfLayer_zlo && k == dom_lo.z);
370  bool SurfLayer_on_zhi = ( SurfLayer_zhi && k == dom_hi.z + 1);
371 
372  if (ext_dir_on_zlo) {
373  // Third order stencil with variable dz
374  Real zm = Compute_Z_AtWFace(i,j,k+1,z_nd);
375  Real dz0 = zm - Compute_Z_AtWFace(i,j,k,z_nd);
376  Real dz1 = Compute_Z_AtWFace(i,j,k+2,z_nd) - zm;
377  Real idz0 = one / dz0;
378  Real f = (dz1 / dz0) + two;
379  Real f2 = f*f;
380  Real c3 = two / (f - f2);
381  Real c2 = -f2*c3;
382  Real c1 = -(one-f2)*c3;
383  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
384  + c2 * cell_prim(i, j, k , prim_index)
385  + c3 * cell_prim(i, j, k+1, prim_index) );
386  } else if (ext_dir_on_zhi) {
387  // Third order stencil with variable dz
388  Real zm = Compute_Z_AtWFace(i,j,k-1,z_nd);
389  Real dz0 = Compute_Z_AtWFace(i,j,k,z_nd) - zm;
390  Real dz1 = zm - Compute_Z_AtWFace(i,j,k-2,z_nd);
391  Real idz0 = one / dz0;
392  Real f = (dz1 / dz0) + two;
393  Real f2 = f*f;
394  Real c3 = two / (f - f2);
395  Real c2 = -f2*c3;
396  Real c1 = -(one-f2)*c3;
397  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
398  + c2 * cell_prim(i, j, k-1, prim_index)
399  + c3 * cell_prim(i, j, k-2, prim_index) ) );
400  } else {
401  Real met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd);
402  GradCz = (dz_inv/met_h_zeta) * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
403  }
404 
405  if (SurfLayer_on_zlo || SurfLayer_on_zhi) {
406  if (qty_index == RhoTheta_comp) {
407  zflux(i,j,k) = hfx_z(i,j,k);
408  } else if (qty_index == RhoQ1_comp) {
409  zflux(i,j,k) = qfx1_z(i,j,k);
410  } else {
411  zflux(i,j,k) = zero;
412  }
413  } else {
414  zflux(i,j,k) = -rhoAlpha * GradCz;
415  }
416 
417  if (qty_index == RhoTheta_comp) {
418  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
419  hfx_z(i,j,k) = zflux(i,j,k);
420  }
421  } else if (qty_index == RhoQ1_comp) {
422  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
423  qfx1_z(i,j,k) = zflux(i,j,k);
424  }
425  } else if (qty_index == RhoQ2_comp) {
426  qfx2_z(i,j,k) = zflux(i,j,k);
427  }
428  });
429  // Constant alpha & no LES/PBL model
430  } else if(l_consA) {
431  ParallelFor(xbx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
432  {
433  const int prim_index = qty_index - 1;
434 
435  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
436  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
437 
438  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
439 
440  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
441  BCVars::RhoScalar_bc_comp : qty_index;
442  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
443  bool SurfLayer_on_xlo = ( SurfLayer_xlo && i == dom_lo.x);
444  bool SurfLayer_on_xhi = ( SurfLayer_xhi && i == dom_hi.x + 1);
445  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
446 
447  Real idz_hi = one / (z_cc(i ,j,k+1) - z_cc(i ,j,k-1));
448  Real idz_lo = one / (z_cc(i-1,j,k+1) - z_cc(i-1,j,k-1));
449  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i-1, j, k+1, prim_index)*idz_lo
450  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i-1, j, k-1, prim_index)*idz_lo );
451  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
452 
453  if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoTheta_comp)) {
454  xflux(i,j,k) = hfx_x(i,j,k);
455  } else if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoQ1_comp)) {
456  xflux(i,j,k) = qfx1_x(i,j,k);
457  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
458  xflux(i,j,k) = hfx_x(i,j,0);
459  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
460  xflux(i,j,k) = qfx1_x(i,j,0);
461  } else {
462  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * ( GradCx - met_h_xi*GradCz );
463  }
464 
465  });
466  ParallelFor(ybx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
467  {
468  const int prim_index = qty_index - 1;
469 
470  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
471  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
472 
473  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
474 
475  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
476  BCVars::RhoScalar_bc_comp : qty_index;
477  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
478  bool SurfLayer_on_ylo = ( SurfLayer_ylo && j == dom_lo.y);
479  bool SurfLayer_on_yhi = ( SurfLayer_yhi && j == dom_hi.y + 1);
480  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
481 
482  Real idz_hi = one / (z_cc(i,j ,k+1) - z_cc(i,j ,k-1));
483  Real idz_lo = one / (z_cc(i,j-1,k+1) - z_cc(i,j-1,k-1));
484  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i, j-1, k+1, prim_index)*idz_lo
485  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i, j-1, k-1, prim_index)*idz_lo );
486  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
487 
488  if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoTheta_comp)) {
489  yflux(i,j,k) = hfx_y(i,j,k);
490  } else if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoQ1_comp)) {
491  yflux(i,j,k) = qfx1_y(i,j,k);
492  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
493  yflux(i,j,k) = hfx_y(i,j,0);
494  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
495  yflux(i,j,k) = qfx1_y(i,j,0);
496  } else {
497  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * ( GradCy - met_h_eta*GradCz );
498  }
499 
500  });
501  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
502  {
503  const int prim_index = qty_index - 1;
504 
505  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
506  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
507 
508  Real GradCz;
509  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
510  BCVars::RhoScalar_bc_comp : qty_index;
511  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
512  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
513  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
514  && k == dom_lo.z);
515  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
516  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
517  && k == dom_hi.z+1);
518  bool SurfLayer_on_zlo = ( SurfLayer_zlo && k == dom_lo.z);
519  bool SurfLayer_on_zhi = ( SurfLayer_zhi && k == dom_hi.z + 1);
520 
521  if (ext_dir_on_zlo) {
522  // Third order stencil with variable dz
523  Real zm = Compute_Z_AtWFace(i,j,k+1,z_nd);
524  Real dz0 = zm - Compute_Z_AtWFace(i,j,k,z_nd);
525  Real dz1 = Compute_Z_AtWFace(i,j,k+2,z_nd) - zm;
526  Real idz0 = one / dz0;
527  Real f = (dz1 / dz0) + two;
528  Real f2 = f*f;
529  Real c3 = two / (f - f2);
530  Real c2 = -f2*c3;
531  Real c1 = -(one-f2)*c3;
532  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
533  + c2 * cell_prim(i, j, k , prim_index)
534  + c3 * cell_prim(i, j, k+1, prim_index) );
535  } else if (ext_dir_on_zhi) {
536  // Third order stencil with variable dz
537  Real zm = Compute_Z_AtWFace(i,j,k-1,z_nd);
538  Real dz0 = Compute_Z_AtWFace(i,j,k,z_nd) - zm;
539  Real dz1 = zm - Compute_Z_AtWFace(i,j,k-2,z_nd);
540  Real idz0 = one / dz0;
541  Real f = (dz1 / dz0) + two;
542  Real f2 = f*f;
543  Real c3 = two / (f - f2);
544  Real c2 = -f2*c3;
545  Real c1 = -(one-f2)*c3;
546  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
547  + c2 * cell_prim(i, j, k-1, prim_index)
548  + c3 * cell_prim(i, j, k-2, prim_index) ) );
549  } else {
550  Real met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd);
551  GradCz = (dz_inv/met_h_zeta) * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
552  }
553 
554  if (SurfLayer_on_zlo || SurfLayer_on_zhi) {
555  if (qty_index == RhoTheta_comp) {
556  zflux(i,j,k) = hfx_z(i,j,k);
557  } else if (qty_index == RhoQ1_comp) {
558  zflux(i,j,k) = qfx1_z(i,j,k);
559  } else {
560  zflux(i,j,k) = zero;
561  }
562  } else {
563  zflux(i,j,k) = -rhoAlpha * GradCz;
564  }
565 
566  if (qty_index == RhoTheta_comp) {
567  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
568  hfx_z(i,j,k) = zflux(i,j,k);
569  }
570  } else if (qty_index == RhoQ1_comp) {
571  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
572  qfx1_z(i,j,k) = zflux(i,j,k);
573  }
574  } else if (qty_index == RhoQ2_comp) {
575  qfx2_z(i,j,k) = zflux(i,j,k);
576  }
577  });
578  // Constant rho*alpha & no LES/PBL model
579  } else {
580  ParallelFor(xbx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
581  {
582  const int prim_index = qty_index - 1;
583 
584  Real rhoAlpha = d_alpha_eff[prim_index];
585 
586  Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd);
587 
588  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
589  BCVars::RhoScalar_bc_comp : qty_index;
590  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
591  bool SurfLayer_on_xlo = ( SurfLayer_xlo && i == dom_lo.x);
592  bool SurfLayer_on_xhi = ( SurfLayer_xhi && i == dom_hi.x + 1);
593  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
594 
595  Real idz_hi = one / (z_cc(i ,j,k+1) - z_cc(i ,j,k-1));
596  Real idz_lo = one / (z_cc(i-1,j,k+1) - z_cc(i-1,j,k-1));
597  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i-1, j, k+1, prim_index)*idz_lo
598  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i-1, j, k-1, prim_index)*idz_lo );
599  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
600 
601  if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoTheta_comp)) {
602  xflux(i,j,k) = hfx_x(i,j,k);
603  } else if ((SurfLayer_on_xlo || SurfLayer_on_xhi) && (qty_index == RhoQ1_comp)) {
604  xflux(i,j,k) = qfx1_x(i,j,k);
605  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
606  xflux(i,j,k) = hfx_x(i,j,0);
607  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
608  xflux(i,j,k) = qfx1_x(i,j,0);
609  } else {
610  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * ( GradCx - met_h_xi*GradCz );
611  }
612 
613  });
614  ParallelFor(ybx_g1, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
615  {
616  const int prim_index = qty_index - 1;
617 
618  Real rhoAlpha = d_alpha_eff[prim_index];
619 
620  Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd);
621 
622  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
623  BCVars::RhoScalar_bc_comp : qty_index;
624  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
625  bool SurfLayer_on_ylo = ( SurfLayer_ylo && j == dom_lo.y);
626  bool SurfLayer_on_yhi = ( SurfLayer_yhi && j == dom_hi.y + 1);
627  bool SurfLayer_on_zlo = ( SurfLayer_zlo && rotate && k == dom_lo.z);
628 
629  Real idz_hi = one / (z_cc(i,j ,k+1) - z_cc(i,j ,k-1));
630  Real idz_lo = one / (z_cc(i,j-1,k+1) - z_cc(i,j-1,k-1));
631  Real GradCz = myhalf * ( cell_prim(i, j, k+1, prim_index)*idz_hi + cell_prim(i, j-1, k+1, prim_index)*idz_lo
632  - cell_prim(i, j, k-1, prim_index)*idz_hi - cell_prim(i, j-1, k-1, prim_index)*idz_lo );
633  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
634 
635  if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoTheta_comp)) {
636  yflux(i,j,k) = hfx_y(i,j,k);
637  } else if ((SurfLayer_on_ylo || SurfLayer_on_yhi) && (qty_index == RhoQ1_comp)) {
638  yflux(i,j,k) = qfx1_y(i,j,k);
639  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
640  yflux(i,j,k) = hfx_y(i,j,0);
641  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
642  yflux(i,j,k) = qfx1_y(i,j,0);
643  } else {
644  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * ( GradCy - met_h_eta*GradCz );
645  }
646 
647  });
648  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
649  {
650  const int prim_index = qty_index - 1;
651 
652  Real rhoAlpha = d_alpha_eff[prim_index];
653 
654 
655  Real GradCz;
656  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
657  BCVars::RhoScalar_bc_comp : qty_index;
658  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
659  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
660  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
661  && k == dom_lo.z);
662  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
663  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
664  && k == dom_hi.z+1);
665  bool SurfLayer_on_zlo = ( SurfLayer_zlo && k == dom_lo.z);
666  bool SurfLayer_on_zhi = ( SurfLayer_zhi && k == dom_hi.z + 1);
667 
668  if (ext_dir_on_zlo) {
669  // Third order stencil with variable dz
670  Real zm = Compute_Z_AtWFace(i,j,k+1,z_nd);
671  Real dz0 = zm - Compute_Z_AtWFace(i,j,k,z_nd);
672  Real dz1 = Compute_Z_AtWFace(i,j,k+2,z_nd) - zm;
673  Real idz0 = one / dz0;
674  Real f = (dz1 / dz0) + two;
675  Real f2 = f*f;
676  Real c3 = two / (f - f2);
677  Real c2 = -f2*c3;
678  Real c1 = -(one-f2)*c3;
679  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
680  + c2 * cell_prim(i, j, k , prim_index)
681  + c3 * cell_prim(i, j, k+1, prim_index) );
682  } else if (ext_dir_on_zhi) {
683  // Third order stencil with variable dz
684  Real zm = Compute_Z_AtWFace(i,j,k-1,z_nd);
685  Real dz0 = Compute_Z_AtWFace(i,j,k,z_nd) - zm;
686  Real dz1 = zm - Compute_Z_AtWFace(i,j,k-2,z_nd);
687  Real idz0 = one / dz0;
688  Real f = (dz1 / dz0) + two;
689  Real f2 = f*f;
690  Real c3 = two / (f - f2);
691  Real c2 = -f2*c3;
692  Real c1 = -(one-f2)*c3;
693  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
694  + c2 * cell_prim(i, j, k-1, prim_index)
695  + c3 * cell_prim(i, j, k-2, prim_index) ) );
696  } else {
697  Real met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd);
698  GradCz = (dz_inv/met_h_zeta) * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
699  }
700 
701  if (SurfLayer_on_zlo || SurfLayer_on_zhi) {
702  if (qty_index == RhoTheta_comp) {
703  zflux(i,j,k) = hfx_z(i,j,k);
704  } else if (qty_index == RhoQ1_comp) {
705  zflux(i,j,k) = qfx1_z(i,j,k);
706  } else {
707  zflux(i,j,k) = zero;
708  }
709  } else {
710  zflux(i,j,k) = -rhoAlpha * GradCz;
711  }
712 
713  if (qty_index == RhoTheta_comp) {
714  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
715  hfx_z(i,j,k) = zflux(i,j,k);
716  }
717  } else if (qty_index == RhoQ1_comp) {
718  if (!(SurfLayer_on_zlo || SurfLayer_on_zhi)) {
719  qfx1_z(i,j,k) = zflux(i,j,k);
720  }
721  } else if (qty_index == RhoQ2_comp) {
722  qfx2_z(i,j,k) = zflux(i,j,k);
723  }
724  });
725  }
726 
727  // NOTE: With terrain, we implicitly treat the leading order vertical gradient (no metric terms)
728  // This allows us to do semi-implicit discretization of the vertical diffusive terms
729  if (qty_index == RhoTheta_comp ||
730  qty_index == RhoKE_comp ||
731  qty_index == RhoQ1_comp) {
732  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
733  {
734  zflux(i,j,k) *= explicit_fac;
735  });
736  }
737 
738  //-----------------------------------------------------------------------------------
739  //
740  // Modify fluxes by terrain and use fluxes to compute RHS
741  //
742  // Note that we combine all of these operations in order to keep this section
743  // of the loop tiling-safe.
744  //-----------------------------------------------------------------------------------
745  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
746  {
747  Real xfluxbar_lo, yfluxbar_lo;
748  Real met_h_xi_lo = Compute_h_xi_AtKface (i,j,k ,cellSizeInv,z_nd);
749  Real met_h_eta_lo = Compute_h_eta_AtKface(i,j,k ,cellSizeInv,z_nd);
750  if (k == dom_lo.z) {
751  Real xfluxlo = myhalf * ( xflux(i,j,k ) + xflux(i+1,j,k ) );
752  Real xfluxhi = myhalf * ( xflux(i,j,k+1) + xflux(i+1,j,k+1) );
753  xfluxbar_lo = Real(1.5)*xfluxlo - myhalf*xfluxhi;
754 
755  Real yfluxlo = myhalf * ( yflux(i,j,k ) + yflux(i,j+1,k ) );
756  Real yfluxhi = myhalf * ( yflux(i,j,k+1) + yflux(i,j+1,k+1) );
757  yfluxbar_lo = Real(1.5)*yfluxlo - myhalf*yfluxhi;
758  } else {
759  xfluxbar_lo = fourth * ( xflux(i,j,k ) + xflux(i+1,j ,k )
760  + xflux(i,j,k-1) + xflux(i+1,j ,k-1) );
761  yfluxbar_lo = fourth * ( yflux(i,j,k ) + yflux(i ,j+1,k )
762  + yflux(i,j,k-1) + yflux(i ,j+1,k-1) );
763  }
764 
765  Real xfluxbar_hi, yfluxbar_hi;
766  Real met_h_xi_hi = Compute_h_xi_AtKface (i,j,k+1,cellSizeInv,z_nd);
767  Real met_h_eta_hi = Compute_h_eta_AtKface(i,j,k+1,cellSizeInv,z_nd);
768  if (k == dom_hi.z) {
769  Real xfluxlo = myhalf * ( xflux(i,j,k-1) + xflux(i+1,j,k-1) );
770  Real xfluxhi = myhalf * ( xflux(i,j,k ) + xflux(i+1,j,k ) );
771  xfluxbar_hi = Real(1.5)*xfluxhi - myhalf*xfluxlo;
772 
773  Real yfluxlo = myhalf * ( yflux(i,j,k-1) + yflux(i,j+1,k-1) );
774  Real yfluxhi = myhalf * ( yflux(i,j,k ) + yflux(i,j+1,k ) );
775  yfluxbar_hi = Real(1.5)*yfluxhi - myhalf*yfluxlo;
776  } else {
777  xfluxbar_hi = fourth * ( xflux(i,j,k+1) + xflux(i+1,j ,k+1)
778  + xflux(i,j,k ) + xflux(i+1,j ,k ) );
779  yfluxbar_hi = fourth * ( yflux(i,j,k+1) + yflux(i ,j+1,k+1)
780  + yflux(i,j,k ) + yflux(i ,j+1,k ) );
781  }
782 
783  // Allow semi-implicit discretization of the vertical diffusive terms
784  Real zflux_lo;
785  if ( SurfLayer_zlo &&
786  k == dom_lo.z &&
787  qty_index != RhoTheta_comp &&
788  qty_index != RhoQ1_comp ) {
789  zflux_lo = zero;
790  } else {
791  zflux_lo = zflux(i,j,k )
792  - met_h_xi_lo * mf_mx(i,j,0) * xfluxbar_lo
793  - met_h_eta_lo * mf_my(i,j,0) * yfluxbar_lo;
794  }
795  Real zflux_hi = zflux(i,j,k+1)
796  - met_h_xi_hi * mf_mx(i,j,0) * xfluxbar_hi
797  - met_h_eta_hi * mf_my(i,j,0) * yfluxbar_hi;
798 
799  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
800  Real stateContrib = ( xflux(i+1,j ,k ) * ax(i+1,j,k) / mf_uy(i+1,j,0)
801  -xflux(i ,j ,k ) * ax(i ,j,k) / mf_uy(i ,j,0) ) * dx_inv * mfsq // Diffusive flux in x-dir
802  +( yflux(i ,j+1,k ) * ay(i,j+1,k) / mf_vx(i,j+1,0)
803  -yflux(i ,j ,k ) * ay(i,j ,k) / mf_vx(i,j ,0) ) * dy_inv * mfsq // Diffusive flux in y-dir
804  +( zflux_hi - zflux_lo) * dz_inv; // Diffusive flux in z-dir
805 
806  stateContrib /= detJ(i,j,k);
807 
808  cell_rhs(i,j,k,qty_index) -= stateContrib;
809  });
810  } // n
811 
812  const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_cc};
813 #include "ERF_AddTKESources.H"
814 #include "ERF_AddQKESources.H"
815 }
ERF_AddQKESources.H
ERF_AddTKESources.H
OpenSide::hi
@ hi
OpenSide::lo
@ lo
DiffusionSrcForState_T
void DiffusionSrcForState_T(const Box &bx, const Box &domain, int start_comp, int num_comp, const bool &rotate, 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 > &z_cc, const Array4< const Real > &ax, const Array4< const Real > &ay, const Array4< const Real > &, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &SmnSmn_a, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_ux, const Array4< const Real > &mf_vx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vy, 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_SurfLayer, const Vector< std::unique_ptr< SurfaceLayer >> &SurfLayer, const Real implicit_fac)
Definition: ERF_DiffusionSrcForState_T.cpp:55
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:43
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:46
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
PrimScalar_comp
#define PrimScalar_comp
Definition: ERF_IndexDefines.H:60
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:41
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:35
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
ERF_SetupDiff.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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:292
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:409
Compute_Z_AtWFace
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Z_AtWFace(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:729
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:433
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:385
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:456
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:93
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:297
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:300
module_model_constants::c2
real(c_double), parameter c2
Definition: ERF_module_model_constants.F90:35
module_mp_morr_two_moment::c1
real(c_double), private c1
Definition: ERF_module_mp_morr_two_moment.F90:212
PBLDerivativeDzInv_T
Functor for inverse vertical spacings for terrain-following grids using cell-center heights.
Definition: ERF_PBLModels.H:461
Here is the call graph for this function: