ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_SlowRhsPost.cpp File Reference
#include <AMReX.H>
#include <ERF_SrcHeaders.H>
#include <ERF_TI_slow_headers.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< ABLMost > &most, 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 > &ax, std::unique_ptr< MultiFab > &ay, std::unique_ptr< MultiFab > &az, std::unique_ptr< MultiFab > &detJ, std::unique_ptr< MultiFab > &detJ_new, std::unique_ptr< MultiFab > &mapfac_m, std::unique_ptr< MultiFab > &mapfac_u, std::unique_ptr< MultiFab > &mapfac_v, 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< ABLMost > &  most,
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 > &  ax,
std::unique_ptr< MultiFab > &  ay,
std::unique_ptr< MultiFab > &  az,
std::unique_ptr< MultiFab > &  detJ,
std::unique_ptr< MultiFab > &  detJ_new,
std::unique_ptr< MultiFab > &  mapfac_m,
std::unique_ptr< MultiFab > &  mapfac_u,
std::unique_ptr< MultiFab > &  mapfac_v,
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]mostPointer to MOST 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]mapfac_mmap factor at cell centers
[in]mapfac_umap factor at x-faces
[in]mapfac_vmap factor at y-faces
[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 (most) t_mean_mf = most->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::OneWay);
114  if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain);
115 
116  const bool l_use_mono_adv = solverChoice.use_mono_adv;
117  const bool l_use_ddorf = (tc.les_type == LESType::Deardorff);
118  const bool l_use_KE = ( (tc.les_type == LESType::Deardorff) ||
119  (tc.pbl_type == PBLType::MYNN25) );
120  const bool l_advect_KE = (tc.use_KE && tc.advect_KE);
121  const bool l_use_diff = ((dc.molec_diff_type != MolecDiffType::None) ||
122  (tc.les_type != LESType::None) ||
123  (tc.pbl_type != PBLType::None) );
124  const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky ||
125  tc.les_type == LESType::Deardorff ||
126  tc.pbl_type == PBLType::MYNN25 ||
127  tc.pbl_type == PBLType::YSU );
128  const bool exp_most = (solverChoice.use_explicit_most);
129  const bool rot_most = (solverChoice.use_rotate_most);
130 
131  const Box& domain = geom.Domain();
132 
133  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
134  const Real* dx = geom.CellSize();
135 
136  // *************************************************************************
137  // Set gravity as a vector
138  // *************************************************************************
139  const Array<Real,AMREX_SPACEDIM> grav{0.0, 0.0, -solverChoice.gravity};
140  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
141 
142  // *************************************************************************
143  // Pre-computed quantities
144  // *************************************************************************
145  int nvars = S_data[IntVars::cons].nComp();
146  const BoxArray& ba = S_data[IntVars::cons].boxArray();
147  const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap();
148 
149  std::unique_ptr<MultiFab> dflux_x;
150  std::unique_ptr<MultiFab> dflux_y;
151  std::unique_ptr<MultiFab> dflux_z;
152 
153  if (l_use_diff) {
154  dflux_x = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)), dm, nvars, 0);
155  dflux_y = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)), dm, nvars, 0);
156  dflux_z = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)), dm, nvars, 0);
157  } else {
158  dflux_x = nullptr;
159  dflux_y = nullptr;
160  dflux_z = nullptr;
161  }
162 
163  // Valid vars
164  Vector<int> is_valid_slow_var; is_valid_slow_var.resize(RhoQ1_comp+1,0);
165  if (l_use_KE) {is_valid_slow_var[ RhoKE_comp] = 1;}
166  is_valid_slow_var[RhoScalar_comp] = 1;
167  if (solverChoice.moisture_type != MoistureType::None) {
168  is_valid_slow_var[RhoQ1_comp] = 1;
169  }
170 
171  // *****************************************************************************
172  // Monotonic advection for scalars
173  // *****************************************************************************
174  int nvar = S_new[IntVars::cons].nComp();
175  Vector<Real> max_scal(nvar, 1.0e34); Gpu::DeviceVector<Real> max_scal_d(nvar);
176  Vector<Real> min_scal(nvar,-1.0e34); Gpu::DeviceVector<Real> min_scal_d(nvar);
177  if (l_use_mono_adv) {
178  auto const& ma_s_arr = S_new[IntVars::cons].const_arrays();
179  for (int ivar(RhoKE_comp); ivar<nvar; ++ivar) {
180  GpuTuple<Real,Real> mm = ParReduce(TypeList<ReduceOpMax,ReduceOpMin>{},
181  TypeList<Real, Real>{},
182  S_new[IntVars::cons], IntVect(0),
183  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
184  -> GpuTuple<Real,Real>
185  {
186  return { ma_s_arr[box_no](i,j,k,ivar), ma_s_arr[box_no](i,j,k,ivar) };
187  });
188  max_scal[ivar] = get<0>(mm);
189  min_scal[ivar] = get<1>(mm);
190  }
191  }
192  Gpu::copy(Gpu::hostToDevice, max_scal.begin(), max_scal.end(), max_scal_d.begin());
193  Gpu::copy(Gpu::hostToDevice, min_scal.begin(), min_scal.end(), min_scal_d.begin());
194  Real* max_s_ptr = max_scal_d.data();
195  Real* min_s_ptr = min_scal_d.data();
196 
197  // *************************************************************************
198  // Calculate cell-centered eddy viscosity & diffusivities
199  //
200  // Notes -- we fill all the data in ghost cells before calling this so
201  // that we can fill the eddy viscosity in the ghost regions and
202  // not have to call a boundary filler on this data itself
203  //
204  // LES - updates both horizontal and vertical eddy viscosityS_tmp components
205  // PBL - only updates vertical eddy viscosity components so horizontal
206  // components come from the LES model or are left as zero.
207  // *************************************************************************
208 
209  // *************************************************************************
210  // Define updates and fluxes in the current RK stage
211  // *************************************************************************
212 #ifdef _OPENMP
213 #pragma omp parallel if (Gpu::notInLaunchRegion())
214 #endif
215  {
216  std::array<FArrayBox,AMREX_SPACEDIM> flux;
217  std::array<FArrayBox,AMREX_SPACEDIM> flux_tmp;
218 
219  int start_comp;
220  int num_comp;
221 
222  for ( MFIter mfi(S_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
223 
224  Box tbx = mfi.tilebox();
225 
226  // *************************************************************************
227  // Define flux arrays for use in advection
228  // *************************************************************************
229  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
230  flux[dir].resize(surroundingNodes(tbx,dir),nvars);
231  flux[dir].setVal<RunOn::Device>(0.);
232  if (l_use_mono_adv) {
233  flux_tmp[dir].resize(surroundingNodes(tbx,dir),nvars);
234  flux_tmp[dir].setVal<RunOn::Device>(0.);
235  }
236  }
237  const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
238  flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
239  Array4<Real> tmpx = (l_use_mono_adv) ? flux_tmp[0].array() : Array4<Real>{};
240  Array4<Real> tmpy = (l_use_mono_adv) ? flux_tmp[1].array() : Array4<Real>{};
241  Array4<Real> tmpz = (l_use_mono_adv) ? flux_tmp[2].array() : Array4<Real>{};
242  const GpuArray<Array4<Real>, AMREX_SPACEDIM> flx_tmp_arr{{AMREX_D_DECL(tmpx,tmpy,tmpz)}};
243 
244  // *************************************************************************
245  // Define Array4's
246  // *************************************************************************
247  const Array4<const Real> & old_cons = S_old[IntVars::cons].array(mfi);
248  const Array4< Real> & cell_rhs = S_rhs[IntVars::cons].array(mfi);
249 
250  const Array4< Real> & new_cons = S_new[IntVars::cons].array(mfi);
251  const Array4< Real> & new_xmom = S_new[IntVars::xmom].array(mfi);
252  const Array4< Real> & new_ymom = S_new[IntVars::ymom].array(mfi);
253  const Array4< Real> & new_zmom = S_new[IntVars::zmom].array(mfi);
254 
255  const Array4< Real> & cur_cons = S_data[IntVars::cons].array(mfi);
256  const Array4<const Real> & cur_prim = S_prim.array(mfi);
257  const Array4< Real> & cur_xmom = S_data[IntVars::xmom].array(mfi);
258  const Array4< Real> & cur_ymom = S_data[IntVars::ymom].array(mfi);
259  const Array4< Real> & cur_zmom = S_data[IntVars::zmom].array(mfi);
260 
261  Array4<Real> avg_xmom = S_scratch[IntVars::xmom].array(mfi);
262  Array4<Real> avg_ymom = S_scratch[IntVars::ymom].array(mfi);
263  Array4<Real> avg_zmom = S_scratch[IntVars::zmom].array(mfi);
264 
265  const Array4<const Real> & u = xvel.array(mfi);
266  const Array4<const Real> & v = yvel.array(mfi);
267 
268  const Array4<Real const>& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4<const Real>{};
269 
270  const Array4<const Real>& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4<const Real>{};
271  const Array4<const Real>& detJ_new_arr = l_moving_terrain ? detJ_new->const_array(mfi) : Array4<const Real>{};
272 
273  // Map factors
274  const Array4<const Real>& mf_m = mapfac_m->const_array(mfi);
275  const Array4<const Real>& mf_u = mapfac_u->const_array(mfi);
276  const Array4<const Real>& mf_v = mapfac_v->const_array(mfi);
277 
278  // SmnSmn for KE src with Deardorff
279  const Array4<const Real>& SmnSmn_a = l_use_ddorf ? SmnSmn->const_array(mfi) : Array4<const Real>{};
280 
281  // **************************************************************************
282  // Here we fill the "current" data with "new" data because that is the result of the previous RK stage
283  // **************************************************************************
284  int nsv = S_old[IntVars::cons].nComp() - 2;
285  const GpuArray<int, IntVars::NumTypes> scomp_slow = { 2,0,0,0};
286  const GpuArray<int, IntVars::NumTypes> ncomp_slow = {nsv,0,0,0};
287 
288  // **************************************************************************
289  // Note that here we do copy only the "slow" variables, not (rho) or (rho theta)
290  // **************************************************************************
291  ParallelFor(tbx, ncomp_slow[IntVars::cons],
292  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
293  const int n = scomp_slow[IntVars::cons] + nn;
294  cur_cons(i,j,k,n) = new_cons(i,j,k,n);
295  });
296 
297  // We have projected the velocities stored in S_data but we will use
298  // the velocities stored in S_scratch to update the scalars, so
299  // we need to copy from S_data (projected) into S_scratch
300  if (solverChoice.anelastic[level]) {
301  Box tbx_inc = mfi.nodaltilebox(0);
302  Box tby_inc = mfi.nodaltilebox(1);
303  Box tbz_inc = mfi.nodaltilebox(2);
304 
305  ParallelFor(tbx_inc, tby_inc, tbz_inc,
306  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
307  avg_xmom(i,j,k) = cur_xmom(i,j,k);
308  },
309  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
310  avg_ymom(i,j,k) = cur_ymom(i,j,k);
311  },
312  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
313  avg_zmom(i,j,k) = cur_zmom(i,j,k);
314  });
315  }
316 
317  // **************************************************************************
318  // Define updates in the RHS of continuity, temperature, and scalar equations
319  // **************************************************************************
320  // Metric terms
321  Array4<const Real> ax_arr;
322  Array4<const Real> ay_arr;
323  Array4<const Real> az_arr;
324  Array4<const Real> detJ_arr;
325  if (solverChoice.terrain_type == TerrainType::EB) {
326  ax_arr = ebfact.getAreaFrac()[0]->const_array(mfi);
327  ay_arr = ebfact.getAreaFrac()[1]->const_array(mfi);
328  az_arr = ebfact.getAreaFrac()[2]->const_array(mfi);
329  detJ_arr = ebfact.getVolFrac().const_array(mfi);
330  } else {
331  ax_arr = ax->const_array(mfi);
332  ay_arr = ay->const_array(mfi);
333  az_arr = az->const_array(mfi);
334  detJ_arr = detJ->const_array(mfi);
335  }
336 
337  AdvType horiz_adv_type, vert_adv_type;
338  Real horiz_upw_frac, vert_upw_frac;
339 
340  Array4<Real> diffflux_x, diffflux_y, diffflux_z;
341  Array4<Real> hfx_x, hfx_y, hfx_z, diss;
342  Array4<Real> q1fx_x, q1fx_y, q1fx_z, q2fx_z;
343  const bool use_most = (most != nullptr);
344 
345  if (l_use_diff) {
346  diffflux_x = dflux_x->array(mfi);
347  diffflux_y = dflux_y->array(mfi);
348  diffflux_z = dflux_z->array(mfi);
349 
350  hfx_x = Hfx1->array(mfi);
351  hfx_y = Hfx2->array(mfi);
352  hfx_z = Hfx3->array(mfi);
353  diss = Diss->array(mfi);
354 
355  if (Q1fx1) q1fx_x = Q1fx1->array(mfi);
356  if (Q1fx2) q1fx_y = Q1fx2->array(mfi);
357  if (Q1fx3) q1fx_z = Q1fx3->array(mfi);
358  if (Q2fx3) q2fx_z = Q2fx3->array(mfi);
359  }
360 
361  //
362  // Note that we either advect and diffuse all or none of the moisture variables
363  //
364  for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar)
365  {
366  if (is_valid_slow_var[ivar])
367  {
368  start_comp = ivar;
369 
370  if (ivar >= RhoQ1_comp) {
371  horiz_adv_type = ac.moistscal_horiz_adv_type;
372  vert_adv_type = ac.moistscal_vert_adv_type;
373  horiz_upw_frac = ac.moistscal_horiz_upw_frac;
374  vert_upw_frac = ac.moistscal_vert_upw_frac;
375 
376  if (ac.use_efficient_advection){
377  horiz_adv_type = EfficientAdvType(nrk,ac.moistscal_horiz_adv_type);
378  vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type);
379  }
380 
381  num_comp = n_qstate;
382 
383  } else {
384  horiz_adv_type = ac.dryscal_horiz_adv_type;
385  vert_adv_type = ac.dryscal_vert_adv_type;
386  horiz_upw_frac = ac.dryscal_horiz_upw_frac;
387  vert_upw_frac = ac.dryscal_vert_upw_frac;
388 
389  if (ac.use_efficient_advection){
390  horiz_adv_type = EfficientAdvType(nrk,ac.dryscal_horiz_adv_type);
391  vert_adv_type = EfficientAdvType(nrk,ac.dryscal_vert_adv_type);
392  }
393  num_comp = 1;
394  }
395 
396  if (( ivar != RhoKE_comp ) ||
397  ((ivar == RhoKE_comp) && l_advect_KE))
398  {
399  AdvectionSrcForScalars(dt, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom,
400  cur_cons, cur_prim, cell_rhs,
401  l_use_mono_adv, max_s_ptr, min_s_ptr,
402  detJ_arr, dxInv, mf_m,
403  horiz_adv_type, vert_adv_type,
404  horiz_upw_frac, vert_upw_frac,
405  flx_arr, flx_tmp_arr, domain, bc_ptr_h);
406  }
407 
408  if (l_use_diff) {
409  const Array4<const Real> tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4<const Real>{};
410  if (l_use_terrain) {
411  DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, exp_most, rot_most, u, v,
412  new_cons, cur_prim, cell_rhs,
413  diffflux_x, diffflux_y, diffflux_z,
414  z_nd, ax_arr, ay_arr, az_arr, detJ_arr,
415  dxInv, SmnSmn_a, mf_m, mf_u, mf_v,
416  hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z,q2fx_z, diss,
417  mu_turb, solverChoice, level,
418  tm_arr, grav_gpu, bc_ptr_d, use_most);
419  } else {
420  DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, exp_most, u, v,
421  new_cons, cur_prim, cell_rhs,
422  diffflux_x, diffflux_y, diffflux_z,
423  dxInv, SmnSmn_a, mf_m, mf_u, mf_v,
424  hfx_z, q1fx_z, q2fx_z, diss,
425  mu_turb, solverChoice, level,
426  tm_arr, grav_gpu, bc_ptr_d, use_most);
427  }
428  } // use_diff
429  } // valid slow var
430  } // loop ivar
431 
432 #if defined(ERF_USE_NETCDF)
433  if (moist_set_rhs_bool)
434  {
435  Box gtbx_moist = mfi.tilebox(IntVect(0),IntVect(2,2,0));
436  const Array4<const Real> & old_cons_const = S_old[IntVars::cons].const_array(mfi);
437  const Array4<const Real> & new_cons_const = S_new[IntVars::cons].const_array(mfi);
438  moist_set_rhs(tbx, gtbx_moist, old_cons_const, new_cons_const, cell_rhs,
439  bdy_time_interval, start_bdy_time, new_stage_time, dt, width, set_width, domain,
440  bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi);
441  }
442 #endif
443 
444  // This updates just the "slow" conserved variables
445  {
446  BL_PROFILE("rhs_post_8");
447 
448  const Real eps = std::numeric_limits<Real>::epsilon();
449 
450  auto const& src_arr = source.const_array(mfi);
451 
452  for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar)
453  {
454  if (is_valid_slow_var[ivar])
455  {
456  start_comp = ivar;
457 
458  if (ivar >= RhoQ1_comp) {
459  num_comp = nvars - RhoQ1_comp;
460  } else {
461  num_comp = 1;
462  }
463 
464  if (l_moving_terrain)
465  {
466  ParallelFor(tbx, num_comp,
467  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) noexcept {
468  const int n = start_comp + nn;
469  cell_rhs(i,j,k,n) += src_arr(i,j,k,n);
470  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);
471  cur_cons(i,j,k,n) = temp_val / detJ_new_arr(i,j,k);
472  if (ivar == RhoKE_comp) {
473  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps);
474  }
475  });
476 
477  } else {
478 
479  ParallelFor(tbx, num_comp,
480  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) noexcept {
481  const int n = start_comp + nn;
482  cell_rhs(i,j,k,n) += src_arr(i,j,k,n);
483  cur_cons(i,j,k,n) = old_cons(i,j,k,n) + dt * cell_rhs(i,j,k,n);
484  if (ivar == RhoKE_comp) {
485  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps);
486  } else if (ivar >= RhoQ1_comp) {
487  cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 0.0);
488  }
489  });
490 
491  } // moving or not?
492 
493  } // is_valid
494  } // ivar
495  } // profile
496 
497  {
498  BL_PROFILE("rhs_post_9");
499  // This updates all the conserved variables (not just the "slow" ones)
500  int num_comp_all = S_data[IntVars::cons].nComp();
501  ParallelFor(tbx, num_comp_all,
502  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept {
503  new_cons(i,j,k,n) = cur_cons(i,j,k,n);
504  });
505  } // end profile
506 
507  Box xtbx = mfi.nodaltilebox(0);
508  Box ytbx = mfi.nodaltilebox(1);
509  Box ztbx = mfi.nodaltilebox(2);
510 
511  {
512  BL_PROFILE("rhs_post_10()");
513  ParallelFor(xtbx, ytbx, ztbx,
514  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
515  new_xmom(i,j,k) = cur_xmom(i,j,k);
516  },
517  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
518  new_ymom(i,j,k) = cur_ymom(i,j,k);
519  },
520  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
521  new_zmom(i,j,k) = cur_zmom(i,j,k);
522  });
523  } // end profile
524 
525  {
526  BL_PROFILE("rhs_post_10");
527  // We only add to the flux registers in the final RK step
528  if (l_reflux && nrk == 2) {
529  int strt_comp_reflux = RhoTheta_comp + 1;
530  int num_comp_reflux = nvars - strt_comp_reflux;
531  if (level < finest_level) {
532  fr_as_crse->CrseAdd(mfi,
533  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
534  dx, dt, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
535  }
536  if (level > 0) {
537  fr_as_fine->FineAdd(mfi,
538  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
539  dx, dt, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
540  }
541 
542  // This is necessary here so we don't go on to the next FArrayBox without
543  // having finished copying the fluxes into the FluxRegisters (since the fluxes
544  // are stored in temporary FArrayBox's)
545  Gpu::streamSynchronize();
546 
547  } // two-way coupling
548  } // end profile
549  } // mfi
550  } // OMP
551 }
void AdvectionSrcForScalars(const amrex::Real &dt, 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 > &cur_cons, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< amrex::Real > &src, const bool &use_mono_adv, amrex::Real *max_s_ptr, amrex::Real *min_s_ptr, const amrex::Array4< const amrex::Real > &vf_arr, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_m, 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::GpuArray< amrex::Array4< amrex::Real >, AMREX_SPACEDIM > &flx_tmp_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:172
@ nvars
Definition: ERF_DataStruct.H:70
void DiffusionSrcForState_N(const amrex::Box &bx, const amrex::Box &domain, int start_comp, int num_comp, const bool &exp_most, 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_m, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, 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_most)
void DiffusionSrcForState_T(const amrex::Box &bx, const amrex::Box &domain, int start_comp, int num_comp, const bool &exp_most, const bool &rot_most, 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 > &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_m, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, 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_most)
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
AdvType
Definition: ERF_IndexDefines.H:191
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
@ ymom
Definition: ERF_IndexDefines.H:141
@ cons
Definition: ERF_IndexDefines.H:139
@ zmom
Definition: ERF_IndexDefines.H:142
@ xmom
Definition: ERF_IndexDefines.H:140
@ xvel
Definition: ERF_IndexDefines.H:130
@ yvel
Definition: ERF_IndexDefines.H:131
Definition: ERF_AdvStruct.H:19
amrex::Real dryscal_vert_upw_frac
Definition: ERF_AdvStruct.H:295
AdvType moistscal_horiz_adv_type
Definition: ERF_AdvStruct.H:286
AdvType moistscal_vert_adv_type
Definition: ERF_AdvStruct.H:287
amrex::Real moistscal_vert_upw_frac
Definition: ERF_AdvStruct.H:297
bool use_efficient_advection
Definition: ERF_AdvStruct.H:281
amrex::Real moistscal_horiz_upw_frac
Definition: ERF_AdvStruct.H:296
AdvType dryscal_horiz_adv_type
Definition: ERF_AdvStruct.H:284
AdvType dryscal_vert_adv_type
Definition: ERF_AdvStruct.H:285
amrex::Real dryscal_horiz_upw_frac
Definition: ERF_AdvStruct.H:294
Definition: ERF_DiffStruct.H:19
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:84
bool use_explicit_most
Definition: ERF_DataStruct.H:654
static MeshType mesh_type
Definition: ERF_DataStruct.H:579
bool use_mono_adv
Definition: ERF_DataStruct.H:670
DiffChoice diffChoice
Definition: ERF_DataStruct.H:588
amrex::Real gravity
Definition: ERF_DataStruct.H:626
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:590
amrex::Vector< int > anelastic
Definition: ERF_DataStruct.H:598
AdvChoice advChoice
Definition: ERF_DataStruct.H:587
MoistureType moisture_type
Definition: ERF_DataStruct.H:673
static TerrainType terrain_type
Definition: ERF_DataStruct.H:576
bool use_rotate_most
Definition: ERF_DataStruct.H:657
CouplingType coupling_type
Definition: ERF_DataStruct.H:672
Definition: ERF_TurbStruct.H:31
PBLType pbl_type
Definition: ERF_TurbStruct.H:205
bool use_KE
Definition: ERF_TurbStruct.H:219
LESType les_type
Definition: ERF_TurbStruct.H:176
bool advect_KE
Definition: ERF_TurbStruct.H:221
Here is the call graph for this function: