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