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_uy = mapfac[
MapFacType::u_y]->const_array(mfi);
274 const Array4<const Real>& mf_vx = mapfac[
MapFacType::v_x]->const_array(mfi);
275 const Array4<const Real>& mf_vy = mapfac[
MapFacType::v_y]->const_array(mfi);
287 BL_PROFILE(
"substep_xymom_T");
289 const auto& bx_lo = lbound(bx);
290 const auto& bx_hi = ubound(bx);
293 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
298 Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi;
299 Real gp_zeta_on_iface = (k == 0) ?
300 myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
301 - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) :
302 fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1)
303 - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) );
304 Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ?
Real(0.) :
305 gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface;
309 Real q = (l_use_moisture) ?
myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) :
zero;
311 Real pi_c =
myhalf * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i ,j,k,0));
314 new_drho_u(i, j, k) = old_drho_u(i,j,k) + dtau * fast_rhs_rho_u
315 + dtau * slow_rhs_rho_u(i,j,k)
316 + dtau * xmom_src_arr(i,j,k);
317 if (k == bx_lo.z && k != domlo.z) {
318 new_drho_u(i,j,k-1) = new_drho_u(i,j,k);
319 }
else if (k == bx_hi.z) {
320 new_drho_u(i,j,k+1) = new_drho_u(i,j,k);
327 avg_xmom_arr(i,j,k) += facinv * new_drho_u(i,j,k) * met_h_zeta / mf_uy(i,j,0);
329 cur_xmom(i,j,k) = stage_xmom(i,j,k) + new_drho_u(i,j,k);
331 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
336 Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi;
337 Real gp_zeta_on_jface = (k == 0) ?
338 myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
339 - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) :
340 fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1)
341 - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) );
342 Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ?
Real(0.) :
343 gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface;
347 Real q = (l_use_moisture) ?
myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) :
zero;
349 Real pi_c =
myhalf * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j ,k,0));
352 new_drho_v(i, j, k) = old_drho_v(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);
356 if (k == bx_lo.z && k != domlo.z) {
357 new_drho_v(i,j,k-1) = new_drho_v(i,j,k);
358 }
else if (k == bx_hi.z) {
359 new_drho_v(i,j,k+1) = new_drho_v(i,j,k);
366 avg_ymom_arr(i,j,k) += facinv * new_drho_v(i,j,k) * met_h_zeta / mf_vx(i,j,0);
368 cur_ymom(i,j,k) = stage_ymom(i,j,k) + new_drho_v(i,j,k);
376 #pragma omp parallel if (Gpu::notInLaunchRegion())
379 std::array<FArrayBox,AMREX_SPACEDIM> flux;
382 Box bx = mfi.tilebox();
383 Box tbz = surroundingNodes(bx,2);
385 Box vbx = mfi.validbox();
386 const auto& vbx_hi = ubound(vbx);
388 const Array4<Real const>& zmom_src_arr = zmom_src.const_array(mfi);
389 const Array4<Real const>& cc_src_arr = cc_src.const_array(mfi);
391 const Array4<const Real> & stage_zmom = S_stage_data[
IntVars::zmom].const_array(mfi);
392 const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
394 const Array4<Real>& old_drho_u = Delta_rho_u.array(mfi);
395 const Array4<Real>& old_drho_v = Delta_rho_v.array(mfi);
396 const Array4<Real>& old_drho_w = Delta_rho_w.array(mfi);
397 const Array4<Real>& old_drho = Delta_rho.array(mfi);
398 const Array4<Real>& old_drho_theta = Delta_rho_theta.array(mfi);
400 const Array4<const Real>& slow_rhs_cons = S_slow_rhs[
IntVars::cons].const_array(mfi);
401 const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[
IntVars::zmom].const_array(mfi);
403 const Array4<Real>& new_drho_u = New_rho_u.array(mfi);
404 const Array4<Real>& new_drho_v = New_rho_v.array(mfi);
406 const Array4<Real>& cur_cons = S_data[
IntVars::cons].array(mfi);
407 const Array4<Real>& cur_zmom = S_data[
IntVars::zmom].array(mfi);
410 const Array4<Real>& avg_zmom_arr = avg_zmom.array(mfi);
412 const Array4<const Real>& z_nd = z_phys_nd->const_array(mfi);
413 const Array4<const Real>& detJ = detJ_cc->const_array(mfi);
415 const Array4< Real>& omega_arr = Omega.array(mfi);
418 const Array4<const Real>& mf_mx = mapfac[
MapFacType::m_x]->const_array(mfi);
419 const Array4<const Real>& mf_my = mapfac[
MapFacType::m_y]->const_array(mfi);
420 const Array4<const Real>& mf_ux = mapfac[
MapFacType::u_x]->const_array(mfi);
421 const Array4<const Real>& mf_uy = mapfac[
MapFacType::u_y]->const_array(mfi);
422 const Array4<const Real>& mf_vx = mapfac[
MapFacType::v_x]->const_array(mfi);
423 const Array4<const Real>& mf_vy = mapfac[
MapFacType::v_y]->const_array(mfi);
431 FArrayBox temp_rhs_fab;
435 RHS_fab.resize (tbz,1,The_Async_Arena());
436 soln_fab.resize (tbz,1,The_Async_Arena());
437 temp_rhs_fab.resize(tbz,2,The_Async_Arena());
439 auto const& RHS_a = RHS_fab.array();
440 auto const& soln_a = soln_fab.array();
441 auto const& temp_rhs_arr = temp_rhs_fab.array();
443 auto const& coeffA_a = coeff_A_mf.array(mfi);
444 auto const& inv_coeffB_a = inv_coeff_B_mf.array(mfi);
445 auto const& coeffC_a = coeff_C_mf.array(mfi);
446 auto const& coeffP_a = coeff_P_mf.array(mfi);
447 auto const& coeffQ_a = coeff_Q_mf.array(mfi);
452 for (
int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
453 flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena());
454 flux[dir].setVal<RunOn::Device>(0);
456 const GpuArray<const Array4<Real>, AMREX_SPACEDIM>
457 flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}};
461 BL_PROFILE(
"fast_T_making_rho_rhs");
462 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
464 ( z_nd(i+1,j ,k+1) + z_nd(i+1,j+1,k+1)
465 -z_nd(i+1,j ,k ) - z_nd(i+1,j+1,k ) );
468 ( z_nd(i ,j ,k+1) + z_nd(i ,j+1,k+1)
469 -z_nd(i ,j ,k ) - z_nd(i ,j+1,k ) );
472 ( z_nd(i ,j+1,k+1) + z_nd(i+1,j+1,k+1)
473 -z_nd(i ,j+1,k ) - z_nd(i+1,j+1,k ) );
476 ( z_nd(i ,j ,k+1) + z_nd(i+1,j ,k+1)
477 -z_nd(i ,j ,k ) - z_nd(i+1,j ,k ) );
479 Real xflux_lo = new_drho_u(i ,j,k)*h_zeta_cc_xface_lo / mf_uy(i ,j,0);
480 Real xflux_hi = new_drho_u(i+1,j,k)*h_zeta_cc_xface_hi / mf_uy(i+1,j,0);
481 Real yflux_lo = new_drho_v(i,j ,k)*h_zeta_cc_yface_lo / mf_vx(i,j ,0);
482 Real yflux_hi = new_drho_v(i,j+1,k)*h_zeta_cc_yface_hi / mf_vx(i,j+1,0);
484 Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
487 temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq +
488 ( yflux_hi - yflux_lo ) * dyi * mfsq;
489 temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) -
490 xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+
491 ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) -
492 yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) *
myhalf;
495 (flx_arr[0])(i,j,k,0) = xflux_lo;
496 (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));
498 (flx_arr[1])(i,j,k,0) = yflux_lo;
499 (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));
502 (flx_arr[0])(i+1,j,k,0) = xflux_hi;
503 (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));
506 (flx_arr[1])(i,j+1,k,0) = yflux_hi;
507 (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));
515 Box gbxo = mfi.nodaltilebox(2);
518 if (gbxo.smallEnd(2) == domlo.z) {
519 Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2));
520 gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1);
521 ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
522 omega_arr(i,j,k) =
zero;
525 if (gbxo.bigEnd(2) == domhi.z+1) {
526 Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2));
527 gbxo_mid.setBig(2,gbxo.bigEnd(2)-1);
528 ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
529 omega_arr(i,j,k) = old_drho_w(i,j,k);
532 ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
533 omega_arr(i,j,k) =
OmegaFromW(i,j,k,old_drho_w(i,j,k),
534 old_drho_u,old_drho_v,
535 mf_ux,mf_vy,z_nd,
dxInv);
540 Box bx_shrunk_in_k = bx;
541 int klo = tbz.smallEnd(2);
542 int khi = tbz.bigEnd(2);
543 bx_shrunk_in_k.setSmall(2,klo+1);
544 bx_shrunk_in_k.setBig(2,
khi-1);
552 BL_PROFILE(
"fast_loop_on_shrunk_t");
554 ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
556 Real coeff_P = coeffP_a(i,j,k);
557 Real coeff_Q = coeffQ_a(i,j,k);
564 Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k )
565 -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1);
573 Real Omega_kp1 = omega_arr(i,j,k+1);
574 Real Omega_k = omega_arr(i,j,k );
575 Real Omega_km1 = omega_arr(i,j,k-1);
577 Real detJdiff = (detJ(i,j,k) - detJ(i,j,k-1)) / (detJ(i,j,k)*detJ(i,j,k-1));
580 R1_tmp += halfg * ( beta_1 * dzi * (Omega_kp1/detJ(i,j,k) + detJdiff*Omega_k - Omega_km1/detJ(i,j,k-1))
581 + 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) );
584 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) )
585 + 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) ) );
588 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);
592 new_drho_u,new_drho_v,
593 mf_ux,mf_vy,z_nd,
dxInv);
600 auto const lo = lbound(bx);
601 auto const hi = ubound(bx);
604 BL_PROFILE(
"substep_b2d_loop_t");
606 ParallelFor(b2d, [=] AMREX_GPU_DEVICE (
int i,
int j,
int)
609 RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z ) + zmom_src_arr(i,j,lo.z ));
610 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));
613 soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
616 for (
int k = lo.z+1; k <= hi.z+1; k++) {
617 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);
620 cur_zmom(i,j,lo.z ) = stage_zmom(i,j,lo.z ) + soln_a(i,j,lo.z );
621 cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
624 for (
int k = hi.z; k >= lo.z; k--) {
625 soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1);
630 for (
int j = lo.y; j <= hi.y; ++j) {
632 for (
int i = lo.x; i <= hi.x; ++i)
634 RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z ) + zmom_src_arr(i,j,lo.z ));
635 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));
638 soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
645 for (
int k = lo.z+1; k <= hi.z+1; ++k) {
646 for (
int j = lo.y; j <= hi.y; ++j) {
648 for (
int i = lo.x; i <= hi.x; ++i) {
649 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);
654 for (
int j = lo.y; j <= hi.y; ++j) {
656 for (
int i = lo.x; i <= hi.x; ++i) {
657 cur_zmom(i,j,lo.z ) = stage_zmom(i,j,lo.z ) + soln_a(i,j,lo.z );
658 cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1);
663 for (
int k = hi.z; k >= lo.z; --k) {
664 for (
int j = lo.y; j <= hi.y; ++j) {
666 for (
int i = lo.x; i <= hi.x; ++i) {
667 soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1);
674 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
676 cur_zmom(i,j,k) = stage_zmom(i,j,k);
679 if (lo.z == domlo.z) {
680 tbz.setSmall(2,domlo.z+1);
682 if (hi.z == domhi.z) {
683 tbz.setBig(2,domhi.z);
685 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
688 new_drho_u,new_drho_v,
689 mf_ux,mf_vy,z_nd,
dxInv);
691 cur_zmom(i,j,k) += wpp;
693 if (l_rayleigh_impl_for_w) {
694 Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k];
695 cur_zmom(i,j,k) /= (
one + damping_coeff);
703 BL_PROFILE(
"fast_rho_final_update");
704 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
706 Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k);
707 Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1);
711 avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
713 (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0));
717 avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
719 (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0));
720 (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));
724 Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k);
726 cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho);
728 Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) +
myhalf *
729 ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) -
730 zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k);
732 cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta);
735 (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) *
myhalf * (prim(i,j,k) + prim(i,j,k-1));
746 int strt_comp_reflux = 0;
748 int num_comp_reflux = 1;
749 if (level < finest_level) {
750 fr_as_crse->CrseAdd(mfi,
751 {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
752 dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
755 fr_as_fine->FineAdd(mfi,
756 {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}},
757 dx, dtau, strt_comp_reflux, strt_comp_reflux, num_comp_reflux, RunOn::Device);
763 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:27
@ u_y
Definition: ERF_DataStruct.H:28
@ v_y
Definition: ERF_DataStruct.H:28
@ m_y
Definition: ERF_DataStruct.H:28
@ u_x
Definition: ERF_DataStruct.H:27
@ m_x
Definition: ERF_DataStruct.H:27
#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(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);})
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:187
@ gpx
Definition: ERF_IndexDefines.H:186
@ ymom
Definition: ERF_IndexDefines.H:196
@ cons
Definition: ERF_IndexDefines.H:194
@ zmom
Definition: ERF_IndexDefines.H:197
@ xmom
Definition: ERF_IndexDefines.H:195
@ qt
Definition: ERF_Kessler.H:29
@ qv
Definition: ERF_Kessler.H:30
@ q
Definition: ERF_WSM6.H:184