ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_DiffusionSrcForState_N.cpp File Reference
#include "ERF_Diffusion.H"
#include "ERF_EddyViscosity.H"
#include "ERF_PBLModels.H"
#include "ERF_DiffSetup.H"
#include "ERF_DiffTKEAdjustment.H"
#include "ERF_DiffQKEAdjustment.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)
 

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 
)

Function for computing the scalar RHS for diffusion operator without terrain.

Parameters
[in]bxcell center box to loop over
[in]domainbox of the whole domain
[in]start_compstarting component index
[in]num_compnumber of components
[in]uvelocity in x-dir
[in]vvelocity in y-dir
[in]cell_dataconserved cell center vars
[in]cell_primprimitive cell center vars
[out]cell_rhsRHS for cell center vars
[in]xfluxflux in x-dir
[in]yfluxflux in y-dir
[in]zfluxflux in z-dir
[in]cellSizeInvinverse cell size array
[in]SmnSmn_astrain rate magnitude
[in]mf_mmap factor at cell center
[in]mf_umap factor at x-face
[in]mf_vmap factor at y-face
[in,out]hfx_zheat flux in z-dir
[in,out]qfx1_zheat flux in z-dir
[out]qfx2_zheat flux in z-dir
[in]dissdissipation of TKE
[in]mu_turbturbulent viscosity
[in]diffChoicecontainer of diffusion parameters
[in]turbChoicecontainer of turbulence parameters
[in]tm_arrtheta mean array
[in]grav_gpugravity vector
[in]bc_ptrcontainer with boundary conditions
[in]use_mostwhether we have turned on MOST BCs
69 {
70  BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N);
71 
72 #include "ERF_DiffSetup.H"
73 
74  const Real dz_inv = cellSizeInv[2];
75 
76  for (int n(0); n<num_comp; ++n) {
77  const int qty_index = start_comp + n;
78 
79  // Compute fluxes at each face
80  if (l_consA && l_turb) {
81  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
82  {
83  const int prim_index = qty_index - 1;
84  // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
85  const int prim_scal_index = prim_index;
86 
87  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
88  BCVars::RhoScalar_bc_comp : qty_index;
89  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
90 
91  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
92  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
93  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
94  ext_dir_on_xlo &= (i == dom_lo.x);
95 
96  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
97  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
98  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
99  ext_dir_on_xlo &= (i == dom_hi.x+1);
100 
101  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
102  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
103  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index])
104  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) );
105 
106  if (ext_dir_on_xlo) {
107  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
108  + 3. * cell_prim(i , j, k, prim_index)
109  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
110  } else if (ext_dir_on_xhi) {
111  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
112  - 3. * cell_prim(i-1, j, k, prim_index)
113  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
114  } else {
115  xflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) -
116  cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
117  }
118  });
119  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
120  {
121  const int prim_index = qty_index - 1;
122  // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
123  const int prim_scal_index = prim_index;
124 
125  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
126  BCVars::RhoScalar_bc_comp : qty_index;
127  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
128  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
129  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
130  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
131  ext_dir_on_ylo &= (j == dom_lo.y);
132  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
133  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
134  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
135  ext_dir_on_yhi &= (j == dom_hi.y+1);
136 
137  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
138  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
139  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index])
140  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) );
141 
142  if (ext_dir_on_ylo) {
143  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
144  + 3. * cell_prim(i, j , k, prim_index)
145  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
146  } else if (ext_dir_on_yhi) {
147  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
148  - 3. * cell_prim(i, j-1, k, prim_index)
149  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
150  } else {
151  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);
152  }
153  });
154  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
155  {
156  const int prim_index = qty_index - 1;
157  // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
158  const int prim_scal_index = prim_index;
159 
160  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
161  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
162  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index])
163  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) );
164 
165  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
166  BCVars::RhoScalar_bc_comp : qty_index;
167  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
168 
169  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
170  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
171  && k == dom_lo.z);
172  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
173  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
174  && k == dom_hi.z+1);
175  bool SurfLayer_on_zlo = ( use_SurfLayer && k==dom_lo.z);
176 
177  if (ext_dir_on_zlo) {
178  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
179  + 3. * cell_prim(i, j, k , prim_index)
180  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
181  } else if (ext_dir_on_zhi) {
182  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
183  - 3. * cell_prim(i, j, k-1, prim_index)
184  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
185  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
186  zflux(i,j,k) = hfx_z(i,j,0);
187  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
188  zflux(i,j,k) = qfx1_z(i,j,0);
189  } else {
190  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
191  }
192 
193  if (qty_index == RhoTheta_comp) {
194  if (!SurfLayer_on_zlo) {
195  hfx_z(i,j,k) = zflux(i,j,k);
196  }
197  } else if (qty_index == RhoQ1_comp) {
198  if (!SurfLayer_on_zlo) {
199  qfx1_z(i,j,k) = zflux(i,j,k);
200  }
201  } else if (qty_index == RhoQ2_comp) {
202  qfx2_z(i,j,k) = zflux(i,j,k);
203  }
204  });
205  } else if (l_turb) {
206  // with MolecDiffType::Constant or None
207  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
208  {
209  const int prim_index = qty_index - 1;
210 
211  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
212  BCVars::RhoScalar_bc_comp : qty_index;
213  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
214 
215  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
216  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
217  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
218  ext_dir_on_xlo &= (i == dom_lo.x);
219 
220  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
221  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
222  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
223  ext_dir_on_xhi &= (i == dom_hi.x+1);
224 
225  Real rhoAlpha = d_alpha_eff[prim_index];
226  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
227  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
228 
229  if (ext_dir_on_xlo) {
230  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
231  + 3. * cell_prim(i , j, k, prim_index)
232  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
233  } else if (ext_dir_on_xhi) {
234  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
235  - 3. * cell_prim(i-1, j, k, prim_index)
236  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
237  } else {
238  xflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
239  }
240  });
241  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
242  {
243  const int prim_index = qty_index - 1;
244 
245  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
246  BCVars::RhoScalar_bc_comp : qty_index;
247  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
248 
249  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
250  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
251  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
252  ext_dir_on_ylo &= (j == dom_lo.y);
253 
254  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
255  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
256  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
257  ext_dir_on_yhi &= (j == dom_hi.y+1);
258 
259  Real rhoAlpha = d_alpha_eff[prim_index];
260  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
261  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
262 
263  if (ext_dir_on_ylo) {
264  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
265  + 3. * cell_prim(i, j , k, prim_index)
266  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
267  } else if (ext_dir_on_yhi) {
268  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
269  - 3. * cell_prim(i, j-1, k, prim_index)
270  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
271  } else {
272  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);
273  }
274  });
275  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
276  {
277  const int prim_index = qty_index - 1;
278 
279  Real rhoAlpha = d_alpha_eff[prim_index];
280  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
281  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
282 
283  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
284  BCVars::RhoScalar_bc_comp : qty_index;
285  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
286  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
287  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
288  && k == dom_lo.z);
289  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
290  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
291  && k == dom_hi.z+1);
292  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
293 
294  if (ext_dir_on_zlo) {
295  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
296  + 3. * cell_prim(i, j, k , prim_index)
297  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
298  } else if (ext_dir_on_zhi) {
299  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
300  - 3. * cell_prim(i, j, k-1, prim_index)
301  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
302  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
303  zflux(i,j,k) = hfx_z(i,j,0);
304  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
305  zflux(i,j,k) = qfx1_z(i,j,0);
306  } else {
307  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
308  }
309 
310  if (qty_index == RhoTheta_comp) {
311  if (!SurfLayer_on_zlo) {
312  hfx_z(i,j,k) = zflux(i,j,k);
313  }
314  } else if (qty_index == RhoQ1_comp) {
315  if (!SurfLayer_on_zlo) {
316  qfx1_z(i,j,k) = zflux(i,j,k);
317  }
318  } else if (qty_index == RhoQ2_comp) {
319  qfx2_z(i,j,k) = zflux(i,j,k);
320  }
321  });
322  } else if(l_consA) {
323  // without an LES/PBL model
324  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
325  {
326  const int prim_index = qty_index - 1;
327 
328  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
329  BCVars::RhoScalar_bc_comp : qty_index;
330  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
331 
332  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
333  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
334  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
335  ext_dir_on_xlo &= (i == dom_lo.x);
336 
337  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
338  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
339  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
340  ext_dir_on_xhi &= (i == dom_hi.x+1);
341 
342  Real rhoFace = 0.5 * ( 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  if (ext_dir_on_xlo) {
346  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
347  + 3. * cell_prim(i , j, k, prim_index)
348  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
349  } else if (ext_dir_on_xhi) {
350  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
351  - 3. * cell_prim(i-1, j, k, prim_index)
352  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
353  } else {
354  xflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
355  }
356  });
357  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
358  {
359  const int prim_index = qty_index - 1;
360 
361  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
362  BCVars::RhoScalar_bc_comp : qty_index;
363  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
364 
365  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
366  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
367  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
368  ext_dir_on_ylo &= (j == dom_lo.y);
369 
370  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
371  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
372  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
373  ext_dir_on_yhi &= (j == dom_hi.y+1);
374 
375  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
376  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
377 
378  if (ext_dir_on_ylo) {
379  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
380  + 3. * cell_prim(i, j , k, prim_index)
381  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
382  } else if (ext_dir_on_yhi) {
383  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
384  - 3. * cell_prim(i, j-1, k, prim_index)
385  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
386  } else {
387  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);
388  }
389  });
390  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
391  {
392  const int prim_index = qty_index - 1;
393 
394  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
395  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
396 
397  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
398  BCVars::RhoScalar_bc_comp : qty_index;
399  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
400  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
401  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
402  && k == dom_lo.z);
403  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
404  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
405  && k == dom_hi.z+1);
406  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
407 
408  if (ext_dir_on_zlo) {
409  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
410  + 3. * cell_prim(i, j, k , prim_index)
411  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
412  } else if (ext_dir_on_zhi) {
413  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
414  - 3. * cell_prim(i, j, k-1, prim_index)
415  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
416  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
417  zflux(i,j,k) = hfx_z(i,j,0);
418  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
419  zflux(i,j,k) = qfx1_z(i,j,0);
420  } else {
421  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
422  }
423 
424  if (qty_index == RhoTheta_comp) {
425  if (!SurfLayer_on_zlo) {
426  hfx_z(i,j,k) = zflux(i,j,k);
427  }
428  } else if (qty_index == RhoQ1_comp) {
429  if (!SurfLayer_on_zlo) {
430  qfx1_z(i,j,k) = zflux(i,j,k);
431  }
432  } else if (qty_index == RhoQ2_comp) {
433  qfx2_z(i,j,k) = zflux(i,j,k);
434  }
435  });
436  } else {
437  // with MolecDiffType::Constant or None
438  // without an LES/PBL model
439  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
440  {
441  const int prim_index = qty_index - 1;
442 
443  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
444  BCVars::RhoScalar_bc_comp : qty_index;
445  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
446 
447  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
448  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
449  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
450  ext_dir_on_xlo &= (i == dom_lo.x);
451 
452  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
453  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
454  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
455  ext_dir_on_xhi &= (i == dom_hi.x+1);
456 
457  Real rhoAlpha = d_alpha_eff[prim_index];
458 
459  if (ext_dir_on_xlo) {
460  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
461  + 3. * cell_prim(i , j, k, prim_index)
462  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
463  } else if (ext_dir_on_xhi) {
464  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
465  - 3. * cell_prim(i-1, j, k, prim_index)
466  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
467  } else {
468  xflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
469  }
470  });
471  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
472  {
473  const int prim_index = qty_index - 1;
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 
479  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
480  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
481  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
482  ext_dir_on_ylo &= (j == dom_lo.y);
483 
484  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
485  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
486  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
487  ext_dir_on_yhi &= (j == dom_hi.y+1);
488 
489  Real rhoAlpha = d_alpha_eff[prim_index];
490 
491  if (ext_dir_on_ylo) {
492  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
493  + 3. * cell_prim(i, j , k, prim_index)
494  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
495  } else if (ext_dir_on_yhi) {
496  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
497  - 3. * cell_prim(i, j-1, k, prim_index)
498  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
499  } else {
500  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);
501  }
502  });
503  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
504  {
505  const int prim_index = qty_index - 1;
506 
507  Real rhoAlpha = d_alpha_eff[prim_index];
508 
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 = ( use_SurfLayer && k == dom_lo.z);
519 
520  if (ext_dir_on_zlo) {
521  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
522  + 3. * cell_prim(i, j, k , prim_index)
523  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
524  } else if (ext_dir_on_zhi) {
525  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
526  - 3. * cell_prim(i, j, k-1, prim_index)
527  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
528  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
529  zflux(i,j,k) = hfx_z(i,j,0);
530  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
531  zflux(i,j,k) = qfx1_z(i,j,0);
532  } else {
533  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
534  }
535 
536  if (qty_index == RhoTheta_comp) {
537  if (!SurfLayer_on_zlo) {
538  hfx_z(i,j,k) = zflux(i,j,k);
539  }
540  } else if (qty_index == RhoQ1_comp) {
541  if (!SurfLayer_on_zlo) {
542  qfx1_z(i,j,k) = zflux(i,j,k);
543  }
544  } else if (qty_index == RhoQ2_comp) {
545  qfx2_z(i,j,k) = zflux(i,j,k);
546  }
547  });
548  }
549 
550  // Use fluxes to compute RHS
551  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
552  {
553  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
554  cell_rhs(i,j,k,qty_index) -= (xflux(i+1,j ,k ) - xflux(i, j, k)) * dx_inv * mfsq // Diffusive flux in x-dir
555  +(yflux(i ,j+1,k ) - yflux(i, j, k)) * dy_inv * mfsq // Diffusive flux in y-dir
556  +(zflux(i ,j ,k+1) - zflux(i, j, k)) * dz_inv; // Diffusive flux in z-dir
557  });
558  } // n
559 
560 #include "ERF_DiffTKEAdjustment.H"
561 #include "ERF_DiffQKEAdjustment.H"
562 }
ERF_DiffQKEAdjustment.H
ERF_DiffSetup.H
dom_hi
const auto & dom_hi
Definition: ERF_DiffSetup.H:10
l_turb
bool l_turb
Definition: ERF_DiffSetup.H:19
dom_lo
const auto & dom_lo
Definition: ERF_DiffSetup.H:9
l_consA
bool l_consA
Definition: ERF_DiffSetup.H:18
d_eddy_diff_idx
int * d_eddy_diff_idx
Definition: ERF_DiffSetup.H:120
zbx
const Box zbx
Definition: ERF_DiffSetup.H:23
d_eddy_diff_idz
int * d_eddy_diff_idz
Definition: ERF_DiffSetup.H:122
dx_inv
const Real dx_inv
Definition: ERF_DiffSetup.H:6
dy_inv
const Real dy_inv
Definition: ERF_DiffSetup.H:7
d_eddy_diff_idy
int * d_eddy_diff_idy
Definition: ERF_DiffSetup.H:121
xbx
const Box xbx
Definition: ERF_DiffSetup.H:21
d_alpha_eff
Real * d_alpha_eff
Definition: ERF_DiffSetup.H:119
ybx
const Box ybx
Definition: ERF_DiffSetup.H:22
ERF_DiffTKEAdjustment.H
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)
Definition: ERF_DiffusionSrcForState_N.cpp:40
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
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:80
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:209
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:211
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:217