76 BL_PROFILE_REGION(
"erf_substep_T()");
78 Real dtau =
static_cast<Real>(dtau_d);
80 const Box& domain = geom.Domain();
81 auto const domlo = lbound(domain);
82 auto const domhi = ubound(domain);
85 int ihi = domhi.x + 1;
87 int jhi = domhi.y + 1;
97 bool l_rayleigh_impl_for_w = (sinesq_stag_d !=
nullptr);
99 const Real*
dx = geom.CellSize();
100 const GpuArray<Real, AMREX_SPACEDIM>
dxInv = geom.InvCellSizeArray();
106 const auto& dm = S_stage_data[
IntVars::cons].DistributionMap();
108 MultiFab Delta_rho_u( convert(ba,IntVect(1,0,0)), dm, 1, 1);
109 MultiFab Delta_rho_v( convert(ba,IntVect(0,1,0)), dm, 1, 1);
110 MultiFab Delta_rho_w( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,0));
111 MultiFab Delta_rho ( ba , dm, 1, 1);
112 MultiFab Delta_rho_theta( ba , dm, 1, 1);
114 MultiFab New_rho_u(convert(ba,IntVect(1,0,0)), dm, 1, 1);
115 MultiFab New_rho_v(convert(ba,IntVect(0,1,0)), dm, 1, 1);
117 MultiFab coeff_A_mf(fast_coeffs, make_alias, 0, 1);
118 MultiFab inv_coeff_B_mf(fast_coeffs, make_alias, 1, 1);
119 MultiFab coeff_C_mf(fast_coeffs, make_alias, 2, 1);
120 MultiFab coeff_P_mf(fast_coeffs, make_alias, 3, 1);
121 MultiFab coeff_Q_mf(fast_coeffs, make_alias, 4, 1);
125 const Array<Real,AMREX_SPACEDIM> grav{
zero,
zero, -gravity};
126 const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
135 #pragma omp parallel if (Gpu::notInLaunchRegion())
137 for ( MFIter mfi(S_stage_data[
IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
139 const Array4<Real> & cur_cons = S_data[
IntVars::cons].array(mfi);
140 const Array4<const Real>& prev_cons = S_prev[
IntVars::cons].const_array(mfi);
141 const Array4<const Real>& stage_cons = S_stage_data[
IntVars::cons].const_array(mfi);
142 const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
144 const Array4<Real>& old_drho = Delta_rho.array(mfi);
145 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
146 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
147 const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
148 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
150 const Array4<const Real>& prev_xmom = S_prev[
IntVars::xmom].const_array(mfi);
151 const Array4<const Real>& prev_ymom = S_prev[
IntVars::ymom].const_array(mfi);
152 const Array4<const Real>& prev_zmom = S_prev[
IntVars::zmom].const_array(mfi);
154 const Array4<const Real>& stage_xmom = S_stage_data[
IntVars::xmom].const_array(mfi);
155 const Array4<const Real>& stage_ymom = S_stage_data[
IntVars::ymom].const_array(mfi);
156 const Array4<const Real>& stage_zmom = S_stage_data[
IntVars::zmom].const_array(mfi);
158 Box bx = mfi.validbox();
159 Box gbx = mfi.tilebox(); gbx.grow(1);
162 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
168 Box gtbx = mfi.nodaltilebox(0); gtbx.grow(IntVect(1,1,0));
169 Box gtby = mfi.nodaltilebox(1); gtby.grow(IntVect(1,1,0));
170 Box gtbz = mfi.nodaltilebox(2); gtbz.grow(IntVect(1,1,0));
172 const auto& bx_lo = lbound(bx);
173 const auto& bx_hi = ubound(bx);
176 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
177 old_drho_u(i,j,k) = prev_xmom(i,j,k) - stage_xmom(i,j,k);
178 if (k == bx_lo.z && k != domlo.z) {
179 old_drho_u(i,j,k-1) = old_drho_u(i,j,k);
180 }
else if (k == bx_hi.z) {
181 old_drho_u(i,j,k+1) = old_drho_u(i,j,k);
184 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
185 old_drho_v(i,j,k) = prev_ymom(i,j,k) - stage_ymom(i,j,k);
186 if (k == bx_lo.z && k != domlo.z) {
187 old_drho_v(i,j,k-1) = old_drho_v(i,j,k);
188 }
else if (k == bx_hi.z) {
189 old_drho_v(i,j,k+1) = old_drho_v(i,j,k);
192 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
193 old_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k);
196 const Array4<Real>& theta_extrap = extrap.array(mfi);
197 const Array4<const Real>& prim = S_stage_prim.const_array(mfi);
199 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
203 theta_extrap(i,j,k) = old_drho_theta(i,j,k);
205 theta_extrap(i,j,k) = old_drho_theta(i,j,k) + beta_d *
206 ( old_drho_theta(i,j,k) - lagged_arr(i,j,k) );
211 theta_extrap(i,j,k) *= (
one + RvOverRd*
qv);
216 #pragma omp parallel if (Gpu::notInLaunchRegion())
218 for ( MFIter mfi(S_stage_data[
IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
221 Box gbx = mfi.tilebox(); gbx.grow(1);
222 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
223 const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
224 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
225 lagged_arr(i,j,k) = old_drho_theta(i,j,k);
234 #pragma omp parallel if (Gpu::notInLaunchRegion())
236 for ( MFIter mfi(S_stage_data[
IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
238 Box bx = mfi.validbox();
239 Box tbx = mfi.nodaltilebox(0);
240 Box tby = mfi.nodaltilebox(1);
242 const Array4<Real const>& xmom_src_arr = xmom_src.const_array(mfi);
243 const Array4<Real const>& ymom_src_arr = ymom_src.const_array(mfi);
245 const Array4<const Real> & stage_xmom = S_stage_data[
IntVars::xmom].const_array(mfi);
246 const Array4<const Real> & stage_ymom = S_stage_data[
IntVars::ymom].const_array(mfi);
247 const Array4<const Real> & qt_arr =
qt.const_array(mfi);
249 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
250 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
252 const Array4<const Real>& slow_rhs_rho_u = S_slow_rhs[
IntVars::xmom].const_array(mfi);
253 const Array4<const Real>& slow_rhs_rho_v = S_slow_rhs[
IntVars::ymom].const_array(mfi);
255 const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
256 const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
258 const Array4<Real>& cur_xmom = S_data[
IntVars::xmom].array(mfi);
259 const Array4<Real>& cur_ymom = S_data[
IntVars::ymom].array(mfi);
262 const Array4<Real>& avg_xmom_arr = avg_xmom.array(mfi);
263 const Array4<Real>& avg_ymom_arr = avg_ymom.array(mfi);
265 const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
267 const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
269 const Array4<Real>& theta_extrap = extrap.array(mfi);
272 const Array4<const Real>& mf_ux = mapfac[
MapFacType::u_x]->const_array(mfi);
273 const Array4<const Real>& mf_vy = mapfac[
MapFacType::v_y]->const_array(mfi);
285 BL_PROFILE(
"substep_xymom_T");
287 const auto& bx_lo = lbound(bx);
288 const auto& bx_hi = ubound(bx);
291 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
296 Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi;
297 Real gp_zeta_on_iface = (k == 0) ?
298 myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
299 - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) :
300 fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
301 - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) );
302 Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ?
Real(0.) :
303 gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface;
307 Real q = (l_use_moisture) ?
myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) :
zero;
309 Real pi_c =
myhalf * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i ,j,k,0));
312 new_drho_u(i, j, k) = old_drho_u(i,j,k) + dtau * fast_rhs_rho_u
313 + dtau * slow_rhs_rho_u(i,j,k)
314 + dtau * xmom_src_arr(i,j,k);
315 if (k == bx_lo.z && k != domlo.z) {
316 new_drho_u(i,j,k-1) = new_drho_u(i,j,k);
317 }
else if (k == bx_hi.z) {
318 new_drho_u(i,j,k+1) = new_drho_u(i,j,k);
321 avg_xmom_arr(i,j,k) += facinv*new_drho_u(i,j,k);
323 cur_xmom(i,j,k) = stage_xmom(i,j,k) + new_drho_u(i,j,k);
325 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
330 Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi;
331 Real gp_zeta_on_jface = (k == 0) ?
332 myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
333 - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) :
334 fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
335 - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) );
336 Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ?
Real(0.) :
337 gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface;
341 Real q = (l_use_moisture) ?
myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) :
zero;
343 Real pi_c =
myhalf * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j ,k,0));
346 new_drho_v(i, j, k) = old_drho_v(i,j,k) + dtau * fast_rhs_rho_v
347 + dtau * slow_rhs_rho_v(i,j,k)
348 + dtau * ymom_src_arr(i,j,k);
350 if (k == bx_lo.z && k != domlo.z) {
351 new_drho_v(i,j,k-1) = new_drho_v(i,j,k);
352 }
else if (k == bx_hi.z) {
353 new_drho_v(i,j,k+1) = new_drho_v(i,j,k);
356 avg_ymom_arr(i,j,k) += facinv*new_drho_v(i,j,k);
358 cur_ymom(i,j,k) = stage_ymom(i,j,k) + new_drho_v(i,j,k);
366 #pragma omp parallel if (Gpu::notInLaunchRegion())
369 std::array<FArrayBox,AMREX_SPACEDIM> flux;
372 Box bx = mfi.tilebox();
373 Box tbz = surroundingNodes(bx,2);
375 Box vbx = mfi.validbox();
376 const auto& vbx_hi = ubound(vbx);
378 const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
379 const Array4<Real const>& cc_src_arr = cc_src.const_array(mfi);
381 const Array4<const Real> & stage_zmom = S_stage_data[
IntVars::zmom].const_array(mfi);
382 const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
384 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
385 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
386 const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
387 const Array4<Real>& old_drho = Delta_rho.array(mfi);
388 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
390 const Array4<const Real>& slow_rhs_cons = S_slow_rhs[
IntVars::cons].const_array(mfi);
391 const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[
IntVars::zmom].const_array(mfi);
393 const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
394 const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
396 const Array4<Real>& cur_cons = S_data[
IntVars::cons].array(mfi);
397 const Array4<Real>& cur_zmom = S_data[
IntVars::zmom].array(mfi);
400 const Array4<Real>& avg_zmom_arr = avg_zmom.array(mfi);
402 const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
403 const Array4<const Real>& detJ = detJ_cc->const_array(mfi);
405 const Array4< Real>& omega_arr = Omega.array(mfi);
408 const Array4<const Real>& mf_mx = mapfac[
MapFacType::m_x]->const_array(mfi);
409 const Array4<const Real>& mf_my = mapfac[
MapFacType::m_y]->const_array(mfi);
410 const Array4<const Real>& mf_ux = mapfac[
MapFacType::u_x]->const_array(mfi);
411 const Array4<const Real>& mf_uy = mapfac[
MapFacType::u_y]->const_array(mfi);
412 const Array4<const Real>& mf_vx = mapfac[
MapFacType::v_x]->const_array(mfi);
413 const Array4<const Real>& mf_vy = mapfac[
MapFacType::v_y]->const_array(mfi);
421 FArrayBox temp_rhs_fab;
425 RHS_fab.resize (tbz,1,The_Async_Arena());
426 soln_fab.resize (tbz,1,The_Async_Arena());
427 temp_rhs_fab.resize(tbz,2,The_Async_Arena());
429 auto const& RHS_a = RHS_fab.array();
430 auto const& soln_a = soln_fab.array();
431 auto const& temp_rhs_arr = temp_rhs_fab.array();
433 auto const& coeffA_a = coeff_A_mf.array(mfi);
434 auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
435 auto const& coeffC_a = coeff_C_mf.array(mfi);
436 auto const& coeffP_a = coeff_P_mf.array(mfi);
437 auto const& coeffQ_a = coeff_Q_mf.array(mfi);
442 for (
int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
443 flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
444 flux[dir].setVal<RunOn::Device>(0);
446 const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
447 flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
451 BL_PROFILE(
"fast_T_making_rho_rhs");
452 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
454 ( z_nd(i+1,j ,k+1) + z_nd(i+1,j+1,k+1)
455 -z_nd(i+1,j ,k ) - z_nd(i+1,j+1,k ) );
458 ( z_nd(i ,j ,k+1) + z_nd(i ,j+1,k+1)
459 -z_nd(i ,j ,k ) - z_nd(i ,j+1,k ) );
462 ( z_nd(i ,j+1,k+1) + z_nd(i+1,j+1,k+1)
463 -z_nd(i ,j+1,k ) - z_nd(i+1,j+1,k ) );
466 ( z_nd(i ,j ,k+1) + z_nd(i+1,j ,k+1)
467 -z_nd(i ,j ,k ) - z_nd(i+1,j ,k ) );
469 Real xflux_lo = new_drho_u(i ,j,k)*h_zeta_cc_xface_lo / mf_uy(i ,j,0);
470 Real xflux_hi = new_drho_u(i+1,j,k)*h_zeta_cc_xface_hi / mf_uy(i+1,j,0);
471 Real yflux_lo = new_drho_v(i,j ,k)*h_zeta_cc_yface_lo / mf_vx(i,j ,0);
472 Real yflux_hi = new_drho_v(i,j+1,k)*h_zeta_cc_yface_hi / mf_vx(i,j+1,0);
474 Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
477 temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq +
478 ( yflux_hi - yflux_lo ) * dyi * mfsq;
479 temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
480 xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+
481 ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
482 yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) *
myhalf;
485 (flx_arr[0])(i,j,k,0) = xflux_lo;
486 (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));
488 (flx_arr[1])(i,j,k,0) = yflux_lo;
489 (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));
492 (flx_arr[0])(i+1,j,k,0) = xflux_hi;
493 (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));
496 (flx_arr[1])(i,j+1,k,0) = yflux_hi;
497 (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));
505 Box gbxo = mfi.nodaltilebox(2);
508 if (gbxo.smallEnd(2) == domlo.z) {
509 Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2));
510 gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1);
511 ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
512 omega_arr(i,j,k) =
zero;
515 if (gbxo.bigEnd(2) == domhi.z+1) {
516 Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
517 gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
518 ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
519 omega_arr(i,j,k) = old_drho_w(i,j,k);
522 ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
523 omega_arr(i,j,k) =
OmegaFromW(i,j,k,old_drho_w(i,j,k),
524 old_drho_u,old_drho_v,
525 mf_ux,mf_vy,z_nd,
dxInv);
530 Box bx_shrunk_in_k = bx;
531 int klo = tbz.smallEnd(2);
532 int khi = tbz.bigEnd(2);
533 bx_shrunk_in_k.setSmall(2,klo+1);
534 bx_shrunk_in_k.setBig(2,
khi-1);
542 BL_PROFILE(
"fast_loop_on_shrunk_t");
544 ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
546 Real coeff_P = coeffP_a(i,j,k);
547 Real coeff_Q = coeffQ_a(i,j,k);
554 Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k )
555 -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1);
563 Real Omega_kp1 = omega_arr(i,j,k+1);
564 Real Omega_k = omega_arr(i,j,k );
565 Real Omega_km1 = omega_arr(i,j,k-1);
567 Real detJdiff = (detJ(i,j,k) - detJ(i,j,k-1)) / (detJ(i,j,k)*detJ(i,j,k-1));
570 R1_tmp += halfg * ( beta_1 * dzi * (Omega_kp1/detJ(i,j,k) + detJdiff*Omega_k - Omega_km1/detJ(i,j,k-1))
571 + 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) );
574 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) )
575 + 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) ) );
578 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);
582 new_drho_u,new_drho_v,
583 mf_ux,mf_vy,z_nd,
dxInv);
590 auto const lo = lbound(bx);
591 auto const hi = ubound(bx);
594 BL_PROFILE(
"substep_b2d_loop_t");
596 ParallelFor(b2d, [=] AMREX_GPU_DEVICE (
int i,
int j,
int)
599 RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z) + zmom_src_arr(i,j,lo.z));
600 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));
603 soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
606 for (
int k = lo.z+1; k <= hi.z+1; k++) {
607 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);
610 cur_zmom(i,j,lo.z ) = stage_zmom(i,j,lo.z ) + soln_a(i,j,lo.z );
611 cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
614 for (
int k = hi.z; k >= lo.z; k--) {
615 soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1);
619 for (
int j = lo.y; j <= hi.y; ++j) {
621 for (
int i = lo.x; i <= hi.x; ++i)
623 RHS_a(i,j,lo.z) = dtau * (slow_rhs_rho_w(i,j,lo.z) + zmom_src_arr(i,j,lo.z));
624 soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
628 for (
int i = lo.x; i <= hi.x; ++i)
630 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));
631 soln_a(i,j,hi.z+1) = RHS_a(i,j,hi.z+1) * inv_coeffB_a(i,j,hi.z+1);
635 for (
int k = lo.z+1; k <= hi.z; ++k) {
636 for (
int j = lo.y; j <= hi.y; ++j) {
638 for (
int i = lo.x; i <= hi.x; ++i) {
639 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);
643 for (
int k = hi.z; k > lo.z; --k) {
644 for (
int j = lo.y; j <= hi.y; ++j) {
646 for (
int i = lo.x; i <= hi.x; ++i) {
647 soln_a(i,j,k) -= (coeffC_a(i,j,k) * inv_coeffB_a(i,j,k)) * soln_a(i,j,k+1);
651 if (hi.z == domhi.z) {
652 for (
int j = lo.y; j <= hi.y; ++j) {
654 for (
int i = lo.x; i <= hi.x; ++i) {
655 cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
662 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
664 cur_zmom(i,j,k) = stage_zmom(i,j,k);
667 if (lo.z == domlo.z) {
668 tbz.setSmall(2,domlo.z+1);
670 if (hi.z == domhi.z) {
671 tbz.setBig(2,domhi.z);
673 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
676 new_drho_u,new_drho_v,
677 mf_ux,mf_vy,z_nd,
dxInv);
679 cur_zmom(i,j,k) += wpp;
681 if (l_rayleigh_impl_for_w) {
682 Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
683 cur_zmom(i,j,k) /= (
one + damping_coeff);
691 BL_PROFILE(
"fast_rho_final_update");
692 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
694 Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k);
695 Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1);
699 avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
701 (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
705 avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
707 (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
708 (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));
712 Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k);
714 cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho);
716 Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) +
myhalf *
717 ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) -
718 zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k);
720 cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta);
723 (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) *
myhalf * (prim(i,j,k) + prim(i,j,k-1));
734 int strt_comp_reflux = 0;
736 int num_comp_reflux = 1;
737 if (level < finest_level) {
738 fr_as_crse->CrseAdd(mfi,
739 {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
740 dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
743 fr_as_fine->FineAdd(mfi,
744 {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
745 dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
751 Gpu::streamSynchronize();
constexpr amrex::Real R_v
Definition: ERF_Constants.H:48
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
constexpr amrex::Real R_d
Definition: ERF_Constants.H:47
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:62
@ v_x
Definition: ERF_DataStruct.H:24
@ u_y
Definition: ERF_DataStruct.H:25
@ v_y
Definition: ERF_DataStruct.H:25
@ m_y
Definition: ERF_DataStruct.H:25
@ u_x
Definition: ERF_DataStruct.H:24
@ m_x
Definition: ERF_DataStruct.H:24
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:58
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:55
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
ParallelFor(grown_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);})
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:414
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
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
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
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
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
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
@ gpy
Definition: ERF_IndexDefines.H:186
@ gpx
Definition: ERF_IndexDefines.H:185
@ ymom
Definition: ERF_IndexDefines.H:195
@ cons
Definition: ERF_IndexDefines.H:193
@ zmom
Definition: ERF_IndexDefines.H:196
@ xmom
Definition: ERF_IndexDefines.H:194
@ qt
Definition: ERF_Kessler.H:29
@ qv
Definition: ERF_Kessler.H:30
@ q
Definition: ERF_WSM6.H:182