From d9500d69960042018ff7bb767a1f1ebdbc0ea42c Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Mon, 20 Jul 2026 14:39:35 -0700 Subject: [PATCH 1/3] rebase without macro coupling. --- Source/Diffusion/ERF_Diffusion.H | 15 +- Source/Diffusion/ERF_ImplicitDiff_N.cpp | 25 +- Source/Diffusion/ERF_ImplicitDiff_S.cpp | 24 +- Source/Diffusion/ERF_ImplicitDiff_T.cpp | 24 +- Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp | 122 +++---- Source/TimeIntegration/ERF_ImplicitPost.H | 58 +-- Source/TimeIntegration/ERF_ImplicitPre.H | 50 ++- Source/TimeIntegration/ERF_MakeFastCoeffs.cpp | 42 ++- Source/TimeIntegration/ERF_SlowRhsPost.cpp | 29 +- Source/TimeIntegration/ERF_SlowRhsPre.cpp | 83 +++++ Source/TimeIntegration/ERF_Substep_T.cpp | 335 ++++++++++-------- Source/TimeIntegration/ERF_TI_slow_rhs_pre.H | 16 +- Source/Utils/ERF_MoistUtils.H | 4 +- 13 files changed, 481 insertions(+), 346 deletions(-) diff --git a/Source/Diffusion/ERF_Diffusion.H b/Source/Diffusion/ERF_Diffusion.H index 9c2cfcd074..3690117b52 100644 --- a/Source/Diffusion/ERF_Diffusion.H +++ b/Source/Diffusion/ERF_Diffusion.H @@ -461,8 +461,9 @@ void ComputeStrain_T (amrex::Box bxcc, amrex::Box tbxxy, void ImplicitDiffForStateLU_N (const amrex::Box& bx, const amrex::Box& domain, const int level, - const int n, - const double dt, + const int tmp_index, + const int qty_index, + const amrex::Real dt, const amrex::GpuArray& bc_neumann_vals, const amrex::Array4< amrex::Real>& cell_data, const amrex::GpuArray& cellSizeInv, @@ -477,8 +478,9 @@ void ImplicitDiffForStateLU_N (const amrex::Box& bx, void ImplicitDiffForStateLU_S (const amrex::Box& bx, const amrex::Box& domain, const int level, - const int n, - const double dt, + const int tmp_index, + const int qty_index, + const amrex::Real dt, const amrex::GpuArray& bc_neumann_vals, const amrex::Array4< amrex::Real>& cell_data, const amrex::Gpu::DeviceVector& stretched_dz_d, @@ -493,8 +495,9 @@ void ImplicitDiffForStateLU_S (const amrex::Box& bx, void ImplicitDiffForStateLU_T (const amrex::Box& bx, const amrex::Box& domain, const int level, - const int n, - const double dt, + const int tmp_index, + const int qty_index, + const amrex::Real dt, const amrex::GpuArray& bc_neumann_vals, const amrex::Array4< amrex::Real>& cell_data, const amrex::Array4& z_nd, diff --git a/Source/Diffusion/ERF_ImplicitDiff_N.cpp b/Source/Diffusion/ERF_ImplicitDiff_N.cpp index d0a91e3372..83be8ab7c2 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_N.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_N.cpp @@ -30,8 +30,9 @@ void ImplicitDiffForStateLU_N (const Box& bx, const Box& domain, const int level, - const int n, - const double dt_d, + const int tmp_index, + const int qty_index, + const Real dt, const GpuArray& bc_neumann_vals, const Array4< Real>& cell_data, const GpuArray& cellSizeInv, @@ -45,11 +46,8 @@ ImplicitDiffForStateLU_N (const Box& bx, { BL_PROFILE_VAR("ImplicitDiffForState_N()",ImplicitDiffForState_N); - Real dt = static_cast(dt_d); - // setup quantities for getRhoAlpha() #include "ERF_SetupVertDiff.H" - const int qty_index = n; const int prim_index = qty_index - 1; const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; @@ -112,7 +110,7 @@ ImplicitDiffForStateLU_N (const Box& bx, b_tmp = cell_data(i,j,klo,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one; - RHS_a(i,j,klo) = cell_data(i,j,klo,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,klo) = cell_data(i,j,klo,tmp_index); // NOTE: this is rho*phi; solution is phi if (use_SurfLayer && scalar_zflux) { RHS_a(i,j,klo) += Fact * scalar_zflux(i,j,klo); // NOTE: scalar_zflux = -K*d_z(\phi) } else if (neumann_on_zlo) { @@ -121,8 +119,8 @@ ImplicitDiffForStateLU_N (const Box& bx, // Add countergradient correction to RHS at bottom boundary // Only upper face contributes (no flux below surface) - if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) { - const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; + if (use_mrf_countergradient && (qty_index == RhoTheta_comp || qty_index == RhoQ1_comp)) { + const int gam_comp = (qty_index == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; const Real gam_hi = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo+1, gam_comp)); RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi; } @@ -143,11 +141,11 @@ ImplicitDiffForStateLU_N (const Box& bx, b_tmp = cell_data(i,j,k,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1)); - RHS_a(i,j,k) = cell_data(i,j,k,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,k) = cell_data(i,j,k,tmp_index); // NOTE: this is rho*phi; solution is phi // Add countergradient correction to RHS in interior - if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) { - const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; + if (use_mrf_countergradient && (qty_index == RhoTheta_comp || qty_index == RhoQ1_comp)) { + const int gam_comp = (qty_index == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; const Real gam_k = mu_turb(i, j, k, gam_comp); const Real gam_km1 = mu_turb(i, j, k-1, gam_comp); const Real gam_kp1 = mu_turb(i, j, k+1, gam_comp); @@ -174,7 +172,7 @@ ImplicitDiffForStateLU_N (const Box& bx, b_tmp = cell_data(i,j,khi,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1)); - RHS_a(i,j,khi) = cell_data(i,j,khi,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,khi) = cell_data(i,j,khi,tmp_index); // NOTE: this is rho*phi; solution is phi if (neumann_on_zhi) { RHS_a(i,j,khi) -= -Fact * rhoAlpha_hi * bc_neumann_vals[5]; // NOTE: N_val = d_z(\phi) } @@ -192,7 +190,7 @@ ImplicitDiffForStateLU_N (const Box& bx, // Convert back to rho*theta //=================================================== for (int k(klo); k<=khi; ++k) { - cell_data(i,j,k,n) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k); + cell_data(i,j,k,tmp_index) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k); } #ifdef AMREX_USE_GPU @@ -496,4 +494,3 @@ INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(0) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(1) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(2) #undef INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU - \ No newline at end of file diff --git a/Source/Diffusion/ERF_ImplicitDiff_S.cpp b/Source/Diffusion/ERF_ImplicitDiff_S.cpp index c27512fecc..ad191080d1 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_S.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_S.cpp @@ -28,8 +28,9 @@ void ImplicitDiffForStateLU_S (const Box& bx, const Box& domain, const int level, - const int n, - const double dt_d, + const int tmp_index, + const int qty_index, + const Real dt, const GpuArray& bc_neumann_vals, const Array4< Real>& cell_data, const Gpu::DeviceVector& stretched_dz_d, @@ -43,11 +44,8 @@ ImplicitDiffForStateLU_S (const Box& bx, { BL_PROFILE_VAR("ImplicitDiffForState_S()",ImplicitDiffForState_S); - Real dt = static_cast(dt_d); - // setup quantities for getRhoAlpha() #include "ERF_SetupVertDiff.H" - const int qty_index = n; const int prim_index = qty_index - 1; const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; @@ -115,7 +113,7 @@ ImplicitDiffForStateLU_S (const Box& bx, b_tmp = cell_data(i,j,klo,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one; - RHS_a(i,j,klo) = cell_data(i,j,klo,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,klo) = cell_data(i,j,klo,tmp_index); // NOTE: this is rho*phi; solution is phi if (use_SurfLayer && scalar_zflux) { RHS_a(i,j,klo) += Fact * dz_inv * scalar_zflux(i,j,klo); // NOTE: scalar_zflux = -K*d_z(\phi) } else if (neumann_on_zlo) { @@ -123,8 +121,8 @@ ImplicitDiffForStateLU_S (const Box& bx, } // Add countergradient correction to RHS at bottom boundary - if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) { - const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; + if (use_mrf_countergradient && (qty_index == RhoTheta_comp || qty_index == RhoQ1_comp)) { + const int gam_comp = (qty_index == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; const Real gam_hi = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo+1, gam_comp)); RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi * dz_inv_hi; } @@ -149,11 +147,11 @@ ImplicitDiffForStateLU_S (const Box& bx, b_tmp = cell_data(i,j,k,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1)); - RHS_a(i,j,k) = cell_data(i,j,k,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,k) = cell_data(i,j,k,tmp_index); // NOTE: this is rho*phi; solution is phi // Add countergradient correction to RHS in interior - if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) { - const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; + if (use_mrf_countergradient && (qty_index == RhoTheta_comp || qty_index == RhoQ1_comp)) { + const int gam_comp = (qty_index == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; const Real gam_k = mu_turb(i, j, k, gam_comp); const Real gam_km1 = mu_turb(i, j, k-1, gam_comp); const Real gam_kp1 = mu_turb(i, j, k+1, gam_comp); @@ -184,7 +182,7 @@ ImplicitDiffForStateLU_S (const Box& bx, b_tmp = cell_data(i,j,khi,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1)); - RHS_a(i,j,khi) = cell_data(i,j,khi,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,khi) = cell_data(i,j,khi,tmp_index); // NOTE: this is rho*phi; solution is phi if (neumann_on_zhi) { RHS_a(i,j,khi) -= -Fact * dz_inv * rhoAlpha_hi * bc_neumann_vals[5]; // NOTE: N_val = d_z(\phi) } @@ -202,7 +200,7 @@ ImplicitDiffForStateLU_S (const Box& bx, // Convert back to rho*theta //=================================================== for (int k(klo); k<=khi; ++k) { - cell_data(i,j,k,n) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k); + cell_data(i,j,k,tmp_index) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k); } #ifdef AMREX_USE_GPU diff --git a/Source/Diffusion/ERF_ImplicitDiff_T.cpp b/Source/Diffusion/ERF_ImplicitDiff_T.cpp index 90c79cc8c7..6030d38039 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_T.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_T.cpp @@ -30,8 +30,9 @@ void ImplicitDiffForStateLU_T (const Box& bx, const Box& domain, const int level, - const int n, - const double dt_d, + const int tmp_index, + const int qty_index, + const Real dt, const GpuArray& bc_neumann_vals, const Array4< Real>& cell_data, const Array4& z_nd, @@ -47,11 +48,8 @@ ImplicitDiffForStateLU_T (const Box& bx, { BL_PROFILE_VAR("ImplicitDiffForState_T()",ImplicitDiffForState_T); - Real dt = static_cast(dt_d); - // setup quantities for getRhoAlpha() #include "ERF_SetupVertDiff.H" - const int qty_index = n; const int prim_index = qty_index - 1; const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; @@ -122,7 +120,7 @@ ImplicitDiffForStateLU_T (const Box& bx, b_tmp = detJ(i,j,klo) * cell_data(i,j,klo,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one; - RHS_a(i,j,klo) = detJ(i,j,klo) * cell_data(i,j,klo,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,klo) = detJ(i,j,klo) * cell_data(i,j,klo,tmp_index); // NOTE: this is rho*phi; solution is phi if (use_SurfLayer && scalar_zflux) { RHS_a(i,j,klo) += Fact * scalar_zflux(i,j,klo); // NOTE: scalar_zflux = -K*d_z(\phi) } else if (neumann_on_zlo) { @@ -130,8 +128,8 @@ ImplicitDiffForStateLU_T (const Box& bx, } // Add countergradient correction to RHS at bottom boundary - if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) { - const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; + if (use_mrf_countergradient && (qty_index == RhoTheta_comp || qty_index == RhoQ1_comp)) { + const int gam_comp = (qty_index == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; const Real gam_hi = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo+1, gam_comp)); RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi / met_h_zeta_hi; } @@ -155,11 +153,11 @@ ImplicitDiffForStateLU_T (const Box& bx, b_tmp = detJ(i,j,k) * cell_data(i,j,k,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1)); - RHS_a(i,j,k) = detJ(i,j,k) * cell_data(i,j,k,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,k) = detJ(i,j,k) * cell_data(i,j,k,tmp_index); // NOTE: this is rho*phi; solution is phi // Add countergradient correction to RHS in interior - if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) { - const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; + if (use_mrf_countergradient && (qty_index == RhoTheta_comp || qty_index == RhoQ1_comp)) { + const int gam_comp = (qty_index == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v; const Real gam_k = mu_turb(i, j, k, gam_comp); const Real gam_km1 = mu_turb(i, j, k-1, gam_comp); const Real gam_kp1 = mu_turb(i, j, k+1, gam_comp); @@ -189,7 +187,7 @@ ImplicitDiffForStateLU_T (const Box& bx, b_tmp = detJ(i,j,khi) * cell_data(i,j,khi,Rho_comp) - a_tmp - c_tmp; inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1)); - RHS_a(i,j,khi) = detJ(i,j,khi) * cell_data(i,j,khi,n); // NOTE: this is rho*phi; solution is phi + RHS_a(i,j,khi) = detJ(i,j,khi) * cell_data(i,j,khi,tmp_index); // NOTE: this is rho*phi; solution is phi if (neumann_on_zhi) { RHS_a(i,j,khi) -= -Fact * rhoAlpha_hi * bc_neumann_vals[5]; // NOTE: N_val = d_z(\phi) } @@ -207,7 +205,7 @@ ImplicitDiffForStateLU_T (const Box& bx, // Convert back to rho*theta //=================================================== for (int k(klo); k<=khi; ++k) { - cell_data(i,j,k,n) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k); + cell_data(i,j,k,tmp_index) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k); } #ifdef AMREX_USE_GPU diff --git a/Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp b/Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp index eac91fcd0e..58789b9bcf 100644 --- a/Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp +++ b/Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp @@ -191,8 +191,8 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, Array4{}; const Array4 z_nd_arr = z_phys_nd->array(mfi); // Get qheating_rates if provided (for LW radiation coupling to top-down mixing) - const Array4 qheat_arr = (qheating_rates != nullptr) - ? qheating_rates->const_array(mfi) + const Array4 qheat_arr = (qheating_rates != nullptr) + ? qheating_rates->const_array(mfi) : Array4{}; const bool has_qheating_rates = (qheating_rates != nullptr); @@ -331,7 +331,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // Rib = (g*z/θv0) * (θv(z) - θv_surf) / ws² (lines 180-200) // References: Hong et al., Mon. Wea. Rev., 134, 2318-2341 (2006) // - + BL_PROFILE_VAR("YSUNew_PBLH_Passes", prof_pblh); ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { @@ -356,13 +356,13 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, int kpbl = klo; Real Rib = rib_base_arr(i,j,klo); // Seed from surface level bool above_critical = (Rib >= Ribcr); // Check if already above critical - + // Scan from klo+1 if not already above critical for (int kk = klo+1; !above_critical && kk <= khi; ++kk) { - if (rib_base_arr(i,j,kk) >= Ribcr) { - kpbl = kk; - above_critical = true; - break; + if (rib_base_arr(i,j,kk) >= Ribcr) { + kpbl = kk; + above_critical = true; + break; } kpbl = kk; // keep updating so kpbl = khi if never exceeded } @@ -491,7 +491,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, : (klo + myhalf) * gdata.CellSize(2); Rib0 = Rib; bool above_critical = (Rib >= Ribcr); // Check if already above critical - + // Scan rib_enhan_arr to find first crossing of Ribcr with linear interpolation for (int kk = klo+1; !above_critical && kk <= khi; ++kk) { if (rib_enhan_arr(i,j,kk) >= Ribcr) { kpbl = kk; above_critical = true; break; } @@ -513,7 +513,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr) : (khi + myhalf) * gdata.CellSize(2); const Real pblh_max = Real(0.9) * z_max; - + amrex::Real pblh_min; if (turbChoice.enable_ysu_terrain_pblh_floor && use_terrain_fitted_coords) { const Real dz_inv = geom.InvCellSize(2); @@ -577,18 +577,18 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // Countergradient: HGAMT = min(CFAC * u* * θ*, GAMCRT), where CFAC=7.8, GAMCRT=3K // WRF Reference: module_bl_ysu.F lines 220-250 const bool enable_ysu_sat_limiter = turbChoice.enable_ysu_sat_limiter; - + // ======================================================================== // Cloud-top detection for top-down mixing (H10 Section 3b) // ======================================================================== const bool enable_ysu_topdown = turbChoice.enable_ysu_topdown; const amrex::Real ysu_qcloud_threshold = turbChoice.ysu_qcloud_threshold; - + if (enable_ysu_topdown && use_moisture) { ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { cloud_top_arr(i, j, 0) = -1; // Initialize: no cloud found - + // Search from PBL top downward for cloud-top cell int kpbl = pbli_arr(i, j, 0); for (int kk = kpbl - 1; kk >= klo; --kk) { @@ -597,7 +597,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, qc_kk = cell_data(i, j, kk, moisture_indices.qc) / cell_data(i, j, kk, Rho_comp); if (moisture_indices.qi >= 0) qi_kk = cell_data(i, j, kk, moisture_indices.qi) / cell_data(i, j, kk, Rho_comp); - + if (qc_kk + qi_kk > ysu_qcloud_threshold) { cloud_top_arr(i, j, 0) = kk; break; @@ -611,7 +611,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, cloud_top_arr(i, j, 0) = -1; }); } - + ParallelFor(xybx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { const Real t_layer = t10av_arr(i, j, 0); @@ -641,17 +641,17 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, Real wscale = u_star_arr(i, j, 0) / phiM_safe; wscale = amrex::max(wscale, u_star_arr(i, j, 0) / Real(5.0)); // Mechanical turbulence floor wscale = amrex::min(wscale, Real(16.0) * u_star_arr(i, j, 0)); // Free convection ceiling - + // Compute top-down convective velocity scale wstar3_down (H10 Eq. 12) // This represents turbulence driven by radiative cooling at cloud top amrex::Real wstar3_down = zero; if (enable_ysu_topdown && cloud_top_arr(i, j, 0) >= klo) { int k_cloud_top = cloud_top_arr(i, j, 0); - + // Longwave cooling flux at cloud top (W/m^2 equivalent in K*m/s) // When RRTMGP is active, read from qheating_rates if available. // Otherwise use zero (feature disabled when radiation is off). - // Compute LRAD by integrating the LW heating rate (component index 1) + // Compute LRAD by integrating the LW heating rate (component index 1) // over the PBL column from surface (klo) to cloud top (k_cloud_top) amrex::Real LRAD = zero; if (has_qheating_rates) { @@ -666,7 +666,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, LRAD += -qheat_arr(i, j, kk, 1) * dz; } } - + // Top-down convective velocity (H10 Eq. 12): // wstar_down^3 = g/theta * LRAD/(rho*cp) * pblh const amrex::Real t_local = cell_data(i, j, k_cloud_top, RhoTheta_comp) @@ -677,13 +677,13 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, * pblh_corr_arr(i, j, 0), zero); } wstar3_down_arr(i, j, 0) = wstar3_down; - + // Recompute wstar3 and wscale with corrected pblh_corr_arr const Real bfx0_corr = amrex::max(sflux_arr(i, j, 0), zero); const Real t_dry = t10av_arr(i, j, 0); const Real wstar3_corr = (CONST_GRAV / t_dry) * bfx0_corr * pblh_corr_arr(i, j, 0); wstar3_arr(i, j, 0) = wstar3_corr; - + // Recompute wscale with corrected pblh: const Real ust3 = u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0); Real wscale_corr = std::cbrt(ust3 + amrex::Real(8.0) * KAPPA * wstar3_corr * amrex::Real(0.5)); @@ -785,7 +785,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, : (klo + myhalf) * gdata.CellSize(2); constexpr Real sfcfrac_h = amrex::Real(0.1); // WRF SFCFRAC const Real height_lim = amrex::min(zl1_col / (sfcfrac_h * pblh), one); - + const Real VPERT_raw = HGAMT + amrex::Real(0.61) * t_layer * HGAMQ; const Real VPERT_capped = enable_ysu_unbounded_vpert ? VPERT_raw @@ -870,7 +870,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, : (klo + myhalf) * gdata.CellSize(2); Rib0 = Rib; bool above_critical = (Rib >= Ribcr); - + for (int kk = klo+1; !above_critical && kk <= khi; ++kk) { if (rib_enhan_arr(i,j,kk) >= Ribcr) { kpbl = kk; above_critical = true; break; } zval0 = (use_terrain_fitted_coords) @@ -891,7 +891,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, ? Compute_Zrel_AtCellCenter(i, j, khi, z_nd_arr) : (khi + myhalf) * gdata.CellSize(2); const Real pblh_max = Real(0.9) * z_max; - + amrex::Real pblh_min; if (turbChoice.enable_ysu_terrain_pblh_floor && use_terrain_fitted_coords) { const Real dz_inv = geom.InvCellSize(2); @@ -932,17 +932,17 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, { const int kpblold = pbli_arr(i,j,0); // Starting PBLH index from Pass 2 bool definebrup = false; - + // Surface liquid-theta virtual potential temperature reference const Real thermalli = GetThetavl(i, j, klo, cell_data, moisture_indices); - + // Upward scan using liquid-theta stability criterion for (int kk = kpblold; kk < khi; ++kk) { // Velocity shear for Richardson number calculation const Real ws2_raw = fourth * ((uvel(i,j,kk)+uvel(i+1,j,kk))*(uvel(i,j,kk)+uvel(i+1,j,kk)) + (vvel(i,j,kk)+vvel(i,j+1,kk))*(vvel(i,j,kk)+vvel(i,j+1,kk))); const Real ws2 = amrex::max(ws2_raw, amrex::Real(1.0)); - + // Elevation at level kk const Real z_sfc = (use_terrain_fitted_coords) ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr) @@ -951,28 +951,28 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, ? Compute_Zrel_AtCellCenter(i, j, kk, z_nd_arr) : (kk + myhalf) * gdata.CellSize(2); const Real zrel_kk = amrex::max(zval_kk - z_sfc, amrex::Real(1.0e-4)); - + // Liquid-theta virtual potential temperature at current level and surface const Real thlix_kk = GetThetavl(i, j, kk, cell_data, moisture_indices); const Real thlix_klo = GetThetavl(i, j, klo, cell_data, moisture_indices); - + // Bulk Richardson number using liquid-theta potential temperature const Real bruptmp = CONST_GRAV * zrel_kk * (thlix_kk - thermalli) / (ws2 * thlix_klo); - + // Stability check: threshold is zero for extension scan const bool stable = (bruptmp >= zero); - + if (definebrup) { pbli_arr(i,j,0) = kk; definebrup = false; } - + if (!stable) { // Continue scanning while Richardson number indicates mixing definebrup = true; } } - + // Bound result to domain top pbli_arr(i,j,0) = amrex::min(pbli_arr(i,j,0), izmax); }); @@ -1008,7 +1008,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, int kpbl_zero = klo; Real Rib = rib_enhan_arr(i,j,klo); // Seed from surface level bool above_critical = (Rib >= Ribcr_zero); // Check if already above critical (Ribcr=0) - + // Scan rib_enhan_arr(i,j,klo..khi) for Ribcr=0 crossing (no interpolation needed) for (int kk = klo+1; !above_critical && kk <= khi; ++kk) { if (rib_enhan_arr(i,j,kk) >= Ribcr_zero) { kpbl_zero = kk; above_critical = true; break; } @@ -1074,14 +1074,14 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // Apply entrainment only when we < 0 (downward entrainment occurring) const Real K_cap = Real(5.0) * rho_kpbl * wscale * KAPPA * pblh * Real(0.01); - + // Cloudy entrainment correction (WRF bl_ysu.F90 lines 840-875) // Applied when qc+qi exceeds threshold at cell kpbl-1 (layer below PBL top) // Uses cloud liquid water content and buoyancy jump for entrainment efficiency. // Radiation coupling is not required; wstar3_2 and hgamt2 remain zero. Real we_final = we; Real K_entr_final = (we < zero) ? rho_kpbl * (-we) * dz_kpbl : zero; - + if (enable_ysu_topdown && use_moisture) { const int k_below = amrex::max(kpbl - 1, klo); const Real rho_k = cell_data(i, j, k_below, Rho_comp); @@ -1130,25 +1130,25 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // Recompute entrainment coefficient with adjusted entrainment velocity K_entr_final = (we_final < zero) ? rho_kpbl * (-we_final) * dz_kpbl : zero; - + // Note: WRF zeroes xkzm at k_below when cloud is detected. // In ERF, entrainment is applied only at kpbl level, making this implicit. } } - + entr_arr(i,j,0) = amrex::min(K_entr_final, K_cap); - + // Cloud PBLH extension: check for cloud at PBL top cell and extend if present if (enable_ysu_cloud_pblh && use_moisture) { const int kpbl_current = pbli_arr(i, j, 0); - + // Check for cloud condensate at PBL top cell amrex::Real qc_kpbl = zero, qi_kpbl = zero; if (moisture_indices.qc >= 0) qc_kpbl = cell_data(i, j, kpbl_current, moisture_indices.qc) / cell_data(i, j, kpbl_current, Rho_comp); if (moisture_indices.qi >= 0) qi_kpbl = cell_data(i, j, kpbl_current, moisture_indices.qi) / cell_data(i, j, kpbl_current, Rho_comp); - + // Extend PBL top by one cell if cloud present exceeds threshold if ((qc_kpbl + qi_kpbl) > ysu_qcloud_threshold) { pbli_arr(i, j, 0) = amrex::min(kpbl_current + 1, izmax); @@ -1161,7 +1161,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, FArrayBox K_down_fab(gbx, 1, The_Async_Arena()); K_down_fab.setVal(zero); const auto& K_down_arr = K_down_fab.array(); - + if (enable_ysu_topdown) { ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -1359,7 +1359,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // Step 2: Compute height-dependent zq_kp1 for prnumfac calculation const Real zq_kp1_prandtl = zval + myhalf * dz_terrain; - + // Step 3: Compute prfac (surface layer Prandtl correction) // WRF bl_ysu.F90 line 951: prfac = conpr if SFCFLG, else 0 const Real prfac = SFCFLG ? conpr : zero; @@ -1372,7 +1372,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // Top-down wstar3 term (wstar3_2 in WRF bl_ysu.F90 line 949) // Non-zero only when top-down convection is active (requires radiation coupling for LRAD) const Real wstar3_2 = wstar3_down_arr(i, j, 0); - + Real prfac2 = zero; if (SFCFLG && ust3 > amrex::Real(1.0e-10)) { const Real wstar_tot3 = wstar3_col + wstar3_2; @@ -1417,11 +1417,11 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // PHASE 13 CHANGES: Use layer-top interface height (zq_kp1) and // level-dependent wscalek for proper K-profile formulation. // WRF Reference: module_bl_ysu.F lines 943-961 - + // Step 1: Compute zq_kp1 (top interface of cell k) // zq(k+1) in WRF = cell center + dz/2 const Real zq_kp1 = zval + myhalf * dz_terrain; - + // Step 2: Get surface height and first-level height (zl1) const Real z_sfc = (use_terrain_fitted_coords) ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr) @@ -1429,29 +1429,29 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, const Real zl1 = (use_terrain_fitted_coords) ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr) : (klo + myhalf) * gdata.CellSize(2); - + // Step 3: Compute zfac using zq_kp1 and zl1 (WRF bl_ysu.F90 line 943) // zfac = min(max((1-(zq(k+1)-zl1)/(hpbl-zl1)), zfmin), 1.) constexpr Real zfacmin = amrex::Real(1.0e-8); const Real pblh_rel = amrex::max(pblh_corr_arr(i, j, 0) - zl1, amrex::Real(1.0e-4)); const Real zfac = amrex::min( amrex::max(one - (zq_kp1 - zl1) / pblh_rel, zfacmin), one); - + // Step 4: Compute level-dependent wscalek (for K-profile, NOT for HGAMT/HGAMQ) // WRF stores wstar3 per column; we now have it from Phase 12 const Real ust3_wscale = u_star_arr(i, j, 0) * u_star_arr(i, j, 0) * u_star_arr(i, j, 0); - + Real wscalek_val; // Unstable/neutral: wscalek = (ust3 + phifac*karman*wstar3*(1-zfac))^(1/3) wscalek_val = std::cbrt(ust3_wscale + amrex::Real(8.0) * KAPPA * wstar3_col * (one - zfac)); - + // Step 5: Compute K_m using zq_kp1 and wscalek_val (WRF line 961) // WRF: xkzm(i,k) = wscalek(k)*karman*zq(k+1)*zfac(k)**pfac constexpr Real ckz_pbl = Real(0.001); const Real K_base = ckz_pbl * dz_terrain * rho; constexpr Real pfac = amrex::Real(2.0); K_turb(i, j, k, EddyDiff::Mom_v) = K_base + rho * wscalek_val * KAPPA * zq_kp1 * std::pow(zfac, pfac); - + // Apply Prandtl number to get heat and moisture diffusivity K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / Prt; if (turbChoice.ysu_moistvars) { @@ -1459,7 +1459,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, } else { K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Theta_v); } - + // Add top-down contribution (H10 Eq. 11) if (enable_ysu_topdown && k < pbli_arr(i, j, 0)) { K_turb(i, j, k, EddyDiff::Mom_v) += K_down_arr(i, j, k); @@ -1470,21 +1470,21 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // WRF uses wscalek = ust/phi_m(zq(k+1)/L), NOT Richardson mixing inside PBL // Richardson mixing is ONLY for k >= pbli_extent (free atmosphere) // WRF Reference: bl_ysu.F90 lines 951-957 - + // Step 1: Compute zq_kp1 (top interface of cell k) const Real zq_kp1_stable = zval + myhalf * dz_terrain; - + // Step 2: Get first-level height (zl1) const Real zl1_stable = (use_terrain_fitted_coords) ? Compute_Zrel_AtCellCenter(i, j, klo, z_nd_arr) : (klo + myhalf) * gdata.CellSize(2); - + // Step 3: Compute zfac for stable PBL constexpr Real zfacmin_stable = amrex::Real(1.0e-8); const Real pblh_rel_stable = amrex::max(pblh_corr_arr(i, j, 0) - zl1_stable, amrex::Real(1.0e-4)); const Real zfac_stable = amrex::min( amrex::max(one - (zq_kp1_stable - zl1_stable) / pblh_rel_stable, zfacmin_stable), one); - + // Step 4: Compute stable wscalek using phi_m // WRF bl_ysu.F90 lines 951-957: wscalek = ust / phi_m(zq(k+1)/L) const Real zol1_stable = zol1_arr_cap(i, j, 0); // stored from Phase 12 @@ -1494,19 +1494,19 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, const Real wscalek_stable = amrex::max( u_star_arr(i, j, 0) / amrex::max(phim_stable, amrex::Real(0.01)), amrex::Real(0.001)); - + // Step 5: Compute K_m for stable PBL using wscalek constexpr Real ckz_pbl_stable = Real(0.001); const Real K_base_stable = ckz_pbl_stable * dz_terrain * rho; constexpr Real pfac_stable = amrex::Real(2.0); K_turb(i, j, k, EddyDiff::Mom_v) = K_base_stable + rho * wscalek_stable * KAPPA * zq_kp1_stable * std::pow(zfac_stable, pfac_stable); - + // Step 6: Apply Prandtl number for stable PBL // For stable, prfac=0 (from step 3 of SECTION A) // prnum = phiH/phiM without height blending adjustment const Real prnum_stable = one + (prnum0 - one) * std::exp(prnumfac); K_turb(i, j, k, EddyDiff::Theta_v) = K_turb(i, j, k, EddyDiff::Mom_v) / prnum_stable; - + if (turbChoice.ysu_moistvars) { K_turb(i, j, k, EddyDiff::Q_v) = K_turb(i, j, k, EddyDiff::Mom_v) / prnum_q; } else { @@ -1559,7 +1559,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, // to prevent extreme floating-point scales from causing numerical instability Real grad_Ri = CONST_GRAV / theta_v * dtheta_v_dz / wind_shear_safe; grad_Ri = std::max(std::min(grad_Ri, Real(100.0)), -Real(100.0)); - + // GAP 10: In-cloud moist Richardson number modification // WRF bl_ysu.F90 lines 992-1000 (imvdif == 1, which is hardcoded in WRF) // Only apply when both the current AND next cell are in-cloud (qc+qi > 0.01e-3 kg/kg) @@ -1601,7 +1601,7 @@ ComputeDiffusivityYSUNew (const MultiFab& xvel, grad_Ri = amrex::max(amrex::min(grad_Ri, amrex::Real(100.0)), amrex::Real(-100.0)); } } - + // YSU stability functions (Hong et al. 2006, MWR, Appendix A) // Reference: https://doi.org/10.1175/MWR3250.1 // See equations A19-A20 diff --git a/Source/TimeIntegration/ERF_ImplicitPost.H b/Source/TimeIntegration/ERF_ImplicitPost.H index 41ce0401a1..ddd018463c 100644 --- a/Source/TimeIntegration/ERF_ImplicitPost.H +++ b/Source/TimeIntegration/ERF_ImplicitPost.H @@ -1,8 +1,7 @@ // ***************************************************************************** - // Do semi-implicit solve for diffusion: q1 + // Do semi-implicit solve for diffusion: TKE // ***************************************************************************** const bool l_do_implicit_ke = solverChoice.implicit_ke_diffusion; - const bool l_do_implicit_moist = solverChoice.implicit_moisture_diffusion; const Real l_vert_implicit_fac = solverChoice.vert_implicit_fac[level][nrk]; if ( l_vert_implicit_fac > zero ) { @@ -34,19 +33,19 @@ Array4{}; if (l_use_stretched_dz) { - ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoKE_comp, + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoKE_comp, RhoKE_comp, slow_dt, l_bc_neumann_vals_d, cell_data, stretched_dz_d[level], Array4{}, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } else if (l_use_terrain_fitted_coords) { - ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoKE_comp, + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoKE_comp, RhoKE_comp, slow_dt, l_bc_neumann_vals_d, cell_data, z_nd_arr, detJ_arr, dxInv, Array4{}, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } else { // no terrain - ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoKE_comp, + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoKE_comp, RhoKE_comp, slow_dt, l_bc_neumann_vals_d, cell_data, dxInv, Array4{}, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); @@ -54,53 +53,4 @@ } // mfi } // do implicit tke and evolve tke - if ( (l_do_implicit_moist) && - (solverChoice.moisture_type != MoistureType::None) ) { - - // BCs for q1 - const BCRec* bc_ptr_h = domain_bcs_type.data(); - GpuArray l_bc_neumann_vals_d; - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { - l_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoQ1_comp][ori]; - } - const bool l_use_SurfLayer = (m_SurfaceLayer != nullptr); - - for ( MFIter mfi(S_new[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - Box bx = mfi.tilebox(); - - // NOTE: We have updated the state, use that directly - const Array4< Real>& cell_data = S_new[IntVars::cons].array(mfi); - - const Array4& z_nd_arr = z_phys_nd[level]->const_array(mfi); - const Array4& detJ_arr = detJ_cc[level]->const_array(mfi); - - const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : - Array4{}; - const Array4& q1fx_z = Q1fx3->const_array(mfi); - const bool l_use_mrf_cg = solverChoice.turbChoice[level].enable_mrf_countergradient - && (solverChoice.turbChoice[level].pbl_type == PBLType::MRF); - if (l_use_stretched_dz) { - ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoQ1_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, - stretched_dz_d[level], q1fx_z, - mu_turb, solverChoice, - bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, - l_use_mrf_cg); - } else if (l_use_terrain_fitted_coords) { - ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoQ1_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, - z_nd_arr, detJ_arr, dxInv, q1fx_z, - mu_turb, solverChoice, - bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, - l_use_mrf_cg); - } else { // no terrain - ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoQ1_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, - dxInv, q1fx_z, mu_turb, solverChoice, - bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, - l_use_mrf_cg); - } - } // mfi - } // if moist model and do implicit diff } // implicit vertical factor > 0 diff --git a/Source/TimeIntegration/ERF_ImplicitPre.H b/Source/TimeIntegration/ERF_ImplicitPre.H index 5644097e0b..1eaf36130f 100644 --- a/Source/TimeIntegration/ERF_ImplicitPre.H +++ b/Source/TimeIntegration/ERF_ImplicitPre.H @@ -6,7 +6,9 @@ if (l_vert_implicit_fac > zero) { + const bool has_moisture = (solverChoice.moisture_type != MoistureType::None); const bool l_do_implicit_theta = solverChoice.implicit_thermal_diffusion; + const bool l_do_implicit_moist = solverChoice.implicit_moisture_diffusion; const bool l_do_implicit_mom = solverChoice.implicit_momentum_diffusion; // If we're doing an implicit solve for momenta (u and v only), @@ -22,11 +24,16 @@ MultiFab* Tau33corr = (l_do_implicit_mom) ? Tau_corr[level][2].get() : nullptr; #endif - // BCs const BCRec* bc_ptr_h = domain_bcs_type.data(); - GpuArray l_bc_neumann_vals_d; + + GpuArray l_th_bc_neumann_vals_d; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_th_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoTheta_comp][ori]; + } + + GpuArray l_qv_bc_neumann_vals_d; for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { - l_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoTheta_comp][ori]; + l_qv_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoQ1_comp][ori]; } const bool l_use_SurfLayer = (m_SurfaceLayer != nullptr); @@ -62,6 +69,8 @@ const Array4 tau23 = Tau[level][TauType::tau23]->array(mfi); [[maybe_unused]] const Array4 tau33 = Tau[level][TauType::tau33]->array(mfi); const Array4& hfx_z = Hfx3->const_array(mfi); + const Array4 qfx_z = Q1fx3 ? Q1fx3->const_array(mfi) : + Array4{}; const bool l_use_mrf_cg = solverChoice.turbChoice[level].enable_mrf_countergradient && (solverChoice.turbChoice[level].pbl_type == PBLType::MRF); const bool l_use_ysu_mom_cg = solverChoice.turbChoice[level].enable_ysu_countergradient @@ -69,13 +78,21 @@ if (l_use_stretched_dz) { if (l_do_implicit_theta) { - ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoTheta_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, + slow_dt, l_th_bc_neumann_vals_d, cell_data, stretched_dz_d[level], hfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, l_use_mrf_cg); } + if (l_do_implicit_moist && has_moisture) { + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, + slow_dt, l_qv_bc_neumann_vals_d, cell_data, + stretched_dz_d[level], qfx_z, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, + l_use_mrf_cg); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_S<0>(tbx, fine_geom.Domain(), level, slow_dt, cell_data, rho_u, tau13, tau13_corr, @@ -101,13 +118,21 @@ } } else if (l_use_terrain_fitted_coords) { if (l_do_implicit_theta) { - ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoTheta_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, + slow_dt, l_th_bc_neumann_vals_d, cell_data, z_nd_arr, detJ_arr, dxInv, hfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, l_use_mrf_cg); } + if (l_do_implicit_moist && has_moisture) { + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, + slow_dt, l_qv_bc_neumann_vals_d, cell_data, + z_nd_arr, detJ_arr, dxInv, qfx_z, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, + l_use_mrf_cg); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_T<0>(tbx, fine_geom.Domain(), level, slow_dt, cell_data, rho_u, tau13, tau13_corr, @@ -133,12 +158,19 @@ } } else { // no terrain if (l_do_implicit_theta) { - ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoTheta_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, + slow_dt, l_th_bc_neumann_vals_d, cell_data, dxInv, hfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, l_use_mrf_cg); } + if (l_do_implicit_moist && has_moisture) { + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, + slow_dt, l_qv_bc_neumann_vals_d, cell_data, + dxInv, qfx_z, mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, + l_use_mrf_cg); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_N<0>(tbx, fine_geom.Domain(), level, slow_dt, cell_data, rho_u, tau13, tau13_corr, diff --git a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp index 1038798e2d..4afd81728c 100644 --- a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp +++ b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp @@ -123,16 +123,32 @@ void make_fast_coeffs (int /*level*/, Real detJ_on_kface = myhalf * (detJ(i,j,k) + detJ(i,j,k-1)); Real inv_detJ_on_kface = one / detJ_on_kface; - Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; - Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; - - Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_p) + Real thd_km2 = prim(i,j,k-2,PrimTheta_comp); + Real thd_km1 = prim(i,j,k-1,PrimTheta_comp); + Real thd_k = prim(i,j,k ,PrimTheta_comp); + Real thd_kp1 = prim(i,j,k+1,PrimTheta_comp); + + Real qv_km2 = (l_use_moisture) ? prim(i,j,k-2,PrimQ1_comp) : zero; + Real qv_km1 = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; + Real qv_k = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; + Real qv_kp1 = (l_use_moisture) ? prim(i,j,k+1,PrimQ1_comp) : zero; + + Real thm_km2 = thd_km2 * (one + RvOverRd*qv_km2); + Real thm_km1 = thd_km1 * (one + RvOverRd*qv_km1); + Real thm_k = thd_k * (one + RvOverRd*qv_k ); + Real thm_kp1 = thd_kp1 * (one + RvOverRd*qv_kp1); + + Real thm_t_lo = myhalf * (thm_km2 + thm_km1); + Real thm_t_mid = myhalf * (thm_km1 + thm_k ); + Real thm_t_hi = myhalf * (thm_k + thm_kp1); + + Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) / - ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) ); + ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) * (one + RvOverRd*qv_k) ); - Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_q) + Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * pi_c + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) / - ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) ); + ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) * (one + RvOverRd*qv_km1) ); coeffP_a(i,j,k) = coeff_P; coeffQ_a(i,j,k) = coeff_Q; @@ -144,17 +160,13 @@ void make_fast_coeffs (int /*level*/, coeff_Q /= (one + q); } - Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) ); - Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) ); - Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) ); - // LHS for tri-diagonal system Real D = beta_2 * beta_2 * dzi * static_cast(dtau * dtau); - coeffA_a(i,j,k) = D * (one/detJ(i,j,k-1)) * ( halfg - coeff_Q * theta_t_lo ); - coeffC_a(i,j,k) = D * (one/detJ(i,j,k )) * (-halfg + coeff_P * theta_t_hi ); + coeffA_a(i,j,k) = D * (one/detJ(i,j,k-1)) * ( halfg - coeff_Q * thm_t_lo ); + coeffC_a(i,j,k) = D * (one/detJ(i,j,k )) * (-halfg + coeff_P * thm_t_hi ); - coeffB_a(i,j,k) = one + D * ( (coeff_Q/detJ(i,j,k-1) - coeff_P/detJ(i,j,k)) * theta_t_mid - + halfg * (Real(1.0)/detJ(i,j,k) - Real(1.0)/detJ(i,j,k-1)) ); + coeffB_a(i,j,k) = one + D * ( (coeff_Q/detJ(i,j,k-1) - coeff_P/detJ(i,j,k)) * thm_t_mid + + halfg * (one/detJ(i,j,k) - one/detJ(i,j,k-1)) ); }); } else { diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index eb6384c7fa..931dbf4dc5 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -171,11 +171,11 @@ void erf_slow_rhs_post (int level, int finest_level, } // Valid vars - Vector is_valid_slow_var; is_valid_slow_var.resize(RhoQ1_comp+1,0); + Vector is_valid_slow_var; is_valid_slow_var.resize(RhoQ2_comp+1,0); if (l_use_KE) { is_valid_slow_var[ RhoKE_comp] = 1; } if (l_do_scalar) { is_valid_slow_var[RhoScalar_comp] = 1; } if (solverChoice.moisture_type != MoistureType::None) { - is_valid_slow_var[RhoQ1_comp] = 1; + is_valid_slow_var[RhoQ2_comp] = 1; } // ************************************************************************* @@ -275,12 +275,12 @@ void erf_slow_rhs_post (int level, int finest_level, const GpuArray ncomp_slow = {nsv,0,0,0}; // ************************************************************************** - // Note that here we do copy only the "slow" variables, not (rho) or (rho theta) + // Note that here we do copy only the "slow" variables, not (Rho, RhoTheta, RhoQv) // ************************************************************************** ParallelFor(tbx, ncomp_slow[IntVars::cons], [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) { const int n = scomp_slow[IntVars::cons] + nn; - cur_cons(i,j,k,n) = new_cons(i,j,k,n); + if (n != RhoQ1_comp) { cur_cons(i,j,k,n) = new_cons(i,j,k,n); } }); // We have projected the velocities stored in S_data but we will use @@ -373,14 +373,14 @@ void erf_slow_rhs_post (int level, int finest_level, // // Note that we either advect and diffuse all or none of the moisture variables // - for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar) + for (int ivar(RhoKE_comp); ivar<= RhoQ2_comp; ++ivar) { if (is_valid_slow_var[ivar]) { start_comp = ivar; - num_comp = 1; + num_comp = 1; - if (ivar == RhoQ1_comp) { + if (ivar == RhoQ2_comp) { horiz_adv_type = ac.moistscal_horiz_adv_type; vert_adv_type = ac.moistscal_vert_adv_type; horiz_upw_frac = ac.moistscal_horiz_upw_frac; @@ -391,7 +391,7 @@ void erf_slow_rhs_post (int level, int finest_level, vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type); } - num_comp = n_qstate; + num_comp = n_qstate - 1; } else { horiz_adv_type = ac.dryscal_horiz_adv_type; @@ -436,10 +436,9 @@ void erf_slow_rhs_post (int level, int finest_level, if (l_use_diff) { - // Allow for implicit moisture diffusion + // Allow for implicit TKE diffusion Real l_vert_implicit_fac = zero; - if ( (ivar == RhoKE_comp && solverChoice.implicit_ke_diffusion ) || - (ivar == RhoQ1_comp && solverChoice.implicit_moisture_diffusion) ) { + if (ivar == RhoKE_comp && solverChoice.implicit_ke_diffusion) { l_vert_implicit_fac = solverChoice.vert_implicit_fac[level][nrk]; } @@ -493,14 +492,14 @@ void erf_slow_rhs_post (int level, int finest_level, auto const& src_arr = source.const_array(mfi); - for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar) + for (int ivar(RhoKE_comp); ivar<= RhoQ2_comp; ++ivar) { if (is_valid_slow_var[ivar]) { start_comp = ivar; - num_comp = 1; - if (ivar == RhoQ1_comp) { - num_comp = nvars - RhoQ1_comp; + num_comp = 1; + if (ivar == RhoQ2_comp) { + num_comp = nvars - RhoQ2_comp; } else if (ivar == RhoScalar_comp) { num_comp = NSCALARS; } diff --git a/Source/TimeIntegration/ERF_SlowRhsPre.cpp b/Source/TimeIntegration/ERF_SlowRhsPre.cpp index 58f89052dd..1fc0165005 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPre.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPre.cpp @@ -133,6 +133,10 @@ void erf_slow_rhs_pre (int level, int finest_level, const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type; const Real l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac; const Real l_vert_upw_frac = solverChoice.advChoice.dycore_vert_upw_frac; + const AdvType l_moist_horiz_adv_type = solverChoice.advChoice.moistscal_horiz_adv_type; + const AdvType l_moist_vert_adv_type = solverChoice.advChoice.moistscal_vert_adv_type; + const Real l_moist_horiz_upw_frac = solverChoice.advChoice.moistscal_horiz_upw_frac; + const Real l_moist_vert_upw_frac = solverChoice.advChoice.moistscal_vert_upw_frac; const bool l_use_stretched_dz = (solverChoice.mesh_type == MeshType::StretchedDz); const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type == MeshType::VariableDz); const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh); @@ -614,6 +618,15 @@ void erf_slow_rhs_pre (int level, int finest_level, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, flx_arr, domain, bc_ptr_h); + if (l_use_moisture) { + AdvectionSrcForScalars(bx, RhoQ1_comp, ncomp, + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + cell_prim, cell_rhs, + detJ_arr, dxInv, mf_mx, mf_my, + l_moist_horiz_adv_type, l_moist_vert_adv_type, + l_moist_horiz_upw_frac, l_moist_vert_upw_frac, + flx_arr, domain, bc_ptr_h); + } } else { EBAdvectionSrcForRho(bx, cell_rhs, rho_u, rho_v, omega_arr, @@ -634,6 +647,18 @@ void erf_slow_rhs_pre (int level, int finest_level, l_horiz_upw_frac, l_vert_upw_frac, flx_arr, domain, bc_ptr_h, already_on_centroids); + if (l_use_moisture) { + EBAdvectionSrcForScalars(bx, RhoQ1_comp, ncomp, + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + cell_prim, cell_rhs, + mask_arr, cfg_arr, ax_arr, ay_arr, az_arr, + fcx_arr, fcy_arr, fcz_arr, + detJ_arr, dxInv, mf_mx, mf_my, + l_moist_horiz_adv_type, l_moist_vert_adv_type, + l_moist_horiz_upw_frac, l_moist_vert_upw_frac, + flx_arr, domain, bc_ptr_h, + already_on_centroids); + } } if (l_use_diff) { @@ -675,6 +700,18 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_apply_surface_layer_fluxes_in_diffusion, l_vert_implicit_fac); + if (l_use_moisture) { + DiffusionSrcForState_S(bx, domain, RhoQ1_comp, n_comp, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + stretched_dz_d, dxInv, SmnSmn_a, + mf_mx, mf_ux, mf_vx, + mf_my, mf_uy, mf_vy, + hfx_z, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, + l_apply_surface_layer_fluxes_in_diffusion, l_vert_implicit_fac); + } } else if (l_use_terrain_fitted_coords) { DiffusionSrcForState_T(bx, domain, n_start, n_comp, l_rotate, u, v, cell_data, cell_prim, cell_rhs, @@ -686,6 +723,19 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_apply_surface_layer_fluxes_in_diffusion, l_vert_implicit_fac); + if (l_use_moisture) { + DiffusionSrcForState_T(bx, domain, RhoQ1_comp, n_comp, l_rotate, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + z_nd, z_cc, ax_arr, ay_arr, az_arr, detJ_arr, + dxInv, SmnSmn_a, + mf_mx, mf_ux, mf_vx, + mf_my, mf_uy, mf_vy, + hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, + l_apply_surface_layer_fluxes_in_diffusion, l_vert_implicit_fac); + } } else if (l_use_eb) { DiffusionSrcForState_EB(bx, domain, n_start, n_comp, u, v, cell_data, cell_prim, cell_rhs, @@ -696,6 +746,17 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, q1fx_z, q2fx_z, hfx_EB, mu_turb, solverChoice, level, bc_ptr_d, l_apply_surface_layer_fluxes_in_diffusion); + if (l_use_moisture) { + DiffusionSrcForState_EB(bx, domain, RhoQ1_comp, n_comp, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + cfg_arr, ax_arr, ay_arr, az_arr, detJ_arr, + barea_arr, bcent_arr, + dx, dxInv, + hfx_z, q1fx_z, q2fx_z, hfx_EB, + mu_turb, solverChoice, level, + bc_ptr_d, l_apply_surface_layer_fluxes_in_diffusion); + } } else { DiffusionSrcForState_N(bx, domain, n_start, n_comp, u, v, cell_data, cell_prim, cell_rhs, @@ -706,6 +767,18 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_apply_surface_layer_fluxes_in_diffusion, l_vert_implicit_fac); + if (l_use_moisture) { + DiffusionSrcForState_N(bx, domain, RhoQ1_comp, n_comp, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + dxInv, SmnSmn_a, + mf_mx, mf_ux, mf_vx, + mf_my, mf_uy, mf_vy, + hfx_z, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, + l_apply_surface_layer_fluxes_in_diffusion, l_vert_implicit_fac); + } } } @@ -714,6 +787,9 @@ void erf_slow_rhs_pre (int level, int finest_level, { cell_rhs(i,j,k,Rho_comp) += source_arr(i,j,k,Rho_comp); cell_rhs(i,j,k,RhoTheta_comp) += source_arr(i,j,k,RhoTheta_comp); + if (l_use_moisture) { + cell_rhs(i,j,k,RhoQ1_comp) += source_arr(i,j,k,RhoQ1_comp); + } }); Real half_dt = static_cast(myhalf/dt); @@ -725,9 +801,16 @@ void erf_slow_rhs_pre (int level, int finest_level, { cell_rhs(i,j,k, Rho_comp) *= myhalf; cell_rhs(i,j,k,RhoTheta_comp) *= myhalf; + if (l_use_moisture) { + cell_rhs(i,j,k,RhoQ1_comp) *= myhalf; + } cell_rhs(i,j,k, Rho_comp) += half_dt * (cell_data(i,j,k, Rho_comp) - cell_old(i,j,k, Rho_comp)); cell_rhs(i,j,k,RhoTheta_comp) += half_dt * (cell_data(i,j,k,RhoTheta_comp) - cell_old(i,j,k,RhoTheta_comp)); + if (l_use_moisture) { + cell_rhs(i,j,k,RhoQ1_comp) += half_dt * + (cell_data(i,j,k,RhoQ1_comp) - cell_old(i,j,k,RhoQ1_comp)); + } }); } diff --git a/Source/TimeIntegration/ERF_Substep_T.cpp b/Source/TimeIntegration/ERF_Substep_T.cpp index a56c726e55..7172dee905 100644 --- a/Source/TimeIntegration/ERF_Substep_T.cpp +++ b/Source/TimeIntegration/ERF_Substep_T.cpp @@ -145,7 +145,7 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& old_drho_u = Delta_rho_u.array(mfi); const Array4& old_drho_v = Delta_rho_v.array(mfi); const Array4& old_drho_w = Delta_rho_w.array(mfi); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); + const Array4& old_drho_thm = Delta_rho_theta.array(mfi); const Array4& prev_xmom = S_prev[IntVars::xmom].const_array(mfi); const Array4& prev_ymom = S_prev[IntVars::ymom].const_array(mfi); @@ -158,13 +158,6 @@ void erf_substep_T (int step, int /*nrk*/, Box bx = mfi.validbox(); Box gbx = mfi.tilebox(); gbx.grow(1); - if (step == 0) { - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); - cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); - }); - } // step = 0 - Box gtbx = mfi.nodaltilebox(0); gtbx.grow(IntVect(1,1,0)); Box gtby = mfi.nodaltilebox(1); gtby.grow(IntVect(1,1,0)); Box gtbz = mfi.nodaltilebox(2); gtbz.grow(IntVect(1,1,0)); @@ -172,6 +165,20 @@ void erf_substep_T (int step, int /*nrk*/, const auto& bx_lo = lbound(bx); const auto& bx_hi = ubound(bx); + // ************************************************************************* + // Transfer from S_old to S_data in first step (update in place in S_data) + // ************************************************************************* + if (step == 0) { + ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); + cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); + if (l_use_moisture) { cur_cons(i,j,k,RhoQ1_comp) = prev_cons(i,j,k,RhoQ1_comp); } + }); + } + + // ************************************************************************* + // Compute perturbational lateral momenta + // ************************************************************************* ParallelFor(gtbx, gtby, gtbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { old_drho_u(i,j,k) = prev_xmom(i,j,k) - stage_xmom(i,j,k); @@ -193,43 +200,30 @@ void erf_substep_T (int step, int /*nrk*/, old_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k); }); - const Array4& theta_extrap = extrap.array(mfi); - const Array4& prim = S_stage_prim.const_array(mfi); - + // ************************************************************************* + // Compute pert density, moist potential temperature, and thm extrap + // ************************************************************************* + const Array4& thm_extrap = extrap.array(mfi); + const Array4& prim = S_stage_prim.const_array(mfi); ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp); - old_drho_theta(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) - stage_cons(i,j,k,RhoTheta_comp); + Real qv_stg = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero; + Real qv_cur = (l_use_moisture) ? cur_cons(i,j,k,RhoQ1_comp ) / cur_cons(i,j,k,Rho_comp) : zero; + old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp); + old_drho_thm(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) * (Real(1.0) + RvOverRd*qv_cur) + - stage_cons(i,j,k,RhoTheta_comp) * (Real(1.0) + RvOverRd*qv_stg); if (step == 0) { - theta_extrap(i,j,k) = old_drho_theta(i,j,k); + thm_extrap(i,j,k) = old_drho_thm(i,j,k); } else { - theta_extrap(i,j,k) = old_drho_theta(i,j,k) + beta_d * - ( old_drho_theta(i,j,k) - lagged_arr(i,j,k) ); + thm_extrap(i,j,k) = old_drho_thm(i,j,k) + beta_d * ( old_drho_thm(i,j,k) - lagged_arr(i,j,k) ); } - - // NOTE: qv is not changing over the fast steps so we use the stage data - Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero; - theta_extrap(i,j,k) *= (one + RvOverRd*qv); + lagged_arr(i,j,k) = old_drho_thm(i,j,k); }); } // mfi -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // We define lagged_delta_rt for our next step as the current delta_rt - Box gbx = mfi.tilebox(); gbx.grow(1); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); - const Array4& lagged_arr = lagged_delta_rt.array(mfi); - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - lagged_arr(i,j,k) = old_drho_theta(i,j,k); - }); - } // mfi - - // ************************************************************************* - // Define updates in the current RK stage - // ************************************************************************* + // ********************************************************************* + // Update lateral momenta (explicit) + // ********************************************************************* #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -266,21 +260,12 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& pi_stage_ca = pi_stage.const_array(mfi); - const Array4& theta_extrap = extrap.array(mfi); + const Array4& thm_extrap = extrap.array(mfi); // Map factors const Array4& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi); const Array4& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi); - // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs - // Note that we do the Copy and Subtract including one ghost cell - // so that we don't have to fill ghost cells of the new MultiFabs - // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that - // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w - - // ********************************************************************* - // Define updates in the RHS of {x, y, z}-momentum equations - // ********************************************************************* { BL_PROFILE("substep_xymom_T"); @@ -293,14 +278,14 @@ void erf_substep_T (int step, int /*nrk*/, // Add (negative) gradient of (rho theta) multiplied by lagged "pi" Real met_h_xi = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); - Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi; + Real gp_xi = (thm_extrap(i,j,k) - thm_extrap(i-1,j,k)) * dxi; Real gp_zeta_on_iface = (k == 0) ? - myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1) - - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) : - fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1) - - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) ); - Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? Real(0.) : - gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface; + myhalf * dzi * ( thm_extrap(i-1,j,k+1) + thm_extrap(i,j,k+1) + - thm_extrap(i-1,j,k ) - thm_extrap(i,j,k ) ) : + fourth * dzi * ( thm_extrap(i-1,j,k+1) + thm_extrap(i,j,k+1) + - thm_extrap(i-1,j,k-1) - thm_extrap(i,j,k-1) ); + Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? zero : + gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface; gpx *= mf_ux(i,j,0); @@ -327,13 +312,13 @@ void erf_substep_T (int step, int /*nrk*/, // Add (negative) gradient of (rho theta) multiplied by lagged "pi" Real met_h_eta = Compute_h_eta_AtJface(i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); - Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi; + Real gp_eta = (thm_extrap(i,j,k) - thm_extrap(i,j-1,k)) * dyi; Real gp_zeta_on_jface = (k == 0) ? - myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1) - - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) : - fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1) - - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) ); - Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? Real(0.) : + myhalf * dzi * ( thm_extrap(i,j,k+1) + thm_extrap(i,j-1,k+1) + - thm_extrap(i,j,k ) - thm_extrap(i,j-1,k ) ) : + fourth * dzi * ( thm_extrap(i,j,k+1) + thm_extrap(i,j-1,k+1) + - thm_extrap(i,j,k-1) - thm_extrap(i,j-1,k-1) ); + Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? zero : gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface; gpy *= mf_vy(i,j,0); @@ -381,11 +366,11 @@ void erf_substep_T (int step, int /*nrk*/, const Array4 & stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi); const Array4 & prim = S_stage_prim.const_array(mfi); - const Array4& old_drho_u = Delta_rho_u.array(mfi); - const Array4& old_drho_v = Delta_rho_v.array(mfi); - const Array4& old_drho_w = Delta_rho_w.array(mfi); - const Array4& old_drho = Delta_rho.array(mfi); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); + const Array4& old_drho_u = Delta_rho_u.array(mfi); + const Array4& old_drho_v = Delta_rho_v.array(mfi); + const Array4& old_drho_w = Delta_rho_w.array(mfi); + const Array4& old_drho = Delta_rho.array(mfi); + const Array4& old_drho_thm = Delta_rho_theta.array(mfi); const Array4& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi); const Array4& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi); @@ -396,7 +381,6 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& cur_cons = S_data[IntVars::cons].array(mfi); const Array4& cur_zmom = S_data[IntVars::zmom].array(mfi); - // These store the advection momenta which we will use to update the slow variables const Array4& avg_zmom_arr = avg_zmom.array(mfi); const Array4& z_nd = z_phys_nd->const_array(mfi); @@ -404,7 +388,6 @@ void erf_substep_T (int step, int /*nrk*/, const Array4< Real>& omega_arr = Omega.array(mfi); - // Map factors const Array4& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi); const Array4& mf_my = mapfac[MapFacType::m_y]->const_array(mfi); const Array4& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi); @@ -412,19 +395,13 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi); const Array4& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi); - // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs - // Note that we do the Copy and Subtract including one ghost cell - // so that we don't have to fill ghost cells of the new MultiFabs - // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that - // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w - FArrayBox temp_rhs_fab; FArrayBox RHS_fab; FArrayBox soln_fab; RHS_fab.resize (tbz,1,The_Async_Arena()); soln_fab.resize (tbz,1,The_Async_Arena()); - temp_rhs_fab.resize(tbz,2,The_Async_Arena()); + temp_rhs_fab.resize(tbz,3,The_Async_Arena()); auto const& RHS_a = RHS_fab.array(); auto const& soln_a = soln_fab.array(); @@ -440,13 +417,43 @@ void erf_substep_T (int step, int /*nrk*/, // Define flux arrays for use in advection // ************************************************************************* for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - flux[dir].resize(surroundingNodes(bx,dir),2,The_Async_Arena()); + flux[dir].resize(surroundingNodes(bx,dir),3,The_Async_Arena()); flux[dir].setVal(0); } const GpuArray, AMREX_SPACEDIM> flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}}; - // ********************************************************************* + // ************************************************************************* + // Define old contravariant velocity (start of sub-step) + // ************************************************************************* + { + Box gbxo = mfi.nodaltilebox(2); + Box gbxo_mid = gbxo; + + if (gbxo.smallEnd(2) == domlo.z) { + Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2)); + gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1); + ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + omega_arr(i,j,k) = zero; + }); + } + if (gbxo.bigEnd(2) == domhi.z+1) { + Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2)); + gbxo_mid.setBig(2,gbxo.bigEnd(2)-1); + ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + omega_arr(i,j,k) = old_drho_w(i,j,k); + }); + } + ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + omega_arr(i,j,k) = OmegaFromW(i,j,k,old_drho_w(i,j,k), + old_drho_u,old_drho_v, + mf_ux,mf_vy,z_nd,dxInv); + }); + } // end profile + + // ************************************************************************* + // Define lateral update to rho & rho theta & qv^{*} + // ************************************************************************* { BL_PROFILE("fast_T_making_rho_rhs"); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -476,57 +483,50 @@ void erf_substep_T (int step, int /*nrk*/, // NOTE: we are saving the (1/J) weighting for later when we add this to rho and theta temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq + ( yflux_hi - yflux_lo ) * dyi * mfsq; - temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) - - xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+ - ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) - - yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) * myhalf; + temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,PrimTheta_comp) + prim(i+1,j,k,PrimTheta_comp)) - + xflux_lo * (prim(i,j,k,PrimTheta_comp) + prim(i-1,j,k,PrimTheta_comp)) ) * dxi * mfsq+ + ( yflux_hi * (prim(i,j,k,PrimTheta_comp) + prim(i,j+1,k,PrimTheta_comp)) - + yflux_lo * (prim(i,j,k,PrimTheta_comp) + prim(i,j-1,k,PrimTheta_comp)) ) * dyi * mfsq) * myhalf; + if (l_use_moisture) { + temp_rhs_arr(i,j,k,2) = (( xflux_hi * (prim(i,j,k,PrimQ1_comp) + prim(i+1,j,k,PrimQ1_comp)) - + xflux_lo * (prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp)) ) * dxi * mfsq+ + ( yflux_hi * (prim(i,j,k,PrimQ1_comp) + prim(i,j+1,k,PrimQ1_comp)) - + yflux_lo * (prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp)) ) * dyi * mfsq) * myhalf; + } else { + temp_rhs_arr(i,j,k,2) = zero; + } if (l_reflux) { (flx_arr[0])(i,j,k,0) = xflux_lo; - (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)); + (flx_arr[0])(i,j,k,1) = (flx_arr[0])(i,j,k,0) * myhalf * (prim(i,j,k,PrimTheta_comp) + prim(i-1,j,k,PrimTheta_comp)); + if (l_use_moisture) { + (flx_arr[0])(i,j,k,2) = (flx_arr[0])(i,j,k,0) * myhalf * (prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp)); + } (flx_arr[1])(i,j,k,0) = yflux_lo; - (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)); + (flx_arr[1])(i,j,k,1) = (flx_arr[1])(i,j,k,0) * myhalf * (prim(i,j,k,PrimTheta_comp) + prim(i,j-1,k,PrimTheta_comp)); + if (l_use_moisture) { + (flx_arr[1])(i,j,k,2) = (flx_arr[1])(i,j,k,0) * myhalf * (prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp)); + } if (i == vbx_hi.x) { (flx_arr[0])(i+1,j,k,0) = xflux_hi; - (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)); + (flx_arr[0])(i+1,j,k,1) = (flx_arr[0])(i+1,j,k,0) * myhalf * (prim(i,j,k,PrimTheta_comp) + prim(i+1,j,k,PrimTheta_comp)); + if (l_use_moisture) { + (flx_arr[0])(i+1,j,k,2) = (flx_arr[0])(i+1,j,k,0) * myhalf * (prim(i,j,k,PrimQ1_comp) + prim(i+1,j,k,PrimQ1_comp)); + } } if (j == vbx_hi.y) { (flx_arr[1])(i,j+1,k,0) = yflux_hi; - (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)); + (flx_arr[1])(i,j+1,k,1) = (flx_arr[1])(i,j+1,k,0) * myhalf * (prim(i,j,k,PrimTheta_comp) + prim(i,j+1,k,PrimTheta_comp)); + if (l_use_moisture) { + (flx_arr[1])(i,j+1,k,2) = (flx_arr[1])(i,j+1,k,0) * myhalf * (prim(i,j,k,PrimQ1_comp) + prim(i,j+1,k,PrimQ1_comp)); + } } } }); } // end profile - // ********************************************************************* - { - Box gbxo = mfi.nodaltilebox(2); - Box gbxo_mid = gbxo; - - if (gbxo.smallEnd(2) == domlo.z) { - Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2)); - gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1); - ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = zero; - }); - } - if (gbxo.bigEnd(2) == domhi.z+1) { - Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2)); - gbxo_mid.setBig(2,gbxo.bigEnd(2)-1); - ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = old_drho_w(i,j,k); - }); - } - ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = OmegaFromW(i,j,k,old_drho_w(i,j,k), - old_drho_u,old_drho_v, - mf_ux,mf_vy,z_nd,dxInv); - }); - } // end profile - // ********************************************************************* - Box bx_shrunk_in_k = bx; int klo = tbz.smallEnd(2); int khi = tbz.bigEnd(2); @@ -538,6 +538,9 @@ void erf_substep_T (int step, int /*nrk*/, // We define halfg to match the notes (which is why we take the absolute value) Real halfg = std::abs(myhalf * grav_gpu[2]); + // ************************************************************************* + // RHS for tridiagonal system + // ************************************************************************* { BL_PROFILE("fast_loop_on_shrunk_t"); //Note we don't act on the bottom or top boundaries of the domain @@ -546,19 +549,60 @@ void erf_substep_T (int step, int /*nrk*/, Real coeff_P = coeffP_a(i,j,k); Real coeff_Q = coeffQ_a(i,j,k); - Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) ); - Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) ); - Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) ); + // Cell center vars for constructing theta_m + Real thd_km2 = prim(i,j,k-2,PrimTheta_comp); + Real thd_km1 = prim(i,j,k-1,PrimTheta_comp); + Real thd_k = prim(i,j,k ,PrimTheta_comp); + Real thd_kp1 = prim(i,j,k+1,PrimTheta_comp); + + Real qv_km2 = (l_use_moisture) ? prim(i,j,k-2,PrimQ1_comp) : zero; + Real qv_km1 = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; + Real qv_k = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; + Real qv_kp1 = (l_use_moisture) ? prim(i,j,k+1,PrimQ1_comp) : zero; + + Real thm_km2 = thd_km2 * (Real(1.0) + RvOverRd*qv_km2); + Real thm_km1 = thd_km1 * (Real(1.0) + RvOverRd*qv_km1); + Real thm_k = thd_k * (Real(1.0) + RvOverRd*qv_k ); + Real thm_kp1 = thd_kp1 * (Real(1.0) + RvOverRd*qv_kp1); + + // Stage theta_m at w-faces + Real thm_t_lo = myhalf * (thm_km2 + thm_km1); + Real thm_t_mid = myhalf * (thm_km1 + thm_k ); + Real thm_t_hi = myhalf * (thm_k + thm_kp1); + + // Stage coefficients for linearization + Real A_T_k = (one + RvOverRd*qv_k ); + Real A_T_km1 = (one + RvOverRd*qv_km1); + Real A_Q_k = RvOverRd*prim(i,j,k ,PrimTheta_comp); + Real A_Q_km1 = RvOverRd*prim(i,j,k-1,PrimTheta_comp); + Real A_D_k = -RvOverRd*prim(i,j,k ,PrimTheta_comp)*qv_k; + Real A_D_km1 = -RvOverRd*prim(i,j,k-1,PrimTheta_comp)*qv_km1; + + // Sums of slow RHS + Real Q_srhs_k = (l_use_moisture) ? slow_rhs_cons(i,j,k ,RhoQ1_comp) : zero; + Real Q_srhs_km1 = (l_use_moisture) ? slow_rhs_cons(i,j,k-1,RhoQ1_comp) : zero; + Real A_sum_srhs_k = A_D_k * slow_rhs_cons(i,j,k ,Rho_comp) + + A_T_k * slow_rhs_cons(i,j,k ,RhoTheta_comp) + + A_Q_k * Q_srhs_k; + Real A_sum_srhs_km1 = A_D_km1 * slow_rhs_cons(i,j,k-1,Rho_comp) + + A_T_km1 * slow_rhs_cons(i,j,k-1,RhoTheta_comp) + + A_Q_km1 * Q_srhs_km1; + + // Sums of the temp RHS + Real A_sum_trhs_k = A_D_k * temp_rhs_arr(i,j,k ,Rho_comp) + + A_T_k * temp_rhs_arr(i,j,k ,RhoTheta_comp) + + A_Q_k * temp_rhs_arr(i,j,k ,RhoTheta_comp+1); + Real A_sum_trhs_km1 = A_D_km1 * temp_rhs_arr(i,j,k-1,Rho_comp) + + A_T_km1 * temp_rhs_arr(i,j,k-1,RhoTheta_comp) + + A_Q_km1 * temp_rhs_arr(i,j,k-1,RhoTheta_comp+1); // line 2 last two terms (order dtau) - Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k ) - -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1); + Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_thm(i,j,k ) + -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_thm(i,j,k-1); // line 3 residuals (order dtau^2) one <-> beta_2 Real R1_tmp = -halfg * ( slow_rhs_cons(i,j,k ,Rho_comp) + slow_rhs_cons(i,j,k-1,Rho_comp) ); - - R1_tmp += coeff_P * slow_rhs_cons(i,j,k ,RhoTheta_comp) - + coeff_Q * slow_rhs_cons(i,j,k-1,RhoTheta_comp); + R1_tmp += coeff_P * A_sum_srhs_k + coeff_Q * A_sum_srhs_km1; Real Omega_kp1 = omega_arr(i,j,k+1); Real Omega_k = omega_arr(i,j,k ); @@ -571,8 +615,8 @@ void erf_substep_T (int step, int /*nrk*/, + 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) ); // consolidate lines 6&7 (order dtau^2) - 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) ) - + 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) ) ); + R1_tmp += -( coeff_P/detJ(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*thm_t_hi - Omega_k*thm_t_mid ) + A_sum_trhs_k ) + + coeff_Q/detJ(i,j,k-1) * ( beta_1 * dzi * (Omega_k*thm_t_mid - Omega_km1*thm_t_lo) + A_sum_trhs_km1) ); // line 1 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); @@ -584,19 +628,22 @@ void erf_substep_T (int step, int /*nrk*/, }); } // end profile - Box b2d = tbz; // Copy constructor + Box b2d = tbz; b2d.setRange(2,0); auto const lo = lbound(bx); auto const hi = ubound(bx); + // ************************************************************************* + // Solve tridiagonal system + // ************************************************************************* { BL_PROFILE("substep_b2d_loop_t"); #ifdef AMREX_USE_GPU ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { // w_klo, w_khi given by specified Dirichlet values - RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z) + zmom_src_arr(i,j,lo.z)); + RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z ) + zmom_src_arr(i,j,lo.z )); 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)); // w = specified Dirichlet value at k = lo.z @@ -685,47 +732,59 @@ void erf_substep_T (int step, int /*nrk*/, }); // ************************************************************************** - // Define updates in the RHS of rho and (rho theta) + // Define updates in the RHS of rho, rho theta, rho qv // ************************************************************************** { BL_PROFILE("fast_rho_final_update"); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=, zero_d=zero] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k); + Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k ); Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1); // Note that in the solve we effectively impose new_drho_w(i,j,vbx_hi.z+1)=0 // so we don't update avg_zmom at k=vbx_hi.z+1 - avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0)); + avg_zmom_arr(i,j,k) += facinv*zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0)); if (l_reflux) { - (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0)); + (flx_arr[2])(i,j,k,0) = zflux_lo / (mf_mx(i,j,0) * mf_my(i,j,0)); } if (k == vbx_hi.z) { - avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0)); + avg_zmom_arr(i,j,k+1) += facinv * zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0)); if (l_reflux) { - (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0)); + (flx_arr[2])(i,j,k+1,0) = zflux_hi / (mf_mx(i,j,0) * mf_my(i,j,0)); (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)); } } - Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k); - - cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho); + Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k); + cur_cons(i,j,k,Rho_comp) += dtau * (slow_rhs_cons(i,j,k,Rho_comp) + fast_rhs_rho); Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) + myhalf * - ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) - - zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k); - - cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta); + ( zflux_hi * (prim(i,j,k,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp)) - + zflux_lo * (prim(i,j,k,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp)) ) * dzi ) / detJ(i,j,k); + cur_cons(i,j,k,RhoTheta_comp) += dtau * (slow_rhs_cons(i,j,k,RhoTheta_comp) + fast_rhs_rhotheta); + + if (l_use_moisture) { + Real fast_rhs_rhoqv = -( temp_rhs_arr(i,j,k,2) + myhalf * + ( zflux_hi * (prim(i,j,k,PrimQ1_comp) + prim(i,j,k+1,PrimQ1_comp)) - + zflux_lo * (prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)) ) * dzi ) / detJ(i,j,k); + cur_cons(i,j,k,RhoQ1_comp) += dtau * (slow_rhs_cons(i,j,k,RhoQ1_comp) + fast_rhs_rhoqv); + cur_cons(i,j,k,RhoQ1_comp) = amrex::max(zero_d,cur_cons(i,j,k,RhoQ1_comp)); + } if (l_reflux) { - (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k) + prim(i,j,k-1)); + (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp)); + if (l_use_moisture) { + (flx_arr[2])(i,j,k,2) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)); + } } // add in source terms for cell-centered conserved variables cur_cons(i,j,k,Rho_comp) += dtau * cc_src_arr(i,j,k,Rho_comp); cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp); + if (l_use_moisture) { + cur_cons(i,j,k,RhoQ1_comp) += dtau * cc_src_arr(i,j,k,RhoQ1_comp); + } }); } // end profile diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H index d37d7f8156..81c986aa66 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H @@ -199,10 +199,12 @@ Real sdt = static_cast(slow_dt); if ((solverChoice.vert_implicit_fac[level][nrk] > zero_d) && solverChoice.implicit_before_substep) { - MultiFab scratch(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(), 2, + MultiFab scratch(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(), 3, S_data[IntVars::cons].nGrowVect()); - MultiFab::Copy(scratch, S_old[IntVars::cons], 0, 0, 2, S_data[IntVars::cons].nGrowVect()); // scratch := S_old (for rho, rhotheta) - MultiFab::Saxpy(scratch, sdt, S_rhs[IntVars::cons], 0, 0, 2, 0); // scratch := S_old + slow_dt*Src (for rho, rhotheta) + MultiFab::Copy(scratch, S_old[IntVars::cons], 0, 0, 2, S_data[IntVars::cons].nGrowVect()); + MultiFab::Copy(scratch, S_old[IntVars::cons], RhoQ1_comp, 2, 1, S_data[IntVars::cons].nGrowVect()); + MultiFab::Saxpy(scratch, sdt, S_rhs[IntVars::cons], 0, 0, 2, 0); + MultiFab::Saxpy(scratch, sdt, S_rhs[IntVars::cons], RhoQ1_comp, 2, 1, 0); scratch.FillBoundary(geom[level].periodicity()); MultiFab scratch_xmom(S_data[IntVars::xmom].boxArray(), @@ -233,9 +235,11 @@ #include "ERF_ImplicitPre.H" - MultiFab::Saxpy(scratch, -one_d, S_old[IntVars::cons], 1, 1, 1, 0); // scratch := (S_new - S_old) (for rhotheta only) - scratch.mult(one_d / sdt); // scratch := (S_new - S_old) / slow_dt - MultiFab::Copy(S_rhs[IntVars::cons], scratch, 1, 1, 1, 0); // slow_rhs := (S_new - S_old) / slow_dt (for rhotheta only) + MultiFab::Saxpy(scratch, -one_d, S_old[IntVars::cons], 1, 1, 1, 0); + MultiFab::Saxpy(scratch, -one_d, S_old[IntVars::cons], RhoQ1_comp, 2, 1, 0); + scratch.mult(one_d / sdt); + MultiFab::Copy(S_rhs[IntVars::cons], scratch, 1, 1, 1, 0); + MultiFab::Copy(S_rhs[IntVars::cons], scratch, 2, RhoQ1_comp, 1, 0); if (solverChoice.implicit_momentum_diffusion) { MultiFab::Saxpy(scratch_xmom, -one_d, S_old[IntVars::xmom], 0, 0, 1, 0); // scratch := (S_new - S_old) diff --git a/Source/Utils/ERF_MoistUtils.H b/Source/Utils/ERF_MoistUtils.H index 093952be66..74d5b2228d 100644 --- a/Source/Utils/ERF_MoistUtils.H +++ b/Source/Utils/ERF_MoistUtils.H @@ -95,7 +95,7 @@ amrex::Real GetThetavl (int i, int j, int k, { const amrex::Real rho = cell_data(i, j, k, Rho_comp); const amrex::Real theta = cell_data(i, j, k, RhoTheta_comp) / rho; - + amrex::Real qv = zero, qc = zero, qi = zero; if (moisture_indices.qv >= 0) qv = cell_data(i, j, k, moisture_indices.qv) / rho; @@ -103,7 +103,7 @@ amrex::Real GetThetavl (int i, int j, int k, qc = cell_data(i, j, k, moisture_indices.qc) / rho; if (moisture_indices.qi >= 0) qi = cell_data(i, j, k, moisture_indices.qi) / rho; - + return theta * (one + amrex::Real(0.61) * qv - qc - qi); } From 7652304d0da5ca7b63b5d2c6740c442c667710d9 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Mon, 20 Jul 2026 15:01:07 -0700 Subject: [PATCH 2/3] Clean up constants. --- Source/TimeIntegration/ERF_MakeFastCoeffs.cpp | 17 +++++------ Source/TimeIntegration/ERF_Substep_T.cpp | 30 +++++++++---------- 2 files changed, 22 insertions(+), 25 deletions(-) diff --git a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp index 4afd81728c..6e6c25043e 100644 --- a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp +++ b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp @@ -42,7 +42,6 @@ void make_fast_coeffs (int /*level*/, Real beta_2 = myhalf * (one + beta_s); // multiplies implicit terms Real c_v = c_p - R_d; - Real RvOverRd = R_v / R_d; const GpuArray dxInv = geom.InvCellSizeArray(); Real dzi = dxInv[2]; @@ -133,10 +132,10 @@ void make_fast_coeffs (int /*level*/, Real qv_k = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; Real qv_kp1 = (l_use_moisture) ? prim(i,j,k+1,PrimQ1_comp) : zero; - Real thm_km2 = thd_km2 * (one + RvOverRd*qv_km2); - Real thm_km1 = thd_km1 * (one + RvOverRd*qv_km1); - Real thm_k = thd_k * (one + RvOverRd*qv_k ); - Real thm_kp1 = thd_kp1 * (one + RvOverRd*qv_kp1); + Real thm_km2 = thd_km2 * (one + RvoRd*qv_km2); + Real thm_km1 = thd_km1 * (one + RvoRd*qv_km1); + Real thm_k = thd_k * (one + RvoRd*qv_k ); + Real thm_kp1 = thd_kp1 * (one + RvoRd*qv_kp1); Real thm_t_lo = myhalf * (thm_km2 + thm_km1); Real thm_t_mid = myhalf * (thm_km1 + thm_k ); @@ -144,11 +143,11 @@ void make_fast_coeffs (int /*level*/, Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) / - ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) * (one + RvOverRd*qv_k) ); + ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) * (one + RvoRd*qv_k) ); Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * pi_c + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) / - ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) * (one + RvOverRd*qv_km1) ); + ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) * (one + RvoRd*qv_km1) ); coeffP_a(i,j,k) = coeff_P; coeffQ_a(i,j,k) = coeff_Q; @@ -184,11 +183,11 @@ void make_fast_coeffs (int /*level*/, Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; - Real coeff_P = -Gamma * R_d * dzi * pi_c * (one + RvOverRd*qv_p) + Real coeff_P = -Gamma * R_d * dzi * pi_c * (one + RvoRd*qv_p) + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) / ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) ); - Real coeff_Q = Gamma * R_d * dzi * pi_c * (one + RvOverRd*qv_q) + Real coeff_Q = Gamma * R_d * dzi * pi_c * (one + RvoRd*qv_q) + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) / ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) ); diff --git a/Source/TimeIntegration/ERF_Substep_T.cpp b/Source/TimeIntegration/ERF_Substep_T.cpp index 7172dee905..934c0800cc 100644 --- a/Source/TimeIntegration/ERF_Substep_T.cpp +++ b/Source/TimeIntegration/ERF_Substep_T.cpp @@ -92,8 +92,6 @@ void erf_substep_T (int step, int /*nrk*/, // How much do we project forward the (rho theta) that is used in the horizontal momentum equations Real beta_d = Real(0.1); - Real RvOverRd = R_v / R_d; - bool l_rayleigh_impl_for_w = (sinesq_stag_d != nullptr); const Real* dx = geom.CellSize(); @@ -209,8 +207,8 @@ void erf_substep_T (int step, int /*nrk*/, Real qv_stg = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero; Real qv_cur = (l_use_moisture) ? cur_cons(i,j,k,RhoQ1_comp ) / cur_cons(i,j,k,Rho_comp) : zero; old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp); - old_drho_thm(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) * (Real(1.0) + RvOverRd*qv_cur) - - stage_cons(i,j,k,RhoTheta_comp) * (Real(1.0) + RvOverRd*qv_stg); + old_drho_thm(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) * (one + RvoRd*qv_cur) + - stage_cons(i,j,k,RhoTheta_comp) * (one + RvoRd*qv_stg); if (step == 0) { thm_extrap(i,j,k) = old_drho_thm(i,j,k); } else { @@ -560,10 +558,10 @@ void erf_substep_T (int step, int /*nrk*/, Real qv_k = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; Real qv_kp1 = (l_use_moisture) ? prim(i,j,k+1,PrimQ1_comp) : zero; - Real thm_km2 = thd_km2 * (Real(1.0) + RvOverRd*qv_km2); - Real thm_km1 = thd_km1 * (Real(1.0) + RvOverRd*qv_km1); - Real thm_k = thd_k * (Real(1.0) + RvOverRd*qv_k ); - Real thm_kp1 = thd_kp1 * (Real(1.0) + RvOverRd*qv_kp1); + Real thm_km2 = thd_km2 * (one + RvoRd*qv_km2); + Real thm_km1 = thd_km1 * (one + RvoRd*qv_km1); + Real thm_k = thd_k * (one + RvoRd*qv_k ); + Real thm_kp1 = thd_kp1 * (one + RvoRd*qv_kp1); // Stage theta_m at w-faces Real thm_t_lo = myhalf * (thm_km2 + thm_km1); @@ -571,12 +569,12 @@ void erf_substep_T (int step, int /*nrk*/, Real thm_t_hi = myhalf * (thm_k + thm_kp1); // Stage coefficients for linearization - Real A_T_k = (one + RvOverRd*qv_k ); - Real A_T_km1 = (one + RvOverRd*qv_km1); - Real A_Q_k = RvOverRd*prim(i,j,k ,PrimTheta_comp); - Real A_Q_km1 = RvOverRd*prim(i,j,k-1,PrimTheta_comp); - Real A_D_k = -RvOverRd*prim(i,j,k ,PrimTheta_comp)*qv_k; - Real A_D_km1 = -RvOverRd*prim(i,j,k-1,PrimTheta_comp)*qv_km1; + Real A_T_k = (one + RvoRd*qv_k ); + Real A_T_km1 = (one + RvoRd*qv_km1); + Real A_Q_k = RvoRd*prim(i,j,k ,PrimTheta_comp); + Real A_Q_km1 = RvoRd*prim(i,j,k-1,PrimTheta_comp); + Real A_D_k = -RvoRd*prim(i,j,k ,PrimTheta_comp)*qv_k; + Real A_D_km1 = -RvoRd*prim(i,j,k-1,PrimTheta_comp)*qv_km1; // Sums of slow RHS Real Q_srhs_k = (l_use_moisture) ? slow_rhs_cons(i,j,k ,RhoQ1_comp) : zero; @@ -726,8 +724,8 @@ void erf_substep_T (int step, int /*nrk*/, cur_zmom(i,j,k) += wpp; if (l_rayleigh_impl_for_w) { - Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k]; - cur_zmom(i,j,k) /= (one + damping_coeff); + Real damping_coeff = l_damp_coef * dtau * sinesq_stag_d[k]; + cur_zmom(i,j,k) /= (one + damping_coeff); } }); From cd6f00924d7b3f38480403d1e48a769de26526a8 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 14 Aug 2026 16:31:30 -0700 Subject: [PATCH 3/3] Fix imp diff of qv. --- Source/TimeIntegration/ERF_ImplicitPre.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/TimeIntegration/ERF_ImplicitPre.H b/Source/TimeIntegration/ERF_ImplicitPre.H index 3310d88c3f..b13c10fd73 100644 --- a/Source/TimeIntegration/ERF_ImplicitPre.H +++ b/Source/TimeIntegration/ERF_ImplicitPre.H @@ -87,7 +87,7 @@ } if (l_do_implicit_moist && has_moisture) { ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, - slow_dt, l_qv_bc_neumann_vals_d, cell_data, + stage_dt, l_qv_bc_neumann_vals_d, cell_data, stretched_dz_d[level], qfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, @@ -127,7 +127,7 @@ } if (l_do_implicit_moist && has_moisture) { ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, - slow_dt, l_qv_bc_neumann_vals_d, cell_data, + stage_dt, l_qv_bc_neumann_vals_d, cell_data, z_nd_arr, detJ_arr, dxInv, qfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, @@ -166,7 +166,7 @@ } if (l_do_implicit_moist && has_moisture) { ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, - slow_dt, l_qv_bc_neumann_vals_d, cell_data, + stage_dt, l_qv_bc_neumann_vals_d, cell_data, dxInv, qfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac, l_use_mrf_cg);