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