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, 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,
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], dy = dx_arr[1], dz = dx_arr[2];
72  const Real vol = dx * dy * dz;
73 
74  for (int n(0); n<num_comp; ++n) {
75  const int qty_index = start_comp + n;
76  const int prim_index = qty_index - 1;
77  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
78  const int eff_index = (l_consA && l_turb) ? prim_scal_index : prim_index;
79  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
80  BCVars::RhoScalar_bc_comp : qty_index;
81  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
82  const Real alpha_mol = d_alpha_eff[eff_index];
83  const int eddy_x = d_eddy_diff_idx[eff_index];
84  const int eddy_y = d_eddy_diff_idy[eff_index];
85  const int eddy_z = d_eddy_diff_idz[eff_index];
86 
87  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
88  {
89  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ) : one;
90  Real rhoAlpha = rhoFace * alpha_mol;
91  if (l_turb) {
92  rhoAlpha += myhalf * ( mu_turb(i , j, k, eddy_x)
93  + mu_turb(i-1, j, k, eddy_x) );
94  }
95 
96  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
97  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
98  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= zero) );
99  ext_dir_on_xlo &= (i == dom_lo.x);
100 
101  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
102  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
103  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= zero) );
104  ext_dir_on_xhi &= (i == dom_hi.x+1);
105 
106  if (ext_dir_on_xlo) {
107  xflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i-1, j, k, prim_index)
108  + three * cell_prim(i , j, k, prim_index)
109  - (one/three) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
110  } else if (ext_dir_on_xhi) {
111  xflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i , j, k, prim_index)
112  - three * cell_prim(i-1, j, k, prim_index)
113  + (one/three) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
114  } else {
115  if (cfg_arr(i,j,k).isCovered()) {
116  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i-3, j, k, prim_index)
117  - three*cell_prim(i-2, j, k, prim_index)
118  + two*cell_prim(i-1, j, k, prim_index) ) * dx_inv;
119  } else if (cfg_arr(i-1,j,k).isCovered()) {
120  xflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i+1, j, k, prim_index)
121  - cell_prim(i+2, j, k, prim_index)
122  - two*cell_prim(i, j, k, prim_index) ) * dx_inv;
123  } else {
124  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
125  - cell_prim(i-1, j, k, prim_index) ) * dx_inv;
126  }
127  }
128  });
129  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
130  {
131  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ) : one;
132  Real rhoAlpha = rhoFace * alpha_mol;
133  if (l_turb) {
134  rhoAlpha += myhalf * ( mu_turb(i, j , k, eddy_y)
135  + mu_turb(i, j-1, k, eddy_y) );
136  }
137 
138  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
139  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
140  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= zero) );
141  ext_dir_on_ylo &= (j == dom_lo.y);
142 
143  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
144  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
145  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= zero) );
146  ext_dir_on_yhi &= (j == dom_hi.y+1);
147 
148  if (ext_dir_on_ylo) {
149  yflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j-1, k, prim_index)
150  + three * cell_prim(i, j , k, prim_index)
151  - (one/three) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
152  } else if (ext_dir_on_yhi) {
153  yflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j , k, prim_index)
154  - three * cell_prim(i, j-1, k, prim_index)
155  + (one/three) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
156  } else {
157  if (cfg_arr(i,j,k).isCovered()) {
158  yflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j-3, k, prim_index)
159  - three*cell_prim(i, j-2, k, prim_index)
160  + two*cell_prim(i, j-1, k, prim_index) ) * dy_inv;
161  } else if (cfg_arr(i,j-1,k).isCovered()) {
162  yflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i, j+1, k, prim_index)
163  - cell_prim(i, j+2, k, prim_index)
164  - two*cell_prim(i, j, k, prim_index) ) * dy_inv;
165  } else {
166  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
167  - cell_prim(i, j-1, k, prim_index)) * dy_inv;
168  }
169  }
170  });
171  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
172  {
173  Real rhoFace = l_consA ? myhalf * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ) : one;
174  Real rhoAlpha = rhoFace * alpha_mol;
175  if (l_turb) {
176  rhoAlpha += myhalf * ( mu_turb(i, j, k , eddy_z)
177  + mu_turb(i, j, k-1, eddy_z) );
178  }
179 
180  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
181  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
182  && k == dom_lo.z);
183  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
184  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
185  && k == dom_hi.z+1);
186  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
187 
188  if (ext_dir_on_zlo) {
189  zflux(i,j,k) = -rhoAlpha * ( -(Real(8.)/three) * cell_prim(i, j, k-1, prim_index)
190  + three * cell_prim(i, j, k , prim_index)
191  - (one/three) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
192  } else if (ext_dir_on_zhi) {
193  zflux(i,j,k) = -rhoAlpha * ( (Real(8.)/three) * cell_prim(i, j, k , prim_index)
194  - three * cell_prim(i, j, k-1, prim_index)
195  + (one/three) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
196  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
197  zflux(i,j,k) = hfx_z(i,j,0);
198  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
199  zflux(i,j,k) = qfx1_z(i,j,0);
200  } else {
201  if (cfg_arr(i,j,k).isCovered()) {
202  zflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j, k-3, prim_index)
203  - three*cell_prim(i, j, k-2, prim_index)
204  + two*cell_prim(i, j, k-1, prim_index) ) * dz_inv;
205  } else if (cfg_arr(i,j,k-1).isCovered()) {
206  zflux(i,j,k) = -rhoAlpha * ( three*cell_prim(i, j, k+1, prim_index)
207  - cell_prim(i, j, k+2, prim_index)
208  - two*cell_prim(i, j, k, prim_index) ) * dz_inv;
209  } else {
210  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
211  - cell_prim(i, j, k-1, prim_index)) * dz_inv;
212  }
213  }
214 
215  // Store z-boundary fluxes.
216  // if (qty_index == RhoTheta_comp) {
217  // if (!SurfLayer_on_zlo) {
218  // hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
219  // }
220  // } else if (qty_index == RhoQ1_comp) {
221  // if (!SurfLayer_on_zlo) {
222  // qfx1_z(i,j,k) = zflux(i,j,k);
223  // }
224  // } else if (qty_index == RhoQ2_comp) {
225  // qfx2_z(i,j,k) = zflux(i,j,k);
226  // }
227  });
228 
229  // Use fluxes to compute RHS
230  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
231  {
232  if (!cfg_arr(i,j,k).isCovered()) {
233  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
234  +(ay_arr(i,j+1,k) * yflux(i ,j+1,k ) - ay_arr(i,j,k) * yflux(i, j, k)) * dy_inv
235  +(az_arr(i,j,k+1) * zflux(i ,j ,k+1) - az_arr(i,j,k) * zflux(i, j, k)) * dz_inv)
236  / detJ(i,j,k);
237  }
238  });
239 
240  // Add EB boundary contributions to fluxes
241  const bool l_rhotheta = (qty_index == RhoTheta_comp);
242  if (l_surface_layer && l_rhotheta) {
243  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
244  {
245  if (cfg_arr(i,j,k).isSingleValued()) {
246 
247  cell_rhs(i,j,k,qty_index) += barea_arr(i,j,k) * hfx_EB(i,j,k) / (vol * detJ(i,j,k));
248  }
249  });
250  }
251 
252  } // n
253 }
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, 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:52
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:24
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
dz
const Real dz
Definition: ERF_InitCustomPert_Bubble.H:25
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+myhalf) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=one) { Real dT=T_pert *(std::cos(PI *L)+one)/two;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
ERF_SetupDiff.H
d_eddy_diff_idx
int * d_eddy_diff_idx
Definition: ERF_SetupDiff.H:47
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:48
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:108
d_alpha_eff
Real * d_alpha_eff
Definition: ERF_SetupVertDiff.H:105
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:227
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:229
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:235
EBChoice
Definition: ERF_EBStruct.H:23
EBChoice::eb_boundary_type
EBBoundaryType eb_boundary_type
Definition: ERF_EBStruct.H:55
SolverChoice::ebChoice
EBChoice ebChoice
Definition: ERF_DataStruct.H:1092
Here is the call graph for this function: