ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_DiffusionSrcForState_EB.cpp File Reference
#include "ERF_Diffusion.H"
#include "ERF_EddyViscosity.H"
#include <ERF_EBStruct.H>
#include "ERF_SetupDiff.H"
Include dependency graph for ERF_DiffusionSrcForState_EB.cpp:

Functions

void DiffusionSrcForState_EB (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 Array4< const EBCellFlag > &cfg_arr, const Array4< const Real > &ax_arr, const Array4< const Real > &ay_arr, const Array4< const Real > &az_arr, const Array4< const Real > &detJ, [[maybe_unused]] const Array4< const Real > &barea_arr, [[maybe_unused]] const Array4< const Real > &bcent_arr, const Real *dx_arr, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, [[maybe_unused]] Array4< Real > &hfx_z, [[maybe_unused]] Array4< Real > &qfx1_z, [[maybe_unused]] Array4< Real > &qfx2_z, Array4< Real > &hfx_EB, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const int level, const BCRec *bc_ptr, const bool use_SurfLayer)
 

Function Documentation

◆ DiffusionSrcForState_EB()

void DiffusionSrcForState_EB ( 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 Array4< const EBCellFlag > &  cfg_arr,
const Array4< const Real > &  ax_arr,
const Array4< const Real > &  ay_arr,
const Array4< const Real > &  az_arr,
const Array4< const Real > &  detJ,
[[maybe_unused] ] const Array4< const Real > &  barea_arr,
[[maybe_unused] ] const Array4< const Real > &  bcent_arr,
const Real *  dx_arr,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
[[maybe_unused] ] Array4< Real > &  hfx_z,
[[maybe_unused] ] Array4< Real > &  qfx1_z,
[[maybe_unused] ] Array4< Real > &  qfx2_z,
Array4< Real > &  hfx_EB,
const Array4< const Real > &  mu_turb,
const SolverChoice &  solverChoice,
const int  level,
const BCRec *  bc_ptr,
const bool  use_SurfLayer 
)

Function for computing the scalar RHS for diffusion operator without terrain.

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]cfg_arrEB cell flags
[in]ax_arrarea fractions on x-faces
[in]ay_arrarea fractions on y-faces
[in]az_arrarea fractions on z-faces
[in]detJJacobian determinant
[in]barea_arrEB boundary area
[in]bcent_arrEB boundary centroid
[in]dx_arrcell size array
[in]cellSizeInvinverse cell size array
[in,out]hfx_zheat flux in z-dir
[in,out]qfx1_zheat flux in z-dir
[out]qfx2_zheat flux in z-dir
[in,out]hfx_EBheat flux through EB faces
[in]mu_turbturbulent viscosity
[in]solverChoicecontainer of solver and diffusion parameters
[in]levelAMR level
[in]bc_ptrcontainer with boundary conditions
[in]use_SurfLayerwhether we have turned on subgrid diffusion
70 {
71  BL_PROFILE_VAR("DiffusionSrcForState_EB()",DiffusionSrcForState_EB);
72 
73 #include "ERF_SetupDiff.H"
74 
75  EBChoice ebChoice = solverChoice.ebChoice;
76  const bool l_surface_layer = (ebChoice.eb_boundary_type == EBBoundaryType::SurfaceLayer);
77 
78  const Real dz_inv = cellSizeInv[2];
79  const Real dx = dx_arr[0];
80  const Real dy = dx_arr[1];
81  const Real dz = dx_arr[2];
82  const Real vol = dx * dy * dz;
83 
84  for (int n(0); n<num_comp; ++n) {
85  const int qty_index = start_comp + n;
86  const int prim_index = qty_index - 1;
87  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
88  const int eff_index = (l_consA && l_turb) ? prim_scal_index : prim_index;
89  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
90  BCVars::RhoScalar_bc_comp : qty_index;
91  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
92  const Real alpha_mol = d_alpha_eff[eff_index];
93  const int eddy_x = d_eddy_diff_idx[eff_index];
94  const int eddy_y = d_eddy_diff_idy[eff_index];
95  const int eddy_z = d_eddy_diff_idz[eff_index];
96 
97  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
98  {
99  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ) : one;
100  Real rhoAlpha = rhoFace * alpha_mol;
101  if (l_turb) {
102  rhoAlpha += myhalf * ( mu_turb(i , j, k, eddy_x)
103  + mu_turb(i-1, j, k, eddy_x) );
104  }
105 
106  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
107  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
108  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= zero) );
109  ext_dir_on_xlo &= (i == dom_lo.x);
110 
111  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
112  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
113  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= zero) );
114  ext_dir_on_xhi &= (i == dom_hi.x+1);
115 
116  if (ext_dir_on_xlo) {
117  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
118  + three * cell_prim(i , j, k, prim_index)
119  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
120  } else if (ext_dir_on_xhi) {
121  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
122  - three * cell_prim(i-1, j, k, prim_index)
123  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
124  } else {
125  if (cfg_arr(i,j,k).isCovered()) {
126  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i-3, j, k, prim_index)
127  - three*cell_prim(i-2, j, k, prim_index)
128  + two*cell_prim(i-1, j, k, prim_index) ) * dx_inv;
129  } else if (cfg_arr(i-1,j,k).isCovered()) {
130  xflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i+1, j, k, prim_index)
131  - cell_prim(i+2, j, k, prim_index)
132  - two*cell_prim(i, j, k, prim_index) ) * dx_inv;
133  } else {
134  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
135  - cell_prim(i-1, j, k, prim_index) ) * dx_inv;
136  }
137  }
138  });
139  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
140  {
141  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ) : one;
142  Real rhoAlpha = rhoFace * alpha_mol;
143  if (l_turb) {
144  rhoAlpha += myhalf * ( mu_turb(i, j , k, eddy_y)
145  + mu_turb(i, j-1, k, eddy_y) );
146  }
147 
148  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
149  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
150  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
151  ext_dir_on_ylo &= (j == dom_lo.y);
152 
153  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
154  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
155  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
156  ext_dir_on_yhi &= (j == dom_hi.y+1);
157 
158  if (ext_dir_on_ylo) {
159  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
160  + three * cell_prim(i, j , k, prim_index)
161  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
162  } else if (ext_dir_on_yhi) {
163  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
164  - three * cell_prim(i, j-1, k, prim_index)
165  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
166  } else {
167  if (cfg_arr(i,j,k).isCovered()) {
168  yflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j-3, k, prim_index)
169  - three*cell_prim(i, j-2, k, prim_index)
170  + two*cell_prim(i, j-1, k, prim_index) ) * dy_inv;
171  } else if (cfg_arr(i,j-1,k).isCovered()) {
172  yflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i, j+1, k, prim_index)
173  - cell_prim(i, j+2, k, prim_index)
174  - two*cell_prim(i, j, k, prim_index) ) * dy_inv;
175  } else {
176  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
177  - cell_prim(i, j-1, k, prim_index)) * dy_inv;
178  }
179  }
180  });
181  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
182  {
183  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ) : one;
184  Real rhoAlpha = rhoFace * alpha_mol;
185  if (l_turb) {
186  rhoAlpha += myhalf * ( mu_turb(i, j, k , eddy_z)
187  + mu_turb(i, j, k-1, eddy_z) );
188  }
189 
190  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
191  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
192  && k == dom_lo.z);
193  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
194  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
195  && k == dom_hi.z+1);
196  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
197 
198  if (ext_dir_on_zlo) {
199  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
200  + three * cell_prim(i, j, k , prim_index)
201  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
202  } else if (ext_dir_on_zhi) {
203  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
204  - three * cell_prim(i, j, k-1, prim_index)
205  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
206  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
207  zflux(i,j,k) = hfx_z(i,j,0);
208  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
209  zflux(i,j,k) = qfx1_z(i,j,0);
210  } else {
211  if (cfg_arr(i,j,k).isCovered()) {
212  zflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j, k-3, prim_index)
213  - three*cell_prim(i, j, k-2, prim_index)
214  + two*cell_prim(i, j, k-1, prim_index) ) * dz_inv;
215  } else if (cfg_arr(i,j,k-1).isCovered()) {
216  zflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i, j, k+1, prim_index)
217  - cell_prim(i, j, k+2, prim_index)
218  - two*cell_prim(i, j, k, prim_index) ) * dz_inv;
219  } else {
220  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
221  - cell_prim(i, j, k-1, prim_index)) * dz_inv;
222  }
223  }
224 
225  // Store z-boundary fluxes.
226  // if (qty_index == RhoTheta_comp) {
227  // if (!SurfLayer_on_zlo) {
228  // hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
229  // }
230  // } else if (qty_index == RhoQ1_comp) {
231  // if (!SurfLayer_on_zlo) {
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 
239  // Use fluxes to compute RHS
240  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
241  {
242  if (!cfg_arr(i,j,k).isCovered()) {
243  cell_rhs(i,j,k,qty_index) -= ((ax_arr(i+1,j,k) * xflux(i+1,j ,k ) - ax_arr(i,j,k) * xflux(i, j, k)) * dx_inv
244  +(ay_arr(i,j+1,k) * yflux(i ,j+1,k ) - ay_arr(i,j,k) * yflux(i, j, k)) * dy_inv
245  +(az_arr(i,j,k+1) * zflux(i ,j ,k+1) - az_arr(i,j,k) * zflux(i, j, k)) * dz_inv)
246  / detJ(i,j,k);
247  }
248  });
249 
250  // Add EB boundary contributions to fluxes
251  const bool l_rhotheta = (qty_index == RhoTheta_comp);
252  if (l_surface_layer && l_rhotheta) {
253  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
254  {
255  if (cfg_arr(i,j,k).isSingleValued()) {
256 
257  Real axm = ax_arr(i ,j ,k );
258  Real axp = ax_arr(i+1,j ,k );
259  Real aym = ay_arr(i ,j ,k );
260  Real ayp = ay_arr(i ,j+1,k );
261  Real azm = az_arr(i ,j ,k );
262  Real azp = az_arr(i ,j ,k+1);
263 
264  Real adx = (axm-axp) * dy * dz;
265  Real ady = (aym-ayp) * dx * dz;
266  Real adz = (azm-azp) * dx * dy;
267 
268  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
269 
270  cell_rhs(i,j,k,qty_index) += barea * hfx_EB(i,j,k) / (vol * detJ(i,j,k));
271  }
272  });
273  }
274 
275  } // n
276 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
DiffusionSrcForState_EB
void DiffusionSrcForState_EB(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 Array4< const EBCellFlag > &cfg_arr, const Array4< const Real > &ax_arr, const Array4< const Real > &ay_arr, const Array4< const Real > &az_arr, const Array4< const Real > &detJ, [[maybe_unused]] const Array4< const Real > &barea_arr, [[maybe_unused]] const Array4< const Real > &bcent_arr, const Real *dx_arr, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, [[maybe_unused]] Array4< Real > &hfx_z, [[maybe_unused]] Array4< Real > &qfx1_z, [[maybe_unused]] Array4< Real > &qfx2_z, Array4< Real > &hfx_EB, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const int level, const BCRec *bc_ptr, const bool use_SurfLayer)
Definition: ERF_DiffusionSrcForState_EB.cpp:42
EBBoundaryType::SurfaceLayer
@ SurfaceLayer
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
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
PrimScalar_comp
#define PrimScalar_comp
Definition: ERF_IndexDefines.H:57
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:24
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
ERF_SetupDiff.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
BCVars::RhoScalar_bc_comp
@ RhoScalar_bc_comp
Definition: ERF_IndexDefines.H:90
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:249
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:252
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:257
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
EBChoice
Definition: ERF_EBStruct.H:36
EBChoice::eb_boundary_type
EBBoundaryType eb_boundary_type
Boundary condition model applied on embedded-boundary surfaces.
Definition: ERF_EBStruct.H:75
SolverChoice::ebChoice
EBChoice ebChoice
Embedded-boundary options.
Definition: ERF_DataStruct.H:1394
Here is the call graph for this function: