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

Functions

void erf_substep_MT (int step, int, int level, int finest_level, Vector< MultiFab > &S_slow_rhs, const Vector< MultiFab > &S_prev, Vector< MultiFab > &S_stg_data, const MultiFab &S_stg_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, const bool use_lagged_delta_rt, std::unique_ptr< MultiFab > &z_t_rk, const MultiFab *z_t_pert, std::unique_ptr< MultiFab > &z_phys_nd_old, std::unique_ptr< MultiFab > &z_phys_nd_new, std::unique_ptr< MultiFab > &z_phys_nd_stg, std::unique_ptr< MultiFab > &detJ_cc_old, std::unique_ptr< MultiFab > &detJ_cc_new, std::unique_ptr< MultiFab > &detJ_cc_stg, 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, const Real *sinesq_stag_d, const Real l_damp_coef)
 

Function Documentation

◆ erf_substep_MT()

void erf_substep_MT ( int  step,
int  ,
int  level,
int  finest_level,
Vector< MultiFab > &  S_slow_rhs,
const Vector< MultiFab > &  S_prev,
Vector< MultiFab > &  S_stg_data,
const MultiFab &  S_stg_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,
const bool  use_lagged_delta_rt,
std::unique_ptr< MultiFab > &  z_t_rk,
const MultiFab *  z_t_pert,
std::unique_ptr< MultiFab > &  z_phys_nd_old,
std::unique_ptr< MultiFab > &  z_phys_nd_new,
std::unique_ptr< MultiFab > &  z_phys_nd_stg,
std::unique_ptr< MultiFab > &  detJ_cc_old,
std::unique_ptr< MultiFab > &  detJ_cc_new,
std::unique_ptr< MultiFab > &  detJ_cc_stg,
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  ,
const Real *  sinesq_stag_d,
const Real  l_damp_coef 
)

