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