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