diff --git a/Source/Diffusion/ERF_Diffusion.H b/Source/Diffusion/ERF_Diffusion.H index 5f2724e507..381a2cdadd 100644 --- a/Source/Diffusion/ERF_Diffusion.H +++ b/Source/Diffusion/ERF_Diffusion.H @@ -482,8 +482,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, @@ -499,8 +500,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, @@ -516,8 +518,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 5156fdd18a..d7782178d3 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_N.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_N.cpp @@ -31,8 +31,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, @@ -46,11 +47,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; @@ -113,7 +111,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) { @@ -122,8 +120,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; } @@ -144,11 +142,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); @@ -175,7 +173,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) } @@ -193,7 +191,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 diff --git a/Source/Diffusion/ERF_ImplicitDiff_S.cpp b/Source/Diffusion/ERF_ImplicitDiff_S.cpp index 5da41244b8..0bd76653df 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_S.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_S.cpp @@ -30,8 +30,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, @@ -45,11 +46,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; @@ -117,7 +115,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) { @@ -125,8 +123,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)); // rhoAlpha*gam is already a flux, so its divergence over cell klo // is scaled by the *cell* spacing, not the face spacing @@ -153,11 +151,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); @@ -193,7 +191,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) } @@ -211,7 +209,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 d27bea0695..bfd24a849b 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_T.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_T.cpp @@ -32,8 +32,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, @@ -49,11 +50,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; @@ -123,7 +121,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) { @@ -131,8 +129,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; } @@ -156,11 +154,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); @@ -190,7 +188,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) } @@ -208,7 +206,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/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 1587147f3e..b13c10fd73 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, - stage_dt, l_bc_neumann_vals_d, cell_data, + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, + stage_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, + 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, + l_use_mrf_cg); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_S<0>(tbx, fine_geom.Domain(), level, stage_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, - stage_dt, l_bc_neumann_vals_d, cell_data, + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, + stage_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, + 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, + l_use_mrf_cg); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_T<0>(tbx, fine_geom.Domain(), level, stage_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, - stage_dt, l_bc_neumann_vals_d, cell_data, + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, + stage_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, + 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); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_N<0>(tbx, fine_geom.Domain(), level, stage_dt, cell_data, rho_u, tau13, tau13_corr, diff --git a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp index 1038798e2d..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]; @@ -123,16 +122,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 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 + 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 ); + Real thm_t_hi = myhalf * (thm_k + thm_kp1); - Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_p) + 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 + RvoRd*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 + RvoRd*qv_km1) ); coeffP_a(i,j,k) = coeff_P; coeffQ_a(i,j,k) = coeff_Q; @@ -144,17 +159,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 { @@ -172,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_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index 07923269de..50abb4e25e 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -215,11 +215,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; } // ************************************************************************* @@ -335,12 +335,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); } }); // Non-EB Anelastic: Per-tile copy of projected momentum (EB done above) @@ -440,15 +440,15 @@ 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; @@ -465,7 +465,7 @@ void erf_slow_rhs_post (int level, int finest_level, // are advanced by the state update below and included in the reflux. // Computing residuals for only the first n_qstate would leave the // rest to be updated with a residual nothing ever wrote. - num_comp = n_qstate_total; + num_comp = n_qstate_total - 1; } else { horiz_adv_type = ac.dryscal_horiz_adv_type; @@ -510,10 +510,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]; } @@ -616,14 +615,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 = n_qstate_total; + if (ivar == RhoQ2_comp) { + num_comp = n_qstate_total - 1; } else if (ivar == RhoScalar_comp) { num_comp = NSCALARS; } diff --git a/Source/TimeIntegration/ERF_SlowRhsPre.cpp b/Source/TimeIntegration/ERF_SlowRhsPre.cpp index 61077abf45..04e5ecb903 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPre.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPre.cpp @@ -144,6 +144,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); @@ -629,6 +633,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, @@ -649,6 +662,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) { @@ -690,6 +715,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, @@ -701,6 +738,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, @@ -711,6 +761,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, @@ -721,6 +782,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); + } } if (use_physical_chamber_wall_flux) { erf_resolved_wall_flux::apply( @@ -737,6 +810,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); @@ -748,9 +824,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 8f90db448d..c8260d08b6 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(); @@ -145,7 +143,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 +156,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 +163,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 +198,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) * (one + RvoRd*qv_cur) + - stage_cons(i,j,k,RhoTheta_comp) * (one + RvoRd*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,7 +258,7 @@ 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); @@ -274,15 +266,6 @@ 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 - - // ********************************************************************* - // Define updates in the RHS of {x, y, z}-momentum equations - // ********************************************************************* { BL_PROFILE("substep_xymom_T"); @@ -295,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); @@ -333,13 +316,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); @@ -391,11 +374,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); @@ -406,7 +389,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); @@ -414,7 +396,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); @@ -422,19 +403,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(); @@ -450,13 +425,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 { @@ -486,57 +491,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); @@ -548,6 +546,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 @@ -556,19 +557,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 * (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); + 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 + 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; + 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 ); @@ -581,8 +623,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); @@ -594,12 +636,15 @@ 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 @@ -691,53 +736,65 @@ 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); } }); // ************************************************************************** - // 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 8cc9b8d087..08df1b2929 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H @@ -209,13 +209,12 @@ if ((solverChoice.vert_implicit_fac[level][nrk] > zero_d) && solverChoice.implicit_before_substep) { const Real stage_dt = static_cast(slow_dt); - /** - * @brief Temporary MultiFab for conserved variable updates. - */ - 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, stage_dt, S_rhs[IntVars::cons], 0, 0, 2, 0); // scratch := S_old + stage_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, stage_dt, S_rhs[IntVars::cons], 0, 0, 2, 0); + MultiFab::Saxpy(scratch, stage_dt, S_rhs[IntVars::cons], RhoQ1_comp, 2, 1, 0); scratch.FillBoundary(geom[level].periodicity()); /** @@ -255,9 +254,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 / stage_dt); // scratch := (S_new - S_old) / stage_dt - MultiFab::Copy(S_rhs[IntVars::cons], scratch, 1, 1, 1, 0); // slow_rhs := (S_new - S_old) / stage_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 / stage_dt); + 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)