ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_DiffusionSrcForState_N.cpp File Reference
#include "ERF_Diffusion.H"
#include "ERF_EddyViscosity.H"
#include "ERF_PBLModels.H"
#include "ERF_SetupDiff.H"
#include "ERF_AddTKESources.H"
#include "ERF_AddQKESources.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 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_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_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_SurfLayer, const Real implicit_fac)
 

Function Documentation

◆ DiffusionSrcForState_N()

void DiffusionSrcForState_N ( const Box &  bx,
const Box &  domain,
int  start_comp,
int  num_comp,
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_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_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_SurfLayer,
const Real  implicit_fac 
)

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_SurfLayerwhether we have turned on subgrid diffusion
[in]implicit_fac– factor of implicitness for vertical differences only
71 {
72  BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N);
73 
74  const Real explicit_fac = one - implicit_fac;
75 
76 #include "ERF_SetupDiff.H"
77  Real l_abs_g = std::abs(grav_gpu[2]);
78 
79  const Real dz_inv = cellSizeInv[2];
80 
81  for (int n(0); n<num_comp; ++n) {
82  const int qty_index = start_comp + n;
83 
84  // Constant alpha & Turb model
85  if (l_consA && l_turb) {
86  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
87  {
88  const int prim_index = qty_index - 1;
89  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
90 
91  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
92  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
93  rhoAlpha += myhalf * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index])
94  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) );
95 
96  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
97  BCVars::RhoScalar_bc_comp : qty_index;
98  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
99 
100  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
101  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
102  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= zero) );
103  ext_dir_on_xlo &= (i == dom_lo.x);
104 
105  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
106  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
107  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= zero) );
108  ext_dir_on_xlo &= (i == dom_hi.x+1);
109 
110  if (ext_dir_on_xlo) {
111  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
112  + three * cell_prim(i , j, k, prim_index)
113  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
114  } else if (ext_dir_on_xhi) {
115  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
116  - three * cell_prim(i-1, j, k, prim_index)
117  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
118  } else {
119  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
120  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
121  }
122  });
123  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
124  {
125  const int prim_index = qty_index - 1;
126  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
127 
128  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
129  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
130  rhoAlpha += myhalf * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index])
131  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) );
132 
133  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
134  BCVars::RhoScalar_bc_comp : qty_index;
135  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
136  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
137  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
138  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
139  ext_dir_on_ylo &= (j == dom_lo.y);
140  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
141  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
142  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
143  ext_dir_on_yhi &= (j == dom_hi.y+1);
144 
145  if (ext_dir_on_ylo) {
146  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
147  + three * cell_prim(i, j , k, prim_index)
148  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
149  } else if (ext_dir_on_yhi) {
150  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
151  - three * cell_prim(i, j-1, k, prim_index)
152  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
153  } else {
154  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
155  }
156  });
157  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
158  {
159  const int prim_index = qty_index - 1;
160  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
161 
162  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
163  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
164  rhoAlpha += myhalf * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index])
165  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) );
166 
167  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
168  BCVars::RhoScalar_bc_comp : qty_index;
169  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
170 
171  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
172  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
173  && k == dom_lo.z);
174  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
175  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
176  && k == dom_hi.z+1);
177  bool SurfLayer_on_zlo = ( use_SurfLayer && k==dom_lo.z);
178 
179  if (ext_dir_on_zlo) {
180  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
181  + three * cell_prim(i, j, k , prim_index)
182  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
183  } else if (ext_dir_on_zhi) {
184  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
185  - three * cell_prim(i, j, k-1, prim_index)
186  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
187  } else if (SurfLayer_on_zlo) {
188  if (qty_index == RhoTheta_comp) {
189  zflux(i,j,k) = hfx_z(i,j,0);
190  } else if (qty_index == RhoQ1_comp) {
191  zflux(i,j,k) = qfx1_z(i,j,0);
192  } else {
193  zflux(i,j,k) = zero;
194  }
195  } else {
196  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
197  }
198 
199  if (qty_index == RhoTheta_comp) {
200  if (!SurfLayer_on_zlo) {
201  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
202  }
203  } else if (qty_index == RhoQ1_comp) {
204  if (!SurfLayer_on_zlo) {
205  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
206  }
207  } else if (qty_index == RhoQ2_comp) {
208  qfx2_z(i,j,k) = zflux(i,j,k);
209  }
210  });
211  // Constant rho*alpha & Turb model
212  } else if (l_turb) {
213  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
214  {
215  const int prim_index = qty_index - 1;
216 
217  Real rhoAlpha = d_alpha_eff[prim_index];
218  rhoAlpha += myhalf * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
219  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
220 
221  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
222  BCVars::RhoScalar_bc_comp : qty_index;
223  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
224 
225  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
226  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
227  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= zero) );
228  ext_dir_on_xlo &= (i == dom_lo.x);
229 
230  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
231  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
232  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= zero) );
233  ext_dir_on_xhi &= (i == dom_hi.x+1);
234 
235  if (ext_dir_on_xlo) {
236  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
237  + three * cell_prim(i , j, k, prim_index)
238  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
239  } else if (ext_dir_on_xhi) {
240  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
241  - three * cell_prim(i-1, j, k, prim_index)
242  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
243  } else {
244  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
245  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
246  }
247  });
248  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
249  {
250  const int prim_index = qty_index - 1;
251 
252  Real rhoAlpha = d_alpha_eff[prim_index];
253  rhoAlpha += myhalf * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
254  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
255 
256  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
257  BCVars::RhoScalar_bc_comp : qty_index;
258  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
259 
260  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
261  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
262  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
263  ext_dir_on_ylo &= (j == dom_lo.y);
264 
265  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
266  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
267  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
268  ext_dir_on_yhi &= (j == dom_hi.y+1);
269 
270  if (ext_dir_on_ylo) {
271  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
272  + three * cell_prim(i, j , k, prim_index)
273  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
274  } else if (ext_dir_on_yhi) {
275  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
276  - three * cell_prim(i, j-1, k, prim_index)
277  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
278  } else {
279  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
280  }
281  });
282  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
283  {
284  const int prim_index = qty_index - 1;
285 
286  Real rhoAlpha = d_alpha_eff[prim_index];
287  rhoAlpha += myhalf * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
288  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
289 
290  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
291  BCVars::RhoScalar_bc_comp : qty_index;
292  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
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 SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
300 
301  if (ext_dir_on_zlo) {
302  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
303  + three * cell_prim(i, j, k , prim_index)
304  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
305  } else if (ext_dir_on_zhi) {
306  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
307  - three * cell_prim(i, j, k-1, prim_index)
308  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
309  } else if (SurfLayer_on_zlo) {
310  if (qty_index == RhoTheta_comp) {
311  zflux(i,j,k) = hfx_z(i,j,0);
312  } else if (qty_index == RhoQ1_comp) {
313  zflux(i,j,k) = qfx1_z(i,j,0);
314  } else {
315  zflux(i,j,k) = zero;
316  }
317  } else {
318  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
319  }
320 
321  if (qty_index == RhoTheta_comp) {
322  if (!SurfLayer_on_zlo) {
323  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
324  }
325  } else if (qty_index == RhoQ1_comp) {
326  if (!SurfLayer_on_zlo) {
327  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
328  }
329  } else if (qty_index == RhoQ2_comp) {
330  qfx2_z(i,j,k) = zflux(i,j,k);
331  }
332  });
333  // Constant alpha & no LES/PBL model
334  } else if(l_consA) {
335  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
336  {
337  const int prim_index = qty_index - 1;
338 
339  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
340  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
341 
342  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
343  BCVars::RhoScalar_bc_comp : qty_index;
344  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
345 
346  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
347  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
348  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= zero) );
349  ext_dir_on_xlo &= (i == dom_lo.x);
350 
351  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
352  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
353  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= zero) );
354  ext_dir_on_xhi &= (i == dom_hi.x+1);
355 
356  if (ext_dir_on_xlo) {
357  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
358  + three * cell_prim(i , j, k, prim_index)
359  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
360  } else if (ext_dir_on_xhi) {
361  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
362  - three * cell_prim(i-1, j, k, prim_index)
363  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
364  } else {
365  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
366  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
367  }
368  });
369  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
370  {
371  const int prim_index = qty_index - 1;
372 
373  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
374  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
375 
376  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
377  BCVars::RhoScalar_bc_comp : qty_index;
378  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
379 
380  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
381  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
382  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
383  ext_dir_on_ylo &= (j == dom_lo.y);
384 
385  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
386  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
387  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
388  ext_dir_on_yhi &= (j == dom_hi.y+1);
389 
390  if (ext_dir_on_ylo) {
391  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
392  + three * cell_prim(i, j , k, prim_index)
393  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
394  } else if (ext_dir_on_yhi) {
395  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
396  - three * cell_prim(i, j-1, k, prim_index)
397  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
398  } else {
399  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
400  }
401  });
402  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
403  {
404  const int prim_index = qty_index - 1;
405 
406  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
407  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
408 
409  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
410  BCVars::RhoScalar_bc_comp : qty_index;
411  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
412  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
413  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
414  && k == dom_lo.z);
415  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
416  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
417  && k == dom_hi.z+1);
418  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
419 
420  if (ext_dir_on_zlo) {
421  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
422  + three * cell_prim(i, j, k , prim_index)
423  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
424  } else if (ext_dir_on_zhi) {
425  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
426  - three * cell_prim(i, j, k-1, prim_index)
427  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
428  } else if (SurfLayer_on_zlo) {
429  if (qty_index == RhoTheta_comp) {
430  zflux(i,j,k) = hfx_z(i,j,0);
431  } else if (qty_index == RhoQ1_comp) {
432  zflux(i,j,k) = qfx1_z(i,j,0);
433  } else {
434  zflux(i,j,k) = zero;
435  }
436  } else {
437  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
438  }
439 
440  if (qty_index == RhoTheta_comp) {
441  if (!SurfLayer_on_zlo) {
442  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
443  }
444  } else if (qty_index == RhoQ1_comp) {
445  if (!SurfLayer_on_zlo) {
446  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
447  }
448  } else if (qty_index == RhoQ2_comp) {
449  qfx2_z(i,j,k) = zflux(i,j,k);
450  }
451  });
452  // Constant rho*alpha & no LES/PBL model
453  } else {
454  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
455  {
456  const int prim_index = qty_index - 1;
457 
458  Real rhoAlpha = d_alpha_eff[prim_index];
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) >= zero) );
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) <= zero) );
472  ext_dir_on_xhi &= (i == dom_hi.x+1);
473 
474  if (ext_dir_on_xlo) {
475  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
476  + three * cell_prim(i , j, k, prim_index)
477  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
478  } else if (ext_dir_on_xhi) {
479  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
480  - three * cell_prim(i-1, j, k, prim_index)
481  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
482  } else {
483  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
484  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
485  }
486  });
487  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
488  {
489  const int prim_index = qty_index - 1;
490 
491  Real rhoAlpha = d_alpha_eff[prim_index];
492 
493  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
494  BCVars::RhoScalar_bc_comp : qty_index;
495  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
496 
497  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
498  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
499  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
500  ext_dir_on_ylo &= (j == dom_lo.y);
501 
502  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
503  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
504  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
505  ext_dir_on_yhi &= (j == dom_hi.y+1);
506 
507  if (ext_dir_on_ylo) {
508  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
509  + three * cell_prim(i, j , k, prim_index)
510  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
511  } else if (ext_dir_on_yhi) {
512  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
513  - three * cell_prim(i, j-1, k, prim_index)
514  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
515  } else {
516  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
517  }
518  });
519  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
520  {
521  const int prim_index = qty_index - 1;
522 
523  Real rhoAlpha = d_alpha_eff[prim_index];
524 
525  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
526  BCVars::RhoScalar_bc_comp : qty_index;
527  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
528  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
529  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
530  && k == dom_lo.z);
531  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
532  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
533  && k == dom_hi.z+1);
534  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
535 
536  if (ext_dir_on_zlo) {
537  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
538  + three * cell_prim(i, j, k , prim_index)
539  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
540  } else if (ext_dir_on_zhi) {
541  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
542  - three * cell_prim(i, j, k-1, prim_index)
543  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
544  } else if (SurfLayer_on_zlo) {
545  if (qty_index == RhoTheta_comp) {
546  zflux(i,j,k) = hfx_z(i,j,0);
547  } else if (qty_index == RhoQ1_comp) {
548  zflux(i,j,k) = qfx1_z(i,j,0);
549  } else {
550  zflux(i,j,k) = zero;
551  }
552  } else {
553  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
554  }
555 
556  if (qty_index == RhoTheta_comp) {
557  if (!SurfLayer_on_zlo) {
558  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
559  }
560  } else if (qty_index == RhoQ1_comp) {
561  if (!SurfLayer_on_zlo) {
562  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
563  }
564  } else if (qty_index == RhoQ2_comp) {
565  qfx2_z(i,j,k) = zflux(i,j,k);
566  }
567  });
568  }
569 
570  // This allows us to do semi-implicit discretization of the vertical diffusive terms
571  if (qty_index == RhoTheta_comp ||
572  qty_index == RhoQ1_comp) {
573  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
574  {
575  zflux(i,j,k) *= explicit_fac;
576  });
577  }
578 
579  // Use fluxes to compute RHS
580  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
581  {
582  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
583  cell_rhs(i,j,k,qty_index) -= (xflux(i+1,j ,k ) - xflux(i, j, k)) * dx_inv * mfsq // Diffusive flux in x-dir
584  +(yflux(i ,j+1,k ) - yflux(i, j, k)) * dy_inv * mfsq // Diffusive flux in y-dir
585  +(zflux(i ,j ,k+1) - zflux(i, j, k)) * dz_inv; // Diffusive flux in z-dir
586  });
587  } // n
588 
589 #include "ERF_AddTKESources.H"
590 #include "ERF_AddQKESources.H"
591 }
ERF_AddQKESources.H
ERF_AddTKESources.H
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
DiffusionSrcForState_N
void DiffusionSrcForState_N(const Box &bx, const Box &domain, int start_comp, int num_comp, 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_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_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_SurfLayer, const Real implicit_fac)
Definition: ERF_DiffusionSrcForState_N.cpp:41
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
PrimScalar_comp
#define PrimScalar_comp
Definition: ERF_IndexDefines.H:57
ParallelFor
ParallelFor(grown_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);})
ERF_SetupDiff.H
d_eddy_diff_idx
int * d_eddy_diff_idx
Definition: ERF_SetupDiff.H:45
zbx
const Box zbx
Definition: ERF_SetupDiff.H:9
dx_inv
const Real dx_inv
Definition: ERF_SetupDiff.H:4
dy_inv
const Real dy_inv
Definition: ERF_SetupDiff.H:5
d_eddy_diff_idy
int * d_eddy_diff_idy
Definition: ERF_SetupDiff.H:46
xbx
const Box xbx
Definition: ERF_SetupDiff.H:7
ybx
const Box ybx
Definition: ERF_SetupDiff.H:8
dom_hi
const auto & dom_hi
Definition: ERF_SetupVertDiff.H:2
l_turb
bool l_turb
Definition: ERF_SetupVertDiff.H:9
dom_lo
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
l_consA
bool l_consA
Definition: ERF_SetupVertDiff.H:8
d_eddy_diff_idz
int * d_eddy_diff_idz
Definition: ERF_SetupVertDiff.H:99
d_alpha_eff
Real * d_alpha_eff
Definition: ERF_SetupVertDiff.H:98
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:90
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:243
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:246
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:251
Here is the call graph for this function: