ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_SlowRhsPost.cpp File Reference
#include <AMReX.H>
#include <ERF_SrcHeaders.H>
#include <ERF_TI_slow_headers.H>
#include <ERF_EBAdvection.H>
#include <ERF_EBRedistribute.H>
Include dependency graph for ERF_SlowRhsPost.cpp:

Functions

void erf_slow_rhs_post (int level, int finest_level, int nrk, Real dt, int n_qstate, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_new, Vector< MultiFab > &S_data, const MultiFab &S_prim, Vector< MultiFab > &S_scratch, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &, const MultiFab &source, const MultiFab *SmnSmn, const 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, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &z_phys_cc, std::unique_ptr< MultiFab > &ax, std::unique_ptr< MultiFab > &ay, std::unique_ptr< MultiFab > &az, std::unique_ptr< MultiFab > &detJ, MultiFab *detJ_new, Gpu::DeviceVector< Real > &stretched_dz_d, Vector< std::unique_ptr< MultiFab >> &mapfac, amrex::EBFArrayBoxFactory const &ebfact, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine)
 

Function Documentation

◆ erf_slow_rhs_post()

void erf_slow_rhs_post ( int  level,
int  finest_level,
int  nrk,
Real  dt,
int  n_qstate,
Vector< MultiFab > &  S_rhs,
Vector< MultiFab > &  S_old,
Vector< MultiFab > &  S_new,
Vector< MultiFab > &  S_data,
const MultiFab &  S_prim,
Vector< MultiFab > &  S_scratch,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  ,
const MultiFab &  source,
const MultiFab *  SmnSmn,
const 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,
std::unique_ptr< MultiFab > &  z_phys_nd,
std::unique_ptr< MultiFab > &  z_phys_cc,
std::unique_ptr< MultiFab > &  ax,
std::unique_ptr< MultiFab > &  ay,
std::unique_ptr< MultiFab > &  az,
std::unique_ptr< MultiFab > &  detJ,
MultiFab *  detJ_new,
Gpu::DeviceVector< Real > &  stretched_dz_d,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
amrex::EBFArrayBoxFactory const &  ebfact,
YAFluxRegister *  fr_as_crse,
YAFluxRegister *  fr_as_fine 
)

Function for computing the slow RHS for the evolution equations for the scalars other than density or potential temperature

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_oldsolution at start of time step
[in]S_newsolution at end of current RK stage
[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]sourcesource terms for conserved variables
[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]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 at start of time step (= 1 if use_terrain is false)
[in]detJ_newJacobian of the metric transformation at new RK stage time (= 1 if use_terrain is false)
[in]mapfacmap factors
[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
98 {
99  BL_PROFILE_REGION("erf_slow_rhs_post()");
100 
101  const BCRec* bc_ptr_d = domain_bcs_type_d.data();
102  const BCRec* bc_ptr_h = domain_bcs_type_h.data();
103 
104  AdvChoice ac = solverChoice.advChoice;
105  DiffChoice dc = solverChoice.diffChoice;
106  TurbChoice tc = solverChoice.turbChoice[level];
107 
108  const MultiFab* t_mean_mf = nullptr;
109  if (SurfLayer) { t_mean_mf = SurfLayer->get_mac_avg(level,2); }
110 
111  const bool l_use_terrain = (solverChoice.mesh_type != MeshType::ConstantDz);
112  const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh);
113  const bool l_reflux = ( (solverChoice.coupling_type == CouplingType::TwoWay) && (nrk == 2) && (finest_level > 0) );
114  if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain);
115 
116  const bool l_anelastic = solverChoice.anelastic[level];
117 
118  const bool l_use_KE = ( tc.use_tke );
119  const bool l_need_SmnSmn = ( tc.les_type == LESType::Deardorff ||
120  tc.rans_type == RANSType::kEqn );
121  const bool l_advect_KE = ( tc.use_tke && tc.advect_tke );
122  const bool l_use_diff = ((dc.molec_diff_type != MolecDiffType::None) ||
123  (tc.les_type != LESType::None) ||
124  (tc.rans_type != RANSType::None) ||
125  (tc.pbl_type != PBLType::None) );
126  const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky ||
127  tc.les_type == LESType::Deardorff ||
128  tc.rans_type == RANSType::kEqn ||
129  tc.pbl_type == PBLType::MYJ ||
130  tc.pbl_type == PBLType::MYNN25 ||
131  tc.pbl_type == PBLType::MYNNEDMF ||
132  tc.pbl_type == PBLType::YSU ||
133  tc.pbl_type == PBLType::MRF );
134  const bool l_rotate = (solverChoice.use_rotate_surface_flux);
135 
136  const Box& domain = geom.Domain();
137 
138  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
139  const Real* dx = geom.CellSize();
140 
141  // *************************************************************************
142  // Set gravity as a vector
143  // *************************************************************************
144  const Array<Real,AMREX_SPACEDIM> grav{0.0, 0.0, -solverChoice.gravity};
145  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
146 
147  // *************************************************************************
148  // Pre-computed quantities
149  // *************************************************************************
150  int nvars = S_data[IntVars::cons].nComp();
151  const BoxArray& ba = S_data[IntVars::cons].boxArray();
152  const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap();
153 
154  std::unique_ptr<MultiFab> dflux_x;
155  std::unique_ptr<MultiFab> dflux_y;
156  std::unique_ptr<MultiFab> dflux_z;
157 
158  if (l_use_diff) {
159  dflux_x = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)), dm, 1, 0);
160  dflux_y = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)), dm, 1, 0);
161  dflux_z = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)), dm, 1, 0);
162  } else {
163  dflux_x = nullptr;
164  dflux_y = nullptr;
165  dflux_z = nullptr;
166  }
167 
168  // Valid vars
169  Vector<int> is_valid_slow_var; is_valid_slow_var.resize(RhoQ1_comp+1,0);
170  if (l_use_KE) {is_valid_slow_var[ RhoKE_comp] = 1;}
171  is_valid_slow_var[RhoScalar_comp] = 1;
172  if (solverChoice.moisture_type != MoistureType::None) {
173  is_valid_slow_var[RhoQ1_comp] = 1;
174  }
175 
176  // *************************************************************************
177  // Calculate cell-centered eddy viscosity & diffusivities
178  //
179  // Notes -- we fill all the data in ghost cells before calling this so
180  // that we can fill the eddy viscosity in the ghost regions and
181  // not have to call a boundary filler on this data itself
182  //
183  // LES - updates both horizontal and vertical eddy viscosityS_tmp components
184  // PBL - only updates vertical eddy viscosity components so horizontal
185  // components come from the LES model or are left as zero.
186  // *************************************************************************
187 
188  // *************************************************************************
189  // Define updates and fluxes in the current RK stage
190  // *************************************************************************
191 #ifdef _OPENMP
192 #pragma omp parallel if (Gpu::notInLaunchRegion())
193 #endif
194  {
195  std::array<FArrayBox,AMREX_SPACEDIM> flux;
196 
197  int start_comp;
198  int num_comp;
199 
200  // Cell-centered masks for EB (used for flux interpolation)
201  iMultiFab physbnd_mask;
202  bool already_on_centroids = false;
203  if (solverChoice.terrain_type == TerrainType::EB) {
204  physbnd_mask.define(S_data[IntVars::cons].boxArray(), S_data[IntVars::cons].DistributionMap(), 1, 1);
205  physbnd_mask.BuildMask(geom.Domain(), geom.periodicity(), 1, 1, 0, 1);
206  }
207 
208  for (MFIter mfi(S_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
209 
210  Box tbx = mfi.tilebox();
211 
212  // *************************************************************************
213  // Define flux arrays for use in advection
214  // *************************************************************************
215  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
216  if (solverChoice.terrain_type != TerrainType::EB) {
217  flux[dir].resize(surroundingNodes(tbx,dir),nvars);
218  } else {
219  flux[dir].resize(surroundingNodes(tbx,dir).grow(1),nvars);
220  }
221  flux[dir].setVal<RunOn::Device>(0.);
222  }
223  const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
224  flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
225 
226  // *************************************************************************
227  // Define Array4's
228  // *************************************************************************
229  const Array4<const Real> & old_cons = S_old[IntVars::cons].array(mfi);
230  const Array4< Real> & cell_rhs = S_rhs[IntVars::cons].array(mfi);
231 
232  const Array4< Real> & new_cons = S_new[IntVars::cons].array(mfi);
233  const Array4< Real> & new_xmom = S_new[IntVars::xmom].array(mfi);
234  const Array4< Real> & new_ymom = S_new[IntVars::ymom].array(mfi);
235  const Array4< Real> & new_zmom = S_new[IntVars::zmom].array(mfi);
236 
237  const Array4< Real> & cur_cons = S_data[IntVars::cons].array(mfi);
238  const Array4<const Real> & cur_prim = S_prim.array(mfi);
239  const Array4< Real> & cur_xmom = S_data[IntVars::xmom].array(mfi);
240  const Array4< Real> & cur_ymom = S_data[IntVars::ymom].array(mfi);
241  const Array4< Real> & cur_zmom = S_data[IntVars::zmom].array(mfi);
242 
243  Array4<Real> avg_xmom = S_scratch[IntVars::xmom].array(mfi);
244  Array4<Real> avg_ymom = S_scratch[IntVars::ymom].array(mfi);
245  Array4<Real> avg_zmom = S_scratch[IntVars::zmom].array(mfi);
246 
247  const Array4<const Real> & u = xvel.array(mfi);
248  const Array4<const Real> & v = yvel.array(mfi);
249 
250  const Array4<Real const>& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4<const Real>{};
251 
252  const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
253  const Array4<const Real>& z_cc = z_phys_cc->const_array(mfi);
254  const Array4<const Real>& detJ_new_arr = l_moving_terrain ? detJ_new->const_array(mfi) : Array4<const Real>{};
255 
256  // Map factors
257  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
258  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
259  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
260  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
261  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
262  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
263 
264  // SmnSmn for KE src with Deardorff or k-eqn RANS
265  const Array4<const Real>& SmnSmn_a = l_need_SmnSmn ? SmnSmn->const_array(mfi) : Array4<const Real>{};
266 
267  // **************************************************************************
268  // Here we fill the "current" data with "new" data because that is the result of the previous RK stage
269  // **************************************************************************
270  int nsv = S_old[IntVars::cons].nComp() - 2;
271  const GpuArray<int, IntVars::NumTypes> scomp_slow = { 2,0,0,0};
272  const GpuArray<int, IntVars::NumTypes> ncomp_slow = {nsv,0,0,0};
273 
274  // **************************************************************************
275  // Note that here we do copy only the "slow" variables, not (rho) or (rho theta)
276  // **************************************************************************
277  ParallelFor(tbx, ncomp_slow[IntVars::cons],
278  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
279  const int n = scomp_slow[IntVars::cons] + nn;
280  cur_cons(i,j,k,n) = new_cons(i,j,k,n);
281  });
282 
283  // We have projected the velocities stored in S_data but we will use
284  // the velocities stored in S_scratch to update the scalars, so
285  // we need to copy from S_data (projected) into S_scratch
286  if (l_anelastic) {
287  Box tbx_inc = mfi.nodaltilebox(0);
288  Box tby_inc = mfi.nodaltilebox(1);
289  Box tbz_inc = mfi.nodaltilebox(2);
290 
291  ParallelFor(tbx_inc, tby_inc, tbz_inc,
292  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
293  avg_xmom(i,j,k) = cur_xmom(i,j,k);
294  },
295  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
296  avg_ymom(i,j,k) = cur_ymom(i,j,k);
297  },
298  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
299  avg_zmom(i,j,k) = cur_zmom(i,j,k);
300  });
301  }
302 
303  // **************************************************************************
304  // Define updates in the RHS of continuity, temperature, and scalar equations
305  // **************************************************************************
306  bool l_eb_terrain_cc = false; // EB terrain on cell-centered grid
307  Array4<const int> mask_arr{};
308  Array4<const EBCellFlag> cfg_arr{};
309  Array4<const Real> ax_arr{};
310  Array4<const Real> ay_arr{};
311  Array4<const Real> az_arr{};
312  Array4<const Real> fcx_arr{};
313  Array4<const Real> fcy_arr{};
314  Array4<const Real> fcz_arr{};
315  Array4<const Real> detJ_arr{};
316  if (solverChoice.terrain_type == TerrainType::EB) {
317  EBCellFlagFab const& cfg = ebfact.getMultiEBCellFlagFab()[mfi];
318  cfg_arr = cfg.const_array();
319  if (cfg.getType(tbx) == FabType::singlevalued) {
320  l_eb_terrain_cc = true;
321  ax_arr = ebfact.getAreaFrac()[0]->const_array(mfi);
322  ay_arr = ebfact.getAreaFrac()[1]->const_array(mfi);
323  az_arr = ebfact.getAreaFrac()[2]->const_array(mfi);
324  fcx_arr = ebfact.getFaceCent()[0]->const_array(mfi);
325  fcy_arr = ebfact.getFaceCent()[1]->const_array(mfi);
326  fcz_arr = ebfact.getFaceCent()[2]->const_array(mfi);
327  detJ_arr = ebfact.getVolFrac().const_array(mfi);
328  // if (!already_on_centroids) {mask_arr = physbnd_mask.const_array(mfi);}
329  mask_arr = physbnd_mask.const_array(mfi);
330  }
331  }
332  if (!l_eb_terrain_cc) {
333  ax_arr = ax->const_array(mfi);
334  ay_arr = ay->const_array(mfi);
335  az_arr = az->const_array(mfi);
336  detJ_arr = detJ->const_array(mfi);
337  }
338 
339  AdvType horiz_adv_type, vert_adv_type;
340  Real horiz_upw_frac, vert_upw_frac;
341 
342  Array4<Real> diffflux_x, diffflux_y, diffflux_z;
343  Array4<Real> hfx_x, hfx_y, hfx_z, diss;
344  Array4<Real> q1fx_x, q1fx_y, q1fx_z, q2fx_z;
345  const bool use_SurfLayer = (SurfLayer != nullptr);
346 
347  if (l_use_diff) {
348  diffflux_x = dflux_x->array(mfi);
349  diffflux_y = dflux_y->array(mfi);
350  diffflux_z = dflux_z->array(mfi);
351 
352  hfx_x = Hfx1->array(mfi);
353  hfx_y = Hfx2->array(mfi);
354  hfx_z = Hfx3->array(mfi);
355  diss = Diss->array(mfi);
356 
357  if (Q1fx1) q1fx_x = Q1fx1->array(mfi);
358  if (Q1fx2) q1fx_y = Q1fx2->array(mfi);
359  if (Q1fx3) q1fx_z = Q1fx3->array(mfi);
360  if (Q2fx3) q2fx_z = Q2fx3->array(mfi);
361  }
362 
363  //
364  // Note that we either advect and diffuse all or none of the moisture variables
365  //
366  for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar)
367  {
368  if (is_valid_slow_var[ivar])
369  {
370  start_comp = ivar;
371  num_comp = 1;
372 
373  if (ivar == RhoQ1_comp) {
374  horiz_adv_type = ac.moistscal_horiz_adv_type;
375  vert_adv_type = ac.moistscal_vert_adv_type;
376  horiz_upw_frac = ac.moistscal_horiz_upw_frac;
377  vert_upw_frac = ac.moistscal_vert_upw_frac;
378 
379  if (ac.use_efficient_advection){
380  horiz_adv_type = EfficientAdvType(nrk,ac.moistscal_horiz_adv_type);
381  vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type);
382  }
383 
384  num_comp = n_qstate;
385 
386  } else {
387  horiz_adv_type = ac.dryscal_horiz_adv_type;
388  vert_adv_type = ac.dryscal_vert_adv_type;
389  horiz_upw_frac = ac.dryscal_horiz_upw_frac;
390  vert_upw_frac = ac.dryscal_vert_upw_frac;
391 
392  if (ac.use_efficient_advection){
393  horiz_adv_type = EfficientAdvType(nrk,ac.dryscal_horiz_adv_type);
394  vert_adv_type = EfficientAdvType(nrk,ac.dryscal_vert_adv_type);
395  }
396 
397  if (ivar == RhoScalar_comp) {
398  num_comp = NSCALARS;
399  }
400  }
401 
402  if (( ivar != RhoKE_comp ) ||
403  ((ivar == RhoKE_comp) && l_advect_KE))
404  {
405  if (!l_eb_terrain_cc){
406  AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom,
407  cur_prim, cell_rhs,
408  detJ_arr, dxInv, mf_mx, mf_my,
409  horiz_adv_type, vert_adv_type,
410  horiz_upw_frac, vert_upw_frac,
411  flx_arr, domain, bc_ptr_h);
412  } else {
413  EBAdvectionSrcForScalars(tbx, start_comp, num_comp,
414  avg_xmom, avg_ymom, avg_zmom,
415  cur_prim, cell_rhs,
416  mask_arr, cfg_arr, ax_arr, ay_arr, az_arr,
417  fcx_arr, fcy_arr, fcz_arr,
418  detJ_arr, dxInv, mf_mx, mf_my,
419  horiz_adv_type, vert_adv_type,
420  horiz_upw_frac, vert_upw_frac,
421  flx_arr, domain, bc_ptr_h,
422  already_on_centroids);
423  }
424  }
425 
426  if (l_use_diff) {
427  const Array4<const Real> tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4<const Real>{};
428  if (solverChoice.mesh_type == MeshType::StretchedDz && solverChoice.terrain_type != TerrainType::EB) {
429  DiffusionSrcForState_S(tbx, domain, start_comp, num_comp, l_rotate, u, v,
430  new_cons, cur_prim, cell_rhs,
431  diffflux_x, diffflux_y, diffflux_z,
432  stretched_dz_d, dxInv, SmnSmn_a,
433  mf_mx, mf_ux, mf_vx,
434  mf_my, mf_uy, mf_vy,
435  hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z,q2fx_z, diss,
436  mu_turb, solverChoice, level,
437  tm_arr, grav_gpu, bc_ptr_d, use_SurfLayer);
438  } else if (l_use_terrain) {
439  DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, l_rotate, u, v,
440  new_cons, cur_prim, cell_rhs,
441  diffflux_x, diffflux_y, diffflux_z,
442  z_nd, z_cc, ax_arr, ay_arr, az_arr,
443  detJ_arr, dxInv, SmnSmn_a,
444  mf_mx, mf_ux, mf_vx,
445  mf_my, mf_uy, mf_vy,
446  hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z,q2fx_z, diss,
447  mu_turb, solverChoice, level,
448  tm_arr, grav_gpu, bc_ptr_d, use_SurfLayer);
449  } else {
450  DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, u, v,
451  new_cons, cur_prim, cell_rhs,
452  diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a,
453  mf_mx, mf_ux, mf_vx,
454  mf_my, mf_uy, mf_vy,
455  hfx_z, q1fx_z, q2fx_z, diss,
456  mu_turb, solverChoice, level,
457  tm_arr, grav_gpu, bc_ptr_d, use_SurfLayer);
458  }
459  } // use_diff
460  } // valid slow var
461  } // loop ivar
462 
463 #if defined(ERF_USE_NETCDF)
464  if (moist_set_rhs_bool)
465  {
466  const Array4<const Real> & old_cons_const = S_old[IntVars::cons].const_array(mfi);
467  const Array4<const Real> & new_cons_const = S_new[IntVars::cons].const_array(mfi);
468  moist_set_rhs(geom, tbx, old_cons_const, new_cons_const, cell_rhs, bdy_time_interval,
469  new_stage_time, dt, stop_time_elapsed, width, set_width, domain,
470  bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi);
471  }
472 #endif
473 
474  // This updates just the "slow" conserved variables
475  {
476  BL_PROFILE("rhs_post_8");
477 
479 
480  auto const& src_arr = source.const_array(mfi);
481 
482  for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar)
483  {
484  if (is_valid_slow_var[ivar])
485  {
486  start_comp = ivar;
487  num_comp = 1;
488  if (ivar == RhoQ1_comp) {
489  num_comp = nvars - RhoQ1_comp;
490  } else if (ivar == RhoScalar_comp) {
491  num_comp = NSCALARS;
492  }
493 
494  if (l_moving_terrain)
495  {
496  ParallelFor(tbx, num_comp,
497  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) noexcept {
498  const int n = start_comp + nn;
499  cell_rhs(i,j,k,n) += src_arr(i,j,k,n);
500  Real temp_val = detJ_arr(i,j,k) * old_cons(i,j,k,n) + dt * detJ_arr(i,j,k) * cell_rhs(i,j,k,n);
501  cur_cons(i,j,k,n) = temp_val / detJ_new_arr(i,j,k);
502  if (ivar == RhoKE_comp) {
503  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps);
504  }
505  });
506 
507  } else if (l_anelastic && (nrk == 1)) { // not moving and ( (anelastic) and second RK stage) )
508 
509  ParallelFor(tbx, num_comp,
510  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) noexcept {
511  const int n = start_comp + nn;
512  cell_rhs(i,j,k,n) += src_arr(i,j,k,n);
513 
514  // Re-construct the cell_rhs used in the first RK stage
515  Real dt_times_old_cell_rhs = cur_cons(i,j,k,n) - old_cons(i,j,k,n);
516 
517  // Add the time-averaged RHS to the old state
518  cur_cons(i,j,k,n) = old_cons(i,j,k,n) + 0.5 * (dt_times_old_cell_rhs + dt * cell_rhs(i,j,k,n));
519 
520  if (ivar == RhoKE_comp) {
521  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps);
522  } else if (ivar >= RhoQ1_comp) {
523  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 0.0);
524  }
525  });
526 
527  } else { // not moving and ( (not anelastic) or (first RK stage) )
528 
529  ParallelFor(tbx, num_comp,
530  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) noexcept {
531  const int n = start_comp + nn;
532  cell_rhs(i,j,k,n) += src_arr(i,j,k,n);
533  cur_cons(i,j,k,n) = old_cons(i,j,k,n) + dt * cell_rhs(i,j,k,n);
534  if (ivar == RhoKE_comp) {
535  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps);
536  } else if (ivar >= RhoQ1_comp) {
537  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 0.0);
538  }
539  });
540 
541  } // moving, anelastic or neither?
542 
543  } // is_valid
544  } // ivar
545  } // profile
546 
547  {
548  BL_PROFILE("rhs_post_9");
549  // This updates all the conserved variables (not just the "slow" ones)
550  int num_comp_all = S_data[IntVars::cons].nComp();
551  ParallelFor(tbx, num_comp_all,
552  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept {
553  new_cons(i,j,k,n) = cur_cons(i,j,k,n);
554  });
555  } // end profile
556 
557  Box xtbx = mfi.nodaltilebox(0);
558  Box ytbx = mfi.nodaltilebox(1);
559  Box ztbx = mfi.nodaltilebox(2);
560 
561  {
562  BL_PROFILE("rhs_post_10()");
563  ParallelFor(xtbx, ytbx, ztbx,
564  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
565  new_xmom(i,j,k) = cur_xmom(i,j,k);
566  },
567  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
568  new_ymom(i,j,k) = cur_ymom(i,j,k);
569  },
570  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
571  new_zmom(i,j,k) = cur_zmom(i,j,k);
572  });
573  } // end profile
574 
575  {
576  BL_PROFILE("rhs_post_10");
577  // We only add to the flux registers in the final RK step
578  if (l_reflux) {
579  int strt_comp_reflux = RhoTheta_comp + 1;
580  int num_comp_reflux = nvars - strt_comp_reflux;
581  if (level < finest_level) {
582  fr_as_crse->CrseAdd(mfi,
583  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
584  dx, dt, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
585  }
586  if (level > 0) {
587  fr_as_fine->FineAdd(mfi,
588  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
589  dx, dt, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
590  }
591 
592  // This is necessary here so we don't go on to the next FArrayBox without
593  // having finished copying the fluxes into the FluxRegisters (since the fluxes
594  // are stored in temporary FArrayBox's)
595  Gpu::streamSynchronize();
596 
597  } // two-way coupling
598  } // end profile
599  } // mfi
600  } // OMP
601 }
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)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE AdvType EfficientAdvType(int nrk, AdvType adv_type)
Definition: ERF_Advection.H:281
@ 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 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)
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define NSCALARS
Definition: ERF_IndexDefines.H:16
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
AdvType
Definition: ERF_IndexDefines.H:221
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
amrex::Real Real
Definition: ERF_ShocInterface.H:19
@ ymom
Definition: ERF_IndexDefines.H:160
@ cons
Definition: ERF_IndexDefines.H:158
@ zmom
Definition: ERF_IndexDefines.H:161
@ xmom
Definition: ERF_IndexDefines.H:159
@ xvel
Definition: ERF_IndexDefines.H:141
@ yvel
Definition: ERF_IndexDefines.H:142
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
real(c_double), private ac
Definition: ERF_module_mp_morr_two_moment.F90:180
Definition: ERF_AdvStruct.H:19
Definition: ERF_DiffStruct.H:19
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:84
static MeshType mesh_type
Definition: ERF_DataStruct.H:815
DiffChoice diffChoice
Definition: ERF_DataStruct.H:824
amrex::Real gravity
Definition: ERF_DataStruct.H:864
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:826
amrex::Vector< int > anelastic
Definition: ERF_DataStruct.H:831
AdvChoice advChoice
Definition: ERF_DataStruct.H:823
MoistureType moisture_type
Definition: ERF_DataStruct.H:911
static TerrainType terrain_type
Definition: ERF_DataStruct.H:806
bool use_rotate_surface_flux
Definition: ERF_DataStruct.H:892
CouplingType coupling_type
Definition: ERF_DataStruct.H:910
Definition: ERF_TurbStruct.H:41
PBLType pbl_type
Definition: ERF_TurbStruct.H:375
RANSType rans_type
Definition: ERF_TurbStruct.H:372
bool advect_tke
Definition: ERF_TurbStruct.H:415
bool use_tke
Definition: ERF_TurbStruct.H:387
LESType les_type
Definition: ERF_TurbStruct.H:330
Here is the call graph for this function: