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