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 = 1.0 - 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 = 0.5 * ( 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 += 0.5 * ( 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) >= 0.) );
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) <= 0.) );
108  ext_dir_on_xlo &= (i == dom_hi.x+1);
109 
110  if (ext_dir_on_xlo) {
111  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
112  + 3. * cell_prim(i , j, k, prim_index)
113  - (1./3.) * 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 * ( (8./3.) * cell_prim(i , j, k, prim_index)
116  - 3. * cell_prim(i-1, j, k, prim_index)
117  + (1./3.) * 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 = 0.5 * ( 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 += 0.5 * ( 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) >= 0.) );
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) <= 0.) );
143  ext_dir_on_yhi &= (j == dom_hi.y+1);
144 
145  if (ext_dir_on_ylo) {
146  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
147  + 3. * cell_prim(i, j , k, prim_index)
148  - (1./3.) * 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 * ( (8./3.) * cell_prim(i, j , k, prim_index)
151  - 3. * cell_prim(i, j-1, k, prim_index)
152  + (1./3.) * 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 = 0.5 * ( 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 += 0.5 * ( 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 * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
181  + 3. * cell_prim(i, j, k , prim_index)
182  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
183  } else if (ext_dir_on_zhi) {
184  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
185  - 3. * cell_prim(i, j, k-1, prim_index)
186  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
187  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
188  zflux(i,j,k) = hfx_z(i,j,0);
189  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
190  zflux(i,j,k) = qfx1_z(i,j,0);
191  } else {
192  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
193  }
194 
195  if (qty_index == RhoTheta_comp) {
196  if (!SurfLayer_on_zlo) {
197  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
198  }
199  } else if (qty_index == RhoQ1_comp) {
200  if (!SurfLayer_on_zlo) {
201  qfx1_z(i,j,k) = zflux(i,j,k);
202  }
203  } else if (qty_index == RhoQ2_comp) {
204  qfx2_z(i,j,k) = zflux(i,j,k);
205  }
206  });
207  // Constant rho*alpha & Turb model
208  } else if (l_turb) {
209  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
210  {
211  const int prim_index = qty_index - 1;
212 
213  Real rhoAlpha = d_alpha_eff[prim_index];
214  rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
215  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
216 
217  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
218  BCVars::RhoScalar_bc_comp : qty_index;
219  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
220 
221  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
222  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
223  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
224  ext_dir_on_xlo &= (i == dom_lo.x);
225 
226  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
227  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
228  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
229  ext_dir_on_xhi &= (i == dom_hi.x+1);
230 
231  if (ext_dir_on_xlo) {
232  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
233  + 3. * cell_prim(i , j, k, prim_index)
234  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
235  } else if (ext_dir_on_xhi) {
236  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
237  - 3. * cell_prim(i-1, j, k, prim_index)
238  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
239  } else {
240  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
241  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
242  }
243  });
244  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
245  {
246  const int prim_index = qty_index - 1;
247 
248  Real rhoAlpha = d_alpha_eff[prim_index];
249  rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
250  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
251 
252  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
253  BCVars::RhoScalar_bc_comp : qty_index;
254  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
255 
256  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
257  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
258  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
259  ext_dir_on_ylo &= (j == dom_lo.y);
260 
261  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
262  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
263  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
264  ext_dir_on_yhi &= (j == dom_hi.y+1);
265 
266  if (ext_dir_on_ylo) {
267  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
268  + 3. * cell_prim(i, j , k, prim_index)
269  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
270  } else if (ext_dir_on_yhi) {
271  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
272  - 3. * cell_prim(i, j-1, k, prim_index)
273  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
274  } else {
275  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);
276  }
277  });
278  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
279  {
280  const int prim_index = qty_index - 1;
281 
282  Real rhoAlpha = d_alpha_eff[prim_index];
283  rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
284  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
285 
286  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
287  BCVars::RhoScalar_bc_comp : qty_index;
288  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
289  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
290  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
291  && k == dom_lo.z);
292  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
293  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
294  && k == dom_hi.z+1);
295  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
296 
297  if (ext_dir_on_zlo) {
298  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
299  + 3. * cell_prim(i, j, k , prim_index)
300  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
301  } else if (ext_dir_on_zhi) {
302  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
303  - 3. * cell_prim(i, j, k-1, prim_index)
304  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
305  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
306  zflux(i,j,k) = hfx_z(i,j,0);
307  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
308  zflux(i,j,k) = qfx1_z(i,j,0);
309  } else {
310  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
311  }
312 
313  if (qty_index == RhoTheta_comp) {
314  if (!SurfLayer_on_zlo) {
315  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
316  }
317  } else if (qty_index == RhoQ1_comp) {
318  if (!SurfLayer_on_zlo) {
319  qfx1_z(i,j,k) = zflux(i,j,k);
320  }
321  } else if (qty_index == RhoQ2_comp) {
322  qfx2_z(i,j,k) = zflux(i,j,k);
323  }
324  });
325  // Constant alpha & no LES/PBL model
326  } else if(l_consA) {
327  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
328  {
329  const int prim_index = qty_index - 1;
330 
331  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
332  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
333 
334  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
335  BCVars::RhoScalar_bc_comp : qty_index;
336  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
337 
338  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
339  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
340  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
341  ext_dir_on_xlo &= (i == dom_lo.x);
342 
343  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
344  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
345  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
346  ext_dir_on_xhi &= (i == dom_hi.x+1);
347 
348  if (ext_dir_on_xlo) {
349  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
350  + 3. * cell_prim(i , j, k, prim_index)
351  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
352  } else if (ext_dir_on_xhi) {
353  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
354  - 3. * cell_prim(i-1, j, k, prim_index)
355  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
356  } else {
357  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
358  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
359  }
360  });
361  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
362  {
363  const int prim_index = qty_index - 1;
364 
365  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
366  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
367 
368  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
369  BCVars::RhoScalar_bc_comp : qty_index;
370  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
371 
372  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
373  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
374  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
375  ext_dir_on_ylo &= (j == dom_lo.y);
376 
377  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
378  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
379  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
380  ext_dir_on_yhi &= (j == dom_hi.y+1);
381 
382  if (ext_dir_on_ylo) {
383  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
384  + 3. * cell_prim(i, j , k, prim_index)
385  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
386  } else if (ext_dir_on_yhi) {
387  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
388  - 3. * cell_prim(i, j-1, k, prim_index)
389  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
390  } else {
391  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);
392  }
393  });
394  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
395  {
396  const int prim_index = qty_index - 1;
397 
398  Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
399  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
400 
401  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
402  BCVars::RhoScalar_bc_comp : qty_index;
403  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
404  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
405  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
406  && k == dom_lo.z);
407  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
408  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
409  && k == dom_hi.z+1);
410  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
411 
412  if (ext_dir_on_zlo) {
413  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
414  + 3. * cell_prim(i, j, k , prim_index)
415  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
416  } else if (ext_dir_on_zhi) {
417  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
418  - 3. * cell_prim(i, j, k-1, prim_index)
419  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
420  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
421  zflux(i,j,k) = hfx_z(i,j,0);
422  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
423  zflux(i,j,k) = qfx1_z(i,j,0);
424  } else {
425  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
426  }
427 
428  if (qty_index == RhoTheta_comp) {
429  if (!SurfLayer_on_zlo) {
430  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
431  }
432  } else if (qty_index == RhoQ1_comp) {
433  if (!SurfLayer_on_zlo) {
434  qfx1_z(i,j,k) = zflux(i,j,k);
435  }
436  } else if (qty_index == RhoQ2_comp) {
437  qfx2_z(i,j,k) = zflux(i,j,k);
438  }
439  });
440  // Constant rho*alpha & no LES/PBL model
441  } else {
442  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
443  {
444  const int prim_index = qty_index - 1;
445 
446  Real rhoAlpha = d_alpha_eff[prim_index];
447 
448  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
449  BCVars::RhoScalar_bc_comp : qty_index;
450  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
451 
452  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
453  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
454  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
455  ext_dir_on_xlo &= (i == dom_lo.x);
456 
457  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
458  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
459  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
460  ext_dir_on_xhi &= (i == dom_hi.x+1);
461 
462  if (ext_dir_on_xlo) {
463  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
464  + 3. * cell_prim(i , j, k, prim_index)
465  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
466  } else if (ext_dir_on_xhi) {
467  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
468  - 3. * cell_prim(i-1, j, k, prim_index)
469  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
470  } else {
471  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
472  - cell_prim(i-1, j, k, prim_index) ) * dx_inv * mf_ux(i,j,0)/mf_uy(i,j,0);
473  }
474  });
475  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
476  {
477  const int prim_index = qty_index - 1;
478 
479  Real rhoAlpha = d_alpha_eff[prim_index];
480 
481  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
482  BCVars::RhoScalar_bc_comp : qty_index;
483  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
484 
485  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
486  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
487  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
488  ext_dir_on_ylo &= (j == dom_lo.y);
489 
490  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
491  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
492  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
493  ext_dir_on_yhi &= (j == dom_hi.y+1);
494 
495  if (ext_dir_on_ylo) {
496  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
497  + 3. * cell_prim(i, j , k, prim_index)
498  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
499  } else if (ext_dir_on_yhi) {
500  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
501  - 3. * cell_prim(i, j-1, k, prim_index)
502  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv * mf_vy(i,j,0)/mf_vx(i,j,0);
503  } else {
504  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);
505  }
506  });
507  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
508  {
509  const int prim_index = qty_index - 1;
510 
511  Real rhoAlpha = d_alpha_eff[prim_index];
512 
513  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
514  BCVars::RhoScalar_bc_comp : qty_index;
515  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
516  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
517  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
518  && k == dom_lo.z);
519  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
520  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
521  && k == dom_hi.z+1);
522  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
523 
524  if (ext_dir_on_zlo) {
525  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
526  + 3. * cell_prim(i, j, k , prim_index)
527  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
528  } else if (ext_dir_on_zhi) {
529  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
530  - 3. * cell_prim(i, j, k-1, prim_index)
531  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
532  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
533  zflux(i,j,k) = hfx_z(i,j,0);
534  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
535  zflux(i,j,k) = qfx1_z(i,j,0);
536  } else {
537  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
538  }
539 
540  if (qty_index == RhoTheta_comp) {
541  if (!SurfLayer_on_zlo) {
542  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
543  }
544  } else if (qty_index == RhoQ1_comp) {
545  if (!SurfLayer_on_zlo) {
546  qfx1_z(i,j,k) = zflux(i,j,k);
547  }
548  } else if (qty_index == RhoQ2_comp) {
549  qfx2_z(i,j,k) = zflux(i,j,k);
550  }
551  });
552  }
553 
554  // This allows us to do semi-implicit discretization of the vertical diffusive terms
555  if (qty_index == RhoTheta_comp) {
556  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
557  {
558  zflux(i,j,k) *= explicit_fac;
559  });
560  }
561 
562  // Use fluxes to compute RHS
563  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
564  {
565  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
566  cell_rhs(i,j,k,qty_index) -= (xflux(i+1,j ,k ) - xflux(i, j, k)) * dx_inv * mfsq // Diffusive flux in x-dir
567  +(yflux(i ,j+1,k ) - yflux(i, j, k)) * dy_inv * mfsq // Diffusive flux in y-dir
568  +(zflux(i ,j ,k+1) - zflux(i, j, k)) * dz_inv; // Diffusive flux in z-dir
569  });
570  } // n
571 
572 #include "ERF_AddTKESources.H"
573 #include "ERF_AddQKESources.H"
574 }
ERF_AddQKESources.H
ERF_AddTKESources.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, const Real implicit_fac)
Definition: ERF_DiffusionSrcForState_N.cpp:41
bc_comp
const int bc_comp
Definition: ERF_Implicit.H:8
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:52
ERF_SetupDiff.H
d_eddy_diff_idx
int * d_eddy_diff_idx
Definition: ERF_SetupDiff.H:47
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:48
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:8
dom_lo
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
l_consA
bool l_consA
Definition: ERF_SetupVertDiff.H:7
d_eddy_diff_idz
int * d_eddy_diff_idz
Definition: ERF_SetupVertDiff.H:107
d_alpha_eff
Real * d_alpha_eff
Definition: ERF_SetupVertDiff.H:104
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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