ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_SlowRhsPre.cpp File Reference
#include "AMReX_MultiFab.H"
#include "AMReX_iMultiFab.H"
#include "AMReX_ArrayLim.H"
#include "AMReX_BCRec.H"
#include "AMReX_GpuContainers.H"
#include "AMReX_GpuPrint.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_EOS.H"
#include "ERF_Utils.H"
#include "ERF_Diffusion.H"
#include "ERF_EBAdvection.H"
#include "ERF_EB.H"
#include "ERF_SurfaceLayer.H"
Include dependency graph for ERF_SlowRhsPre.cpp:

Functions

void erf_slow_rhs_pre (int level, int finest_level, int nrk, Real dt, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_data, const MultiFab &S_prim, const MultiFab &qt, Vector< MultiFab > &S_scratch, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, std::unique_ptr< MultiFab > &z_t_mf, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const MultiFab &buoyancy, const MultiFab *zmom_crse_rhs, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab *SmnSmn, MultiFab *eddyDiffs, MultiFab *Hfx1, MultiFab *Hfx2, MultiFab *Hfx3, MultiFab *Q1fx1, MultiFab *Q1fx2, MultiFab *Q1fx3, MultiFab *Q2fx3, MultiFab *Diss, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, const Gpu::DeviceVector< BCRec > &domain_bcs_type_d, const Vector< BCRec > &domain_bcs_type_h, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, const MultiFab &ax, const MultiFab &ay, const MultiFab &az, const MultiFab &detJ, Gpu::DeviceVector< Real > &stretched_dz_d, Vector< MultiFab > &gradp, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine)
 

Function Documentation

◆ erf_slow_rhs_pre()

void erf_slow_rhs_pre ( int  level,
int  finest_level,
int  nrk,
Real  dt,
Vector< MultiFab > &  S_rhs,
Vector< MultiFab > &  S_old,
Vector< MultiFab > &  S_data,
const MultiFab &  S_prim,
const MultiFab &  qt,
Vector< MultiFab > &  S_scratch,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  zvel,
std::unique_ptr< MultiFab > &  z_t_mf,
const MultiFab &  cc_src,
const MultiFab &  xmom_src,
const MultiFab &  ymom_src,
const MultiFab &  zmom_src,
const MultiFab &  buoyancy,
const MultiFab *  zmom_crse_rhs,
Vector< std::unique_ptr< MultiFab >> &  Tau_lev,
MultiFab *  SmnSmn,
MultiFab *  eddyDiffs,
MultiFab *  Hfx1,
MultiFab *  Hfx2,
MultiFab *  Hfx3,
MultiFab *  Q1fx1,
MultiFab *  Q1fx2,
MultiFab *  Q1fx3,
MultiFab *  Q2fx3,
MultiFab *  Diss,
const Geometry  geom,
const SolverChoice solverChoice,
std::unique_ptr< SurfaceLayer > &  SurfLayer,
const Gpu::DeviceVector< BCRec > &  domain_bcs_type_d,
const Vector< BCRec > &  domain_bcs_type_h,
const MultiFab &  z_phys_nd,
const MultiFab &  z_phys_cc,
const MultiFab &  ax,
const MultiFab &  ay,
const MultiFab &  az,
const MultiFab &  detJ,
Gpu::DeviceVector< Real > &  stretched_dz_d,
Vector< MultiFab > &  gradp,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
const eb_ ebfact,
YAFluxRegister *  fr_as_crse,
YAFluxRegister *  fr_as_fine 
)

Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum.

Parameters
[in]levellevel of resolution
[in]finest_levelfinest level of resolution
[in]nrkwhich RK stage
[in]dtslow time step
[out]S_rhsRHS computed here
[in]S_oldold-time solution – used only for anelastic
[in]S_datacurrent solution
[in]S_primprimitive variables (i.e. conserved variables divided by density)
[in]S_scratchscratch space
[in]xvelx-component of velocity
[in]yvely-component of velocity
[in]zvelz-component of velocity
[in]z_t_mf rate of change of grid height – only relevant for moving terrain
[in]cc_srcsource terms for conserved variables
[in]xmom_srcsource terms for x-momentum
[in]ymom_srcsource terms for y-momentum
[in]zmom_srcsource terms for z-momentum
[in]buoyancybuoyancy source term for z-momentum
[in]zmom_crse_rhsupdate term from coarser level for z-momentum; non-zero on c/f boundary only
[in]Tau_levcomponents of stress tensor
[in]SmnSmnstrain rate magnitude
[in]eddyDiffsdiffusion coefficients for LES turbulence models
[in]Hfx3heat flux in z-dir
[in]Dissdissipation of turbulent kinetic energy
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]SurfLayerPointer to SurfaceLayer class for Monin-Obukhov Similarity Theory boundary condition
[in]domain_bcs_type_ddevice vector for domain boundary conditions
[in]domain_bcs_type_hhost vector for domain boundary conditions
[in]z_phys_ndheight coordinate at nodes
[in]axarea fractions on x-faces
[in]ayarea fractions on y-faces
[in]azarea fractions on z-faces
[in]detJJacobian of the metric transformation (= 1 if use_terrain_fitted_coords is false)
[in]gradppressure gradient
[in]mapfacmap factors
[in]ebfactEB factories for cell- and face-centered variables
[in,out]fr_as_crseYAFluxRegister at level l at level l / l+1 interface
[in,out]fr_as_fineYAFluxRegister at level l at level l-1 / l interface
108 {
109  BL_PROFILE_REGION("erf_slow_rhs_pre()");
110 
111  const BCRec* bc_ptr_d = domain_bcs_type_d.data();
112  const BCRec* bc_ptr_h = domain_bcs_type_h.data();
113 
114  DiffChoice dc = solverChoice.diffChoice;
115  TurbChoice tc = solverChoice.turbChoice[level];
116 
117  const MultiFab* t_mean_mf = nullptr;
118  if (SurfLayer) { t_mean_mf = SurfLayer->get_mac_avg(level,2); }
119 
120  const Box& domain = geom.Domain();
121  int klo = domain.smallEnd(2);
122  int khi = domain.bigEnd(2);
123 
124  const AdvType l_horiz_adv_type = solverChoice.advChoice.dycore_horiz_adv_type;
125  const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type;
126  const Real l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac;
127  const Real l_vert_upw_frac = solverChoice.advChoice.dycore_vert_upw_frac;
128  const bool l_use_stretched_dz = (solverChoice.mesh_type == MeshType::StretchedDz);
129  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type == MeshType::VariableDz);
130  const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh);
131  if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_stretched_dz || l_use_terrain_fitted_coords);
132 
133  const bool l_reflux = ( (solverChoice.coupling_type == CouplingType::TwoWay) && (nrk == 2) && (finest_level > 0) );
134 
135  const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) ||
136  (tc.les_type != LESType::None) ||
137  (tc.rans_type != RANSType::None) ||
138  (tc.pbl_type != PBLType::None) );
139  const bool l_use_turb = tc.use_kturb;
140  const bool l_need_SmnSmn = tc.use_keqn;
141 
142  const bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None);
143  const bool l_use_SurfLayer = (SurfLayer != nullptr);
144  const bool l_rotate = (solverChoice.use_rotate_surface_flux);
145 
146  const bool l_anelastic = solverChoice.anelastic[level];
147  const bool l_fixed_rho = solverChoice.fixed_density;
148 
149  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
150  const Real* dx = geom.CellSize();
151 
152  // *****************************************************************************
153  // Combine external forcing terms
154  // *****************************************************************************
155  const Array<Real,AMREX_SPACEDIM> grav{0.0, 0.0, -solverChoice.gravity};
156  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
157 
158  // *****************************************************************************
159  // Pre-computed quantities
160  // *****************************************************************************
161  int nvars = S_data[IntVars::cons].nComp();
162  const BoxArray& ba = S_data[IntVars::cons].boxArray();
163  const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap();
164 
165  int nGhost = (solverChoice.terrain_type == TerrainType::EB) ? 2 : 1;
166  MultiFab Omega(convert(ba,IntVect(0,0,1)), dm, 1, nGhost);
167 
168  std::unique_ptr<MultiFab> expr;
169  std::unique_ptr<MultiFab> dflux_x;
170  std::unique_ptr<MultiFab> dflux_y;
171  std::unique_ptr<MultiFab> dflux_z;
172 
173  if (l_use_diff) {
174  erf_make_tau_terms(level,nrk,domain_bcs_type_h,z_phys_nd,
175  S_data,xvel,yvel,zvel,Tau_lev,
176  SmnSmn,eddyDiffs,geom,solverChoice,SurfLayer,
177  stretched_dz_d, detJ,mapfac);
178 
179  dflux_x = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)), dm, nvars, 0);
180  dflux_y = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)), dm, nvars, 0);
181  dflux_z = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)), dm, nvars, 0);
182 
183  if (l_use_SurfLayer) {
184  Vector<const MultiFab*> mfs = {&S_data[IntVars::cons], &xvel, &yvel, &zvel};
185  SurfLayer->impose_SurfaceLayer_bcs(level, mfs, Tau_lev,
186  Hfx1, Hfx2, Hfx3,
187  Q1fx1, Q1fx2, Q1fx3,
188  &z_phys_nd);
189  }
190  } // l_use_diff
191 
192  // This is just cautionary to deal with grid boundaries that aren't domain boundaries
193  S_rhs[IntVars::zmom].setVal(0.0);
194 
195  // *****************************************************************************
196  // Define updates and fluxes in the current RK stage
197  // *****************************************************************************
198 #ifdef _OPENMP
199 #pragma omp parallel if (Gpu::notInLaunchRegion())
200 #endif
201  {
202  std::array<FArrayBox,AMREX_SPACEDIM> flux;
203  std::array<FArrayBox,AMREX_SPACEDIM> flux_u;
204  std::array<FArrayBox,AMREX_SPACEDIM> flux_v;
205  std::array<FArrayBox,AMREX_SPACEDIM> flux_w;
206 
207  // Cell-centered masks for EB (used for flux interpolation)
208  bool already_on_centroids = false;
209  Vector<iMultiFab> physbnd_mask;
210  physbnd_mask.resize(IntVars::NumTypes);
211  if (solverChoice.terrain_type == TerrainType::EB) {
212  physbnd_mask[IntVars::cons].define(S_data[IntVars::cons].boxArray(), S_data[IntVars::cons].DistributionMap(), 1, 1);
213  physbnd_mask[IntVars::cons].BuildMask(geom.Domain(), geom.periodicity(), 1, 1, 0, 1);
214  // physbnd_mask[IntVars::cons].FillBoundary(geom.periodicity());
215  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
216  physbnd_mask[1+dir].define(S_data[1+dir].boxArray(), S_data[1+dir].DistributionMap(), 1, 1);
217  physbnd_mask[1+dir].BuildMask(geom.Domain(), geom.periodicity(), 1, 1, 0, 1);
218  // physbnd_mask[1+dir].FillBoundary(geom.periodicity());
219  }
220  }
221 
222  for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
223  {
224  Box bx = mfi.tilebox();
225  Box tbx = mfi.nodaltilebox(0);
226  Box tby = mfi.nodaltilebox(1);
227  Box tbz = mfi.nodaltilebox(2);
228 
229  // Boxes for momentum fluxes
230  Vector<Box> tbx_grown(AMREX_SPACEDIM);
231  Vector<Box> tby_grown(AMREX_SPACEDIM);
232  Vector<Box> tbz_grown(AMREX_SPACEDIM);
233  if (solverChoice.terrain_type == TerrainType::EB) {
234  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
235  tbx_grown[dir] = tbx;
236  tby_grown[dir] = tby;
237  tbz_grown[dir] = tbz;
238  IntVect iv(1, 1, 1);
239  iv[dir] = 0;
240  tbx_grown[dir] = (tbx_grown[dir].growHi(dir,1)).grow(iv);
241  tby_grown[dir] = (tby_grown[dir].growHi(dir,1)).grow(iv);
242  tbz_grown[dir] = (tbz_grown[dir].growHi(dir,1)).grow(iv);
243  }
244  }
245 
246  // We don't compute a source term for z-momentum on the bottom or top domain boundary
247  if (tbz.smallEnd(2) == domain.smallEnd(2)) {
248  tbz.growLo(2,-1);
249  }
250  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) {
251  tbz.growHi(2,-1);
252  }
253 
254  const Array4<const Real> & cell_data = S_data[IntVars::cons].array(mfi);
255  const Array4<const Real> & cell_prim = S_prim.array(mfi);
256  const Array4<Real> & cell_rhs = S_rhs[IntVars::cons].array(mfi);
257 
258  const Array4<const Real> & cell_old = S_old[IntVars::cons].array(mfi);
259 
260  const Array4<Real const>& xmom_src_arr = xmom_src.const_array(mfi);
261  const Array4<Real const>& ymom_src_arr = ymom_src.const_array(mfi);
262  const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
263  const Array4<Real const>& buoyancy_arr = buoyancy.const_array(mfi);
264 
265  const Array4<Real const>& gpx_arr = gradp[GpVars::gpx].const_array(mfi);
266  const Array4<Real const>& gpy_arr = gradp[GpVars::gpy].const_array(mfi);
267  const Array4<Real const>& gpz_arr = gradp[GpVars::gpz].const_array(mfi);
268 
269  const Array4<Real const>& qt_arr = qt.const_array(mfi);
270 
271  const Array4<Real>& rho_u_old = S_old[IntVars::xmom].array(mfi);
272  const Array4<Real>& rho_v_old = S_old[IntVars::ymom].array(mfi);
273 
274  if (l_anelastic) {
275  // When anelastic we must reset these to 0 each RK step
276  S_scratch[IntVars::xmom][mfi].template setVal<RunOn::Device>(0.0,tbx);
277  S_scratch[IntVars::ymom][mfi].template setVal<RunOn::Device>(0.0,tby);
278  S_scratch[IntVars::zmom][mfi].template setVal<RunOn::Device>(0.0,tbz);
279  }
280 
281  Array4<Real> avg_xmom = S_scratch[IntVars::xmom].array(mfi);
282  Array4<Real> avg_ymom = S_scratch[IntVars::ymom].array(mfi);
283  Array4<Real> avg_zmom = S_scratch[IntVars::zmom].array(mfi);
284 
285  const Array4<const Real> & u = xvel.array(mfi);
286  const Array4<const Real> & v = yvel.array(mfi);
287  const Array4<const Real> & w = zvel.array(mfi);
288 
289  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
290  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
291  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
292 
293  // Map factors
294  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
295  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
296  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
297  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
298  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
299  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
300 
301  const Array4< Real>& omega_arr = Omega.array(mfi);
302 
303  Array4<const Real> z_t;
304  if (z_t_mf) {
305  z_t = z_t_mf->array(mfi);
306  } else {
307  z_t = Array4<const Real>{};
308  }
309 
310  const Array4<Real>& rho_u_rhs = S_rhs[IntVars::xmom].array(mfi);
311  const Array4<Real>& rho_v_rhs = S_rhs[IntVars::ymom].array(mfi);
312  const Array4<Real>& rho_w_rhs = S_rhs[IntVars::zmom].array(mfi);
313 
314  const Array4<Real const>& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4<const Real>{};
315 
316  // Terrain metrics
317  const Array4<const Real>& z_nd = z_phys_nd.const_array(mfi);
318  const Array4<const Real>& z_cc = z_phys_cc.const_array(mfi);
319 
320  // *****************************************************************************
321  // Define flux arrays for use in advection
322  // *****************************************************************************
323  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
324  if (solverChoice.terrain_type != TerrainType::EB) {
325  flux[dir].resize(surroundingNodes(bx,dir),2);
326  } else {
327  flux[dir].resize(surroundingNodes(bx,dir).grow(1),2);
328  }
329  flux[dir].setVal<RunOn::Device>(0.);
330  }
331  const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
332  flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
333 
334  // Define flux arrays for momentum variables (used only for EB now)
335  GpuArray<Array4<Real>, AMREX_SPACEDIM> flx_u_arr{};
336  GpuArray<Array4<Real>, AMREX_SPACEDIM> flx_v_arr{};
337  GpuArray<Array4<Real>, AMREX_SPACEDIM> flx_w_arr{};
338  if (solverChoice.terrain_type == TerrainType::EB) {
339  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
340  flux_u[dir].resize(tbx_grown[dir],1);
341  flux_v[dir].resize(tby_grown[dir],1);
342  flux_w[dir].resize(tbz_grown[dir],1);
343  flux_u[dir].setVal<RunOn::Device>(0.);
344  flux_v[dir].setVal<RunOn::Device>(0.);
345  flux_w[dir].setVal<RunOn::Device>(0.);
346  flx_u_arr[dir] = flux_u[dir].array();
347  flx_v_arr[dir] = flux_v[dir].array();
348  flx_w_arr[dir] = flux_w[dir].array();
349  }
350  }
351 
352  // *****************************************************************************
353  // Contravariant flux field
354  // *****************************************************************************
355  {
356  BL_PROFILE("slow_rhs_making_omega");
357  IntVect nGrowVect = (solverChoice.terrain_type == TerrainType::EB)
358  ? IntVect(AMREX_D_DECL(2, 2, 2)) : IntVect(AMREX_D_DECL(1, 1, 1));
359  Box gbxo = surroundingNodes(bx,2); gbxo.grow(nGrowVect);
360  //
361  // Now create Omega with momentum (not velocity) with z_t subtracted if moving terrain
362  // ONLY if not doing anelastic + terrain -- in that case Omega will be defined coming
363  // out of the projection
364  //
365  if (!l_use_terrain_fitted_coords) {
366  ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
367  omega_arr(i,j,k) = rho_w(i,j,k);
368  });
369 
370  } else {
371 
372  Box gbxo_lo = gbxo; gbxo_lo.setBig(2,domain.smallEnd(2));
373  int lo_z_face = domain.smallEnd(2);
374  if (gbxo_lo.smallEnd(2) <= lo_z_face) {
375  ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
376  omega_arr(i,j,k) = 0.;
377  });
378  }
379  Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
380  int hi_z_face = domain.bigEnd(2)+1;
381  if (gbxo_hi.bigEnd(2) >= hi_z_face) {
382  ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
383  omega_arr(i,j,k) = rho_w(i,j,k);
384  });
385  }
386 
387  if (z_t) { // Note we never do anelastic with moving terrain
388  Box gbxo_mid = gbxo; gbxo_mid.setSmall(2,1); gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
389  ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
390  // We define rho on the z-face the same way as in MomentumToVelocity/VelocityToMomentum
391  Real rho_at_face = 0.5 * (cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp));
392  omega_arr(i,j,k) = OmegaFromW(i,j,k,rho_w(i,j,k),
393  rho_u,rho_v,mf_ux,mf_vy,z_nd,dxInv) -
394  rho_at_face * z_t(i,j,k);
395  });
396  } else {
397  Box gbxo_mid = gbxo;
398  if (gbxo_mid.smallEnd(2) <= domain.smallEnd(2)) {
399  gbxo_mid.setSmall(2,1);
400  }
401  if (gbxo_mid.bigEnd(2) >= domain.bigEnd(2)+1) {
402  gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
403  }
404  ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
405  omega_arr(i,j,k) = OmegaFromW(i,j,k,rho_w(i,j,k),
406  rho_u,rho_v,mf_ux,mf_vy,z_nd,dxInv);
407  });
408  }
409  }
410  } // end profile
411 
412  // *****************************************************************************
413  // Diffusive terms (pre-computed above)
414  // *****************************************************************************
415  // No terrain diffusion
416  Array4<Real> tau11,tau22,tau33;
417  Array4<Real> tau12,tau13,tau23;
418  if (Tau_lev[TauType::tau11]) {
419  tau11 = Tau_lev[TauType::tau11]->array(mfi); tau22 = Tau_lev[TauType::tau22]->array(mfi);
420  tau33 = Tau_lev[TauType::tau33]->array(mfi); tau12 = Tau_lev[TauType::tau12]->array(mfi);
421  tau13 = Tau_lev[TauType::tau13]->array(mfi); tau23 = Tau_lev[TauType::tau23]->array(mfi);
422  } else {
423  tau11 = Array4<Real>{}; tau22 = Array4<Real>{}; tau33 = Array4<Real>{};
424  tau12 = Array4<Real>{}; tau13 = Array4<Real>{}; tau23 = Array4<Real>{};
425  }
426  // Terrain diffusion
427  Array4<Real> tau21,tau31,tau32;
428  if (Tau_lev[TauType::tau21]) {
429  tau21 = Tau_lev[TauType::tau21]->array(mfi);
430  tau31 = Tau_lev[TauType::tau31]->array(mfi);
431  tau32 = Tau_lev[TauType::tau32]->array(mfi);
432  } else {
433  tau21 = Array4<Real>{}; tau31 = Array4<Real>{}; tau32 = Array4<Real>{};
434  }
435 
436  // Strain magnitude
437  Array4<Real> SmnSmn_a;
438  if (l_need_SmnSmn) {
439  SmnSmn_a = SmnSmn->array(mfi);
440  } else {
441  SmnSmn_a = Array4<Real>{};
442  }
443 
444  // *****************************************************************************
445  // Define updates in the RHS of continuity and potential temperature equations
446  // *****************************************************************************
447  bool l_eb_terrain_cc = false; // EB terrain on cell-centered grid
448  Array4<const int> mask_arr{};
449  Array4<const EBCellFlag> cfg_arr{};
450  Array4<const Real> ax_arr{};
451  Array4<const Real> ay_arr{};
452  Array4<const Real> az_arr{};
453  Array4<const Real> fcx_arr{};
454  Array4<const Real> fcy_arr{};
455  Array4<const Real> fcz_arr{};
456  Array4<const Real> detJ_arr{};
457  Array4<const Real> u_detJ_arr{};
458  Array4<const Real> v_detJ_arr{};
459  Array4<const Real> w_detJ_arr{};
460 
461  if (solverChoice.terrain_type == TerrainType::EB)
462  {
463  EBCellFlagFab const& cfg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi];
464  cfg_arr = cfg.const_array();
465  if (cfg.getType(bx) == FabType::singlevalued) {
466  l_eb_terrain_cc = true;
467  ax_arr = (ebfact.get_const_factory())->getAreaFrac()[0]->const_array(mfi);
468  ay_arr = (ebfact.get_const_factory())->getAreaFrac()[1]->const_array(mfi);
469  az_arr = (ebfact.get_const_factory())->getAreaFrac()[2]->const_array(mfi);
470  fcx_arr = (ebfact.get_const_factory())->getFaceCent()[0]->const_array(mfi);
471  fcy_arr = (ebfact.get_const_factory())->getFaceCent()[1]->const_array(mfi);
472  fcz_arr = (ebfact.get_const_factory())->getFaceCent()[2]->const_array(mfi);
473  detJ_arr = (ebfact.get_const_factory())->getVolFrac().const_array(mfi);
474  // if (!already_on_centroids) {mask_arr = physbnd_mask[IntVars::cons].const_array(mfi);}
475  mask_arr = physbnd_mask[IntVars::cons].const_array(mfi);
476  } else {
477  ax_arr = ax.const_array(mfi);
478  ay_arr = ay.const_array(mfi);
479  az_arr = az.const_array(mfi);
480  detJ_arr = detJ.const_array(mfi);
481  }
482  u_detJ_arr = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
483  v_detJ_arr = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
484  w_detJ_arr = (ebfact.get_w_const_factory())->getVolFrac().const_array(mfi);
485  } else {
486  ax_arr = ax.const_array(mfi);
487  ay_arr = ay.const_array(mfi);
488  az_arr = az.const_array(mfi);
489  detJ_arr = detJ.const_array(mfi);
490  }
491 
492  int icomp = RhoTheta_comp; int ncomp = 1;
493  if (!l_eb_terrain_cc){
494  AdvectionSrcForRho( bx, cell_rhs,
495  rho_u, rho_v, omega_arr, // these are being used to build the fluxes
496  avg_xmom, avg_ymom, avg_zmom, // these are being defined from the fluxes
497  ax_arr, ay_arr, az_arr, detJ_arr,
498  dxInv, mf_mx, mf_my, mf_uy, mf_vx,
499  flx_arr, l_fixed_rho);
500  AdvectionSrcForScalars(bx, icomp, ncomp,
501  avg_xmom, avg_ymom, avg_zmom,
502  cell_prim, cell_rhs,
503  detJ_arr, dxInv, mf_mx, mf_my,
504  l_horiz_adv_type, l_vert_adv_type,
505  l_horiz_upw_frac, l_vert_upw_frac,
506  flx_arr, domain, bc_ptr_h);
507  } else {
508  EBAdvectionSrcForRho(bx, cell_rhs,
509  rho_u, rho_v, omega_arr,
510  avg_xmom, avg_ymom, avg_zmom,
511  mask_arr, cfg_arr,
512  ax_arr, ay_arr, az_arr,
513  fcx_arr, fcy_arr, fcz_arr, detJ_arr,
514  dxInv, mf_mx, mf_my, mf_uy, mf_vx,
515  flx_arr, l_fixed_rho,
516  already_on_centroids);
517  EBAdvectionSrcForScalars(bx, icomp, ncomp,
518  avg_xmom, avg_ymom, avg_zmom,
519  cell_prim, cell_rhs,
520  mask_arr, cfg_arr, ax_arr, ay_arr, az_arr,
521  fcx_arr, fcy_arr, fcz_arr,
522  detJ_arr, dxInv, mf_mx, mf_my,
523  l_horiz_adv_type, l_vert_adv_type,
524  l_horiz_upw_frac, l_vert_upw_frac,
525  flx_arr, domain, bc_ptr_h,
526  already_on_centroids);
527  }
528 
529  if (l_use_diff) {
530  Array4<Real> diffflux_x = dflux_x->array(mfi);
531  Array4<Real> diffflux_y = dflux_y->array(mfi);
532  Array4<Real> diffflux_z = dflux_z->array(mfi);
533 
534  Array4<Real> hfx_x = Hfx1->array(mfi);
535  Array4<Real> hfx_y = Hfx2->array(mfi);
536  Array4<Real> hfx_z = Hfx3->array(mfi);
537 
538  Array4<Real> q1fx_x = (Q1fx1) ? Q1fx1->array(mfi) : Array4<Real>{};
539  Array4<Real> q1fx_y = (Q1fx2) ? Q1fx2->array(mfi) : Array4<Real>{};
540  Array4<Real> q1fx_z = (Q1fx3) ? Q1fx3->array(mfi) : Array4<Real>{};
541 
542  Array4<Real> q2fx_z = (Q2fx3) ? Q2fx3->array(mfi) : Array4<Real>{};
543  Array4<Real> diss = Diss->array(mfi);
544 
545  const Array4<const Real> tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4<const Real>{};
546 
547  // NOTE: No diffusion for continuity, so n starts at 1.
548  int n_start = RhoTheta_comp;
549  int n_comp = 1;
550 
551  if (l_use_stretched_dz) {
552  DiffusionSrcForState_S(bx, domain, n_start, n_comp, l_rotate, u, v,
553  cell_data, cell_prim, cell_rhs,
554  diffflux_x, diffflux_y, diffflux_z,
555  stretched_dz_d, dxInv, SmnSmn_a,
556  mf_mx, mf_ux, mf_vx,
557  mf_my, mf_uy, mf_vy,
558  hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss,
559  mu_turb, solverChoice, level,
560  tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer);
561  } else if (l_use_terrain_fitted_coords) {
562  DiffusionSrcForState_T(bx, domain, n_start, n_comp, l_rotate, u, v,
563  cell_data, cell_prim, cell_rhs,
564  diffflux_x, diffflux_y, diffflux_z,
565  z_nd, z_cc, ax_arr, ay_arr, az_arr, detJ_arr,
566  dxInv, SmnSmn_a,
567  mf_mx, mf_ux, mf_vx,
568  mf_my, mf_uy, mf_vy,
569  hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss,
570  mu_turb, solverChoice, level,
571  tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer);
572  } else {
573  DiffusionSrcForState_N(bx, domain, n_start, n_comp, u, v,
574  cell_data, cell_prim, cell_rhs,
575  diffflux_x, diffflux_y, diffflux_z,
576  dxInv, SmnSmn_a,
577  mf_mx, mf_ux, mf_vx,
578  mf_my, mf_uy, mf_vy,
579  hfx_z, q1fx_z, q2fx_z, diss,
580  mu_turb, solverChoice, level,
581  tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer);
582  }
583  }
584 
585  const Array4<Real const>& source_arr = cc_src.const_array(mfi);
586  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
587  {
588  cell_rhs(i,j,k,Rho_comp) += source_arr(i,j,k,Rho_comp);
589  cell_rhs(i,j,k,RhoTheta_comp) += source_arr(i,j,k,RhoTheta_comp);
590  });
591 
592  // Multiply the slow RHS for rho and rhotheta by detJ here so we don't have to later
593  if (l_moving_terrain) {
594  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
595  {
596  cell_rhs(i,j,k,Rho_comp) *= detJ_arr(i,j,k);
597  cell_rhs(i,j,k,RhoTheta_comp) *= detJ_arr(i,j,k);
598  });
599  }
600 
601  // If anelastic and in second RK stage, take average of old-time and new-time source
602  if ( l_anelastic && (nrk == 1) )
603  {
604  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
605  {
606  cell_rhs(i,j,k, Rho_comp) *= 0.5;
607  cell_rhs(i,j,k,RhoTheta_comp) *= 0.5;
608 
609  cell_rhs(i,j,k, Rho_comp) += 0.5 / dt * (cell_data(i,j,k, Rho_comp) - cell_old(i,j,k, Rho_comp));
610  cell_rhs(i,j,k,RhoTheta_comp) += 0.5 / dt * (cell_data(i,j,k,RhoTheta_comp) - cell_old(i,j,k,RhoTheta_comp));
611  });
612  }
613 
614  // *****************************************************************************
615  // Define updates in the RHS of {x, y, z}-momentum equations
616  // *****************************************************************************
617  int lo_z_face = domain.smallEnd(2);
618  int hi_z_face = domain.bigEnd(2)+1;
619 
620  AdvectionSrcForMom(mfi, bx, tbx, tby, tbz, tbx_grown, tby_grown, tbz_grown,
621  rho_u_rhs, rho_v_rhs, rho_w_rhs,
622  cell_data, u, v, w,
623  rho_u, rho_v, omega_arr,
624  z_nd, ax_arr, ay_arr, az_arr,
625  detJ_arr, stretched_dz_d,
626  dxInv, mf_mx, mf_ux, mf_vx, mf_my, mf_uy, mf_vy,
627  l_horiz_adv_type, l_vert_adv_type,
628  l_horiz_upw_frac, l_vert_upw_frac,
629  solverChoice.mesh_type, solverChoice.terrain_type,
630  ebfact, flx_u_arr, flx_v_arr, flx_w_arr,
631  physbnd_mask, already_on_centroids,
632  lo_z_face, hi_z_face, domain, bc_ptr_h);
633 
634  if (l_use_diff) {
635  // Note: tau** were calculated with calls to
636  // ComputeStress[Cons|Var]Visc_[N|S|T] in which ConsVisc ("constant
637  // viscosity") means that there is no contribution from a
638  // turbulence model. However, whether this field truly is constant
639  // depends on whether MolecDiffType is Constant or ConstantAlpha.
640  if (solverChoice.terrain_type != TerrainType::EB) {
641  DiffusionSrcForMom(tbx, tby, tbz,
642  rho_u_rhs, rho_v_rhs, rho_w_rhs,
643  tau11, tau22, tau33,
645  detJ_arr, stretched_dz_d, dxInv,
646  mf_mx, mf_ux, mf_vx,
647  mf_my, mf_uy, mf_vy,
648  l_use_stretched_dz,
649  l_use_terrain_fitted_coords);
650  } else {
651  DiffusionSrcForMom_EB(mfi, domain, tbx, tby, tbz,
652  rho_u_rhs, rho_v_rhs, rho_w_rhs,
653  u, v, w,
654  tau11, tau22, tau33,
655  tau12, tau13, tau23,
656  dx, dxInv,
657  mf_mx, mf_ux, mf_vx,
658  mf_my, mf_uy, mf_vy,
659  solverChoice, ebfact, bc_ptr_d);
660  }
661  }
662 
663  auto abl_pressure_grad = solverChoice.abl_pressure_grad;
664 
665  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
666  { // x-momentum equation
667 
668  Real gpx = gpx_arr(i,j,k) * mf_ux(i,j,0);
669 
670  Real q = (l_use_moisture) ? 0.5 * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) : 0.0;
671 
672  rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) + xmom_src_arr(i,j,k);
673 
674  if (l_moving_terrain) {
675  Real h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd);
676  rho_u_rhs(i, j, k) *= h_zeta;
677  }
678 
679  if ( l_anelastic && (nrk == 1) ) {
680  rho_u_rhs(i,j,k) *= 0.5;
681  rho_u_rhs(i,j,k) += 0.5 / dt * (rho_u(i,j,k) - rho_u_old(i,j,k));
682  }
683  });
684 
685  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
686  { // y-momentum equation
687 
688  Real gpy = gpy_arr(i,j,k) * mf_vy(i,j,0);
689 
690  Real q = (l_use_moisture) ? 0.5 * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) : 0.0;
691 
692  rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0 + q) + ymom_src_arr(i,j,k);
693 
694  if (l_moving_terrain) {
695  Real h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd);
696  rho_v_rhs(i, j, k) *= h_zeta;
697  }
698 
699  if ( l_anelastic && (nrk == 1) ) {
700  rho_v_rhs(i,j,k) *= 0.5;
701  rho_v_rhs(i,j,k) += 0.5 / dt * (rho_v(i,j,k) - rho_v_old(i,j,k));
702  }
703  });
704 
705  // *****************************************************************************
706  // Zero out source terms for x- and y- momenta if at walls or inflow
707  // We need to do this -- even though we call the boundary conditions later --
708  // because the slow source is used to update the state in the fast interpolater.
709  // *****************************************************************************
710  if (bx.smallEnd(0) == domain.smallEnd(0)) {
711  Box lo_x_dom_face(bx); lo_x_dom_face.setBig(0,bx.smallEnd(0));
712  if (bc_ptr_h[BCVars::xvel_bc].lo(0) == ERFBCType::ext_dir) {
713  ParallelFor(lo_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
714  rho_u_rhs(i,j,k) = 0.;
715  });
716  } else if (bc_ptr_h[BCVars::xvel_bc].lo(0) == ERFBCType::ext_dir_upwind) {
717  ParallelFor(lo_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
718  if (u(i,j,k) >= 0.) {
719  rho_u_rhs(i,j,k) = 0.;
720  }
721  });
722  }
723  }
724  if (bx.bigEnd(0) == domain.bigEnd(0)) {
725  Box hi_x_dom_face(bx); hi_x_dom_face.setSmall(0,bx.bigEnd(0)+1); hi_x_dom_face.setBig(0,bx.bigEnd(0)+1);
726  if (bc_ptr_h[BCVars::xvel_bc].hi(0) == ERFBCType::ext_dir) {
727  ParallelFor(hi_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
728  rho_u_rhs(i,j,k) = 0.;
729  });
730  } else if (bc_ptr_h[BCVars::xvel_bc].hi(0) == ERFBCType::ext_dir_upwind) {
731  ParallelFor(hi_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
732  if (u(i,j,k) <= 0.) {
733  rho_u_rhs(i,j,k) = 0.;
734  }
735  });
736  }
737  }
738  if (bx.smallEnd(1) == domain.smallEnd(1)) {
739  Box lo_y_dom_face(bx); lo_y_dom_face.setBig(1,bx.smallEnd(1));
740  if (bc_ptr_h[BCVars::yvel_bc].lo(1) == ERFBCType::ext_dir) {
741  ParallelFor(lo_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
742  rho_v_rhs(i,j,k) = 0.;
743  });
744  } else if (bc_ptr_h[BCVars::yvel_bc].lo(1) == ERFBCType::ext_dir_upwind) {
745  ParallelFor(lo_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
746  if (v(i,j,k) >= 0.) {
747  rho_v_rhs(i,j,k) = 0.;
748  }
749  });
750  }
751  }
752  if (bx.bigEnd(1) == domain.bigEnd(1)) {
753  Box hi_y_dom_face(bx); hi_y_dom_face.setSmall(1,bx.bigEnd(1)+1); hi_y_dom_face.setBig(1,bx.bigEnd(1)+1);
754  if (bc_ptr_h[BCVars::yvel_bc].hi(1) == ERFBCType::ext_dir) {
755  ParallelFor(hi_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
756  rho_v_rhs(i,j,k) = 0.;
757  });
758  } else if (bc_ptr_h[BCVars::yvel_bc].hi(1) == ERFBCType::ext_dir_upwind) {
759  ParallelFor(hi_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
760  if (v(i,j,k) <= 0.) {
761  rho_v_rhs(i,j,k) = 0.;
762  }
763  });
764  }
765  }
766 
767  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
768  { // z-momentum equation
769 
770  Real gpz = gpz_arr(i,j,k);
771 
772  Real q = (l_use_moisture) ? 0.5 * (qt_arr(i,j,k) + qt_arr(i,j,k-1)) : 0.0;
773 
774  rho_w_rhs(i, j, k) += (-gpz - abl_pressure_grad[2] + buoyancy_arr(i,j,k)) / (1.0 + q) + zmom_src_arr(i,j,k);
775 
776  if (l_moving_terrain) {
777  rho_w_rhs(i, j, k) *= 0.5 * (detJ_arr(i,j,k) + detJ_arr(i,j,k-1));
778  }
779  });
780 
781  auto const lo = lbound(bx);
782  auto const hi = ubound(bx);
783 
784  // Note: the logic below assumes no tiling in z!
785  if (level > 0) {
786 
787  const Array4<const Real>& rho_w_rhs_crse = zmom_crse_rhs->const_array(mfi);
788 
789  Box b2d = bx; b2d.setRange(2,0);
790 
791  if (lo.z > klo) {
792  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) // bottom of box but not of domain
793  {
794  rho_w_rhs(i,j,lo.z) = rho_w_rhs_crse(i,j,lo.z);
795  });
796  }
797 
798  if (hi.z < khi+1) {
799  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) // top of box but not of domain
800  {
801  rho_w_rhs(i,j,hi.z+1) = rho_w_rhs_crse(i,j,hi.z+1);
802  });
803  }
804  }
805 
806  {
807  BL_PROFILE("slow_rhs_pre_fluxreg");
808  // We only add to the flux registers in the final RK step
809  // NOTE: for now we are only refluxing density not (rho theta) since the latter seems to introduce
810  // a problem at top and bottom boundaries
811  if (l_reflux) {
812  int strt_comp_reflux = (l_fixed_rho) ? 1 : 0;
813  int num_comp_reflux = 1;
814  if (level < finest_level) {
815  fr_as_crse->CrseAdd(mfi,
816  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
817  dx, dt, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
818  }
819  if (level > 0) {
820  fr_as_fine->FineAdd(mfi,
821  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
822  dx, dt, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
823  }
824 
825  // This is necessary here so we don't go on to the next FArrayBox without
826  // having finished copying the fluxes into the FluxRegisters (since the fluxes
827  // are stored in temporary FArrayBox's)
828  Gpu::streamSynchronize();
829 
830  } // two-way coupling
831  } // end profile
832  } // mfi
833  } // OMP
834 }
void AdvectionSrcForRho(const amrex::Box &bx, const amrex::Array4< amrex::Real > &src, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &omega, const amrex::Array4< amrex::Real > &avg_xmom, const amrex::Array4< amrex::Real > &avg_ymom, const amrex::Array4< amrex::Real > &avg_zmom, const amrex::Array4< const amrex::Real > &ax_arr, const amrex::Array4< const amrex::Real > &ay_arr, const amrex::Array4< const amrex::Real > &az_arr, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::GpuArray< const amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_arr, const bool fixed_rho)
void AdvectionSrcForScalars(const amrex::Box &bx, const int icomp, const int ncomp, const amrex::Array4< const amrex::Real > &avg_xmom, const amrex::Array4< const amrex::Real > &avg_ymom, const amrex::Array4< const amrex::Real > &avg_zmom, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< amrex::Real > &src, const amrex::Array4< const amrex::Real > &vf_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const AdvType horiz_adv_type, const AdvType vert_adv_type, const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac, const amrex::GpuArray< const amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_arr, const amrex::Box &domain, const amrex::BCRec *bc_ptr_h)
void AdvectionSrcForMom(const amrex::MFIter &mfi, const amrex::Box &bx, const amrex::Box &bxx, const amrex::Box &bxy, const amrex::Box &bxz, const amrex::Vector< amrex::Box > &bxx_grown, const amrex::Vector< amrex::Box > &bxy_grown, const amrex::Vector< amrex::Box > &bxz_grown, const amrex::Array4< amrex::Real > &rho_u_rhs, const amrex::Array4< amrex::Real > &rho_v_rhs, const amrex::Array4< amrex::Real > &rho_w_rhs, const amrex::Array4< const amrex::Real > &rho, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &Omega, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, amrex::Gpu::DeviceVector< amrex::Real > &stretched_dz_d, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vy, const AdvType horiz_adv_type, const AdvType vert_adv_type, const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac, MeshType &mesh_type, TerrainType &terrain_type, const eb_ &ebfact, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_u_arr, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_v_arr, amrex::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_w_arr, const amrex::Vector< amrex::iMultiFab > &physbnd_mask, const bool already_on_centroids, const int lo_z_face, const int hi_z_face, const amrex::Box &domain, const amrex::BCRec *bc_ptr_h)
@ tau12
Definition: ERF_DataStruct.H:30
@ tau23
Definition: ERF_DataStruct.H:30
@ tau33
Definition: ERF_DataStruct.H:30
@ tau22
Definition: ERF_DataStruct.H:30
@ tau11
Definition: ERF_DataStruct.H:30
@ tau32
Definition: ERF_DataStruct.H:30
@ tau31
Definition: ERF_DataStruct.H:30
@ tau21
Definition: ERF_DataStruct.H:30
@ tau13
Definition: ERF_DataStruct.H:30
@ nvars
Definition: ERF_DataStruct.H:91
@ v_x
Definition: ERF_DataStruct.H:22
@ u_y
Definition: ERF_DataStruct.H:23
@ v_y
Definition: ERF_DataStruct.H:23
@ m_y
Definition: ERF_DataStruct.H:23
@ u_x
Definition: ERF_DataStruct.H:22
@ m_x
Definition: ERF_DataStruct.H:22
void DiffusionSrcForState_S(const amrex::Box &bx, const amrex::Box &domain, int start_comp, int num_comp, const bool &rotate, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< amrex::Real > &cell_rhs, const amrex::Array4< amrex::Real > &xflux, const amrex::Array4< amrex::Real > &yflux, const amrex::Array4< amrex::Real > &zflux, const amrex::Gpu::DeviceVector< amrex::Real > &stretched_dz_d, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &SmnSmn_a, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vy, amrex::Array4< amrex::Real > &hfx_x, amrex::Array4< amrex::Real > &hfx_y, amrex::Array4< amrex::Real > &hfx_z, amrex::Array4< amrex::Real > &qfx1_x, amrex::Array4< amrex::Real > &qfx1_y, amrex::Array4< amrex::Real > &qfx1_z, amrex::Array4< amrex::Real > &qfx2_z, amrex::Array4< amrex::Real > &diss, const amrex::Array4< const amrex::Real > &mu_turb, const SolverChoice &solverChoice, const int level, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > grav_gpu, const amrex::BCRec *bc_ptr, const bool use_SurfLayer)
void DiffusionSrcForMom(const amrex::Box &bxx, const amrex::Box &bxy, const amrex::Box &bxz, const amrex::Array4< amrex::Real > &rho_u_rhs, const amrex::Array4< amrex::Real > &rho_v_rhs, const amrex::Array4< amrex::Real > &rho_w_rhs, const amrex::Array4< const amrex::Real > &tau11, const amrex::Array4< const amrex::Real > &tau22, const amrex::Array4< const amrex::Real > &tau33, const amrex::Array4< const amrex::Real > &tau12, const amrex::Array4< const amrex::Real > &tau21, const amrex::Array4< const amrex::Real > &tau13, const amrex::Array4< const amrex::Real > &tau31, const amrex::Array4< const amrex::Real > &tau23, const amrex::Array4< const amrex::Real > &tau32, const amrex::Array4< const amrex::Real > &detJ_arr, const amrex::Gpu::DeviceVector< amrex::Real > &stretched_dz_d, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vy, const bool use_stretched_dz, const bool use_variable_dz)
void DiffusionSrcForMom_EB(const amrex::MFIter &mfi, [[maybe_unused]] const amrex::Box &domain, const amrex::Box &bxx, const amrex::Box &bxy, const amrex::Box &bxz, const amrex::Array4< amrex::Real > &rho_u_rhs, const amrex::Array4< amrex::Real > &rho_v_rhs, const amrex::Array4< amrex::Real > &rho_w_rhs, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &w_arr, const amrex::Array4< const amrex::Real > &tau11, const amrex::Array4< const amrex::Real > &tau22, const amrex::Array4< const amrex::Real > &tau33, const amrex::Array4< const amrex::Real > &tau12, const amrex::Array4< const amrex::Real > &tau13, const amrex::Array4< const amrex::Real > &tau23, const amrex::Real *dx_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vy, const SolverChoice &solverChoice, const eb_ &ebfact, [[maybe_unused]] const amrex::BCRec *bc_ptr)
void DiffusionSrcForState_N(const amrex::Box &bx, const amrex::Box &domain, int start_comp, int num_comp, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< amrex::Real > &cell_rhs, const amrex::Array4< amrex::Real > &xflux, const amrex::Array4< amrex::Real > &yflux, const amrex::Array4< amrex::Real > &zflux, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &SmnSmn_a, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vy, amrex::Array4< amrex::Real > &hfx_z, amrex::Array4< amrex::Real > &qfx1_z, amrex::Array4< amrex::Real > &qfx2_z, amrex::Array4< amrex::Real > &diss, const amrex::Array4< const amrex::Real > &mu_turb, const SolverChoice &solverChoice, const int level, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > grav_gpu, const amrex::BCRec *bc_ptr, const bool use_SurfLayer)
void DiffusionSrcForState_T(const amrex::Box &bx, const amrex::Box &domain, int start_comp, int num_comp, const bool &rotate, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< amrex::Real > &cell_rhs, const amrex::Array4< amrex::Real > &xflux, const amrex::Array4< amrex::Real > &yflux, const amrex::Array4< amrex::Real > &zflux, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &z_cc, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &SmnSmn_a, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vy, amrex::Array4< amrex::Real > &hfx_x, amrex::Array4< amrex::Real > &hfx_y, amrex::Array4< amrex::Real > &hfx_z, amrex::Array4< amrex::Real > &qfx1_x, amrex::Array4< amrex::Real > &qfx1_y, amrex::Array4< amrex::Real > &qfx1_z, amrex::Array4< amrex::Real > &qfx2_z, amrex::Array4< amrex::Real > &diss, const amrex::Array4< const amrex::Real > &mu_turb, const SolverChoice &solverChoice, const int level, const amrex::Array4< const amrex::Real > &tm_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > grav_gpu, const amrex::BCRec *bc_ptr, const bool use_SurfLayer)
void EBAdvectionSrcForScalars(const amrex::Box &bx, const int icomp, const int ncomp, const amrex::Array4< const amrex::Real > &avg_xmom, const amrex::Array4< const amrex::Real > &avg_ymom, const amrex::Array4< const amrex::Real > &avg_zmom, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< amrex::Real > &src, const amrex::Array4< const int > &mask_arr, const amrex::Array4< const amrex::EBCellFlag > &cfg_arr, const amrex::Array4< const amrex::Real > &ax_arr, const amrex::Array4< const amrex::Real > &ay_arr, const amrex::Array4< const amrex::Real > &az_arr, const amrex::Array4< const amrex::Real > &fcx_arr, const amrex::Array4< const amrex::Real > &fcy_arr, const amrex::Array4< const amrex::Real > &fcz_arr, const amrex::Array4< const amrex::Real > &vf_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const AdvType horiz_adv_type, const AdvType vert_adv_type, const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac, const amrex::GpuArray< const amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_arr, const amrex::Box &domain, const amrex::BCRec *bc_ptr_h, bool already_on_centroids)
void EBAdvectionSrcForRho(const amrex::Box &bx, const amrex::Array4< amrex::Real > &src, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &omega, const amrex::Array4< amrex::Real > &avg_xmom, const amrex::Array4< amrex::Real > &avg_ymom, const amrex::Array4< amrex::Real > &avg_zmom, const amrex::Array4< const int > &mask_arr, const amrex::Array4< const amrex::EBCellFlag > &cfg_arr, const amrex::Array4< const amrex::Real > &ax_arr, const amrex::Array4< const amrex::Real > &ay_arr, const amrex::Array4< const amrex::Real > &az_arr, const amrex::Array4< const amrex::Real > &fcx_arr, const amrex::Array4< const amrex::Real > &fcy_arr, const amrex::Array4< const amrex::Real > &fcz_arr, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::GpuArray< const amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_arr, const bool fixed_rho, bool already_on_centroids)
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
AdvType
Definition: ERF_IndexDefines.H:221
void erf_make_tau_terms(int level, int nrk, const Vector< BCRec > &domain_bcs_type_h, const MultiFab &z_phys_nd, Vector< MultiFab > &S_data, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab *SmnSmn, MultiFab *eddyDiffs, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &, Gpu::DeviceVector< Real > &stretched_dz_d, const MultiFab &detJ, Vector< std::unique_ptr< MultiFab >> &mapfac)
Definition: ERF_MakeTauTerms.cpp:12
amrex::Real Real
Definition: ERF_ShocInterface.H:16
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real OmegaFromW(int &i, int &j, int &k, amrex::Real w, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:415
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:96
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:139
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
eb_aux_ const * get_w_const_factory() const noexcept
Definition: ERF_EB.H:56
const std::unique_ptr< amrex::EBFArrayBoxFactory > & get_const_factory() const noexcept
Definition: ERF_EB.H:46
eb_aux_ const * get_v_const_factory() const noexcept
Definition: ERF_EB.H:55
eb_aux_ const * get_u_const_factory() const noexcept
Definition: ERF_EB.H:54
@ yvel_bc
Definition: ERF_IndexDefines.H:88
@ xvel_bc
Definition: ERF_IndexDefines.H:87
@ ext_dir
Definition: ERF_IndexDefines.H:209
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:217
@ gpz
Definition: ERF_IndexDefines.H:152
@ gpy
Definition: ERF_IndexDefines.H:151
@ gpx
Definition: ERF_IndexDefines.H:150
@ NumTypes
Definition: ERF_IndexDefines.H:162
@ ymom
Definition: ERF_IndexDefines.H:160
@ cons
Definition: ERF_IndexDefines.H:158
@ zmom
Definition: ERF_IndexDefines.H:161
@ xmom
Definition: ERF_IndexDefines.H:159
@ qt
Definition: ERF_Kessler.H:27
@ xvel
Definition: ERF_IndexDefines.H:141
@ zvel
Definition: ERF_IndexDefines.H:143
@ yvel
Definition: ERF_IndexDefines.H:142
AdvType dycore_vert_adv_type
Definition: ERF_AdvStruct.H:396
amrex::Real dycore_vert_upw_frac
Definition: ERF_AdvStruct.H:406
AdvType dycore_horiz_adv_type
Definition: ERF_AdvStruct.H:395
amrex::Real dycore_horiz_upw_frac
Definition: ERF_AdvStruct.H:405
Definition: ERF_DiffStruct.H:19
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:84
static MeshType mesh_type
Definition: ERF_DataStruct.H:767
DiffChoice diffChoice
Definition: ERF_DataStruct.H:776
amrex::Real gravity
Definition: ERF_DataStruct.H:816
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_pressure_grad
Definition: ERF_DataStruct.H:870
bool fixed_density
Definition: ERF_DataStruct.H:785
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:778
amrex::Vector< int > anelastic
Definition: ERF_DataStruct.H:783
AdvChoice advChoice
Definition: ERF_DataStruct.H:775
MoistureType moisture_type
Definition: ERF_DataStruct.H:863
static TerrainType terrain_type
Definition: ERF_DataStruct.H:761
bool use_rotate_surface_flux
Definition: ERF_DataStruct.H:844
CouplingType coupling_type
Definition: ERF_DataStruct.H:862
Definition: ERF_TurbStruct.H:41
PBLType pbl_type
Definition: ERF_TurbStruct.H:370
bool use_keqn
Definition: ERF_TurbStruct.H:377
RANSType rans_type
Definition: ERF_TurbStruct.H:367
LESType les_type
Definition: ERF_TurbStruct.H:325
bool use_kturb
Definition: ERF_TurbStruct.H:376
Here is the call graph for this function: