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

Macros

#define CF_OMEGA_STAGE(I, J, K)   cf_omega_stage(I,J,K,dklo,dkhi,stage_zmom,stage_xmom,stage_ymom,mf_ux,mf_vy,z_nd,dxInv)
 

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 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, const Real beta_d, const CFFaceFluxes &cf_fluxes)
 

Macro Definition Documentation

◆ CF_OMEGA_STAGE

#define CF_OMEGA_STAGE (   I,
  J,
  K 
)    cf_omega_stage(I,J,K,dklo,dkhi,stage_zmom,stage_xmom,stage_ymom,mf_ux,mf_vy,z_nd,dxInv)

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  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,
const Real  beta_d,
const CFFaceFluxes &  cf_fluxes 
)

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, and beta_s on the w-faces (component 5)
[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]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_real_bc
[in]l_damp_coef
[in]beta_ddivergence damping coefficient: how far (rho theta) is projected forward
[in,out]cf_fluxescoarse/fine face fluxes: recorded here on a coarse level, imposed on a fine one
102 {
103  BL_PROFILE_REGION("erf_substep_T()");
104 
105  Real dtau = static_cast<Real>(dtau_d);
106 
107  const Box& domain = geom.Domain();
108  auto const domlo = lbound(domain);
109  auto const domhi = ubound(domain);
110 
111  int ilo = domlo.x;
112  int ihi = domhi.x + 1;
113  int jlo = domlo.y;
114  int jhi = domhi.y + 1;
115 
116  Real RvOverRd = R_v / R_d;
117 
118  bool l_rayleigh_impl_for_w = (sinesq_stag_d != nullptr);
119 
120  const Real* dx = geom.CellSize();
121  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
122 
123  // Set gravity as a vector
124  const Array<Real,AMREX_SPACEDIM> grav{zero, zero, -gravity};
125  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
126 
127  Real dxi = dxInv[0];
128  Real dyi = dxInv[1];
129  Real dzi = dxInv[2];
130  const auto& ba = S_stage_data[IntVars::cons].boxArray();
131  const auto& dm = S_stage_data[IntVars::cons].DistributionMap();
132 
133  MultiFab Delta_rho_u( convert(ba,IntVect(1,0,0)), dm, 1, 1);
134  MultiFab Delta_rho_v( convert(ba,IntVect(0,1,0)), dm, 1, 1);
135  MultiFab Delta_rho_w( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,0));
136  MultiFab Delta_rho ( ba , dm, 1, 1);
137  MultiFab Delta_rho_theta( ba , dm, 1, 1);
138 
139  MultiFab New_rho_u(convert(ba,IntVect(1,0,0)), dm, 1, 1);
140  MultiFab New_rho_v(convert(ba,IntVect(0,1,0)), dm, 1, 1);
141 
142  MultiFab coeff_A_mf(fast_coeffs, make_alias, 0, 1);
143  MultiFab inv_coeff_B_mf(fast_coeffs, make_alias, 1, 1);
144  MultiFab coeff_C_mf(fast_coeffs, make_alias, 2, 1);
145  MultiFab coeff_P_mf(fast_coeffs, make_alias, 3, 1);
146  MultiFab coeff_Q_mf(fast_coeffs, make_alias, 4, 1);
147 
148  // beta_s on the w-faces, from which
149  // beta_1 = myhalf * (one - beta_s) multiplies explicit terms
150  // beta_2 = myhalf * (one + beta_s) multiplies implicit terms
151  MultiFab beta_s_mf(fast_coeffs, make_alias, 5, 1);
152 
153  MultiFab Omega(S_data[IntVars::zmom].boxArray(), dm, 1, 1);
154 
155  MultiFab extrap(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(),1,1);
156 
157  // *************************************************************************
158  // Set up data that depends upon the "old" state
159  // *************************************************************************
160 #ifdef _OPENMP
161 #pragma omp parallel if (Gpu::notInLaunchRegion())
162 #endif
163  for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
164  {
165  const Array4<Real> & cur_cons = S_data[IntVars::cons].array(mfi);
166  const Array4<const Real>& prev_cons = S_prev[IntVars::cons].const_array(mfi);
167  const Array4<const Real>& stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
168  const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
169 
170  const Array4<Real>& old_drho = Delta_rho.array(mfi);
171  const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
172  const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
173  const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
174  const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
175 
176  const Array4<const Real>& prev_xmom = S_prev[IntVars::xmom].const_array(mfi);
177  const Array4<const Real>& prev_ymom = S_prev[IntVars::ymom].const_array(mfi);
178  const Array4<const Real>& prev_zmom = S_prev[IntVars::zmom].const_array(mfi);
179 
180  const Array4<const Real>& stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
181  const Array4<const Real>& stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
182  const Array4<const Real>& stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi);
183 
184  const Array4<Real>& theta_extrap = extrap.array(mfi);
185  const Array4<const Real>& prim = S_stage_prim.const_array(mfi);
186 
187  Box bx = mfi.validbox();
188 
189  // NOTE: growntilebox / grownnodaltilebox, not the tilebox grown by hand. Growing by
190  // hand pushes every tile one cell past its own share of the grid, so once the
191  // grid is tiled two tiles store to the same cells -- a data race under OpenMP,
192  // even though both store the same value. The union of the tiles, and the values
193  // stored, are unchanged, and this is identical to the old expression when there
194  // is one tile per grid.
195  Box gbx = mfi.growntilebox(1);
196 
197  if (step == 0) {
198  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
199  cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp);
200  cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp);
201  });
202  } // step = 0
203 
204  Box gtbx = mfi.grownnodaltilebox(0,IntVect(1,1,0));
205  Box gtby = mfi.grownnodaltilebox(1,IntVect(1,1,0));
206  Box gtbz = mfi.grownnodaltilebox(2,IntVect(1,1,0));
207 
208  const auto& bx_lo = lbound(bx);
209  const auto& bx_hi = ubound(bx);
210 
211  // *********************************************************************
212  // Populate "old" perturbations in momenta
213  // *********************************************************************
214  ParallelFor(gtbx, gtby, gtbz,
215  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
216  old_drho_u(i,j,k) = prev_xmom(i,j,k) - stage_xmom(i,j,k);
217  if (k == bx_lo.z && k != domlo.z) {
218  old_drho_u(i,j,k-1) = old_drho_u(i,j,k);
219  }
220  if (k == bx_hi.z) {
221  old_drho_u(i,j,k+1) = old_drho_u(i,j,k);
222  }
223  },
224  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
225  old_drho_v(i,j,k) = prev_ymom(i,j,k) - stage_ymom(i,j,k);
226  if (k == bx_lo.z && k != domlo.z) {
227  old_drho_v(i,j,k-1) = old_drho_v(i,j,k);
228  }
229  if (k == bx_hi.z) {
230  old_drho_v(i,j,k+1) = old_drho_v(i,j,k);
231  }
232  },
233  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
234  old_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k);
235  });
236 
237  // *********************************************************************
238  // Populate "old" perturbations in Rho & RhoTheta and lagged RhoTheta
239  // *********************************************************************
240  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
241  old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp);
242  old_drho_theta(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) - stage_cons(i,j,k,RhoTheta_comp);
243 
244  // Compute theta extrap for lateral momenta pressure gradient
245  if (step == 0) {
246  theta_extrap(i,j,k) = old_drho_theta(i,j,k);
247  } else {
248  theta_extrap(i,j,k) = old_drho_theta(i,j,k) + beta_d *
249  ( old_drho_theta(i,j,k) - lagged_arr(i,j,k) );
250  }
251 
252  // Convert to Theta_m
253  // NOTE: qv is not changing over the fast steps so we use the stage data
254  Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero;
255  theta_extrap(i,j,k) *= (one + RvOverRd*qv);
256 
257  // Save Theta^{'',n-1} for theta extrap next substep
258  lagged_arr(i,j,k) = old_drho_theta(i,j,k);
259  });
260  } // mfi
261 
262  // *************************************************************************
263  // Define updates in the current RK stage
264  // *************************************************************************
265 #ifdef _OPENMP
266 #pragma omp parallel if (Gpu::notInLaunchRegion())
267 #endif
268  for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
269  {
270  Box bx = mfi.validbox();
271  Box tbx = mfi.nodaltilebox(0);
272  Box tby = mfi.nodaltilebox(1);
273 
274  const Array4<Real const>& xmom_src_arr = xmom_src.const_array(mfi);
275  const Array4<Real const>& ymom_src_arr = ymom_src.const_array(mfi);
276 
277  const Array4<const Real> & stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
278  const Array4<const Real> & stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
279  const Array4<const Real> & qt_arr = qt.const_array(mfi);
280 
281  const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
282  const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
283 
284  const Array4<const Real>& slow_rhs_rho_u = S_slow_rhs[IntVars::xmom].const_array(mfi);
285  const Array4<const Real>& slow_rhs_rho_v = S_slow_rhs[IntVars::ymom].const_array(mfi);
286 
287  const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
288  const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
289 
290  const Array4<Real>& cur_xmom = S_data[IntVars::xmom].array(mfi);
291  const Array4<Real>& cur_ymom = S_data[IntVars::ymom].array(mfi);
292 
293  // These store the advection momenta which we will use to update the slow variables
294  const Array4<Real>& avg_xmom_arr = avg_xmom.array(mfi);
295  const Array4<Real>& avg_ymom_arr = avg_ymom.array(mfi);
296 
297  const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
298 
299  const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
300 
301  const Array4<Real>& theta_extrap = extrap.array(mfi);
302 
303  // Map factors
304  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
305  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
306  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
307  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
308 
309  // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs
310  // Note that we do the Copy and Subtract including one ghost cell
311  // so that we don't have to fill ghost cells of the new MultiFabs
312  // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that
313  // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w
314 
315  // *********************************************************************
316  // Explicit update for lateral momenta
317  // *********************************************************************
318  {
319  BL_PROFILE("substep_xymom_T");
320 
321  const auto& bx_lo = lbound(bx);
322  const auto& bx_hi = ubound(bx);
323 
324  ParallelFor(tbx, tby,
325  [=] AMREX_GPU_DEVICE (int i, int j, int k)
326  {
327  // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
328  Real met_h_xi = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd);
329  Real met_h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd);
330  Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi;
331  Real gp_zeta_on_iface = (k == 0) ?
332  myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
333  - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) :
334  fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
335  - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) );
336  Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? Real(0.) :
337  gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface;
338 
339  gpx *= mf_ux(i,j,0);
340 
341  Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) : zero;
342 
343  Real pi_c = myhalf * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i ,j,k,0));
344  Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx / (one + q);
345 
346  new_drho_u(i, j, k) = old_drho_u(i,j,k) + dtau * fast_rhs_rho_u
347  + dtau * slow_rhs_rho_u(i,j,k)
348  + dtau * xmom_src_arr(i,j,k);
349  if (k == bx_lo.z && k != domlo.z) {
350  new_drho_u(i,j,k-1) = new_drho_u(i,j,k);
351  }
352  if (k == bx_hi.z) {
353  new_drho_u(i,j,k+1) = new_drho_u(i,j,k);
354  }
355 
356  // NOTE: met_h_zeta here is identically the x-face area ax computed by
357  // make_areas, so this matches the base value of avg_xmom defined in
358  // AdvectionSrcForRho (ax*rho_u/mf_uy) as well as the density flux
359  // (new_drho_u*h_zeta_cc_xface/mf_uy) formed below.
360  avg_xmom_arr(i,j,k) += facinv * new_drho_u(i,j,k) * met_h_zeta / mf_uy(i,j,0);
361 
362  cur_xmom(i,j,k) = stage_xmom(i,j,k) + new_drho_u(i,j,k);
363  },
364  [=] AMREX_GPU_DEVICE (int i, int j, int k)
365  {
366  // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
367  Real met_h_eta = Compute_h_eta_AtJface(i, j, k, dxInv, z_nd);
368  Real met_h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd);
369  Real gp_eta = (theta_extrap(i,j,k) - theta_extrap(i,j-1,k)) * dyi;
370  Real gp_zeta_on_jface = (k == 0) ?
371  myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
372  - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) :
373  fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
374  - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) );
375  Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? Real(0.) :
376  gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface;
377 
378  gpy *= mf_vy(i,j,0);
379 
380  Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) : zero;
381 
382  Real pi_c = myhalf * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j ,k,0));
383  Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy / (one + q);
384 
385  new_drho_v(i, j, k) = old_drho_v(i,j,k) + dtau * fast_rhs_rho_v
386  + dtau * slow_rhs_rho_v(i,j,k)
387  + dtau * ymom_src_arr(i,j,k);
388 
389  if (k == bx_lo.z && k != domlo.z) {
390  new_drho_v(i,j,k-1) = new_drho_v(i,j,k);
391  }
392  if (k == bx_hi.z) {
393  new_drho_v(i,j,k+1) = new_drho_v(i,j,k);
394  }
395 
396  // NOTE: met_h_zeta here is identically the y-face area ay computed by
397  // make_areas, so this matches the base value of avg_ymom defined in
398  // AdvectionSrcForRho (ay*rho_v/mf_vx) as well as the density flux
399  // (new_drho_v*h_zeta_cc_yface/mf_vx) formed below.
400  avg_ymom_arr(i,j,k) += facinv * new_drho_v(i,j,k) * met_h_zeta / mf_vx(i,j,0);
401 
402  cur_ymom(i,j,k) = stage_ymom(i,j,k) + new_drho_v(i,j,k);
403  });
404 
405  // On the coarse/fine faces of a fine level, the coarse level's flux over this
406  // substep (CFFaceFluxes) in place of the level's own; on a coarse level, record
407  // the flux each face applied. In both, the flux is h_zeta (rho u)/m_uy, the
408  // density flux formed below with h_zeta the x-face area.
409  for (int dir = 0; dir < 2; ++dir) {
410  if (!cf_fluxes.tgt[dir] && !cf_fluxes.rec[dir]) { continue; }
411  const bool has_tgt = (cf_fluxes.tgt[dir] != nullptr);
412  const bool has_rec = (cf_fluxes.rec[dir] != nullptr);
413  const Box& fbx = (dir == 0) ? tbx : tby;
414  const Array4<const Real> tgt = (has_tgt) ? cf_fluxes.tgt[dir]->const_array(mfi) : Array4<const Real>{};
415  const Array4< Real> rec = (has_rec) ? cf_fluxes.rec[dir]->array(mfi) : Array4<Real>{};
416  const Array4<Real>& pert = (dir == 0) ? new_drho_u : new_drho_v;
417  const Array4<Real>& cur = (dir == 0) ? cur_xmom : cur_ymom;
418  const Array4<Real>& avg = (dir == 0) ? avg_xmom_arr : avg_ymom_arr;
419  const Array4<const Real>& stage = (dir == 0) ? stage_xmom : stage_ymom;
420  const Array4<const Real>& mf = (dir == 0) ? mf_uy : mf_vx;
421  ParallelFor(fbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
422  {
423  Real h_zeta = (dir == 0) ? Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd)
424  : Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd);
425  if (has_tgt && tgt(i,j,k,1) > myhalf) {
426  Real new_pert = tgt(i,j,k,0) * mf(i,j,0) / h_zeta - stage(i,j,k);
427  avg(i,j,k) += facinv * (new_pert - pert(i,j,k)) * h_zeta / mf(i,j,0);
428  pert(i,j,k) = new_pert;
429  if (k == bx_lo.z && k != domlo.z) {
430  pert(i,j,k-1) = new_pert;
431  } else if (k == bx_hi.z) {
432  pert(i,j,k+1) = new_pert;
433  }
434  cur(i,j,k) = stage(i,j,k) + new_pert;
435  }
436  if (has_rec) {
437  rec(i,j,k) = (stage(i,j,k) + pert(i,j,k)) * h_zeta / mf(i,j,0);
438  }
439  });
440  }
441  } // end profile
442  }
443 
444 #ifdef _OPENMP
445 #pragma omp parallel if (Gpu::notInLaunchRegion())
446 #endif
447  {
448  std::array<FArrayBox,AMREX_SPACEDIM> flux;
449  for ( MFIter mfi(S_stage_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
450  {
451  Box bx = mfi.tilebox();
452  Box tbz = surroundingNodes(bx,2);
453 
454  Box vbx = mfi.validbox();
455  const auto& vbx_hi = ubound(vbx);
456 
457  const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
458  const Array4<Real const>& cc_src_arr = cc_src.const_array(mfi);
459 
460  const Array4<const Real> & stage_xmom = S_stage_data[IntVars::xmom].const_array(mfi);
461  const Array4<const Real> & stage_ymom = S_stage_data[IntVars::ymom].const_array(mfi);
462  const Array4<const Real> & stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi);
463  const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
464 
465  const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
466  const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
467  const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
468  const Array4<Real>& old_drho = Delta_rho.array(mfi);
469  const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
470 
471  const Array4<const Real>& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi);
472  const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi);
473 
474  const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
475  const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
476 
477  const Array4<Real>& cur_cons = S_data[IntVars::cons].array(mfi);
478  const Array4<Real>& cur_zmom = S_data[IntVars::zmom].array(mfi);
479 
480  // These store the advection momenta which we will use to update the slow variables
481  const Array4<Real>& avg_zmom_arr = avg_zmom.array(mfi);
482 
483  const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
484  const Array4<const Real>& detJ = detJ_cc->const_array(mfi);
485 
486  const Array4< Real>& omega_arr = Omega.array(mfi);
487 
488  // Map factors
489  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
490  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
491  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
492  const Array4<const Real>& mf_uy = mapfac[MapFacType::u_y]->const_array(mfi);
493  const Array4<const Real>& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi);
494  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
495 
496  // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs
497  // Note that we do the Copy and Subtract including one ghost cell
498  // so that we don't have to fill ghost cells of the new MultiFabs
499  // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that
500  // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w
501 
502  FArrayBox temp_rhs_fab;
503  FArrayBox RHS_fab;
504  FArrayBox soln_fab;
505 
506  RHS_fab.resize (tbz,1,The_Async_Arena());
507  soln_fab.resize (tbz,1,The_Async_Arena());
508  temp_rhs_fab.resize(tbz,2,The_Async_Arena());
509 
510  auto const& RHS_a = RHS_fab.array();
511  auto const& soln_a = soln_fab.array();
512  auto const& temp_rhs_arr = temp_rhs_fab.array();
513 
514  auto const& coeffA_a = coeff_A_mf.array(mfi);
515  auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
516  auto const& coeffC_a = coeff_C_mf.array(mfi);
517  auto const& coeffP_a = coeff_P_mf.array(mfi);
518  auto const& coeffQ_a = coeff_Q_mf.array(mfi);
519  auto const& beta_s_a = beta_s_mf.const_array(mfi);
520 
521  // *************************************************************************
522  // Define flux arrays for use in advection
523  // *************************************************************************
524  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
525  flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
526  flux[dir].setVal<RunOn::Device>(0);
527  }
528  const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
529  flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
530 
531  // *********************************************************************
532  // Lateral Advection RHS for Rho & RhoTheta
533  // *********************************************************************
534  {
535  BL_PROFILE("fast_T_making_rho_rhs");
536  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
537  Real h_zeta_cc_xface_hi = myhalf * dzi *
538  ( z_nd(i+1,j ,k+1) + z_nd(i+1,j+1,k+1)
539  -z_nd(i+1,j ,k ) - z_nd(i+1,j+1,k ) );
540 
541  Real h_zeta_cc_xface_lo = myhalf * dzi *
542  ( z_nd(i ,j ,k+1) + z_nd(i ,j+1,k+1)
543  -z_nd(i ,j ,k ) - z_nd(i ,j+1,k ) );
544 
545  Real h_zeta_cc_yface_hi = myhalf * dzi *
546  ( z_nd(i ,j+1,k+1) + z_nd(i+1,j+1,k+1)
547  -z_nd(i ,j+1,k ) - z_nd(i+1,j+1,k ) );
548 
549  Real h_zeta_cc_yface_lo = myhalf * dzi *
550  ( z_nd(i ,j ,k+1) + z_nd(i+1,j ,k+1)
551  -z_nd(i ,j ,k ) - z_nd(i+1,j ,k ) );
552 
553  Real xflux_lo = new_drho_u(i ,j,k)*h_zeta_cc_xface_lo / mf_uy(i ,j,0);
554  Real xflux_hi = new_drho_u(i+1,j,k)*h_zeta_cc_xface_hi / mf_uy(i+1,j,0);
555  Real yflux_lo = new_drho_v(i,j ,k)*h_zeta_cc_yface_lo / mf_vx(i,j ,0);
556  Real yflux_hi = new_drho_v(i,j+1,k)*h_zeta_cc_yface_hi / mf_vx(i,j+1,0);
557 
558  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
559 
560  // NOTE: we are saving the (1/J) weighting for later when we add this to rho and theta
561  temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq +
562  ( yflux_hi - yflux_lo ) * dyi * mfsq;
563  temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
564  xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+
565  ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
566  yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) * myhalf;
567 
568  if (l_reflux) {
569  (flx_arr[0])(i,j,k,0) = xflux_lo;
570  (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));
571 
572  (flx_arr[1])(i,j,k,0) = yflux_lo;
573  (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));
574 
575  if (i == vbx_hi.x) {
576  (flx_arr[0])(i+1,j,k,0) = xflux_hi;
577  (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));
578  }
579  if (j == vbx_hi.y) {
580  (flx_arr[1])(i,j+1,k,0) = yflux_hi;
581  (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));
582  }
583  }
584  });
585  }
586 
587  // *********************************************************************
588  // Populate "old" Omega
589  // *********************************************************************
590  // NOTE: this must stay in a separate MFIter loop from the one that fills
591  // old_drho_u/v, since OmegaFromW reads old_drho_u/v at (i+1, j+1, k-1),
592  // which another tile fills when the grid is tiled under OpenMP
593  {
594  Box gbxo = mfi.nodaltilebox(2);
595  Box gbxo_mid = gbxo;
596 
597  if (gbxo.smallEnd(2) == domlo.z) {
598  Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2));
599  gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1);
600  ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
601  omega_arr(i,j,k) = zero;
602  });
603  }
604  if (gbxo.bigEnd(2) == domhi.z+1) {
605  Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
606  gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
607  ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
608  omega_arr(i,j,k) = old_drho_w(i,j,k);
609  });
610  }
611  ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
612  omega_arr(i,j,k) = OmegaFromW(i,j,k,old_drho_w(i,j,k),
613  old_drho_u,old_drho_v,
614  mf_ux,mf_vy,z_nd,dxInv);
615  });
616  }
617 
618  // *********************************************************************
619  // Build RHS for tridiagonal system
620  // *********************************************************************
621  // NOTE: we don't act on the bottom or top boundaries of the domain
622  Box bx_shrunk_in_k = bx;
623  int klo = tbz.smallEnd(2);
624  int khi = tbz.bigEnd(2);
625  bx_shrunk_in_k.setSmall(2,klo+1);
626  bx_shrunk_in_k.setBig(2,khi-1);
627 
628  // Note that the notes use "g" to mean the magnitude of gravity, so it is positive
629  // We set grav_gpu[2] to be the vector component which is negative
630  // We define halfg to match the notes (which is why we take the absolute value)
631  Real halfg = std::abs(myhalf * grav_gpu[2]);
632 
633  {
634  BL_PROFILE("fast_loop_on_shrunk_t");
635  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
636  {
637  Real coeff_P = coeffP_a(i,j,k);
638  Real coeff_Q = coeffQ_a(i,j,k);
639 
640  Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
641  Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
642  Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
643 
644  // line 2 last two terms (order dtau)
645  Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k )
646  -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1);
647 
648  // line 3 residuals (order dtau^2) one <-> beta_2
649  Real R1_tmp = -halfg * ( slow_rhs_cons(i,j,k ,Rho_comp) + slow_rhs_cons(i,j,k-1,Rho_comp) );
650 
651  R1_tmp += coeff_P * slow_rhs_cons(i,j,k ,RhoTheta_comp)
652  + coeff_Q * slow_rhs_cons(i,j,k-1,RhoTheta_comp);
653 
654  Real Omega_kp1 = omega_arr(i,j,k+1);
655  Real Omega_k = omega_arr(i,j,k );
656  Real Omega_km1 = omega_arr(i,j,k-1);
657 
658  // Explicit part of the vertical flux on each face is beta_1 * Omega. beta_1 on face k
659  // is factored out, and dbeta_1 is how much larger beta_1 is on the neighbouring faces,
660  // which is exactly zero (so the roundoff is that of a scalar beta_1) for uniform beta_s
661  Real beta_1 = myhalf * (one - beta_s_a(i,j,k));
662  Real dbeta_1_kp1 = myhalf * (one - beta_s_a(i,j,k+1)) - beta_1;
663  Real dbeta_1_km1 = myhalf * (one - beta_s_a(i,j,k-1)) - beta_1;
664 
665  Real detJdiff = (detJ(i,j,k) - detJ(i,j,k-1)) / (detJ(i,j,k)*detJ(i,j,k-1));
666 
667  // consolidate lines 4&5 (order dtau^2)
668  R1_tmp += halfg * ( beta_1 * dzi * (Omega_kp1/detJ(i,j,k) + detJdiff*Omega_k - Omega_km1/detJ(i,j,k-1))
669  + dzi * (dbeta_1_kp1*Omega_kp1/detJ(i,j,k) - dbeta_1_km1*Omega_km1/detJ(i,j,k-1))
670  + 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) );
671 
672  // consolidate lines 6&7 (order dtau^2)
673  R1_tmp += -( coeff_P/detJ(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*theta_t_hi - Omega_k*theta_t_mid) + dzi * dbeta_1_kp1*Omega_kp1*theta_t_hi + temp_rhs_arr(i,j,k ,RhoTheta_comp) )
674  + coeff_Q/detJ(i,j,k-1) * ( beta_1 * dzi * (Omega_k*theta_t_mid - Omega_km1*theta_t_lo) - dzi * dbeta_1_km1*Omega_km1*theta_t_lo + temp_rhs_arr(i,j,k-1,RhoTheta_comp) ) );
675 
676  // consolidate terms proportional to dtau
677  Real beta_2 = myhalf * (one + beta_s_a(i,j,k));
678  RHS_a(i,j,k) = dtau * (slow_rhs_rho_w(i,j,k) + zmom_src_arr(i,j,k) + R0_tmp + dtau*beta_2*R1_tmp);
679 
680  // We cannot use omega_arr here since that was built with old_rho_u and old_rho_v ...
681  RHS_a(i,j,k) += old_drho_w(i,j,k) + OmegaFromW(i,j,k,zero,
682  new_drho_u,new_drho_v,
683  mf_ux,mf_vy,z_nd,dxInv);
684  });
685  } // end profile
686 
687  // *********************************************************************
688  // Solve tridiagonal system
689  // *********************************************************************
690  Box b2d = tbz;
691  b2d.setRange(2,0);
692  auto const lo = lbound(bx);
693  auto const hi = ubound(bx);
694 
695  // The coarse/fine face data (CFFaceFluxes) are in Omega, which is the full Omega of the
696  // stage -- built as in erf_slow_rhs_pre, zero at the bottom of the domain and (rho w) at
697  // its top -- plus the Omega'' of the solve
698  const bool has_tgt_z = (cf_fluxes.tgt[2] != nullptr);
699  const bool has_rec_z = (cf_fluxes.rec[2] != nullptr);
700  const Array4<const Real> tgt_z = (has_tgt_z) ? cf_fluxes.tgt[2]->const_array(mfi) : Array4<const Real>{};
701  const Array4<const Real> tgt_om = (has_tgt_z) ? cf_fluxes.tgt_w->const_array(mfi) : Array4<const Real>{};
702  const Array4< Real> rec_z = (has_rec_z) ? cf_fluxes.rec[2]->array(mfi) : Array4<Real>{};
703  const Array4< Real> rec_om = (has_rec_z) ? cf_fluxes.rec_w->array(mfi) : Array4<Real>{};
704  const bool has_rec_z_old = (cf_fluxes.rec_w_old != nullptr);
705  const Array4< Real> rec_om_old = (has_rec_z_old) ? cf_fluxes.rec_w_old->array(mfi) : Array4<Real>{};
706  const int dklo = domlo.z;
707  const int dkhi = domhi.z;
708 #define CF_OMEGA_STAGE(I,J,K) cf_omega_stage(I,J,K,dklo,dkhi,stage_zmom,stage_xmom,stage_ymom,mf_ux,mf_vy,z_nd,dxInv)
709 
710  {
711  BL_PROFILE("substep_b2d_loop_t");
712 
713  // NOTE: RHS is built on a shrunk box in z, explicitly set
714  // the lo and hi faces here
715  // On a box face inside the domain (a c/f face) the row carries the accumulated
716  // perturbation converted to Omega'', as in the interior rows. At the domain bottom
717  // and top the row is zero: Omega'' = 0 at the ground and at a wall, and
718  // Omega''_top = Omega''_top-1 for outflow, the matrix row set in make_fast_coeffs
719  // deciding which
720  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
721  {
722  // OmegaFromW takes its indices by non-const reference, hence the copies
723  int kf_lo = lo.z;
724  int kf_hi = hi.z+1;
725  RHS_a(i,j,lo.z ) = (lo.z == domlo.z)
726  ? zero
727  : old_drho_w(i,j,lo.z ) + OmegaFromW(i,j,kf_lo,zero,
728  new_drho_u,new_drho_v,
729  mf_ux,mf_vy,z_nd,dxInv)
730  + dtau * slow_rhs_rho_w(i,j,lo.z )
731  + dtau * zmom_src_arr(i,j,lo.z );
732  RHS_a(i,j,hi.z+1) = (hi.z == domhi.z)
733  ? zero
734  : old_drho_w(i,j,hi.z+1) + OmegaFromW(i,j,kf_hi,zero,
735  new_drho_u,new_drho_v,
736  mf_ux,mf_vy,z_nd,dxInv)
737  + dtau * slow_rhs_rho_w(i,j,hi.z+1)
738  + dtau * zmom_src_arr(i,j,hi.z+1);
739 
740  // On the coarse/fine faces of a fine level, Omega'' from the coarse level's record
741  if (has_tgt_z) {
742  if (tgt_z(i,j,lo.z ,1) > myhalf) { RHS_a(i,j,lo.z ) = tgt_om(i,j,lo.z ) - CF_OMEGA_STAGE(i,j,lo.z ); }
743  if (tgt_z(i,j,hi.z+1,1) > myhalf) { RHS_a(i,j,hi.z+1) = tgt_om(i,j,hi.z+1) - CF_OMEGA_STAGE(i,j,hi.z+1); }
744  }
745  });
746 
747 #ifdef AMREX_USE_GPU
748  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
749  {
750  // First row (coeffA = coeffC = 0) at k = lo.z
751  soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
752 
753  // Transform the RHS from r_i -> rho_i
754  for (int k = lo.z+1; k <= hi.z+1; k++) {
755  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);
756  }
757 
758  // Back sweep to obtain the solution
759  for (int k = hi.z; k >= lo.z; k--) {
760  soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1);
761  }
762  });
763 #else
764  // First row (coeffA = coeffC = 0) at k = lo.z
765  for (int j = lo.y; j <= hi.y; ++j) {
766  AMREX_PRAGMA_SIMD
767  for (int i = lo.x; i <= hi.x; ++i) {
768  soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
769  }
770  }
771 
772  // Transform the RHS from r_i -> rho_i
773  // NOTE: this must include k = hi.z+1 so that the top row picks up the sub-diagonal
774  // term, which is non-zero when the top of the domain is an outflow boundary
775  for (int k = lo.z+1; k <= hi.z+1; ++k) {
776  for (int j = lo.y; j <= hi.y; ++j) {
777  AMREX_PRAGMA_SIMD
778  for (int i = lo.x; i <= hi.x; ++i) {
779  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);
780  }
781  }
782  }
783 
784  // Back sweep to obtain the solution
785  for (int k = hi.z; k >= lo.z; --k) {
786  for (int j = lo.y; j <= hi.y; ++j) {
787  AMREX_PRAGMA_SIMD
788  for (int i = lo.x; i <= hi.x; ++i) {
789  soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
790  }
791  }
792  }
793 #endif
794  } // end profile
795 
796  // *********************************************************************
797  // Convert from Omega to W
798  // *********************************************************************
799  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
800  {
801  // NOTE: WfromOmega uses an interior stencil at domlo.z to touch valid data
802  // WfromOmega cannot be used at domhi.z+1 for outflow. We enforce
803  // zero metric terms at the domain high (W'' = Omega'')
804  Real wpp;
805  if (k == domhi.z+1) {
806  wpp = soln_a(i,j,k);
807  } else {
808  wpp = WFromOmega(i,j,k,soln_a(i,j,k),
809  new_drho_u,new_drho_v,
810  mf_ux,mf_vy,z_nd,dxInv);
811  }
812 
813  cur_zmom(i,j,k) = stage_zmom(i,j,k) + wpp;
814 
815  // W at the ground is kinematic and is not damped; the top face is, since with
816  // outflow it is no longer pinned to the stage value. A coarse/fine face whose
817  // value is imposed from the coarse level's record is not damped either: the coarse
818  // level has damped it already, and it must match the record
819  const bool imposed = has_tgt_z && (tgt_z(i,j,k,1) > myhalf);
820  if (k > domlo.z && !imposed && l_rayleigh_impl_for_w) {
821  Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
822  cur_zmom(i,j,k) /= (one + damping_coeff);
823  }
824  });
825 
826  // **************************************************************************
827  // Update Rho & RhoTheta with lateral and vertical advection
828  // **************************************************************************
829  {
830  BL_PROFILE("fast_rho_final_update");
831  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
832  {
833  Real zflux_lo = myhalf * (one + beta_s_a(i,j,k )) * soln_a(i,j,k )
834  + myhalf * (one - beta_s_a(i,j,k )) * omega_arr(i,j,k );
835  Real zflux_hi = myhalf * (one + beta_s_a(i,j,k+1)) * soln_a(i,j,k+1)
836  + myhalf * (one - beta_s_a(i,j,k+1)) * omega_arr(i,j,k+1);
837 
838  // On the coarse/fine faces of a fine level, the coarse level's flux over this substep
839  if (has_tgt_z) {
840  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
841  if (tgt_z(i,j,k ,1) > myhalf) { zflux_lo = tgt_z(i,j,k ,0) * mfsq - CF_OMEGA_STAGE(i,j,k ); }
842  if (tgt_z(i,j,k+1,1) > myhalf) { zflux_hi = tgt_z(i,j,k+1,0) * mfsq - CF_OMEGA_STAGE(i,j,k+1); }
843  }
844 
845  // On a coarse level, the flux each face applied over this substep, and Omega at
846  // its end (and at its start when asked)
847  if (has_rec_z) {
848  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
849  Real om_k = CF_OMEGA_STAGE(i,j,k);
850  rec_z (i,j,k) = (zflux_lo + om_k) / mfsq;
851  rec_om(i,j,k) = om_k + soln_a(i,j,k);
852  if (has_rec_z_old) { rec_om_old(i,j,k) = om_k + omega_arr(i,j,k); }
853  if (k == vbx_hi.z) {
854  Real om_kp1 = CF_OMEGA_STAGE(i,j,k+1);
855  rec_z (i,j,k+1) = (zflux_hi + om_kp1) / mfsq;
856  rec_om(i,j,k+1) = om_kp1 + soln_a(i,j,k+1);
857  if (has_rec_z_old) { rec_om_old(i,j,k+1) = om_kp1 + omega_arr(i,j,k+1); }
858  }
859  }
860 
861  // avg_zmom at k=vbx_hi.z+1 is updated by the top cell of the box; the
862  // perturbation there is zero at a wall but not at an outflow top or a c/f face
863  avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
864  if (l_reflux) {
865  (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
866  }
867 
868  if (k == vbx_hi.z) {
869  avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
870  if (l_reflux) {
871  (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
872  (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));
873  }
874  }
875 
876  Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k);
877 
878  cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho);
879 
880  Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) + myhalf *
881  ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) -
882  zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k);
883 
884  cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta);
885 
886  if (l_reflux) {
887  (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k) + prim(i,j,k-1));
888  }
889 
890  // add in source terms for cell-centered conserved variables
891  cur_cons(i,j,k,Rho_comp) += dtau * cc_src_arr(i,j,k,Rho_comp);
892  cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp);
893  });
894  } // end profile
895 
896  // We only add to the flux registers in the final RK step
897  if (l_reflux) {
898  int strt_comp_reflux = 0;
899  // For now we don't reflux (rho theta) because it seems to create issues at c/f boundaries
900  int num_comp_reflux = 1;
901  if (level < finest_level) {
902  fr_as_crse->CrseAdd(mfi,
903  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
904  dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
905  }
906  if (level > 0) {
907  fr_as_fine->FineAdd(mfi,
908  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
909  dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
910  }
911 
912  // This is necessary here so we don't go on to the next FArrayBox without
913  // having finished copying the fluxes into the FluxRegisters (since the fluxes
914  // are stored in temporary FArrayBox's)
915  Gpu::streamSynchronize();
916 
917  } // two-way coupling
918 #undef CF_OMEGA_STAGE
919  } // mfi
920  } // OMP
921 }
OpenSide::hi
@ hi
OpenSide::lo
@ lo
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:50
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:49
Gamma
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:70
v_x
@ v_x
Definition: ERF_DataStruct.H:38
u_y
@ u_y
Definition: ERF_DataStruct.H:39
v_y
@ v_y
Definition: ERF_DataStruct.H:39
m_y
@ m_y
Definition: ERF_DataStruct.H:39
u_x
@ u_x
Definition: ERF_DataStruct.H:38
m_x
@ m_x
Definition: ERF_DataStruct.H:38
PrimQ1_comp
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:82
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:53
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:54
PrimTheta_comp
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:79
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:25
klo
const int klo
Definition: ERF_InitCustomPert_KE.H:9
ParallelFor
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);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
CF_OMEGA_STAGE
#define CF_OMEGA_STAGE(I, J, K)
OmegaFromW
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:802
Compute_h_eta_AtJface
AMREX_GPU_HOST_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:396
Compute_h_zeta_AtIface
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:280
Compute_h_zeta_AtJface
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:350
Compute_h_xi_AtIface
AMREX_GPU_HOST_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:303
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
GpVars::gpy
@ gpy
Definition: ERF_IndexDefines.H:265
GpVars::gpx
@ gpx
Definition: ERF_IndexDefines.H:264
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:274
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:272
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:275
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:273
MicVar_Kess::qt
@ qt
Definition: ERF_Kessler.H:30
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:31
WSM6Ind::q
@ q
Definition: ERF_WSM6.H:294
CFFaceFluxes::rec
std::array< amrex::MultiFab *, AMREX_SPACEDIM > rec
Definition: ERF_TI_fast_headers.H:29
CFFaceFluxes::rec_w
amrex::MultiFab * rec_w
Definition: ERF_TI_fast_headers.H:30
CFFaceFluxes::tgt_w
const amrex::MultiFab * tgt_w
Definition: ERF_TI_fast_headers.H:33
CFFaceFluxes::tgt
std::array< const amrex::MultiFab *, AMREX_SPACEDIM > tgt
Definition: ERF_TI_fast_headers.H:32
CFFaceFluxes::rec_w_old
amrex::MultiFab * rec_w_old
Definition: ERF_TI_fast_headers.H:31
Here is the call graph for this function: