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