ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Substep_T.cpp File Reference
#include <ERF_TI_fast_headers.H>
#include "ERF_Constants.H"
Include dependency graph for ERF_Substep_T.cpp:

Functions

void erf_substep_T (int step, int, int level, int finest_level, Vector< MultiFab > &S_slow_rhs, const Vector< MultiFab > &S_prev, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &qt, const MultiFab &pi_stage, const MultiFab &fast_coeffs, Vector< MultiFab > &S_data, MultiFab &lagged_delta_rt, MultiFab &avg_xmom, MultiFab &avg_ymom, MultiFab &avg_zmom, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const Geometry geom, const Real gravity, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &detJ_cc, const double dtau_d, const Real beta_s, const Real facinv, Vector< std::unique_ptr< MultiFab >> &mapfac, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine, bool l_use_moisture, bool l_reflux, bool l_real_bc, const Real *sinesq_stag_d, const Real l_damp_coef)
 

Function Documentation

◆ erf_substep_T()

void erf_substep_T ( int  step,
int  ,
int  level,
int  finest_level,
Vector< MultiFab > &  S_slow_rhs,
const Vector< MultiFab > &  S_prev,
Vector< MultiFab > &  S_stage_data,
const MultiFab &  S_stage_prim,
const MultiFab &  qt,
const MultiFab &  pi_stage,
const MultiFab &  fast_coeffs,
Vector< MultiFab > &  S_data,
MultiFab &  lagged_delta_rt,
MultiFab &  avg_xmom,
MultiFab &  avg_ymom,
MultiFab &  avg_zmom,
const MultiFab &  cc_src,
const MultiFab &  xmom_src,
const MultiFab &  ymom_src,
const MultiFab &  zmom_src,
const Geometry  geom,
const Real  gravity,
std::unique_ptr< MultiFab > &  z_phys_nd,
std::unique_ptr< MultiFab > &  detJ_cc,
const double  dtau_d,
const Real  beta_s,
const Real  facinv,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
YAFluxRegister *  fr_as_crse,
YAFluxRegister *  fr_as_fine,
bool  l_use_moisture,
bool  l_reflux,
bool  l_real_bc,
const Real sinesq_stag_d,
const Real  l_damp_coef 
)

Function for computing the fast RHS with fixed-in-time terrain

Parameters
[in]stepwhich fast time step within each Runge-Kutta step
[in]nrkwhich Runge-Kutta step
[in]levellevel of resolution
[in]finest_levelfinest level of resolution
[in]S_slow_rhsslow RHS computed in erf_slow_rhs_pre
[in]S_prevprevious solution
[in]S_stage_datasolution at previous RK stage
[in]S_stage_primprimitive variables at previous RK stage
[in]pi_stageExner function at previous RK stage
[in]fast_coeffscoefficients for the tridiagonal solve used in the fast integrator
[out]S_datacurrent solution
[in,out]lagged_delta_rt
[in,out]avg_xmomtime-averaged x-momentum to be used for updating slow variables
[in,out]avg_ymomtime-averaged y-momentum to be used for updating slow variables
[in,out]avg_zmomtime-averaged z-momentum to be used for updating slow variables
[in]cc_srcsource terms for conserved variables
[in]xmom_srcsource terms for x-momentum
[in]ymom_srcsource terms for y-momentum
[in]zmom_srcsource terms for z-momentum
[in]geomcontainer for geometric information
[in]gravitymagnitude of gravity
[in]z_phys_ndheight coordinate at nodes
[in]detJ_ccJacobian of the metric transformation
[in]dtaufast time step
[in]beta_sCoefficient which determines how implicit vs explicit the solve is
[in]facinvinverse factor for time-averaging the momenta
[in]mapfacvector of map 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
[in]l_use_moisture
[in]l_refluxshould we add fluxes to the FluxRegisters?
[in]l_damp_coef
76 {
77  BL_PROFILE_REGION("erf_substep_T()");
78 
79  Real dtau = static_cast<Real>(dtau_d);
80 
81  const Box& domain = geom.Domain();
82  auto const domlo = lbound(domain);
83  auto const domhi = ubound(domain);
84 
85  int ilo = domlo.x;
86  int ihi = domhi.x + 1;
87  int jlo = domlo.y;
88  int jhi = domhi.y + 1;
89 
90  Real beta_1 = myhalf * (one - beta_s); // multiplies explicit terms
91  Real beta_2 = myhalf * (one + beta_s); // multiplies implicit terms
92 
93  // How much do we project forward the (rho theta) that is used in the horizontal momentum equations
94  Real beta_d = Real(0.1);
95 
96  Real RvOverRd = R_v / R_d;
97 
98  bool l_rayleigh_impl_for_w = (sinesq_stag_d != nullptr);
99 
100  const Real* dx = geom.CellSize();
101  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
102 
103  Real dxi = dxInv[0];
104  Real dyi = dxInv[1];
105  Real dzi = dxInv[2];
106  const auto& ba = S_stage_data[IntVars::cons].boxArray();
107  const auto& dm = S_stage_data[IntVars::cons].DistributionMap();
108 
109  MultiFab Delta_rho_u( convert(ba,IntVect(1,0,0)), dm, 1, 1);
110  MultiFab Delta_rho_v( convert(ba,IntVect(0,1,0)), dm, 1, 1);
111  MultiFab Delta_rho_w( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,0));
112  MultiFab Delta_rho ( ba , dm, 1, 1);
113  MultiFab Delta_rho_theta( ba , dm, 1, 1);
114 
115  MultiFab New_rho_u(convert(ba,IntVect(1,0,0)), dm, 1, 1);
116  MultiFab New_rho_v(convert(ba,IntVect(0,1,0)), dm, 1, 1);
117 
118  MultiFab coeff_A_mf(fast_coeffs, make_alias, 0, 1);
119  MultiFab inv_coeff_B_mf(fast_coeffs, make_alias, 1, 1);
120  MultiFab coeff_C_mf(fast_coeffs, make_alias, 2, 1);
121  MultiFab coeff_P_mf(fast_coeffs, make_alias, 3, 1);
122  MultiFab coeff_Q_mf(fast_coeffs, make_alias, 4, 1);
123 
124  // *************************************************************************
125  // Set gravity as a vector
126  const Array<Real,AMREX_SPACEDIM> grav{zero, zero, -gravity};
127  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
128 
129  MultiFab extrap(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(),1,1);
130 
131  // *************************************************************************
132  // First set up some arrays we'll need
133  // *************************************************************************
134 
135 #ifdef _OPENMP
136 #pragma omp parallel if (Gpu::notInLaunchRegion())
137 #endif
138  for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
139  {
140  const Array4<Real> & cur_cons = S_data[IntVars::cons].array(mfi);
141  const Array4<const Real>& prev_cons = S_prev[IntVars::cons].const_array(mfi);
142  const Array4<const Real>& stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
143  const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
144 
145  const Array4<Real>& old_drho = Delta_rho.array(mfi);
146  const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
147  const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
148  const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
149  const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
150 
151  const Array4<const Real>& prev_xmom = S_prev[IntVars::xmom].const_array(mfi);
152  const Array4<const Real>& prev_ymom = S_prev[IntVars::ymom].const_array(mfi);
153  const Array4<const Real>& prev_zmom = S_prev[IntVars::zmom].const_array(mfi);
154 
155  const Array4<const Real>& stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
156  const Array4<const Real>& stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
157  const Array4<const Real>& stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi);
158 
159  Box bx = mfi.validbox();
160  //
161  // NOTE: growntilebox / grownnodaltilebox, not the tilebox grown by hand. Growing by
162  // hand pushes every tile one cell past its own share of the grid, so once the
163  // grid is tiled two tiles store to the same cells -- a data race under OpenMP,
164  // even though both store the same value. The union of the tiles, and the values
165  // stored, are unchanged, and this is identical to the old expression when there
166  // is one tile per grid.
167  //
168  Box gbx = mfi.growntilebox(1);
169 
170  if (step == 0) {
171  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
172  cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp);
173  cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp);
174  });
175  } // step = 0
176 
177  Box gtbx = mfi.grownnodaltilebox(0,IntVect(1,1,0));
178  Box gtby = mfi.grownnodaltilebox(1,IntVect(1,1,0));
179  Box gtbz = mfi.grownnodaltilebox(2,IntVect(1,1,0));
180 
181  const auto& bx_lo = lbound(bx);
182  const auto& bx_hi = ubound(bx);
183 
184  ParallelFor(gtbx, gtby, gtbz,
185  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
186  old_drho_u(i,j,k) = prev_xmom(i,j,k) - stage_xmom(i,j,k);
187  if (k == bx_lo.z && k != domlo.z) {
188  old_drho_u(i,j,k-1) = old_drho_u(i,j,k);
189  } else if (k == bx_hi.z) {
190  old_drho_u(i,j,k+1) = old_drho_u(i,j,k);
191  }
192  },
193  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
194  old_drho_v(i,j,k) = prev_ymom(i,j,k) - stage_ymom(i,j,k);
195  if (k == bx_lo.z && k != domlo.z) {
196  old_drho_v(i,j,k-1) = old_drho_v(i,j,k);
197  } else if (k == bx_hi.z) {
198  old_drho_v(i,j,k+1) = old_drho_v(i,j,k);
199  }
200  },
201  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
202  old_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k);
203  });
204 
205  const Array4<Real>& theta_extrap = extrap.array(mfi);
206  const Array4<const Real>& prim = S_stage_prim.const_array(mfi);
207 
208  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
209  old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp);
210  old_drho_theta(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) - stage_cons(i,j,k,RhoTheta_comp);
211  if (step == 0) {
212  theta_extrap(i,j,k) = old_drho_theta(i,j,k);
213  } else {
214  theta_extrap(i,j,k) = old_drho_theta(i,j,k) + beta_d *
215  ( old_drho_theta(i,j,k) - lagged_arr(i,j,k) );
216  }
217 
218  // NOTE: qv is not changing over the fast steps so we use the stage data
219  Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero;
220  theta_extrap(i,j,k) *= (one + RvOverRd*qv);
221  });
222  } // mfi
223 
224 #ifdef _OPENMP
225 #pragma omp parallel if (Gpu::notInLaunchRegion())
226 #endif
227  for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
228  {
229  // We define lagged_delta_rt for our next step as the current delta_rt
230  // (growntilebox, not a hand-grown tilebox -- see the note above)
231  Box gbx = mfi.growntilebox(1);
232  const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
233  const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
234  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
235  lagged_arr(i,j,k) = old_drho_theta(i,j,k);
236  });
237  } // mfi
238 
239  // *************************************************************************
240  // Define updates in the current RK stage
241  // *************************************************************************
242 
243 #ifdef _OPENMP
244 #pragma omp parallel if (Gpu::notInLaunchRegion())
245 #endif
246  for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
247  {
248  Box bx = mfi.validbox();
249  Box tbx = mfi.nodaltilebox(0);
250  Box tby = mfi.nodaltilebox(1);
251 
252  const Array4<Real const>& xmom_src_arr = xmom_src.const_array(mfi);
253  const Array4<Real const>& ymom_src_arr = ymom_src.const_array(mfi);
254 
255  const Array4<const Real> & stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
256  const Array4<const Real> & stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
257  const Array4<const Real> & qt_arr = qt.const_array(mfi);
258 
259  const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
260  const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
261 
262  const Array4<const Real>& slow_rhs_rho_u = S_slow_rhs[IntVars::xmom].const_array(mfi);
263  const Array4<const Real>& slow_rhs_rho_v = S_slow_rhs[IntVars::ymom].const_array(mfi);
264 
265  const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
266  const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
267 
268  const Array4<Real>& cur_xmom = S_data[IntVars::xmom].array(mfi);
269  const Array4<Real>& cur_ymom = S_data[IntVars::ymom].array(mfi);
270 
271  // These store the advection momenta which we will use to update the slow variables
272  const Array4<Real>& avg_xmom_arr = avg_xmom.array(mfi);
273  const Array4<Real>& avg_ymom_arr = avg_ymom.array(mfi);
274 
275  const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
276 
277  const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
278 
279  const Array4<Real>& theta_extrap = extrap.array(mfi);
280 
281  // Map factors
282  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
283  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
284  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
285  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
286 
287  // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs
288  // Note that we do the Copy and Subtract including one ghost cell
289  // so that we don't have to fill ghost cells of the new MultiFabs
290  // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that
291  // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w
292 
293  // *********************************************************************
294  // Define updates in the RHS of {x, y, z}-momentum equations
295  // *********************************************************************
296  {
297  BL_PROFILE("substep_xymom_T");
298 
299  const auto& bx_lo = lbound(bx);
300  const auto& bx_hi = ubound(bx);
301 
302  ParallelFor(tbx, tby,
303  [=] AMREX_GPU_DEVICE (int i, int j, int k)
304  {
305  // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
306  Real met_h_xi = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd);
307  Real met_h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd);
308  Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi;
309  Real gp_zeta_on_iface = (k == 0) ?
310  myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
311  - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) :
312  fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
313  - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) );
314  Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? Real(0.) :
315  gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface;
316 
317  gpx *= mf_ux(i,j,0);
318 
319  Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) : zero;
320 
321  Real pi_c = myhalf * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i ,j,k,0));
322  Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx / (one + q);
323 
324  new_drho_u(i, j, k) = old_drho_u(i,j,k) + dtau * fast_rhs_rho_u
325  + dtau * slow_rhs_rho_u(i,j,k)
326  + dtau * xmom_src_arr(i,j,k);
327  if (k == bx_lo.z && k != domlo.z) {
328  new_drho_u(i,j,k-1) = new_drho_u(i,j,k);
329  } else if (k == bx_hi.z) {
330  new_drho_u(i,j,k+1) = new_drho_u(i,j,k);
331  }
332 
333  // NOTE: met_h_zeta here is identically the x-face area ax computed by
334  // make_areas, so this matches the base value of avg_xmom defined in
335  // AdvectionSrcForRho (ax*rho_u/mf_uy) as well as the density flux
336  // (new_drho_u*h_zeta_cc_xface/mf_uy) formed below.
337  avg_xmom_arr(i,j,k) += facinv * new_drho_u(i,j,k) * met_h_zeta / mf_uy(i,j,0);
338 
339  cur_xmom(i,j,k) = stage_xmom(i,j,k) + new_drho_u(i,j,k);
340  },
341  [=] AMREX_GPU_DEVICE (int i, int j, int k)
342  {
343  // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
344  Real met_h_eta = Compute_h_eta_AtJface(i, j, k, dxInv, z_nd);
345  Real met_h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd);
346  Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi;
347  Real gp_zeta_on_jface = (k == 0) ?
348  myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
349  - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) :
350  fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
351  - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) );
352  Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? Real(0.) :
353  gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface;
354 
355  gpy *= mf_vy(i,j,0);
356 
357  Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) : zero;
358 
359  Real pi_c = myhalf * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j ,k,0));
360  Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy / (one + q);
361 
362  new_drho_v(i, j, k) = old_drho_v(i,j,k) + dtau * fast_rhs_rho_v
363  + dtau * slow_rhs_rho_v(i,j,k)
364  + dtau * ymom_src_arr(i,j,k);
365 
366  if (k == bx_lo.z && k != domlo.z) {
367  new_drho_v(i,j,k-1) = new_drho_v(i,j,k);
368  } else if (k == bx_hi.z) {
369  new_drho_v(i,j,k+1) = new_drho_v(i,j,k);
370  }
371 
372  // NOTE: met_h_zeta here is identically the y-face area ay computed by
373  // make_areas, so this matches the base value of avg_ymom defined in
374  // AdvectionSrcForRho (ay*rho_v/mf_vx) as well as the density flux
375  // (new_drho_v*h_zeta_cc_yface/mf_vx) formed below.
376  avg_ymom_arr(i,j,k) += facinv * new_drho_v(i,j,k) * met_h_zeta / mf_vx(i,j,0);
377 
378  cur_ymom(i,j,k) = stage_ymom(i,j,k) + new_drho_v(i,j,k);
379  });
380  } // end profile
381  }
382 
383  MultiFab Omega(S_data[IntVars::zmom].boxArray(), dm, 1, 1);
384 
385 #ifdef _OPENMP
386 #pragma omp parallel if (Gpu::notInLaunchRegion())
387 #endif
388  {
389  std::array<FArrayBox,AMREX_SPACEDIM> flux;
390  for ( MFIter mfi(S_stage_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
391  {
392  Box bx = mfi.tilebox();
393  Box tbz = surroundingNodes(bx,2);
394 
395  Box vbx = mfi.validbox();
396  const auto& vbx_hi = ubound(vbx);
397 
398  const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
399  const Array4<Real const>& cc_src_arr = cc_src.const_array(mfi);
400 
401  const Array4<const Real> & stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi);
402  const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
403 
404  const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
405  const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
406  const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
407  const Array4<Real>& old_drho = Delta_rho.array(mfi);
408  const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
409 
410  const Array4<const Real>& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi);
411  const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi);
412 
413  const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
414  const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
415 
416  const Array4<Real>& cur_cons = S_data[IntVars::cons].array(mfi);
417  const Array4<Real>& cur_zmom = S_data[IntVars::zmom].array(mfi);
418 
419  // These store the advection momenta which we will use to update the slow variables
420  const Array4<Real>& avg_zmom_arr = avg_zmom.array(mfi);
421 
422  const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
423  const Array4<const Real>& detJ = detJ_cc->const_array(mfi);
424 
425  const Array4< Real>& omega_arr = Omega.array(mfi);
426 
427  // Map factors
428  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
429  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
430  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
431  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
432  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
433  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
434 
435  // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs
436  // Note that we do the Copy and Subtract including one ghost cell
437  // so that we don't have to fill ghost cells of the new MultiFabs
438  // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that
439  // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w
440 
441  FArrayBox temp_rhs_fab;
442  FArrayBox RHS_fab;
443  FArrayBox soln_fab;
444 
445  RHS_fab.resize (tbz,1,The_Async_Arena());
446  soln_fab.resize (tbz,1,The_Async_Arena());
447  temp_rhs_fab.resize(tbz,2,The_Async_Arena());
448 
449  auto const& RHS_a = RHS_fab.array();
450  auto const& soln_a = soln_fab.array();
451  auto const& temp_rhs_arr = temp_rhs_fab.array();
452 
453  auto const& coeffA_a = coeff_A_mf.array(mfi);
454  auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
455  auto const& coeffC_a = coeff_C_mf.array(mfi);
456  auto const& coeffP_a = coeff_P_mf.array(mfi);
457  auto const& coeffQ_a = coeff_Q_mf.array(mfi);
458 
459  // *************************************************************************
460  // Define flux arrays for use in advection
461  // *************************************************************************
462  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
463  flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
464  flux[dir].setVal<RunOn::Device>(0);
465  }
466  const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
467  flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
468 
469  // *********************************************************************
470  {
471  BL_PROFILE("fast_T_making_rho_rhs");
472  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
473  Real h_zeta_cc_xface_hi = myhalf * dzi *
474  ( z_nd(i+1,j ,k+1) + z_nd(i+1,j+1,k+1)
475  -z_nd(i+1,j ,k ) - z_nd(i+1,j+1,k ) );
476 
477  Real h_zeta_cc_xface_lo = myhalf * dzi *
478  ( z_nd(i ,j ,k+1) + z_nd(i ,j+1,k+1)
479  -z_nd(i ,j ,k ) - z_nd(i ,j+1,k ) );
480 
481  Real h_zeta_cc_yface_hi = myhalf * dzi *
482  ( z_nd(i ,j+1,k+1) + z_nd(i+1,j+1,k+1)
483  -z_nd(i ,j+1,k ) - z_nd(i+1,j+1,k ) );
484 
485  Real h_zeta_cc_yface_lo = myhalf * dzi *
486  ( z_nd(i ,j ,k+1) + z_nd(i+1,j ,k+1)
487  -z_nd(i ,j ,k ) - z_nd(i+1,j ,k ) );
488 
489  Real xflux_lo = new_drho_u(i ,j,k)*h_zeta_cc_xface_lo / mf_uy(i ,j,0);
490  Real xflux_hi = new_drho_u(i+1,j,k)*h_zeta_cc_xface_hi / mf_uy(i+1,j,0);
491  Real yflux_lo = new_drho_v(i,j ,k)*h_zeta_cc_yface_lo / mf_vx(i,j ,0);
492  Real yflux_hi = new_drho_v(i,j+1,k)*h_zeta_cc_yface_hi / mf_vx(i,j+1,0);
493 
494  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
495 
496  // NOTE: we are saving the (1/J) weighting for later when we add this to rho and theta
497  temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq +
498  ( yflux_hi - yflux_lo ) * dyi * mfsq;
499  temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
500  xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+
501  ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
502  yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) * myhalf;
503 
504  if (l_reflux) {
505  (flx_arr[0])(i,j,k,0) = xflux_lo;
506  (flx_arr[0])(i,j,k,1) = (flx_arr[0])(i ,j,k,0) * myhalf * (prim(i,j,k,0) + prim(i-1,j,k,0));
507 
508  (flx_arr[1])(i,j,k,0) = yflux_lo;
509  (flx_arr[1])(i,j,k,1) = (flx_arr[1])(i,j ,k,0) * myhalf * (prim(i,j,k,0) + prim(i,j-1,k,0));
510 
511  if (i == vbx_hi.x) {
512  (flx_arr[0])(i+1,j,k,0) = xflux_hi;
513  (flx_arr[0])(i+1,j,k,1) = (flx_arr[0])(i+1,j,k,0) * myhalf * (prim(i,j,k,0) + prim(i+1,j,k,0));
514  }
515  if (j == vbx_hi.y) {
516  (flx_arr[1])(i,j+1,k,0) = yflux_hi;
517  (flx_arr[1])(i,j+1,k,1) = (flx_arr[1])(i,j+1,k,0) * myhalf * (prim(i,j,k,0) + prim(i,j+1,k,0));
518  }
519  }
520  });
521  } // end profile
522 
523  // *********************************************************************
524  {
525  Box gbxo = mfi.nodaltilebox(2);
526  Box gbxo_mid = gbxo;
527 
528  if (gbxo.smallEnd(2) == domlo.z) {
529  Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2));
530  gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1);
531  ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
532  omega_arr(i,j,k) = zero;
533  });
534  }
535  if (gbxo.bigEnd(2) == domhi.z+1) {
536  Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
537  gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
538  ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
539  omega_arr(i,j,k) = old_drho_w(i,j,k);
540  });
541  }
542  ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
543  omega_arr(i,j,k) = OmegaFromW(i,j,k,old_drho_w(i,j,k),
544  old_drho_u,old_drho_v,
545  mf_ux,mf_vy,z_nd,dxInv);
546  });
547  } // end profile
548  // *********************************************************************
549 
550  Box bx_shrunk_in_k = bx;
551  int klo = tbz.smallEnd(2);
552  int khi = tbz.bigEnd(2);
553  bx_shrunk_in_k.setSmall(2,klo+1);
554  bx_shrunk_in_k.setBig(2,khi-1);
555 
556  // Note that the notes use "g" to mean the magnitude of gravity, so it is positive
557  // We set grav_gpu[2] to be the vector component which is negative
558  // We define halfg to match the notes (which is why we take the absolute value)
559  Real halfg = std::abs(myhalf * grav_gpu[2]);
560 
561  {
562  BL_PROFILE("fast_loop_on_shrunk_t");
563  //Note we don't act on the bottom or top boundaries of the domain
564  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
565  {
566  Real coeff_P = coeffP_a(i,j,k);
567  Real coeff_Q = coeffQ_a(i,j,k);
568 
569  Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
570  Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
571  Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
572 
573  // line 2 last two terms (order dtau)
574  Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k )
575  -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1);
576 
577  // line 3 residuals (order dtau^2) one <-> beta_2
578  Real R1_tmp = -halfg * ( slow_rhs_cons(i,j,k ,Rho_comp) + slow_rhs_cons(i,j,k-1,Rho_comp) );
579 
580  R1_tmp += coeff_P * slow_rhs_cons(i,j,k ,RhoTheta_comp)
581  + coeff_Q * slow_rhs_cons(i,j,k-1,RhoTheta_comp);
582 
583  Real Omega_kp1 = omega_arr(i,j,k+1);
584  Real Omega_k = omega_arr(i,j,k );
585  Real Omega_km1 = omega_arr(i,j,k-1);
586 
587  Real detJdiff = (detJ(i,j,k) - detJ(i,j,k-1)) / (detJ(i,j,k)*detJ(i,j,k-1));
588 
589  // consolidate lines 4&5 (order dtau^2)
590  R1_tmp += halfg * ( beta_1 * dzi * (Omega_kp1/detJ(i,j,k) + detJdiff*Omega_k - Omega_km1/detJ(i,j,k-1))
591  + temp_rhs_arr(i,j,k,Rho_comp)/detJ(i,j,k) + temp_rhs_arr(i,j,k-1,Rho_comp)/detJ(i,j,k-1) );
592 
593  // consolidate lines 6&7 (order dtau^2)
594  R1_tmp += -( coeff_P/detJ(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*theta_t_hi - Omega_k*theta_t_mid) + temp_rhs_arr(i,j,k ,RhoTheta_comp) )
595  + coeff_Q/detJ(i,j,k-1) * ( beta_1 * dzi * (Omega_k*theta_t_mid - Omega_km1*theta_t_lo) + temp_rhs_arr(i,j,k-1,RhoTheta_comp) ) );
596 
597  // line 1
598  RHS_a(i,j,k) = old_drho_w(i,j,k) + dtau * (slow_rhs_rho_w(i,j,k) + zmom_src_arr(i,j,k) + R0_tmp + dtau*beta_2*R1_tmp);
599 
600  // We cannot use omega_arr here since that was built with old_rho_u and old_rho_v ...
601  RHS_a(i,j,k) += OmegaFromW(i,j,k,zero,
602  new_drho_u,new_drho_v,
603  mf_ux,mf_vy,z_nd,dxInv);
604  });
605  } // end profile
606 
607  Box b2d = tbz; // Copy constructor
608  b2d.setRange(2,0);
609 
610  auto const lo = lbound(bx);
611  auto const hi = ubound(bx);
612 
613  {
614  BL_PROFILE("substep_b2d_loop_t");
615 #ifdef AMREX_USE_GPU
616  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
617  {
618  // w_klo, w_khi given by specified Dirichlet values
619  RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z ) + zmom_src_arr(i,j,lo.z ));
620  RHS_a(i,j,hi.z+1) = dtau * (slow_rhs_rho_w(i,j,hi.z+1) + zmom_src_arr(i,j,hi.z+1));
621 
622  // w = specified Dirichlet value at k = lo.z
623  soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
624 
625  // Transform the RHS from r_i -> rho_i
626  for (int k = lo.z+1; k <= hi.z+1; k++) {
627  soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
628  }
629 
630  cur_zmom(i,j,lo.z ) = stage_zmom(i,j,lo.z ) + soln_a(i,j,lo.z );
631  cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
632 
633  // Back sweep to obtain the solution
634  for (int k = hi.z; k >= lo.z; k--) {
635  soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1);
636  }
637  });
638 #else
639  // w_klo, w_khi given by specified Dirichlet values
640  for (int j = lo.y; j <= hi.y; ++j) {
641  AMREX_PRAGMA_SIMD
642  for (int i = lo.x; i <= hi.x; ++i)
643  {
644  RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z ) + zmom_src_arr(i,j,lo.z ));
645  RHS_a(i,j,hi.z+1) = dtau * (slow_rhs_rho_w(i,j,hi.z+1) + zmom_src_arr(i,j,hi.z+1));
646 
647  // w = specified Dirichlet value at k = lo.z
648  soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
649  }
650  }
651 
652  // Transform the RHS from r_i -> rho_i
653  // NOTE: this must include k = hi.z+1 so that the top row picks up the sub-diagonal
654  // term, which is non-zero when the top of the domain is an outflow boundary
655  for (int k = lo.z+1; k <= hi.z+1; ++k) {
656  for (int j = lo.y; j <= hi.y; ++j) {
657  AMREX_PRAGMA_SIMD
658  for (int i = lo.x; i <= hi.x; ++i) {
659  soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
660  }
661  }
662  }
663 
664  for (int j = lo.y; j <= hi.y; ++j) {
665  AMREX_PRAGMA_SIMD
666  for (int i = lo.x; i <= hi.x; ++i) {
667  cur_zmom(i,j,lo.z ) = stage_zmom(i,j,lo.z ) + soln_a(i,j,lo.z );
668  cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
669  }
670  }
671 
672  // Back sweep to obtain the solution
673  for (int k = hi.z; k >= lo.z; --k) {
674  for (int j = lo.y; j <= hi.y; ++j) {
675  AMREX_PRAGMA_SIMD
676  for (int i = lo.x; i <= hi.x; ++i) {
677  soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
678  }
679  }
680  }
681 #endif
682  } // end profile
683 
684  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
685  {
686  cur_zmom(i,j,k) = stage_zmom(i,j,k);
687  });
688 
689  if (lo.z == domlo.z) {
690  tbz.setSmall(2,domlo.z+1);
691  }
692  if (hi.z == domhi.z) {
693  tbz.setBig(2,domhi.z);
694  }
695  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
696  {
697  Real wpp = WFromOmega(i,j,k,soln_a(i,j,k),
698  new_drho_u,new_drho_v,
699  mf_ux,mf_vy,z_nd,dxInv);
700 
701  cur_zmom(i,j,k) += wpp;
702 
703  if (l_rayleigh_impl_for_w) {
704  Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
705  cur_zmom(i,j,k) /= (one + damping_coeff);
706  }
707  });
708 
709  // **************************************************************************
710  // Define updates in the RHS of rho and (rho theta)
711  // **************************************************************************
712  {
713  BL_PROFILE("fast_rho_final_update");
714  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
715  {
716  Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k);
717  Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1);
718 
719  // Note that in the solve we effectively impose new_drho_w(i,j,vbx_hi.z+1)=0
720  // so we don't update avg_zmom at k=vbx_hi.z+1
721  avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
722  if (l_reflux) {
723  (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
724  }
725 
726  if (k == vbx_hi.z) {
727  avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
728  if (l_reflux) {
729  (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
730  (flx_arr[2])(i,j,k+1,1) = (flx_arr[2])(i,j,k+1,0) * myhalf * (prim(i,j,k) + prim(i,j,k+1));
731  }
732  }
733 
734  Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k);
735 
736  cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho);
737 
738  Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) + myhalf *
739  ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) -
740  zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k);
741 
742  cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta);
743 
744  if (l_reflux) {
745  (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k) + prim(i,j,k-1));
746  }
747 
748  // add in source terms for cell-centered conserved variables
749  cur_cons(i,j,k,Rho_comp) += dtau * cc_src_arr(i,j,k,Rho_comp);
750  cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp);
751  });
752  } // end profile
753 
754  // We only add to the flux registers in the final RK step
755  if (l_reflux) {
756  int strt_comp_reflux = 0;
757  // For now we don't reflux (rho theta) because it seems to create issues at c/f boundaries
758  int num_comp_reflux = 1;
759  if (level < finest_level) {
760  fr_as_crse->CrseAdd(mfi,
761  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
762  dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
763  }
764  if (level > 0) {
765  fr_as_fine->FineAdd(mfi,
766  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
767  dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
768  }
769 
770  // This is necessary here so we don't go on to the next FArrayBox without
771  // having finished copying the fluxes into the FluxRegisters (since the fluxes
772  // are stored in temporary FArrayBox's)
773  Gpu::streamSynchronize();
774 
775  } // two-way coupling
776  } // mfi
777  } // OMP
778 }
constexpr amrex::Real R_v
Definition: ERF_Constants.H:35
constexpr amrex::Real R_d
Definition: ERF_Constants.H:34
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:54
@ v_x
Definition: ERF_DataStruct.H:29
@ u_y
Definition: ERF_DataStruct.H:30
@ v_y
Definition: ERF_DataStruct.H:30
@ m_y
Definition: ERF_DataStruct.H:30
@ u_x
Definition: ERF_DataStruct.H:29
@ m_x
Definition: ERF_DataStruct.H:29
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:61
#define Rho_comp
Definition: ERF_IndexDefines.H:39
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:58
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:35
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real OmegaFromW(int &i, int &j, int &k, amrex::Real w, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:791
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:292
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:269
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:339
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:385
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real WFromOmega(int &i, int &j, int &k, amrex::Real omega, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:856
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
@ gpy
Definition: ERF_IndexDefines.H:225
@ gpx
Definition: ERF_IndexDefines.H:224
@ ymom
Definition: ERF_IndexDefines.H:234
@ cons
Definition: ERF_IndexDefines.H:232
@ zmom
Definition: ERF_IndexDefines.H:235
@ xmom
Definition: ERF_IndexDefines.H:233
@ qt
Definition: ERF_Kessler.H:30
@ qv
Definition: ERF_Kessler.H:31
@ q
Definition: ERF_WSM6.H:273
Here is the call graph for this function: