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, const Vector< std::unique_ptr< SurfaceLayer >> &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,
const Vector< std::unique_ptr< SurfaceLayer >> &  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
71 {
72  BL_PROFILE_VAR("DiffusionSrcForState_EB()",DiffusionSrcForState_EB);
73 
74 #include "ERF_SetupDiff.H"
75 
76  EBChoice ebChoice = solverChoice.ebChoice;
77  const bool l_surface_layer = (ebChoice.eb_boundary_type == EBBoundaryType::SurfaceLayer);
78 
79  const Real dz_inv = cellSizeInv[2];
80  const Real dx = dx_arr[0];
81  const Real dy = dx_arr[1];
82  const Real dz = dx_arr[2];
83  const Real vol = dx * dy * dz;
84 
85  for (int n(0); n<num_comp; ++n) {
86  const int qty_index = start_comp + n;
87  const int prim_index = qty_index - 1;
88  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
89  const int eff_index = (l_consA && l_turb) ? prim_scal_index : prim_index;
90  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
91  BCVars::RhoScalar_bc_comp : qty_index;
92  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
93  const Real alpha_mol = d_alpha_eff[eff_index];
94  const int eddy_x = d_eddy_diff_idx[eff_index];
95  const int eddy_y = d_eddy_diff_idy[eff_index];
96  const int eddy_z = d_eddy_diff_idz[eff_index];
97 
98  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
99  {
100  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ) : one;
101  Real rhoAlpha = rhoFace * alpha_mol;
102  if (l_turb) {
103  rhoAlpha += myhalf * ( mu_turb(i , j, k, eddy_x)
104  + mu_turb(i-1, j, k, eddy_x) );
105  }
106 
107  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
108  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
109  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= zero) );
110  ext_dir_on_xlo &= (i == dom_lo.x);
111 
112  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
113  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
114  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= zero) );
115  ext_dir_on_xhi &= (i == dom_hi.x+1);
116 
117  if (ext_dir_on_xlo) {
118  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
119  + three * cell_prim(i , j, k, prim_index)
120  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
121  } else if (ext_dir_on_xhi) {
122  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
123  - three * cell_prim(i-1, j, k, prim_index)
124  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
125  } else {
126  if (cfg_arr(i,j,k).isCovered()) {
127  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i-3, j, k, prim_index)
128  - three*cell_prim(i-2, j, k, prim_index)
129  + two*cell_prim(i-1, j, k, prim_index) ) * dx_inv;
130  } else if (cfg_arr(i-1,j,k).isCovered()) {
131  xflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i+1, j, k, prim_index)
132  - cell_prim(i+2, j, k, prim_index)
133  - two*cell_prim(i, j, k, prim_index) ) * dx_inv;
134  } else {
135  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
136  - cell_prim(i-1, j, k, prim_index) ) * dx_inv;
137  }
138  }
139  });
140  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
141  {
142  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ) : one;
143  Real rhoAlpha = rhoFace * alpha_mol;
144  if (l_turb) {
145  rhoAlpha += myhalf * ( mu_turb(i, j , k, eddy_y)
146  + mu_turb(i, j-1, k, eddy_y) );
147  }
148 
149  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
150  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
151  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
152  ext_dir_on_ylo &= (j == dom_lo.y);
153 
154  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
155  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
156  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
157  ext_dir_on_yhi &= (j == dom_hi.y+1);
158 
159  if (ext_dir_on_ylo) {
160  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
161  + three * cell_prim(i, j , k, prim_index)
162  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
163  } else if (ext_dir_on_yhi) {
164  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
165  - three * cell_prim(i, j-1, k, prim_index)
166  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
167  } else {
168  if (cfg_arr(i,j,k).isCovered()) {
169  yflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j-3, k, prim_index)
170  - three*cell_prim(i, j-2, k, prim_index)
171  + two*cell_prim(i, j-1, k, prim_index) ) * dy_inv;
172  } else if (cfg_arr(i,j-1,k).isCovered()) {
173  yflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i, j+1, k, prim_index)
174  - cell_prim(i, j+2, k, prim_index)
175  - two*cell_prim(i, j, k, prim_index) ) * dy_inv;
176  } else {
177  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
178  - cell_prim(i, j-1, k, prim_index)) * dy_inv;
179  }
180  }
181  });
182  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
183  {
184  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ) : one;
185  Real rhoAlpha = rhoFace * alpha_mol;
186  if (l_turb) {
187  rhoAlpha += myhalf * ( mu_turb(i, j, k , eddy_z)
188  + mu_turb(i, j, k-1, eddy_z) );
189  }
190 
191  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
192  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
193  && k == dom_lo.z);
194  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
195  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
196  && k == dom_hi.z+1);
197  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
198 
199  if (ext_dir_on_zlo) {
200  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
201  + three * cell_prim(i, j, k , prim_index)
202  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
203  } else if (ext_dir_on_zhi) {
204  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
205  - three * cell_prim(i, j, k-1, prim_index)
206  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
207  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
208  zflux(i,j,k) = hfx_z(i,j,0);
209  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
210  zflux(i,j,k) = qfx1_z(i,j,0);
211  } else {
212  if (cfg_arr(i,j,k).isCovered()) {
213  zflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j, k-3, prim_index)
214  - three*cell_prim(i, j, k-2, prim_index)
215  + two*cell_prim(i, j, k-1, prim_index) ) * dz_inv;
216  } else if (cfg_arr(i,j,k-1).isCovered()) {
217  zflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i, j, k+1, prim_index)
218  - cell_prim(i, j, k+2, prim_index)
219  - two*cell_prim(i, j, k, prim_index) ) * dz_inv;
220  } else {
221  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
222  - cell_prim(i, j, k-1, prim_index)) * dz_inv;
223  }
224  }
225 
226  // Store z-boundary fluxes.
227  // if (qty_index == RhoTheta_comp) {
228  // if (!SurfLayer_on_zlo) {
229  // hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
230  // }
231  // } else if (qty_index == RhoQ1_comp) {
232  // if (!SurfLayer_on_zlo) {
233  // qfx1_z(i,j,k) = zflux(i,j,k);
234  // }
235  // } else if (qty_index == RhoQ2_comp) {
236  // qfx2_z(i,j,k) = zflux(i,j,k);
237  // }
238  });
239 
240  // Use fluxes to compute RHS
241  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
242  {
243  if (!cfg_arr(i,j,k).isCovered()) {
244  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
245  +(ay_arr(i,j+1,k) * yflux(i ,j+1,k ) - ay_arr(i,j,k) * yflux(i, j, k)) * dy_inv
246  +(az_arr(i,j,k+1) * zflux(i ,j ,k+1) - az_arr(i,j,k) * zflux(i, j, k)) * dz_inv)
247  / detJ(i,j,k);
248  }
249  });
250 
251  // Add EB boundary contributions to fluxes
252  const bool l_rhotheta = (qty_index == RhoTheta_comp);
253  if (l_surface_layer && l_rhotheta) {
254  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
255  {
256  if (cfg_arr(i,j,k).isSingleValued()) {
257 
258  Real axm = ax_arr(i ,j ,k );
259  Real axp = ax_arr(i+1,j ,k );
260  Real aym = ay_arr(i ,j ,k );
261  Real ayp = ay_arr(i ,j+1,k );
262  Real azm = az_arr(i ,j ,k );
263  Real azp = az_arr(i ,j ,k+1);
264 
265  Real adx = (axm-axp) * dy * dz;
266  Real ady = (aym-ayp) * dx * dz;
267  Real adz = (azm-azp) * dx * dy;
268 
269  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
270 
271  cell_rhs(i,j,k,qty_index) += barea * hfx_EB(i,j,k) / (vol * detJ(i,j,k));
272  }
273  });
274  }
275 
276  } // n
277 }
OpenSide::hi
@ hi
OpenSide::lo
@ lo
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, const Vector< std::unique_ptr< SurfaceLayer >> &SurfLayer)
Definition: ERF_DiffusionSrcForState_EB.cpp:42
EBBoundaryType::SurfaceLayer
@ SurfaceLayer
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
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
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:45
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
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);})
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:32
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
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:305
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
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:1975
Here is the call graph for this function: