From 7849906ab44592a6436ffc145d2b6582dc2389f1 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Mon, 9 Dec 2024 15:30:42 -0800 Subject: [PATCH] Fix thetav for PBL schemes. (#2016) --- Source/BoundaryConditions/ERF_ABLMost.H | 8 ++++-- Source/BoundaryConditions/ERF_ABLMost.cpp | 12 +++++--- Source/DataStructs/ERF_DataStruct.H | 6 ++++ .../ERF_ComputeTurbulentViscosity.cpp | 3 +- .../Diffusion/ERF_DiffusionSrcForState_N.cpp | 3 +- .../Diffusion/ERF_DiffusionSrcForState_T.cpp | 3 +- Source/ERF.cpp | 4 ++- Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp | 3 +- Source/PBL/ERF_ComputeDiffusivityYSU.cpp | 3 +- Source/PBL/ERF_PBLHeight.H | 15 +++++----- Source/PBL/ERF_PBLModels.H | 28 +++++++++++-------- Source/TimeIntegration/ERF_Advance.cpp | 4 ++- Source/Utils/ERF_Thetav.H | 5 +++- 13 files changed, 64 insertions(+), 33 deletions(-) diff --git a/Source/BoundaryConditions/ERF_ABLMost.H b/Source/BoundaryConditions/ERF_ABLMost.H index 177d90c5a..f652b3d69 100644 --- a/Source/BoundaryConditions/ERF_ABLMost.H +++ b/Source/BoundaryConditions/ERF_ABLMost.H @@ -368,7 +368,9 @@ public: update_pblh (const int& lev, amrex::Vector>& vars, amrex::MultiFab* z_phys_cc, - const int RhoQv_comp, const int RhoQr_comp); + const int RhoQv_comp, + const int RhoQc_comp, + const int RhoQr_comp); template void @@ -376,7 +378,9 @@ public: amrex::Vector>& vars, amrex::MultiFab* z_phys_cc, const PBLHeightEstimator& est, - const int RhoQv_comp, const int RhoQr_comp); + const int RhoQv_comp, + const int RhoQc_comp, + const int RhoQr_comp); void read_custom_roughness (const int& lev, diff --git a/Source/BoundaryConditions/ERF_ABLMost.cpp b/Source/BoundaryConditions/ERF_ABLMost.cpp index f95ba33d4..65b3a43b4 100644 --- a/Source/BoundaryConditions/ERF_ABLMost.cpp +++ b/Source/BoundaryConditions/ERF_ABLMost.cpp @@ -561,11 +561,13 @@ void ABLMost::update_pblh (const int& lev, Vector>& vars, MultiFab* z_phys_cc, - int RhoQv_comp, int RhoQr_comp) + int RhoQv_comp, + int RhoQc_comp, + int RhoQr_comp) { if (pblh_type == PBLHeightCalcType::MYNN25) { MYNNPBLH estimator; - compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQr_comp); + compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQc_comp, RhoQr_comp); } else if (pblh_type == PBLHeightCalcType::YSU) { amrex::Error("YSU PBLH calc not implemented yet"); } @@ -577,11 +579,13 @@ ABLMost::compute_pblh (const int& lev, Vector>& vars, MultiFab* z_phys_cc, const PBLHeightEstimator& est, - int RhoQv_comp, int RhoQr_comp) + int RhoQv_comp, + int RhoQc_comp, + int RhoQr_comp) { est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(), vars[lev][Vars::cons],m_lmask_lev[lev][0], - RhoQv_comp, RhoQr_comp); + RhoQv_comp, RhoQc_comp, RhoQr_comp); } void diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 872025b54..3e7349b88 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -107,17 +107,22 @@ struct SolverChoice { pp.query_enum_case_insensitive("moisture_model",moisture_type); if (moisture_type == MoistureType::SAM) { RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; RhoQr_comp = RhoQ4_comp; } else if (moisture_type == MoistureType::SAM_NoIce) { RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; RhoQr_comp = RhoQ4_comp; } else if (moisture_type == MoistureType::SAM_NoPrecip_NoIce) { RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; } else if (moisture_type == MoistureType::Kessler) { RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; RhoQr_comp = RhoQ3_comp; } else if (moisture_type == MoistureType::Kessler_NoRain) { RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; } // TODO: should we set default for dry?? @@ -637,6 +642,7 @@ struct SolverChoice { bool do_precip {true}; bool use_moist_background {false}; int RhoQv_comp {-1}; + int RhoQc_comp {-1}; // This component will be model-dependent: // if a model with no rain, this will stay -1 diff --git a/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp b/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp index a18337a49..0a20802b7 100644 --- a/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp @@ -489,7 +489,8 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel , ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity, geom, turbChoice, most, use_moisture, level, bc_ptr, vert_only, z_phys_nd, - solverChoice.RhoQv_comp, solverChoice.RhoQr_comp); + solverChoice.RhoQv_comp, solverChoice.RhoQc_comp, + solverChoice.RhoQr_comp); } else if (turbChoice.pbl_type == PBLType::YSU) { ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity, geom, turbChoice, most, use_moisture, diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp index e4429425b..5fb22b6dd 100644 --- a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp @@ -713,6 +713,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1; const int rhoqv_comp = solverChoice.RhoQv_comp; + const int rhoqc_comp = solverChoice.RhoQc_comp; const int rhoqr_comp = solverChoice.RhoQr_comp; ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -728,7 +729,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim, mu_turb,cellSizeInv,domain, pbl_mynn_B1_l,tm_arr(i,j,0), - rhoqv_comp, rhoqr_comp, + rhoqv_comp, rhoqc_comp, rhoqr_comp, c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp index 8b4e1bcc5..b588f8117 100644 --- a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp @@ -824,6 +824,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1; const int rhoqv_comp = solverChoice.RhoQv_comp; + const int rhoqc_comp = solverChoice.RhoQc_comp; const int rhoqr_comp = solverChoice.RhoQr_comp; ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -839,7 +840,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim, mu_turb,cellSizeInv,domain, pbl_mynn_B1_l,tm_arr(i,j,0), - rhoqv_comp, rhoqr_comp, + rhoqv_comp, rhoqc_comp, rhoqr_comp, c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 306eaf40f..3c759adcf 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1061,7 +1061,9 @@ ERF::InitData_post () // we don't want to call update_fluxes multiple times because // it will change u* and theta* from their previous values m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get(), - solverChoice.RhoQv_comp, solverChoice.RhoQr_comp); + solverChoice.RhoQv_comp, + solverChoice.RhoQc_comp, + solverChoice.RhoQr_comp); m_most->update_fluxes(lev, time); } } diff --git a/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp b/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp index f1608f078..b8024729d 100644 --- a/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp +++ b/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp @@ -21,6 +21,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, bool /*vert_only*/, const std::unique_ptr& z_phys_nd, const int RhoQv_comp, + const int RhoQc_comp, const int RhoQr_comp) { const bool use_terrain = (z_phys_nd != nullptr); @@ -140,7 +141,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz, dudz, dvdz, - RhoQv_comp, RhoQr_comp, use_most); + RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most); // Spatially varying MOST Real theta0 = tm_arr(i,j,0); diff --git a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp index f197fdf33..36ff5fedc 100644 --- a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp +++ b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp @@ -204,13 +204,14 @@ ComputeDiffusivityYSU (const MultiFab& xvel, constexpr Real prandtl_max = 4.0; Real dthetadz, dudz, dvdz; const int RhoQv_comp = -1; + const int RhoQc_comp = -1; const int RhoQr_comp = -1; ComputeVerticalDerivativesPBL(i, j, k, uvel, vvel, cell_data, izmin, izmax, 1.0/dz_terrain, c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, - dthetadz, dudz, dvdz, RhoQv_comp, RhoQr_comp); + dthetadz, dudz, dvdz, RhoQv_comp, RhoQc_comp, RhoQr_comp); const Real shear_squared = dudz*dudz + dvdz*dvdz + 1.0e-9; // 1.0e-9 from WRF to avoid divide by zero const Real theta = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp); Real richardson = CONST_GRAV / theta * dthetadz / shear_squared; diff --git a/Source/PBL/ERF_PBLHeight.H b/Source/PBL/ERF_PBLHeight.H index 1ca84631f..c5dbcce34 100644 --- a/Source/PBL/ERF_PBLHeight.H +++ b/Source/PBL/ERF_PBLHeight.H @@ -25,6 +25,7 @@ struct MYNNPBLH { const amrex::MultiFab& cons, const amrex::iMultiFab* lmask, const int RhoQv_comp, + const int RhoQc_comp, const int RhoQr_comp) const { #if 0 @@ -35,7 +36,7 @@ struct MYNNPBLH { auto thetav_min = amrex::ReduceToPlane(dir, bxlow, cons, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real { - return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQr_comp); + return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQc_comp,RhoQr_comp); }); #endif @@ -88,7 +89,7 @@ struct MYNNPBLH { int jj = amrex::max(amrex::min(j,jmax),jmin); if (zphys_arr(ii,jj,k) < thetamin_height) { - amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQr_comp); + amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp); if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv; } }); @@ -104,8 +105,8 @@ struct MYNNPBLH { // // Find PBL height based on thetav increase (best for CBLs) // - amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQr_comp); - amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQr_comp); + amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp); + amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp); int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1; if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land) @@ -164,7 +165,7 @@ struct MYNNPBLH { ParallelFor(gtbxlow, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQr_comp); + amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp); if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv; }); @@ -176,8 +177,8 @@ struct MYNNPBLH { // // Find PBL height based on thetav increase (best for CBLs) // - amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQr_comp); - amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQr_comp); + amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp); + amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp); int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1; if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land) diff --git a/Source/PBL/ERF_PBLModels.H b/Source/PBL/ERF_PBLModels.H index 9c0ee5086..41965ba6c 100644 --- a/Source/PBL/ERF_PBLModels.H +++ b/Source/PBL/ERF_PBLModels.H @@ -33,7 +33,9 @@ ComputeDiffusivityMYNN25 (const amrex::MultiFab& xvel, const amrex::BCRec* bc_ptr, bool /*vert_only*/, const std::unique_ptr& z_phys_nd, - const int RhoQv_comp, const int RhoQr_comp); + const int RhoQv_comp, + const int RhoQc_comp, + const int RhoQr_comp); /** * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the @@ -94,23 +96,24 @@ ComputeVerticalDerivativesPBL (int i, int j, int k, amrex::Real& dudz, amrex::Real& dvdz, const int RhoQv_comp, + const int RhoQc_comp, const int RhoQr_comp, const bool use_most=true) { if ( k==izmax && c_ext_dir_on_zhi ) { - dthetadz = (1.0/3.0)*( -Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) - - 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQr_comp) - + 4.0 * Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv; + dthetadz = (1.0/3.0)*( -Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) + - 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) + + 4.0 * Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv; } else if ( k==izmin && c_ext_dir_on_zlo ) { - dthetadz = (1.0/3.0)*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp) - + 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQr_comp) - - 4.0 * Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv; + dthetadz = (1.0/3.0)*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) + + 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) + - 4.0 * Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv; } else if ( k==izmin && use_most ) { - dthetadz = ( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp) - - Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv; + dthetadz = ( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) + - Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv; } else { - dthetadz = 0.5*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp) - - Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv; + dthetadz = 0.5*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) + - Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv; } if ( k==izmax && u_ext_dir_on_zhi ) { @@ -165,6 +168,7 @@ ComputeQKESourceTerms (int i, int j, int k, amrex::Real pbl_mynn_B1_l, const amrex::Real theta_mean, const int RhoQv_comp, + const int RhoQc_comp, const int RhoQr_comp, bool c_ext_dir_on_zlo, bool c_ext_dir_on_zhi, @@ -195,7 +199,7 @@ ComputeQKESourceTerms (int i, int j, int k, u_ext_dir_on_zlo, u_ext_dir_on_zhi, v_ext_dir_on_zlo, v_ext_dir_on_zhi, dthetadz, dudz, dvdz, - RhoQv_comp, RhoQr_comp, use_most); + RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most); // Note: Transport terms due to turbulence and pressure are included when // DiffusionSrcForState_* is called from ERF_slow_rhs_post. diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 6f99b3a0b..52abe93be 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -110,7 +110,9 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) // Reassign the field ptrs for MAC avg computation. m_most->update_mac_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim); m_most->update_pblh(lev, vars_old, z_phys_cc[lev].get(), - solverChoice.RhoQv_comp, solverChoice.RhoQr_comp); + solverChoice.RhoQv_comp, + solverChoice.RhoQc_comp, + solverChoice.RhoQr_comp); m_most->update_fluxes(lev, time); } } diff --git a/Source/Utils/ERF_Thetav.H b/Source/Utils/ERF_Thetav.H index 36d63cc50..cb2efc9a1 100644 --- a/Source/Utils/ERF_Thetav.H +++ b/Source/Utils/ERF_Thetav.H @@ -10,15 +10,18 @@ amrex::Real Thetav (int i, int j, int k, const amrex::Array4& cell_data, const int RhoQv_comp, + const int RhoQc_comp, const int RhoQr_comp) { amrex::Real thetav = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp); if (RhoQr_comp > 0) { thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQv_comp) / cell_data(i,j,k,Rho_comp) + - cell_data(i,j,k,RhoQc_comp) / cell_data(i,j,k,Rho_comp) - cell_data(i,j,k,RhoQr_comp) / cell_data(i,j,k,Rho_comp)); } else if (RhoQv_comp > 0) { - thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQ1_comp) / cell_data(i,j,k,Rho_comp)); + thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQ1_comp) / cell_data(i,j,k,Rho_comp) + - cell_data(i,j,k,RhoQc_comp) / cell_data(i,j,k,Rho_comp)); } return thetav;