ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_DiffusionSrcForState_S.cpp File Reference
#include "ERF_Diffusion.H"
#include "ERF_EddyViscosity.H"
#include "ERF_PBLModels.H"
#include "ERF_DiffSetup.H"
#include "ERF_DiffTKEAdjustment.H"
#include "ERF_DiffQKEAdjustment.H"
Include dependency graph for ERF_DiffusionSrcForState_S.cpp:

Functions

void DiffusionSrcForState_S (const Box &bx, const Box &domain, int start_comp, int num_comp, const bool &rotate, 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_x, Array4< Real > &hfx_y, Array4< Real > &hfx_z, Array4< Real > &qfx1_x, Array4< Real > &qfx1_y, Array4< Real > &qfx1_z, Array4< Real > &qfx2_z, Array4< Real > &diss, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const int level, const Array4< const Real > &tm_arr, const GpuArray< Real, AMREX_SPACEDIM > grav_gpu, const BCRec *bc_ptr, const bool use_SurfLayer)
 

Function Documentation

◆ DiffusionSrcForState_S()

void DiffusionSrcForState_S ( const Box &  bx,
const Box &  domain,
int  start_comp,
int  num_comp,
const bool &  rotate,
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_x,
Array4< Real > &  hfx_y,
Array4< Real > &  hfx_z,
Array4< Real > &  qfx1_x,
Array4< Real > &  qfx1_y,
Array4< Real > &  qfx1_z,
Array4< Real > &  qfx2_z,
Array4< Real > &  diss,
const Array4< const Real > &  mu_turb,
const SolverChoice &  solverChoice,
const int  level,
const Array4< const Real > &  tm_arr,
const GpuArray< Real, AMREX_SPACEDIM >  grav_gpu,
const BCRec *  bc_ptr,
const bool  use_SurfLayer 
)

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