Function for computing the fast RHS with moving 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_stg_datasolution at previous RK stage
[in]S_stg_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]use_lagged_delta_rtdefine lagged_delta_rt for our next step
[in]z_t_rkrate of change of grid height – only relevant for moving terrain
[in]z_t_pertrate of change of grid height – interpolated between RK stages
[in]z_phys_nd_oldheight coordinate at nodes at old time
[in]z_phys_nd_newheight coordinate at nodes at new time
[in]z_phys_nd_stgheight coordinate at nodes at previous stage
[in]detJ_cc_oldJacobian of the metric transformation at old time
[in]detJ_cc_newJacobian of the metric transformation at new time
[in]detJ_cc_stgJacobian of the metric transformation at previous stage
[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
89 {
90  BL_PROFILE_REGION("erf_substep_MT()");
91 
92  Real dtau = static_cast<Real>(dtau_d);
93 
94  Real beta_1 = myhalf * (one - beta_s); // multiplies explicit terms
95  Real beta_2 = myhalf * (one + beta_s); // multiplies implicit terms
96 
97  // How much do we project forward the (rho theta) that is used in the horizontal momentum equations
98  Real beta_d = Real(0.1);
99 
100  Real RvOverRd = R_v / R_d;
101 
102  bool l_rayleigh_impl_for_w = (sinesq_stag_d != nullptr);
103 
104  const Real* dx = geom.CellSize();
105  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
106 
107  Real dxi = dxInv[0];
108  Real dyi = dxInv[1];
109  Real dzi = dxInv[2];
110 
111  MultiFab coeff_A_mf(fast_coeffs, make_alias, 0, 1);
112  MultiFab inv_coeff_B_mf(fast_coeffs, make_alias, 1, 1);
113  MultiFab coeff_C_mf(fast_coeffs, make_alias, 2, 1);
114  MultiFab coeff_P_mf(fast_coeffs, make_alias, 3, 1);
115  MultiFab coeff_Q_mf(fast_coeffs, make_alias, 4, 1);
116 
117  // *************************************************************************
118  // Set gravity as a vector
119  const Array<Real,AMREX_SPACEDIM> grav{zero, zero, -gravity};
120  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
121 
122  MultiFab extrap(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(),1,1);
123 
124  MultiFab Omega(S_data[IntVars::zmom].boxArray(), S_data[IntVars::zmom].DistributionMap(), 1, 1);
125 
126  // *************************************************************************
127  // Define updates in the current RK stg
128  // *************************************************************************
129 #ifdef _OPENMP
130 #pragma omp parallel if (Gpu::notInLaunchRegion())
131 #endif
132  {
133  FArrayBox temp_rhs_fab;
134 
135  FArrayBox RHS_fab;
136  FArrayBox soln_fab;
137 
138  std::array<FArrayBox,AMREX_SPACEDIM> flux;
139 
140  // NOTE: we leave tiling off here for efficiency -- to make this loop work with tiling
141  // will require additional changes
142  for ( MFIter mfi(S_stg_data[IntVars::cons],false); mfi.isValid(); ++mfi)
143  {
144  Box bx = mfi.tilebox();
145  Box tbx = surroundingNodes(bx,0);
146  Box tby = surroundingNodes(bx,1);
147  Box tbz = surroundingNodes(bx,2);
148 
149  Box vbx = mfi.validbox();
150  const auto& vbx_hi = ubound(vbx);
151 
152  const Array4<Real const>& xmom_src_arr = xmom_src.const_array(mfi);
153  const Array4<Real const>& ymom_src_arr = ymom_src.const_array(mfi);
154  const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
155  const Array4<Real const>& cc_src_arr = cc_src.const_array(mfi);
156 
157  const Array4<const Real> & stg_cons = S_stg_data[IntVars::cons].const_array(mfi);
158  const Array4<const Real> & stg_xmom = S_stg_data[IntVars::xmom].const_array(mfi);
159  const Array4<const Real> & stg_ymom = S_stg_data[IntVars::ymom].const_array(mfi);
160  const Array4<const Real> & stg_zmom = S_stg_data[IntVars::zmom].const_array(mfi);
161  const Array4<const Real> & prim = S_stg_prim.const_array(mfi);
162  const Array4<const Real> & qt_arr = qt.const_array(mfi);
163 
164  const Array4<const Real>& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi);
165  const Array4<const Real>& slow_rhs_rho_u = S_slow_rhs[IntVars::xmom].const_array(mfi);
166  const Array4<const Real>& slow_rhs_rho_v = S_slow_rhs[IntVars::ymom].const_array(mfi);
167  const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi);
168 
169  const Array4<Real>& cur_cons = S_data[IntVars::cons].array(mfi);
170  const Array4<Real>& cur_xmom = S_data[IntVars::xmom].array(mfi);
171  const Array4<Real>& cur_ymom = S_data[IntVars::ymom].array(mfi);
172  const Array4<Real>& cur_zmom = S_data[IntVars::zmom].array(mfi);
173 
174  const Array4<Real>& lagged = lagged_delta_rt.array(mfi);
175 
176  const Array4<const Real>& prev_cons = S_prev[IntVars::cons].const_array(mfi);
177  const Array4<const Real>& prev_xmom = S_prev[IntVars::xmom].const_array(mfi);
178  const Array4<const Real>& prev_ymom = S_prev[IntVars::ymom].const_array(mfi);
179  const Array4<const Real>& prev_zmom = S_prev[IntVars::zmom].const_array(mfi);
180 
181  // These store the advection momenta which we will use to update the slow variables
182  const Array4<Real>& avg_xmom_arr = avg_xmom.array(mfi);
183  const Array4<Real>& avg_ymom_arr = avg_ymom.array(mfi);
184  const Array4<Real>& avg_zmom_arr = avg_zmom.array(mfi);
185 
186  const Array4<const Real>& z_nd_old = z_phys_nd_old->const_array(mfi);
187  const Array4<const Real>& z_nd_new = z_phys_nd_new->const_array(mfi);
188  const Array4<const Real>& z_nd_stg = z_phys_nd_stg->const_array(mfi);
189  const Array4<const Real>& detJ_old = detJ_cc_old->const_array(mfi);
190  const Array4<const Real>& detJ_new = detJ_cc_new->const_array(mfi);
191  const Array4<const Real>& detJ_stg = detJ_cc_stg->const_array(mfi);
192 
193  const Array4<const Real>& z_t_arr = z_t_rk->const_array(mfi);
194  const Array4<const Real>& zp_t_arr = z_t_pert->const_array(mfi);
195 
196  const Array4< Real>& omega_arr = Omega.array(mfi);
197 
198  // Map factors
199  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
200  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
201  const Array4<const Real>& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi);
202  const Array4<const Real>& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi);
203 
204  // *********************************************************************
205  // This must be done before we set cur_xmom and cur_ymom, since those
206  // in fact point to the same array as prev_xmom and prev_ymom
207  // *********************************************************************
208  Box gbxo = mfi.nodaltilebox(2);
209  {
210  BL_PROFILE("fast_MT_making_omega");
211  Box gbxo_lo = gbxo; gbxo_lo.setBig(2,0);
212  ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
213  omega_arr(i,j,k) = zero;
214  });
215  Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
216  ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
217  omega_arr(i,j,k) = prev_zmom(i,j,k) - stg_zmom(i,j,k) - zp_t_arr(i,j,k);
218  });
219 
220  // gbxo_lo covers k<=0 and gbxo_hi covers the box's top face, so the mid box is
221  // [max(smallEnd,1), bigEnd-1]. Using max() never EXPANDS the box past the FAB,
222  // which is what setSmall(2,1) does for a box that starts above k=1.
223  Box gbxo_mid = gbxo;
224  gbxo_mid.setSmall(2, std::max(gbxo.smallEnd(2), 1));
225  gbxo_mid.setBig (2, gbxo.bigEnd(2)-1);
226 
227  ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
228  omega_arr(i,j,k) =
229  ( OmegaFromW(i,j,k,prev_zmom(i,j,k),prev_xmom,prev_ymom,mf_ux,mf_vy,z_nd_old,dxInv)
230  -OmegaFromW(i,j,k, stg_zmom(i,j,k), stg_xmom, stg_ymom,mf_ux,mf_vy,z_nd_old,dxInv) )
231  - zp_t_arr(i,j,k);
232  });
233  } // end profile
234  // *********************************************************************
235 
236  const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
237 
238  const Array4<Real>& theta_extrap = extrap.array(mfi);
239 
240  // Note: it is important to grow the tilebox rather than use growntilebox because
241  // we need to fill the ghost cells of the tilebox so we can use them below
242  Box gbx = mfi.tilebox(); gbx.grow(1);
243  Box gtbx = mfi.nodaltilebox(0); gtbx.grow(1); gtbx.setSmall(2,0);
244  Box gtby = mfi.nodaltilebox(1); gtby.grow(1); gtby.setSmall(2,0);
245 
246  if (step == 0) {
247  ParallelFor(gbx,
248  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
249  cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp);
250  cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp);
251 
252  Real delta_rt = cur_cons(i,j,k,RhoTheta_comp) - stg_cons(i,j,k,RhoTheta_comp);
253  theta_extrap(i,j,k) = delta_rt;
254 
255  // NOTE: qv is not changing over the fast steps so we use the stage data
256  Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero;
257  theta_extrap(i,j,k) *= (one + RvOverRd*qv);
258 
259  // We define lagged_delta_rt for our next step as the current delta_rt
260  lagged(i,j,k) = delta_rt;
261  });
262  } else if (use_lagged_delta_rt) {
263  // This is the default for cases with no or static terrain
264  ParallelFor(gbx,
265  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
266  Real delta_rt = cur_cons(i,j,k,RhoTheta_comp) - stg_cons(i,j,k,RhoTheta_comp);
267  theta_extrap(i,j,k) = delta_rt + beta_d * (delta_rt - lagged(i,j,k));
268 
269  // NOTE: qv is not changing over the fast steps so we use the stage data
270  Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero;
271  theta_extrap(i,j,k) *= (one + RvOverRd*qv);
272 
273  // We define lagged_delta_rt for our next step as the current delta_rt
274  lagged(i,j,k) = delta_rt;
275  });
276  } else {
277  // For the moving wave problem, this choice seems more robust
278  ParallelFor(gbx,
279  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
280  theta_extrap(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) - stg_cons(i,j,k,RhoTheta_comp);
281 
282  // NOTE: qv is not changing over the fast steps so we use the stage data
283  Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero;
284  theta_extrap(i,j,k) *= (one + RvOverRd*qv);
285  });
286  } // if step
287 
288  RHS_fab.resize (tbz,1, The_Async_Arena());
289  soln_fab.resize (tbz,1, The_Async_Arena());
290  temp_rhs_fab.resize(tbz,2, The_Async_Arena());
291 
292  auto const& RHS_a = RHS_fab.array();
293  auto const& soln_a = soln_fab.array();
294  auto const& temp_rhs_arr = temp_rhs_fab.array();
295 
296  auto const& coeffA_a = coeff_A_mf.array(mfi);
297  auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
298  auto const& coeffC_a = coeff_C_mf.array(mfi);
299  auto const& coeffP_a = coeff_P_mf.array(mfi);
300  auto const& coeffQ_a = coeff_Q_mf.array(mfi);
301 
302  // *********************************************************************
303  // Define updates in the RHS of {x, y, z}-momentum equations
304  // *********************************************************************
305  {
306  BL_PROFILE("substep_xymom_T");
307  ParallelFor(tbx, tby,
308  [=] AMREX_GPU_DEVICE (int i, int j, int k)
309  {
310  // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
311  Real h_xi_old = Compute_h_xi_AtIface(i, j, k, dxInv, z_nd_old);
312  Real h_zeta_old = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd_old);
313  Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi;
314  Real gp_zeta_on_iface = (k == 0) ?
315  myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
316  - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) :
317  fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
318  - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) );
319  Real gpx = h_zeta_old * gp_xi - h_xi_old * gp_zeta_on_iface;
320  gpx *= mf_ux(i,j,0);
321 
322  Real q = (l_use_moisture) ? myhalf * (qt_arr(i-1,j,k) + qt_arr(i,j,k)) : zero;
323 
324  Real pi_c = myhalf * (pi_stage_ca(i-1,j,k) + pi_stage_ca(i ,j,k));
325  Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx / (one + q);
326 
327  // We have already scaled the source terms to have the extra factor of dJ
328  cur_xmom(i,j,k) = h_zeta_old * prev_xmom(i,j,k) + dtau * fast_rhs_rho_u
329  + dtau * slow_rhs_rho_u(i,j,k)
330  + dtau * xmom_src_arr(i,j,k);
331  },
332  [=] AMREX_GPU_DEVICE (int i, int j, int k)
333  {
334  // Add (negative) gradient of (rho theta) multiplied by lagged "pi"
335  Real h_eta_old = Compute_h_eta_AtJface(i, j, k, dxInv, z_nd_old);
336  Real h_zeta_old = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd_old);
337  Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi;
338  Real gp_zeta_on_jface = (k == 0) ?
339  myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
340  - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) :
341  fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
342  - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) );
343  Real gpy = h_zeta_old * gp_eta - h_eta_old * gp_zeta_on_jface;
344  gpy *= mf_vy(i,j,0);
345 
346  Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j-1,k) + qt_arr(i,j,k)) : zero;
347 
348  Real pi_c = myhalf * (pi_stage_ca(i,j-1,k) + pi_stage_ca(i,j ,k));
349  Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy / (one + q);
350 
351  // We have already scaled the source terms to have the extra factor of dJ
352  cur_ymom(i, j, k) = h_zeta_old * prev_ymom(i,j,k) + dtau * fast_rhs_rho_v
353  + dtau * slow_rhs_rho_v(i,j,k)
354  + dtau * ymom_src_arr(i,j,k);
355  });
356  } // end profile
357 
358  // *************************************************************************
359  // Define flux arrays for use in advection
360  // *************************************************************************
361  for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
362  flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
363  flux[dir].setVal<RunOn::Device>(0);
364  }
365  const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
366  flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
367 
368  // *********************************************************************
369  {
370  BL_PROFILE("fast_T_making_rho_rhs");
371  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
372  {
373  Real h_zeta_stg_xlo = Compute_h_zeta_AtIface(i, j , k, dxInv, z_nd_stg);
374  Real h_zeta_stg_xhi = Compute_h_zeta_AtIface(i+1,j , k, dxInv, z_nd_stg);
375  Real xflux_lo = cur_xmom(i ,j,k) - stg_xmom(i ,j,k)*h_zeta_stg_xlo;
376  Real xflux_hi = cur_xmom(i+1,j,k) - stg_xmom(i+1,j,k)*h_zeta_stg_xhi;
377 
378  Real h_zeta_stg_yhi = Compute_h_zeta_AtJface(i, j+1, k, dxInv, z_nd_stg);
379  Real h_zeta_stg_ylo = Compute_h_zeta_AtJface(i, j , k, dxInv, z_nd_stg);
380  Real yflux_lo = cur_ymom(i,j ,k) - stg_ymom(i,j ,k)*h_zeta_stg_ylo;
381  Real yflux_hi = cur_ymom(i,j+1,k) - stg_ymom(i,j+1,k)*h_zeta_stg_yhi;
382 
383  // NOTE: we are saving the (1/J) weighting for later when we add this to rho and theta
384  temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi + ( yflux_hi - yflux_lo ) * dyi;
385  temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
386  xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi +
387  ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
388  yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi) * myhalf;
389 
390  if (l_reflux) {
391  (flx_arr[0])(i,j,k,0) = xflux_lo;
392  (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));
393 
394  (flx_arr[1])(i,j,k,0) = yflux_lo;
395  (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));
396 
397  if (i == vbx_hi.x) {
398  (flx_arr[0])(i+1,j,k,0) = xflux_hi;
399  (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));
400  }
401  if (j == vbx_hi.y) {
402  (flx_arr[1])(i,j+1,k,0) = yflux_hi;
403  (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));
404  }
405  }
406  });
407  } // end profile
408 
409  ParallelFor(tbx, tby,
410  [=] AMREX_GPU_DEVICE (int i, int j, int k)
411  {
412  Real h_zeta_new = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd_new);
413  cur_xmom(i, j, k) /= h_zeta_new;
414  avg_xmom_arr(i,j,k) += facinv*(cur_xmom(i,j,k) - stg_xmom(i,j,k));
415  },
416  [=] AMREX_GPU_DEVICE (int i, int j, int k)
417  {
418  Real h_zeta_new = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd_new);
419  cur_ymom(i, j, k) /= h_zeta_new;
420  avg_ymom_arr(i,j,k) += facinv*(cur_ymom(i,j,k) - stg_ymom(i,j,k));
421  });
422 
423  Box bx_shrunk_in_k = bx;
424  int klo = tbz.smallEnd(2);
425  int khi = tbz.bigEnd(2);
426  bx_shrunk_in_k.setSmall(2,klo+1);
427  bx_shrunk_in_k.setBig(2,khi-1);
428 
429  // Note that the notes use "g" to mean the magnitude of gravity, so it is positive
430  // We set grav_gpu[2] to be the vector component which is negative
431  // We define halfg to match the notes (which is why we take the absolute value)
432  Real halfg = std::abs(myhalf * grav_gpu[2]);
433 
434  {
435  BL_PROFILE("fast_loop_on_shrunk_t");
436  //Note we don't act on the bottom or top boundaries of the domain
437  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
438  {
439  Real dJ_old_kface = myhalf * (detJ_old(i,j,k) + detJ_old(i,j,k-1));
440  Real dJ_new_kface = myhalf * (detJ_new(i,j,k) + detJ_new(i,j,k-1));
441  Real dJ_stg_kface = myhalf * (detJ_stg(i,j,k) + detJ_stg(i,j,k-1));
442 
443  Real coeff_P = coeffP_a(i,j,k);
444  Real coeff_Q = coeffQ_a(i,j,k);
445 
446  Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
447  Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
448  Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
449 
450  // line 2 last two terms (order dtau)
451  Real R0_tmp = coeff_P * cur_cons(i,j,k ,RhoTheta_comp)
452  + coeff_Q * cur_cons(i,j,k-1,RhoTheta_comp)
453  - coeff_P * stg_cons(i,j,k ,RhoTheta_comp) * (dJ_stg_kface/dJ_old_kface)
454  - coeff_Q * stg_cons(i,j,k-1,RhoTheta_comp) * (dJ_stg_kface/dJ_old_kface)
455  - halfg * ( cur_cons(i,j,k,Rho_comp) + cur_cons(i,j,k-1,Rho_comp) )
456  + halfg * ( stg_cons(i,j,k,Rho_comp) + stg_cons(i,j,k-1,Rho_comp) ) * (dJ_stg_kface/dJ_old_kface);
457 
458  // line 3 residuals (order dtau^2) one <-> beta_2
459  Real R1_tmp = - halfg * ( slow_rhs_cons(i,j,k ,Rho_comp) + slow_rhs_cons(i,j,k-1,Rho_comp))
460  + coeff_P * slow_rhs_cons(i,j,k ,RhoTheta_comp) + coeff_Q * slow_rhs_cons(i,j,k-1,RhoTheta_comp);
461 
462  Real Omega_kp1 = omega_arr(i,j,k+1);
463  Real Omega_k = omega_arr(i,j,k );
464  Real Omega_km1 = omega_arr(i,j,k-1);
465 
466  Real detJdiff = (detJ_old(i,j,k) - detJ_old(i,j,k-1)) / (detJ_old(i,j,k)*detJ_old(i,j,k-1));
467 
468  // consolidate lines 4&5 (order dtau^2)
469  R1_tmp += halfg * ( beta_1 * dzi * (Omega_kp1/detJ_old(i,j,k) + detJdiff*Omega_k - Omega_km1/detJ_old(i,j,k-1))
470  + temp_rhs_arr(i,j,k,Rho_comp)/detJ_old(i,j,k) + temp_rhs_arr(i,j,k-1,Rho_comp)/detJ_old(i,j,k-1) );
471 
472  // consolidate lines 6&7 (order dtau^2)
473  R1_tmp += -(
474  coeff_P/detJ_old(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*theta_t_hi - Omega_k*theta_t_mid)
475  +temp_rhs_arr(i,j,k ,RhoTheta_comp) ) +
476  coeff_Q/detJ_old(i,j,k-1) * ( beta_1 * dzi * (Omega_k*theta_t_mid - Omega_km1*theta_t_lo)
477  +temp_rhs_arr(i,j,k-1,RhoTheta_comp) ) );
478 
479  // line 1
480  RHS_a(i,j,k) = prev_zmom(i,j,k) - (dJ_stg_kface/dJ_old_kface) * stg_zmom(i,j,k)
481  + dtau * slow_rhs_rho_w(i,j,k) / dJ_stg_kface
482  + dtau * zmom_src_arr(i,j,k);
483 
484  RHS_a(i,j,k) += dtau * R0_tmp;
485 
486  RHS_a(i,j,k) += dtau * dtau*beta_2*R1_tmp;
487 
488  // We cannot use omega_arr here since that was built with old_rho_u and old_rho_v ...
489  Real UppVpp = (dJ_new_kface/dJ_old_kface) * OmegaFromW(i,j,k,0.,cur_xmom,cur_ymom,mf_ux,mf_vy,z_nd_new,dxInv)
490  -(dJ_stg_kface/dJ_old_kface) * OmegaFromW(i,j,k,0.,stg_xmom,stg_ymom,mf_ux,mf_vy,z_nd_stg,dxInv);
491  RHS_a(i,j,k) += UppVpp;
492  });
493  } // end profile
494 
495  Box b2d = tbz; // Copy constructor
496  b2d.setRange(2,0);
497 
498  auto const lo = lbound(bx);
499  auto const hi = ubound(bx);
500 
501  {
502  BL_PROFILE("substep_b2d_loop_t");
503 
504 #ifdef AMREX_USE_GPU
505  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
506  {
507  // Moving terrain
508  Real rho_on_bdy = myhalf * ( prev_cons(i,j,lo.z) + prev_cons(i,j,lo.z-1) );
509  RHS_a(i,j,lo.z) = rho_on_bdy * zp_t_arr(i,j,lo.z);
510 
511  soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
512 
513  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));
514 
515  for (int k = lo.z+1; k <= hi.z+1; k++) {
516  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);
517  }
518 
519  for (int k = hi.z; k >= lo.z; k--) {
520  soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
521  }
522 
523  // We assume that Omega == w at the top boundary and that changes in J there are irrelevant
524  cur_zmom(i,j,hi.z+1) = stg_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
525  });
526 #else
527  for (int j = lo.y; j <= hi.y; ++j) {
528  AMREX_PRAGMA_SIMD
529  for (int i = lo.x; i <= hi.x; ++i) {
530 
531  Real rho_on_bdy = myhalf * ( prev_cons(i,j,lo.z) + prev_cons(i,j,lo.z-1) );
532  RHS_a(i,j,lo.z) = rho_on_bdy * zp_t_arr(i,j,lo.z);
533 
534  soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
535  }
536  }
537 
538  for (int j = lo.y; j <= hi.y; ++j) {
539  AMREX_PRAGMA_SIMD
540  for (int i = lo.x; i <= hi.x; ++i) {
541  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));
542  }
543  }
544  for (int k = lo.z+1; k <= hi.z+1; ++k) {
545  for (int j = lo.y; j <= hi.y; ++j) {
546  AMREX_PRAGMA_SIMD
547  for (int i = lo.x; i <= hi.x; ++i) {
548  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);
549  }
550  }
551  }
552  for (int k = hi.z; k >= lo.z; --k) {
553  for (int j = lo.y; j <= hi.y; ++j) {
554  AMREX_PRAGMA_SIMD
555  for (int i = lo.x; i <= hi.x; ++i) {
556  soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
557  }
558  }
559  }
560 
561  // We assume that Omega == w at the top boundary and that changes in J there are irrelevant
562  for (int j = lo.y; j <= hi.y; ++j) {
563  AMREX_PRAGMA_SIMD
564  for (int i = lo.x; i <= hi.x; ++i) {
565  cur_zmom(i,j,hi.z+1) = stg_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
566  }
567  }
568 #endif
569  } // end profile
570 
571  {
572  BL_PROFILE("substep_new_drhow");
573  tbz.setBig(2,hi.z);
574  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
575  {
576  Real rho_on_face = myhalf * (cur_cons(i,j,k,Rho_comp) + cur_cons(i,j,k-1,Rho_comp));
577 
578  if (k == lo.z) {
579  cur_zmom(i,j,k) = WFromOmega(i,j,k,rho_on_face*(z_t_arr(i,j,k)+zp_t_arr(i,j,k)),
580  cur_xmom,cur_ymom,mf_ux,mf_vy,z_nd_new,dxInv);
581 
582  // We need to set this here because it is used to define zflux_lo below
583  soln_a(i,j,k) = zero;
584 
585  } else {
586 
587  Real UppVpp = WFromOmega(i,j,k,zero,cur_xmom,cur_ymom,mf_ux,mf_vy,z_nd_new,dxInv)
588  - WFromOmega(i,j,k,zero,stg_xmom,stg_ymom,mf_ux,mf_vy,z_nd_stg,dxInv);
589  Real wpp = soln_a(i,j,k) + UppVpp;
590  Real dJ_old_kface = myhalf * (detJ_old(i,j,k) + detJ_old(i,j,k-1));
591  Real dJ_new_kface = myhalf * (detJ_new(i,j,k) + detJ_new(i,j,k-1));
592 
593  cur_zmom(i,j,k) = dJ_old_kface * (stg_zmom(i,j,k) + wpp);
594  cur_zmom(i,j,k) /= dJ_new_kface;
595 
596  soln_a(i,j,k) = OmegaFromW(i,j,k,cur_zmom(i,j,k),cur_xmom,cur_ymom,mf_ux,mf_vy,z_nd_new,dxInv)
597  - OmegaFromW(i,j,k,stg_zmom(i,j,k),stg_xmom,stg_ymom,mf_ux,mf_vy,z_nd_stg,dxInv);
598  soln_a(i,j,k) -= rho_on_face * zp_t_arr(i,j,k);
599  }
600 
601  if (l_rayleigh_impl_for_w && k > 0) {
602  Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
603  cur_zmom(i,j,k) /= (one + damping_coeff);
604  }
605  });
606  } // end profile
607 
608  // **************************************************************************
609  // Define updates in the RHS of rho and (rho theta)
610  // **************************************************************************
611  {
612  BL_PROFILE("fast_rho_final_update");
613  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
614  {
615  Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k);
616  Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1);
617 
618  // Note that in the solve we effectively impose new_drho_w(i,j,vbx_hi.z+1)=0
619  // so we don't update avg_zmom at k=vbx_hi.z+1
620  avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
621  if (l_reflux) {
622  (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
623  }
624 
625  // Note that the factor of (1/J) in the fast source term is canceled
626  // when we multiply old and new by detJ_old and detJ_new , respectively
627  // We have already scaled the slow source term to have the extra factor of dJ
628  Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi);
629  Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) + myhalf *
630  ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1))
631  - zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi );
632 
633  cur_cons(i,j,k,0) *= (detJ_old(i,j,k)/detJ_new(i,j,k));
634  cur_cons(i,j,k,1) *= (detJ_old(i,j,k)/detJ_new(i,j,k));
635 
636  cur_cons(i,j,k,0) += dtau * ( slow_rhs_cons(i,j,k,0) + fast_rhs_rho / detJ_new(i,j,k));
637  cur_cons(i,j,k,1) += dtau * ( slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta / detJ_new(i,j,k));
638 
639  if (l_reflux) {
640  (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k) + prim(i,j,k-1));
641  }
642 
643  if (k == vbx_hi.z) {
644  avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
645  if (l_reflux) {
646  (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
647  (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));
648  }
649  }
650 
651  // add in source terms for cell-centered conserved variables
652  cur_cons(i,j,k,Rho_comp) += dtau * cc_src_arr(i,j,k,Rho_comp);
653  cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp);
654  });
655  } // end profile
656 
657  // We only add to the flux registers in the final RK step
658  if (l_reflux) {
659  int strt_comp_reflux = 0;
660  // For now we don't reflux (rho theta) because it seems to create issues at c/f boundaries
661  int num_comp_reflux = 1;
662  if (level < finest_level) {
663  fr_as_crse->CrseAdd(mfi,
664  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
665  dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
666  }
667  if (level > 0) {
668  fr_as_fine->FineAdd(mfi,
669  {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
670  dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
671  }
672 
673  // This is necessary here so we don't go on to the next FArrayBox without
674  // having finished copying the fluxes into the FluxRegisters (since the fluxes
675  // are stored in temporary FArrayBox's)
676  Gpu::streamSynchronize();
677 
678  } // two-way coupling
679 
680  } // mfi
681  } // OMP
682 }
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:48
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:47
Gamma
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:62
v_y
@ v_y
Definition: ERF_DataStruct.H:28
m_y
@ m_y
Definition: ERF_DataStruct.H:28
u_x
@ u_x
Definition: ERF_DataStruct.H:27
m_x
@ m_x
Definition: ERF_DataStruct.H:27
PrimQ1_comp
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:58
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
PrimTheta_comp
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:55
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
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);})
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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:414
Compute_h_xi_AtIface
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:117
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:104
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:144
Compute_h_eta_AtJface
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:170
WFromOmega
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:464
GpVars::gpy
@ gpy
Definition: ERF_IndexDefines.H:187
GpVars::gpx
@ gpx
Definition: ERF_IndexDefines.H:186
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:196
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:194
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:197
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:195
MicVar_Kess::qt
@ qt
Definition: ERF_Kessler.H:29
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:30
WSM6Ind::q
@ q
Definition: ERF_WSM6.H:184
Here is the call graph for this function: