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