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_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 GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, [[maybe_unused]] Array4< Real > &hfx_z, [[maybe_unused]] Array4< Real > &qfx1_z, [[maybe_unused]] Array4< Real > &qfx2_z, 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 GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
[[maybe_unused] ] Array4< Real > &  hfx_z,
[[maybe_unused] ] Array4< Real > &  qfx1_z,
[[maybe_unused] ] Array4< Real > &  qfx2_z,
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
57 {
58  BL_PROFILE_VAR("DiffusionSrcForState_EB()",DiffusionSrcForState_EB);
59 
60 #include "ERF_SetupDiff.H"
61 
62  const Real dz_inv = cellSizeInv[2];
63 
64  for (int n(0); n<num_comp; ++n) {
65  const int qty_index = start_comp + n;
66  const int prim_index = qty_index - 1;
67  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
68  const int eff_index = (l_consA && l_turb) ? prim_scal_index : prim_index;
69  int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ?
70  BCVars::RhoScalar_bc_comp : qty_index;
71  if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1);
72  const Real alpha_mol = d_alpha_eff[eff_index];
73  const int eddy_x = d_eddy_diff_idx[eff_index];
74  const int eddy_y = d_eddy_diff_idy[eff_index];
75  const int eddy_z = d_eddy_diff_idz[eff_index];
76 
77  ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
78  {
79  Real rhoFace = l_consA ? 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ) : 1.0;
80  Real rhoAlpha = rhoFace * alpha_mol;
81  if (l_turb) {
82  rhoAlpha += 0.5 * ( mu_turb(i , j, k, eddy_x)
83  + mu_turb(i-1, j, k, eddy_x) );
84  }
85 
86  bool ext_dir_on_xlo = ( (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
87  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
88  (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind && u(dom_lo.x,j,k) >= 0.) );
89  ext_dir_on_xlo &= (i == dom_lo.x);
90 
91  bool ext_dir_on_xhi = ( (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
92  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
93  (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind && u(dom_hi.x+1,j,k) <= 0.) );
94  ext_dir_on_xhi &= (i == dom_hi.x+1);
95 
96  if (ext_dir_on_xlo) {
97  xflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index)
98  + 3. * cell_prim(i , j, k, prim_index)
99  - (1./3.) * cell_prim(i+1, j, k, prim_index) ) * dx_inv;
100  } else if (ext_dir_on_xhi) {
101  xflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i , j, k, prim_index)
102  - 3. * cell_prim(i-1, j, k, prim_index)
103  + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv;
104  } else {
105  if (cfg_arr(i,j,k).isCovered()) {
106  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i-3, j, k, prim_index)
107  - 3.*cell_prim(i-2, j, k, prim_index)
108  + 2.*cell_prim(i-1, j, k, prim_index) ) * dx_inv;
109  } else if (cfg_arr(i-1,j,k).isCovered()) {
110  xflux(i,j,k) = -rhoAlpha * ( 3.*cell_prim(i+1, j, k, prim_index)
111  - cell_prim(i+2, j, k, prim_index)
112  - 2.*cell_prim(i, j, k, prim_index) ) * dx_inv;
113  } else {
114  xflux(i,j,k) = -rhoAlpha * ( cell_prim(i , j, k, prim_index)
115  - cell_prim(i-1, j, k, prim_index) ) * dx_inv;
116  }
117  }
118  });
119  ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
120  {
121  Real rhoFace = l_consA ? 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ) : 1.0;
122  Real rhoAlpha = rhoFace * alpha_mol;
123  if (l_turb) {
124  rhoAlpha += 0.5 * ( mu_turb(i, j , k, eddy_y)
125  + mu_turb(i, j-1, k, eddy_y) );
126  }
127 
128  bool ext_dir_on_ylo = ( (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
129  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
130  (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind && v(i,dom_lo.y,k) >= 0.) );
131  ext_dir_on_ylo &= (j == dom_lo.y);
132 
133  bool ext_dir_on_yhi = ( (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
134  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
135  (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind && v(i,dom_hi.y+1,k) <= 0.) );
136  ext_dir_on_yhi &= (j == dom_hi.y+1);
137 
138  if (ext_dir_on_ylo) {
139  yflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index)
140  + 3. * cell_prim(i, j , k, prim_index)
141  - (1./3.) * cell_prim(i, j+1, k, prim_index) ) * dy_inv;
142  } else if (ext_dir_on_yhi) {
143  yflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j , k, prim_index)
144  - 3. * cell_prim(i, j-1, k, prim_index)
145  + (1./3.) * cell_prim(i, j-2, k, prim_index) ) * dy_inv;
146  } else {
147  if (cfg_arr(i,j,k).isCovered()) {
148  yflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j-3, k, prim_index)
149  - 3.*cell_prim(i, j-2, k, prim_index)
150  + 2.*cell_prim(i, j-1, k, prim_index) ) * dy_inv;
151  } else if (cfg_arr(i,j-1,k).isCovered()) {
152  yflux(i,j,k) = -rhoAlpha * ( 3.*cell_prim(i, j+1, k, prim_index)
153  - cell_prim(i, j+2, k, prim_index)
154  - 2.*cell_prim(i, j, k, prim_index) ) * dy_inv;
155  } else {
156  yflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
157  - cell_prim(i, j-1, k, prim_index)) * dy_inv;
158  }
159  }
160  });
161  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
162  {
163  Real rhoFace = l_consA ? 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ) : 1.0;
164  Real rhoAlpha = rhoFace * alpha_mol;
165  if (l_turb) {
166  rhoAlpha += 0.5 * ( mu_turb(i, j, k , eddy_z)
167  + mu_turb(i, j, k-1, eddy_z) );
168  }
169 
170  bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
171  (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim))
172  && k == dom_lo.z);
173  bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
174  (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim))
175  && k == dom_hi.z+1);
176  bool SurfLayer_on_zlo = ( use_SurfLayer && k == dom_lo.z);
177 
178  if (ext_dir_on_zlo) {
179  zflux(i,j,k) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index)
180  + 3. * cell_prim(i, j, k , prim_index)
181  - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
182  } else if (ext_dir_on_zhi) {
183  zflux(i,j,k) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index)
184  - 3. * cell_prim(i, j, k-1, prim_index)
185  + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv;
186  } else if (SurfLayer_on_zlo && (qty_index == RhoTheta_comp)) {
187  zflux(i,j,k) = hfx_z(i,j,0);
188  } else if (SurfLayer_on_zlo && (qty_index == RhoQ1_comp)) {
189  zflux(i,j,k) = qfx1_z(i,j,0);
190  } else {
191  if (cfg_arr(i,j,k).isCovered()) {
192  zflux(i,j,k) = -rhoAlpha * ( cell_prim(i, j, k-3, prim_index)
193  - 3.*cell_prim(i, j, k-2, prim_index)
194  + 2.*cell_prim(i, j, k-1, prim_index) ) * dz_inv;
195  } else if (cfg_arr(i,j,k-1).isCovered()) {
196  zflux(i,j,k) = -rhoAlpha * ( 3.*cell_prim(i, j, k+1, prim_index)
197  - cell_prim(i, j, k+2, prim_index)
198  - 2.*cell_prim(i, j, k, prim_index) ) * dz_inv;
199  } else {
200  zflux(i,j,k) = -rhoAlpha * (cell_prim(i, j, k, prim_index)
201  - cell_prim(i, j, k-1, prim_index)) * dz_inv;
202  }
203  }
204 
205  // Store z-boundary fluxes.
206  // if (qty_index == RhoTheta_comp) {
207  // if (!SurfLayer_on_zlo) {
208  // hfx_z(i,j,k) = zflux(i,j,k) * explicit_fac;
209  // }
210  // } else if (qty_index == RhoQ1_comp) {
211  // if (!SurfLayer_on_zlo) {
212  // qfx1_z(i,j,k) = zflux(i,j,k);
213  // }
214  // } else if (qty_index == RhoQ2_comp) {
215  // qfx2_z(i,j,k) = zflux(i,j,k);
216  // }
217  });
218 
219  // Use fluxes to compute RHS
220  ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
221  {
222  if (!cfg_arr(i,j,k).isCovered()) {
223  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
224  +(ay_arr(i,j+1,k) * yflux(i ,j+1,k ) - ay_arr(i,j,k) * yflux(i, j, k)) * dy_inv
225  +(az_arr(i,j,k+1) * zflux(i ,j ,k+1) - az_arr(i,j,k) * zflux(i, j, k)) * dz_inv)
226  / detJ(i,j,k);
227  }
228  });
229  } // n
230 }
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 GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, [[maybe_unused]] Array4< Real > &hfx_z, [[maybe_unused]] Array4< Real > &qfx1_z, [[maybe_unused]] Array4< Real > &qfx2_z, 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:33
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
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+0.5) *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<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;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:8
dom_lo
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
l_consA
bool l_consA
Definition: ERF_SetupVertDiff.H:7
d_eddy_diff_idz
int * d_eddy_diff_idz
Definition: ERF_SetupVertDiff.H:107
d_alpha_eff
Real * d_alpha_eff
Definition: ERF_SetupVertDiff.H:104
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:225
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:227
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:233
Here is the call graph for this function: