77 BL_PROFILE_REGION(
"erf_substep_T()");
79 Real dtau =
static_cast<Real>(dtau_d);
81 const Box& domain = geom.Domain();
82 auto const domlo = lbound(domain);
83 auto const domhi = ubound(domain);
86 int ihi = domhi.x + 1;
88 int jhi = domhi.y + 1;
98 bool l_rayleigh_impl_for_w = (sinesq_stag_d !=
nullptr);
100 const Real*
dx = geom.CellSize();
101 const GpuArray<Real, AMREX_SPACEDIM>
dxInv = geom.InvCellSizeArray();
107 const auto& dm = S_stage_data[
IntVars::cons].DistributionMap();
109 MultiFab Delta_rho_u( convert(ba,IntVect(1,0,0)), dm, 1, 1);
110 MultiFab Delta_rho_v( convert(ba,IntVect(0,1,0)), dm, 1, 1);
111 MultiFab Delta_rho_w( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,0));
112 MultiFab Delta_rho ( ba , dm, 1, 1);
113 MultiFab Delta_rho_theta( ba , dm, 1, 1);
115 MultiFab New_rho_u(convert(ba,IntVect(1,0,0)), dm, 1, 1);
116 MultiFab New_rho_v(convert(ba,IntVect(0,1,0)), dm, 1, 1);
118 MultiFab coeff_A_mf(fast_coeffs, make_alias, 0, 1);
119 MultiFab inv_coeff_B_mf(fast_coeffs, make_alias, 1, 1);
120 MultiFab coeff_C_mf(fast_coeffs, make_alias, 2, 1);
121 MultiFab coeff_P_mf(fast_coeffs, make_alias, 3, 1);
122 MultiFab coeff_Q_mf(fast_coeffs, make_alias, 4, 1);
126 const Array<Real,AMREX_SPACEDIM> grav{
zero,
zero, -gravity};
127 const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
136 #pragma omp parallel if (Gpu::notInLaunchRegion())
138 for ( MFIter mfi(S_stage_data[
IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
140 const Array4<Real> & cur_cons = S_data[
IntVars::cons].array(mfi);
141 const Array4<const Real>& prev_cons = S_prev[
IntVars::cons].const_array(mfi);
142 const Array4<const Real>& stage_cons = S_stage_data[
IntVars::cons].const_array(mfi);
143 const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
145 const Array4<Real>& old_drho = Delta_rho.array(mfi);
146 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
147 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
148 const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
149 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
151 const Array4<const Real>& prev_xmom = S_prev[
IntVars::xmom].const_array(mfi);
152 const Array4<const Real>& prev_ymom = S_prev[
IntVars::ymom].const_array(mfi);
153 const Array4<const Real>& prev_zmom = S_prev[
IntVars::zmom].const_array(mfi);
155 const Array4<const Real>& stage_xmom = S_stage_data[
IntVars::xmom].const_array(mfi);
156 const Array4<const Real>& stage_ymom = S_stage_data[
IntVars::ymom].const_array(mfi);
157 const Array4<const Real>& stage_zmom = S_stage_data[
IntVars::zmom].const_array(mfi);
159 Box bx = mfi.validbox();
168 Box gbx = mfi.growntilebox(1);
171 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
177 Box gtbx = mfi.grownnodaltilebox(0,IntVect(1,1,0));
178 Box gtby = mfi.grownnodaltilebox(1,IntVect(1,1,0));
179 Box gtbz = mfi.grownnodaltilebox(2,IntVect(1,1,0));
181 const auto& bx_lo = lbound(bx);
182 const auto& bx_hi = ubound(bx);
185 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
186 old_drho_u(i,j,k) = prev_xmom(i,j,k) - stage_xmom(i,j,k);
187 if (k == bx_lo.z && k != domlo.z) {
188 old_drho_u(i,j,k-1) = old_drho_u(i,j,k);
189 }
else if (k == bx_hi.z) {
190 old_drho_u(i,j,k+1) = old_drho_u(i,j,k);
193 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
194 old_drho_v(i,j,k) = prev_ymom(i,j,k) - stage_ymom(i,j,k);
195 if (k == bx_lo.z && k != domlo.z) {
196 old_drho_v(i,j,k-1) = old_drho_v(i,j,k);
197 }
else if (k == bx_hi.z) {
198 old_drho_v(i,j,k+1) = old_drho_v(i,j,k);
201 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
202 old_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k);
205 const Array4<Real>& theta_extrap = extrap.array(mfi);
206 const Array4<const Real>& prim = S_stage_prim.const_array(mfi);
208 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
212 theta_extrap(i,j,k) = old_drho_theta(i,j,k);
214 theta_extrap(i,j,k) = old_drho_theta(i,j,k) + beta_d *
215 ( old_drho_theta(i,j,k) - lagged_arr(i,j,k) );
220 theta_extrap(i,j,k) *= (
one + RvOverRd*
qv);
225 #pragma omp parallel if (Gpu::notInLaunchRegion())
227 for ( MFIter mfi(S_stage_data[
IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
231 Box gbx = mfi.growntilebox(1);
232 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
233 const Array4<Real>& lagged_arr = lagged_delta_rt.array(mfi);
234 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
235 lagged_arr(i,j,k) = old_drho_theta(i,j,k);
244 #pragma omp parallel if (Gpu::notInLaunchRegion())
246 for ( MFIter mfi(S_stage_data[
IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
248 Box bx = mfi.validbox();
249 Box tbx = mfi.nodaltilebox(0);
250 Box tby = mfi.nodaltilebox(1);
252 const Array4<Real const>& xmom_src_arr = xmom_src.const_array(mfi);
253 const Array4<Real const>& ymom_src_arr = ymom_src.const_array(mfi);
255 const Array4<const Real> & stage_xmom = S_stage_data[
IntVars::xmom].const_array(mfi);
256 const Array4<const Real> & stage_ymom = S_stage_data[
IntVars::ymom].const_array(mfi);
257 const Array4<const Real> & qt_arr =
qt.const_array(mfi);
259 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
260 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
262 const Array4<const Real>& slow_rhs_rho_u = S_slow_rhs[
IntVars::xmom].const_array(mfi);
263 const Array4<const Real>& slow_rhs_rho_v = S_slow_rhs[
IntVars::ymom].const_array(mfi);
265 const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
266 const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
268 const Array4<Real>& cur_xmom = S_data[
IntVars::xmom].array(mfi);
269 const Array4<Real>& cur_ymom = S_data[
IntVars::ymom].array(mfi);
272 const Array4<Real>& avg_xmom_arr = avg_xmom.array(mfi);
273 const Array4<Real>& avg_ymom_arr = avg_ymom.array(mfi);
275 const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
277 const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
279 const Array4<Real>& theta_extrap = extrap.array(mfi);
282 const Array4<const Real>& mf_ux = mapfac[
MapFacType::u_x]->const_array(mfi);
283 const Array4<const Real>& mf_uy = mapfac[
MapFacType::u_y]->const_array(mfi);
284 const Array4<const Real>& mf_vx = mapfac[
MapFacType::v_x]->const_array(mfi);
285 const Array4<const Real>& mf_vy = mapfac[
MapFacType::v_y]->const_array(mfi);
297 BL_PROFILE(
"substep_xymom_T");
299 const auto& bx_lo = lbound(bx);
300 const auto& bx_hi = ubound(bx);
303 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
308 Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi;
309 Real gp_zeta_on_iface = (k == 0) ?
310 myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
311 - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) :
312 fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
313 - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) );
314 Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ?
Real(0.) :
315 gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface;
319 Real q = (l_use_moisture) ?
myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) :
zero;
321 Real pi_c =
myhalf * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i ,j,k,0));
324 new_drho_u(i, j, k) = old_drho_u(i,j,k) + dtau * fast_rhs_rho_u
325 + dtau * slow_rhs_rho_u(i,j,k)
326 + dtau * xmom_src_arr(i,j,k);
327 if (k == bx_lo.z && k != domlo.z) {
328 new_drho_u(i,j,k-1) = new_drho_u(i,j,k);
329 }
else if (k == bx_hi.z) {
330 new_drho_u(i,j,k+1) = new_drho_u(i,j,k);
337 avg_xmom_arr(i,j,k) += facinv * new_drho_u(i,j,k) * met_h_zeta / mf_uy(i,j,0);
339 cur_xmom(i,j,k) = stage_xmom(i,j,k) + new_drho_u(i,j,k);
341 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
346 Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi;
347 Real gp_zeta_on_jface = (k == 0) ?
348 myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
349 - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) :
350 fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
351 - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) );
352 Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ?
Real(0.) :
353 gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface;
357 Real q = (l_use_moisture) ?
myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) :
zero;
359 Real pi_c =
myhalf * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j ,k,0));
362 new_drho_v(i, j, k) = old_drho_v(i,j,k) + dtau * fast_rhs_rho_v
363 + dtau * slow_rhs_rho_v(i,j,k)
364 + dtau * ymom_src_arr(i,j,k);
366 if (k == bx_lo.z && k != domlo.z) {
367 new_drho_v(i,j,k-1) = new_drho_v(i,j,k);
368 }
else if (k == bx_hi.z) {
369 new_drho_v(i,j,k+1) = new_drho_v(i,j,k);
376 avg_ymom_arr(i,j,k) += facinv * new_drho_v(i,j,k) * met_h_zeta / mf_vx(i,j,0);
378 cur_ymom(i,j,k) = stage_ymom(i,j,k) + new_drho_v(i,j,k);
386 #pragma omp parallel if (Gpu::notInLaunchRegion())
389 std::array<FArrayBox,AMREX_SPACEDIM> flux;
392 Box bx = mfi.tilebox();
393 Box tbz = surroundingNodes(bx,2);
395 Box vbx = mfi.validbox();
396 const auto& vbx_hi = ubound(vbx);
398 const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
399 const Array4<Real const>& cc_src_arr = cc_src.const_array(mfi);
401 const Array4<const Real> & stage_zmom = S_stage_data[
IntVars::zmom].const_array(mfi);
402 const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
404 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
405 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
406 const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
407 const Array4<Real>& old_drho = Delta_rho.array(mfi);
408 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
410 const Array4<const Real>& slow_rhs_cons = S_slow_rhs[
IntVars::cons].const_array(mfi);
411 const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[
IntVars::zmom].const_array(mfi);
413 const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
414 const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
416 const Array4<Real>& cur_cons = S_data[
IntVars::cons].array(mfi);
417 const Array4<Real>& cur_zmom = S_data[
IntVars::zmom].array(mfi);
420 const Array4<Real>& avg_zmom_arr = avg_zmom.array(mfi);
422 const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
423 const Array4<const Real>& detJ = detJ_cc->const_array(mfi);
425 const Array4< Real>& omega_arr = Omega.array(mfi);
428 const Array4<const Real>& mf_mx = mapfac[
MapFacType::m_x]->const_array(mfi);
429 const Array4<const Real>& mf_my = mapfac[
MapFacType::m_y]->const_array(mfi);
430 const Array4<const Real>& mf_ux = mapfac[
MapFacType::u_x]->const_array(mfi);
431 const Array4<const Real>& mf_uy = mapfac[
MapFacType::u_y]->const_array(mfi);
432 const Array4<const Real>& mf_vx = mapfac[
MapFacType::v_x]->const_array(mfi);
433 const Array4<const Real>& mf_vy = mapfac[
MapFacType::v_y]->const_array(mfi);
441 FArrayBox temp_rhs_fab;
445 RHS_fab.resize (tbz,1,The_Async_Arena());
446 soln_fab.resize (tbz,1,The_Async_Arena());
447 temp_rhs_fab.resize(tbz,2,The_Async_Arena());
449 auto const& RHS_a = RHS_fab.array();
450 auto const& soln_a = soln_fab.array();
451 auto const& temp_rhs_arr = temp_rhs_fab.array();
453 auto const& coeffA_a = coeff_A_mf.array(mfi);
454 auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
455 auto const& coeffC_a = coeff_C_mf.array(mfi);
456 auto const& coeffP_a = coeff_P_mf.array(mfi);
457 auto const& coeffQ_a = coeff_Q_mf.array(mfi);
462 for (
int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
463 flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
464 flux[dir].setVal<RunOn::Device>(0);
466 const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
467 flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
471 BL_PROFILE(
"fast_T_making_rho_rhs");
472 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
474 ( z_nd(i+1,j ,k+1) + z_nd(i+1,j+1,k+1)
475 -z_nd(i+1,j ,k ) - z_nd(i+1,j+1,k ) );
478 ( z_nd(i ,j ,k+1) + z_nd(i ,j+1,k+1)
479 -z_nd(i ,j ,k ) - z_nd(i ,j+1,k ) );
482 ( z_nd(i ,j+1,k+1) + z_nd(i+1,j+1,k+1)
483 -z_nd(i ,j+1,k ) - z_nd(i+1,j+1,k ) );
486 ( z_nd(i ,j ,k+1) + z_nd(i+1,j ,k+1)
487 -z_nd(i ,j ,k ) - z_nd(i+1,j ,k ) );
489 Real xflux_lo = new_drho_u(i ,j,k)*h_zeta_cc_xface_lo / mf_uy(i ,j,0);
490 Real xflux_hi = new_drho_u(i+1,j,k)*h_zeta_cc_xface_hi / mf_uy(i+1,j,0);
491 Real yflux_lo = new_drho_v(i,j ,k)*h_zeta_cc_yface_lo / mf_vx(i,j ,0);
492 Real yflux_hi = new_drho_v(i,j+1,k)*h_zeta_cc_yface_hi / mf_vx(i,j+1,0);
494 Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
497 temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq +
498 ( yflux_hi - yflux_lo ) * dyi * mfsq;
499 temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
500 xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+
501 ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
502 yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) *
myhalf;
505 (flx_arr[0])(i,j,k,0) = xflux_lo;
506 (flx_arr[0])(i,j,k,1) = (flx_arr[0])(i ,j,k,0) *
myhalf * (prim(i,j,k,0) + prim(i-1,j,k,0));
508 (flx_arr[1])(i,j,k,0) = yflux_lo;
509 (flx_arr[1])(i,j,k,1) = (flx_arr[1])(i,j ,k,0) *
myhalf * (prim(i,j,k,0) + prim(i,j-1,k,0));
512 (flx_arr[0])(i+1,j,k,0) = xflux_hi;
513 (flx_arr[0])(i+1,j,k,1) = (flx_arr[0])(i+1,j,k,0) *
myhalf * (prim(i,j,k,0) + prim(i+1,j,k,0));
516 (flx_arr[1])(i,j+1,k,0) = yflux_hi;
517 (flx_arr[1])(i,j+1,k,1) = (flx_arr[1])(i,j+1,k,0) *
myhalf * (prim(i,j,k,0) + prim(i,j+1,k,0));
525 Box gbxo = mfi.nodaltilebox(2);
528 if (gbxo.smallEnd(2) == domlo.z) {
529 Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2));
530 gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1);
531 ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
532 omega_arr(i,j,k) =
zero;
535 if (gbxo.bigEnd(2) == domhi.z+1) {
536 Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
537 gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
538 ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
539 omega_arr(i,j,k) = old_drho_w(i,j,k);
542 ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
543 omega_arr(i,j,k) =
OmegaFromW(i,j,k,old_drho_w(i,j,k),
544 old_drho_u,old_drho_v,
545 mf_ux,mf_vy,z_nd,
dxInv);
550 Box bx_shrunk_in_k = bx;
551 int klo = tbz.smallEnd(2);
552 int khi = tbz.bigEnd(2);
553 bx_shrunk_in_k.setSmall(2,
klo+1);
554 bx_shrunk_in_k.setBig(2,
khi-1);
562 BL_PROFILE(
"fast_loop_on_shrunk_t");
564 ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
566 Real coeff_P = coeffP_a(i,j,k);
567 Real coeff_Q = coeffQ_a(i,j,k);
574 Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k )
575 -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1);
583 Real Omega_kp1 = omega_arr(i,j,k+1);
584 Real Omega_k = omega_arr(i,j,k );
585 Real Omega_km1 = omega_arr(i,j,k-1);
587 Real detJdiff = (detJ(i,j,k) - detJ(i,j,k-1)) / (detJ(i,j,k)*detJ(i,j,k-1));
590 R1_tmp += halfg * ( beta_1 * dzi * (Omega_kp1/detJ(i,j,k) + detJdiff*Omega_k - Omega_km1/detJ(i,j,k-1))
591 + temp_rhs_arr(i,j,k,
Rho_comp)/detJ(i,j,k) + temp_rhs_arr(i,j,k-1,
Rho_comp)/detJ(i,j,k-1) );
594 R1_tmp += -( coeff_P/detJ(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*theta_t_hi - Omega_k*theta_t_mid) + temp_rhs_arr(i,j,k ,
RhoTheta_comp) )
595 + coeff_Q/detJ(i,j,k-1) * ( beta_1 * dzi * (Omega_k*theta_t_mid - Omega_km1*theta_t_lo) + temp_rhs_arr(i,j,k-1,
RhoTheta_comp) ) );
598 RHS_a(i,j,k) = old_drho_w(i,j,k) + dtau * (slow_rhs_rho_w(i,j,k) + zmom_src_arr(i,j,k) + R0_tmp + dtau*beta_2*R1_tmp);
602 new_drho_u,new_drho_v,
603 mf_ux,mf_vy,z_nd,
dxInv);
610 auto const lo = lbound(bx);
611 auto const hi = ubound(bx);
614 BL_PROFILE(
"substep_b2d_loop_t");
616 ParallelFor(b2d, [=] AMREX_GPU_DEVICE (
int i,
int j,
int)
619 RHS_a(i,j,
lo.z ) = dtau * (slow_rhs_rho_w(i,j,
lo.z ) + zmom_src_arr(i,j,
lo.z ));
620 RHS_a(i,j,
hi.z+1) = dtau * (slow_rhs_rho_w(i,j,
hi.z+1) + zmom_src_arr(i,j,
hi.z+1));
623 soln_a(i,j,
lo.z) = RHS_a(i,j,
lo.z) * inv_coeffB_a(i,j,
lo.z);
626 for (
int k =
lo.z+1; k <=
hi.z+1; k++) {
627 soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
630 cur_zmom(i,j,
lo.z ) = stage_zmom(i,j,
lo.z ) + soln_a(i,j,
lo.z );
631 cur_zmom(i,j,
hi.z+1) = stage_zmom(i,j,
hi.z+1) + soln_a(i,j,
hi.z+1);
634 for (
int k =
hi.z; k >=
lo.z; k--) {
635 soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1);
640 for (
int j =
lo.y; j <=
hi.y; ++j) {
642 for (
int i =
lo.x; i <=
hi.x; ++i)
644 RHS_a(i,j,
lo.z ) = dtau * (slow_rhs_rho_w(i,j,
lo.z ) + zmom_src_arr(i,j,
lo.z ));
645 RHS_a(i,j,
hi.z+1) = dtau * (slow_rhs_rho_w(i,j,
hi.z+1) + zmom_src_arr(i,j,
hi.z+1));
648 soln_a(i,j,
lo.z) = RHS_a(i,j,
lo.z) * inv_coeffB_a(i,j,
lo.z);
655 for (
int k =
lo.z+1; k <=
hi.z+1; ++k) {
656 for (
int j =
lo.y; j <=
hi.y; ++j) {
658 for (
int i =
lo.x; i <=
hi.x; ++i) {
659 soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
664 for (
int j =
lo.y; j <=
hi.y; ++j) {
666 for (
int i =
lo.x; i <=
hi.x; ++i) {
667 cur_zmom(i,j,
lo.z ) = stage_zmom(i,j,
lo.z ) + soln_a(i,j,
lo.z );
668 cur_zmom(i,j,
hi.z+1) = stage_zmom(i,j,
hi.z+1) + soln_a(i,j,
hi.z+1);
673 for (
int k =
hi.z; k >=
lo.z; --k) {
674 for (
int j =
lo.y; j <=
hi.y; ++j) {
676 for (
int i =
lo.x; i <=
hi.x; ++i) {
677 soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
684 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
686 cur_zmom(i,j,k) = stage_zmom(i,j,k);
689 if (
lo.z == domlo.z) {
690 tbz.setSmall(2,domlo.z+1);
692 if (
hi.z == domhi.z) {
693 tbz.setBig(2,domhi.z);
695 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
698 new_drho_u,new_drho_v,
699 mf_ux,mf_vy,z_nd,
dxInv);
701 cur_zmom(i,j,k) += wpp;
703 if (l_rayleigh_impl_for_w) {
704 Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
705 cur_zmom(i,j,k) /= (
one + damping_coeff);
713 BL_PROFILE(
"fast_rho_final_update");
714 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
716 Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k);
717 Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1);
721 avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
723 (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
727 avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
729 (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
730 (flx_arr[2])(i,j,k+1,1) = (flx_arr[2])(i,j,k+1,0) *
myhalf * (prim(i,j,k) + prim(i,j,k+1));
734 Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k);
736 cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho);
738 Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) +
myhalf *
739 ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) -
740 zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k);
742 cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta);
745 (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) *
myhalf * (prim(i,j,k) + prim(i,j,k-1));
756 int strt_comp_reflux = 0;
758 int num_comp_reflux = 1;
759 if (level < finest_level) {
760 fr_as_crse->CrseAdd(mfi,
761 {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
762 dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
765 fr_as_fine->FineAdd(mfi,
766 {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
767 dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
773 Gpu::streamSynchronize();
constexpr amrex::Real R_v
Definition: ERF_Constants.H:35
constexpr amrex::Real R_d
Definition: ERF_Constants.H:34
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:54
@ v_x
Definition: ERF_DataStruct.H:29
@ u_y
Definition: ERF_DataStruct.H:30
@ v_y
Definition: ERF_DataStruct.H:30
@ m_y
Definition: ERF_DataStruct.H:30
@ u_x
Definition: ERF_DataStruct.H:29
@ m_x
Definition: ERF_DataStruct.H:29
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:61
#define Rho_comp
Definition: ERF_IndexDefines.H:39
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:58
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:35
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real OmegaFromW(int &i, int &j, int &k, amrex::Real w, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:791
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:292
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:269
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:339
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:385
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real WFromOmega(int &i, int &j, int &k, amrex::Real omega, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:856
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
@ gpy
Definition: ERF_IndexDefines.H:225
@ gpx
Definition: ERF_IndexDefines.H:224
@ ymom
Definition: ERF_IndexDefines.H:234
@ cons
Definition: ERF_IndexDefines.H:232
@ zmom
Definition: ERF_IndexDefines.H:235
@ xmom
Definition: ERF_IndexDefines.H:233
@ qt
Definition: ERF_Kessler.H:30
@ qv
Definition: ERF_Kessler.H:31
@ q
Definition: ERF_WSM6.H:273