ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_DiffusionSrcForState_S.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_S.cpp:

Functions

void DiffusionSrcForState_S (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 Gpu::DeviceVector< Real > &stretched_dz_d, 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_S()

void DiffusionSrcForState_S ( 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 Gpu::DeviceVector< Real > &  stretched_dz_d,
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 with stretched dz

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
72 {
73  BL_PROFILE_VAR("DiffusionSrcForState_S()",DiffusionSrcForState_S);
74 
75  const Real explicit_fac = one - implicit_fac;
76 
77 #include "ERF_SetupDiff.H"
78  Real l_abs_g = std::abs(grav_gpu[2]);
79 
80  int klo = domain.smallEnd(2);
81  int khi = domain.bigEnd(2);
82  auto dz_ptr = stretched_dz_d.data();
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  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
104 
105  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * GradCx;
106  });
107  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
108  {
109  const int prim_index = qty_index - 1;
110  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
111 
112  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
113  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
114  rhoAlpha += myhalf * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index])
115  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) );
116 
117  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
118  BCVars::RhoScalar_bc_comp : qty_index;
119  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
120 
121  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
122 
123  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * GradCy;
124  });
125  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
126  {
127  const int prim_index = qty_index - 1;
128  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
129 
130  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
131  Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index];
132  rhoAlpha += myhalf * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index])
133  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) );
134 
135  Real GradCz;
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_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
140  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
141  && k == dom_lo.z);
142  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
143  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim) )
144  && k == dom_hi.z+1);
145  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
146 
147  if (ext_dir_on_zlo) {
148  // Third order stencil with variable dz
149  Real dz0 = dz_ptr[k];
150  Real dz1 = dz_ptr[k+1];
151  Real idz0 = one / dz0;
152  Real f = (dz1 / dz0) + two;
153  Real f2 = f*f;
154  Real c3 = two / (f - f2);
155  Real c2 = -f2*c3;
156  Real c1 = -(one-f2)*c3;
157  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
158  + c2 * cell_prim(i, j, k , prim_index)
159  + c3 * cell_prim(i, j, k+1, prim_index) );
160  } else if (ext_dir_on_zhi) {
161  // Third order stencil with variable dz
162  Real dz0 = dz_ptr[k-1];
163  Real dz1 = dz_ptr[k-2];
164  Real idz0 = one / dz0;
165  Real f = (dz1 / dz0) + two;
166  Real f2 = f*f;
167  Real c3 = two / (f - f2);
168  Real c2 = -f2*c3;
169  Real c1 = -(one-f2)*c3;
170  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
171  + c2 * cell_prim(i, j, k-1, prim_index)
172  + c3 * cell_prim(i, j, k-2, prim_index) ) );
173  } else {
174  Real dzk_inv;
175  if (k==klo) {
176  dzk_inv = one / dz_ptr[k];
177  } else if (k==(khi+1)) {
178  dzk_inv = one / dz_ptr[k-1];
179  } else {
180  dzk_inv = two / (dz_ptr[k] + dz_ptr[k-1]);
181  }
182  GradCz = dzk_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
183  }
184 
185  if (SurfLayer_on_zlo) {
186  if (qty_index == RhoTheta_comp) {
187  zflux(i,j,k) = hfx_z(i,j,0);
188  } else if (qty_index == RhoQ1_comp) {
189  zflux(i,j,k) = qfx1_z(i,j,0);
190  } else {
191  zflux(i,j,k) = zero;
192  }
193  } else {
194  zflux(i,j,k) = -rhoAlpha * GradCz;
195  }
196 
197  if (qty_index == RhoTheta_comp) {
198  if (!SurfLayer_on_zlo) {
199  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
200  }
201  } else if (qty_index == RhoQ1_comp) {
202  if (!SurfLayer_on_zlo) {
203  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
204  }
205  } else if (qty_index == RhoQ2_comp) {
206  qfx2_z(i,j,k) = zflux(i,j,k);
207  }
208  });
209  // Constant rho*alpha & Turb model
210  } else if (l_turb) {
211  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
212  {
213  const int prim_index = qty_index - 1;
214 
215  Real rhoAlpha = d_alpha_eff[prim_index];
216  rhoAlpha += myhalf * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index])
217  + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) );
218 
219  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
220  BCVars::RhoScalar_bc_comp : qty_index;
221  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
222 
223  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
224 
225  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * GradCx;
226  });
227  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
228  {
229  const int prim_index = qty_index - 1;
230 
231  Real rhoAlpha = d_alpha_eff[prim_index];
232  rhoAlpha += myhalf * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index])
233  + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) );
234 
235  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
236  BCVars::RhoScalar_bc_comp : qty_index;
237  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
238 
239  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
240 
241  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * GradCy;
242  });
243  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
244  {
245  const int prim_index = qty_index - 1;
246 
247  Real rhoAlpha = d_alpha_eff[prim_index];
248  rhoAlpha += myhalf * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
249  + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
250 
251  Real GradCz;
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  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
256  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
257  && k == dom_lo.z);
258  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
259  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
260  && k == dom_hi.z+1);
261 
262  if (ext_dir_on_zlo) {
263  // Third order stencil with variable dz
264  Real dz0 = dz_ptr[k];
265  Real dz1 = dz_ptr[k+1];
266  Real idz0 = one / dz0;
267  Real f = (dz1 / dz0) + two;
268  Real f2 = f*f;
269  Real c3 = two / (f - f2);
270  Real c2 = -f2*c3;
271  Real c1 = -(one-f2)*c3;
272  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
273  + c2 * cell_prim(i, j, k , prim_index)
274  + c3 * cell_prim(i, j, k+1, prim_index) );
275  } else if (ext_dir_on_zhi) {
276  // Third order stencil with variable dz
277  Real dz0 = dz_ptr[k-1];
278  Real dz1 = dz_ptr[k-2];
279  Real idz0 = one / dz0;
280  Real f = (dz1 / dz0) + two;
281  Real f2 = f*f;
282  Real c3 = two / (f - f2);
283  Real c2 = -f2*c3;
284  Real c1 = -(one-f2)*c3;
285  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
286  + c2 * cell_prim(i, j, k-1, prim_index)
287  + c3 * cell_prim(i, j, k-2, prim_index) ) );
288  } else {
289  Real dzk_inv;
290  if (k==klo) {
291  dzk_inv = one / dz_ptr[k];
292  } else if (k==(khi+1)) {
293  dzk_inv = one / dz_ptr[k-1];
294  } else {
295  dzk_inv = two / (dz_ptr[k] + dz_ptr[k-1]);
296  }
297  GradCz = dzk_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
298  }
299 
300  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
301 
302  if (SurfLayer_on_zlo) {
303  if (qty_index == RhoTheta_comp) {
304  zflux(i,j,k) = hfx_z(i,j,0);
305  } else if (qty_index == RhoQ1_comp) {
306  zflux(i,j,k) = qfx1_z(i,j,0);
307  } else {
308  zflux(i,j,k) = zero;
309  }
310  } else {
311  zflux(i,j,k) = -rhoAlpha * GradCz;
312  }
313 
314  if (qty_index == RhoTheta_comp) {
315  if (!SurfLayer_on_zlo) {
316  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
317  }
318  } else if (qty_index == RhoQ1_comp) {
319  if (!SurfLayer_on_zlo) {
320  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
321  }
322  } else if (qty_index == RhoQ2_comp) {
323  qfx2_z(i,j,k) = zflux(i,j,k);
324  }
325  });
326  // Constant alpha & no LES/PBL model
327  } else if(l_consA) {
328  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
329  {
330  const int prim_index = qty_index - 1;
331 
332  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
333  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
334 
335  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
336  BCVars::RhoScalar_bc_comp : qty_index;
337  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
338 
339  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
340 
341  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * GradCx;
342  });
343  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
344  {
345  const int prim_index = qty_index - 1;
346 
347  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
348  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
349 
350  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
351  BCVars::RhoScalar_bc_comp : qty_index;
352  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
353 
354  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
355 
356  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * GradCy;
357  });
358  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
359  {
360  const int prim_index = qty_index - 1;
361 
362  Real rhoFace = myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
363  Real rhoAlpha = rhoFace * d_alpha_eff[prim_index];
364 
365  Real GradCz;
366  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
367  BCVars::RhoScalar_bc_comp : qty_index;
368  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
369  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
370  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
371  && k == dom_lo.z);
372  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
373  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
374  && k == dom_hi.z+1);
375 
376  if (ext_dir_on_zlo) {
377  // Third order stencil with variable dz
378  Real dz0 = dz_ptr[k];
379  Real dz1 = dz_ptr[k+1];
380  Real idz0 = one / dz0;
381  Real f = (dz1 / dz0) + two;
382  Real f2 = f*f;
383  Real c3 = two / (f - f2);
384  Real c2 = -f2*c3;
385  Real c1 = -(one-f2)*c3;
386  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
387  + c2 * cell_prim(i, j, k , prim_index)
388  + c3 * cell_prim(i, j, k+1, prim_index) );
389  } else if (ext_dir_on_zhi) {
390  // Third order stencil with variable dz
391  Real dz0 = dz_ptr[k-1];
392  Real dz1 = dz_ptr[k-2];
393  Real idz0 = one / dz0;
394  Real f = (dz1 / dz0) + two;
395  Real f2 = f*f;
396  Real c3 = two / (f - f2);
397  Real c2 = -f2*c3;
398  Real c1 = -(one-f2)*c3;
399  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
400  + c2 * cell_prim(i, j, k-1, prim_index)
401  + c3 * cell_prim(i, j, k-2, prim_index) ) );
402  } else {
403  Real dzk_inv;
404  if (k==klo) {
405  dzk_inv = one / dz_ptr[k];
406  } else if (k==(khi+1)) {
407  dzk_inv = one / dz_ptr[k-1];
408  } else {
409  dzk_inv = two / (dz_ptr[k] + dz_ptr[k-1]);
410  }
411  GradCz = dzk_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
412  }
413 
414  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
415 
416  if (SurfLayer_on_zlo) {
417  if (qty_index == RhoTheta_comp) {
418  zflux(i,j,k) = hfx_z(i,j,0);
419  } else if (qty_index == RhoQ1_comp) {
420  zflux(i,j,k) = qfx1_z(i,j,0);
421  } else {
422  zflux(i,j,k) = zero;
423  }
424  } else {
425  zflux(i,j,k) = -rhoAlpha * GradCz;
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) * explicit_fac;
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  Real GradCx = dx_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i-1, j, k , prim_index) );
453 
454  xflux(i,j,k) = -rhoAlpha * mf_ux(i,j,0) * GradCx;
455  });
456  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
457  {
458  const int prim_index = qty_index - 1;
459 
460  Real rhoAlpha = d_alpha_eff[prim_index];
461 
462  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
463  BCVars::RhoScalar_bc_comp : qty_index;
464  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
465 
466  Real GradCy = dy_inv * ( cell_prim(i, j, k , prim_index) - cell_prim(i, j-1, k , prim_index) );
467 
468  yflux(i,j,k) = -rhoAlpha * mf_vy(i,j,0) * GradCy;
469  });
470  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
471  {
472  const int prim_index = qty_index - 1;
473 
474  Real rhoAlpha = d_alpha_eff[prim_index];
475 
476  Real GradCz;
477  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
478  BCVars::RhoScalar_bc_comp : qty_index;
479  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
480  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
481  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
482  && k == dom_lo.z);
483  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
484  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
485  && k == dom_hi.z+1);
486 
487  if (ext_dir_on_zlo) {
488  // Third order stencil with variable dz
489  Real dz0 = dz_ptr[k];
490  Real dz1 = dz_ptr[k+1];
491  Real idz0 = one / dz0;
492  Real f = (dz1 / dz0) + two;
493  Real f2 = f*f;
494  Real c3 = two / (f - f2);
495  Real c2 = -f2*c3;
496  Real c1 = -(one-f2)*c3;
497  GradCz = idz0 * ( c1 * cell_prim(i, j, k-1, prim_index)
498  + c2 * cell_prim(i, j, k , prim_index)
499  + c3 * cell_prim(i, j, k+1, prim_index) );
500  } else if (ext_dir_on_zhi) {
501  // Third order stencil with variable dz
502  Real dz0 = dz_ptr[k-1];
503  Real dz1 = dz_ptr[k-2];
504  Real idz0 = one / dz0;
505  Real f = (dz1 / dz0) + two;
506  Real f2 = f*f;
507  Real c3 = two / (f - f2);
508  Real c2 = -f2*c3;
509  Real c1 = -(one-f2)*c3;
510  GradCz = idz0 * ( -( c1 * cell_prim(i, j, k , prim_index)
511  + c2 * cell_prim(i, j, k-1, prim_index)
512  + c3 * cell_prim(i, j, k-2, prim_index) ) );
513  } else {
514  Real dzk_inv;
515  if (k==klo) {
516  dzk_inv = one / dz_ptr[k];
517  } else if (k==(khi+1)) {
518  dzk_inv = one / dz_ptr[k-1];
519  } else {
520  dzk_inv = two / (dz_ptr[k] + dz_ptr[k-1]);
521  }
522  GradCz = dzk_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) );
523  }
524 
525  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
526 
527  if (SurfLayer_on_zlo) {
528  if (qty_index == RhoTheta_comp) {
529  zflux(i,j,k) = hfx_z(i,j,0);
530  } else if (qty_index == RhoQ1_comp) {
531  zflux(i,j,k) = qfx1_z(i,j,0);
532  } else {
533  zflux(i,j,k) = zero;
534  }
535  } else {
536  zflux(i,j,k) = -rhoAlpha * GradCz;
537  }
538 
539  if (qty_index == RhoTheta_comp) {
540  if (!SurfLayer_on_zlo) {
541  hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
542  }
543  } else if (qty_index == RhoQ1_comp) {
544  if (!SurfLayer_on_zlo) {
545  qfx1_z(i,j,k) = zflux(i,j,k) * explicit_fac;
546  }
547  } else if (qty_index == RhoQ2_comp) {
548  qfx2_z(i,j,k) = zflux(i,j,k);
549  }
550  });
551  }
552 
553  // Adjust with map factors
554  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
555  {
556  xflux(i,j,k) /= mf_uy(i,j,0);
557  });
558  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
559  {
560  yflux(i,j,k) /= mf_vx(i,j,0);
561  });
562 
563  // This allows us to do semi-implicit discretization of the vertical diffusive terms
564  if (qty_index == RhoTheta_comp ||
565  qty_index == RhoQ1_comp) {
566  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
567  {
568  zflux(i,j,k) *= explicit_fac;
569  });
570  }
571 
572  // Use fluxes to compute RHS
573  //-----------------------------------------------------------------------------------
574  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
575  {
576  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
577  Real dzk_inv = one / dz_ptr[k];
578  Real stateContrib = (xflux(i+1,j ,k ) - xflux(i, j, k)) * dx_inv * mfsq // Diffusive flux in x-dir
579  +(yflux(i ,j+1,k ) - yflux(i, j, k)) * dy_inv * mfsq // Diffusive flux in y-dir
580  +(zflux(i ,j ,k+1) - zflux(i, j, k)) * dzk_inv; // Diffusive flux in z-dir
581 
582  cell_rhs(i,j,k,qty_index) -= stateContrib;
583  });
584 
585  } // n
586 
587 #include "ERF_AddTKESources.H"
588 #include "ERF_AddQKESources.H"
589 }
ERF_AddQKESources.H
ERF_AddTKESources.H
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
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_S
void DiffusionSrcForState_S(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 Gpu::DeviceVector< Real > &stretched_dz_d, 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_S.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
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
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
module_model_constants::c2
real(c_double), parameter c2
Definition: ERF_module_model_constants.F90:35
module_mp_morr_two_moment::c1
real(c_double), private c1
Definition: ERF_module_mp_morr_two_moment.F90:212
Here is the call graph for this function: