diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_realbdy.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_realbdy.cpp index 5d6be8b7d..65cf8f039 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_realbdy.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_realbdy.cpp @@ -40,7 +40,7 @@ ERF::fill_from_realbdy (const Vector& mfs, 0, 0, 0, 0, 0}; Vector cons_map = {Rho_comp, RealBdyVars::T, RhoKE_comp, - RhoQKE_comp, RhoScalar_comp, RealBdyVars::QV, + RhoScalar_comp, RealBdyVars::QV, RhoQ2_comp, RhoQ3_comp, RhoQ4_comp, RhoQ5_comp, RhoQ6_comp}; diff --git a/Source/DataStructs/ERF_TurbStruct.H b/Source/DataStructs/ERF_TurbStruct.H index a8a10f72a..1468f26b7 100644 --- a/Source/DataStructs/ERF_TurbStruct.H +++ b/Source/DataStructs/ERF_TurbStruct.H @@ -76,13 +76,6 @@ struct TurbChoice { } } - // Right now, solving the QKE equation is only supported when MYNN PBL is turned on - if (pbl_type == PBLType::MYNN25) { - use_QKE = true; - query_one_or_per_level(pp, "advect_QKE" , advect_QKE, lev, max_level); - query_one_or_per_level(pp, "diffuse_QKE_3D", diffuse_QKE_3D, lev, max_level); - } - // LES constants... query_one_or_per_level(pp, "Cs" ,Cs, lev, max_level); query_one_or_per_level(pp, "CI" ,CI, lev, max_level); @@ -201,10 +194,5 @@ struct TurbChoice { bool pbl_ysu_force_over_water = false; // Force YSU to act as if it is over water regardless of other inputs (for testing) amrex::Real pbl_ysu_land_Ribcr = 0.25; // Critical Bulk Richardson number of Land for stable conditions amrex::Real pbl_ysu_unst_Ribcr = 0.0; // Critical Bulk Richardson number for unstable conditions - - // QKE stuff - default is not to use it, if MYNN2.5 PBL is used default is turb transport in Z-direction only - bool use_QKE = false; - bool diffuse_QKE_3D = false; - bool advect_QKE = true; }; #endif diff --git a/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp b/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp index 5e9bc6e5c..d75ffb663 100644 --- a/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp @@ -84,7 +84,8 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most,exp_most); + bool only_pbl = false; + Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most,exp_most,only_pbl); Real dxInv = cellSizeInv[0]; Real dyInv = cellSizeInv[1]; Real dzInv = cellSizeInv[2]; @@ -207,6 +208,8 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, mu_turb(i,j,k,EddyDiff::Mom_v) = mu_turb(i,j,k,EddyDiff::Mom_h); // KH = (1 + 2*l/delta) * mu_turb mu_turb(i,j,k,EddyDiff::Theta_v) = (1.+2.*length/DeltaMsf) * mu_turb(i,j,k,EddyDiff::Mom_v); + // Store lengthscale for TKE source terms + mu_turb(i,j,k,EddyDiff::Turb_lengthscale) = length; // Calculate SFS quantities // - dissipation @@ -217,6 +220,7 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf; } diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length; + // - heat flux // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will // be overwritten when BCs are applied) @@ -230,14 +234,14 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0) //*********************************************************************************** int ngc(1); - // EddyDiff mapping : Theta_h KE_h QKE_h Scalar_h Q_h - Vector Factors = {inv_Pr_t, inv_sigma_k, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr + // EddyDiff mapping : Theta_h KE_h Scalar_h Q_h + Vector Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr Gpu::AsyncVector d_Factors; d_Factors.resize(Factors.size()); Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin()); Real* fac_ptr = d_Factors.data(); - bool use_KE = (turbChoice.les_type == LESType::Deardorff); - bool use_QKE = turbChoice.use_QKE; + const bool use_KE = ( (turbChoice.les_type == LESType::Deardorff) || + (turbChoice.pbl_type == PBLType::MYNN25) ); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -324,18 +328,6 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, int offset = (EddyDiff::NumDiffs-1)/2; switch (n) { - case EddyDiff::QKE_h: - // Populate element other than mom_h/v on the whole grid - if(use_QKE) { - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int indx = n; - int indx_v = indx + offset; - mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; - mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); - }); - } - break; case EddyDiff::KE_h: if (use_KE) { ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -389,20 +381,6 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, int offset = (EddyDiff::NumDiffs-1)/2; switch (n) { - case EddyDiff::QKE_h: - // Extrap all components at top & bottom - if(use_QKE) { - ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int indx = n; - int indx_v = indx + offset; - mu_turb(i, j, k_lo-k, indx ) = mu_turb(i, j, k_lo, indx ); - mu_turb(i, j, k_hi+k, indx ) = mu_turb(i, j, k_hi, indx ); - mu_turb(i, j, k_lo-k, indx_v) = mu_turb(i, j, k_lo, indx_v); - mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); - }); - } - break; case EddyDiff::KE_h: if (use_KE) { ParallelFor(planez, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -509,7 +487,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel , } if (turbChoice.pbl_type == PBLType::MYNN25) { - ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity, + ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity, Diss, geom, turbChoice, most, use_moisture, level, bc_ptr, vert_only, z_phys_nd, solverChoice.RhoQv_comp, solverChoice.RhoQr_comp); diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp index 58fde87ac..10b27ea95 100644 --- a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp @@ -79,11 +79,11 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, const auto& dom_lo = lbound(domain); const auto& dom_hi = ubound(domain); - bool l_use_QKE = turbChoice.use_QKE; - bool l_use_deardorff = (turbChoice.les_type == LESType::Deardorff); Real l_inv_theta0 = 1.0 / turbChoice.theta_ref; Real l_abs_g = std::abs(grav_gpu[2]); + bool l_use_KE = ( (turbChoice.les_type == LESType::Deardorff) || + (turbChoice.pbl_type == PBLType::MYNN) ); bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha); bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) || (turbChoice.les_type == LESType::Deardorff ) || @@ -96,7 +96,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, const int end_comp = start_comp + num_comp - 1; - // Theta, KE, QKE, Scalar + // Theta, KE, Scalar Vector alpha_eff(NPRIMVAR_max, 0.0); if (l_consA) { for (int i = 0; i < NPRIMVAR_max; ++i) { @@ -164,15 +164,15 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, } } - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, - EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v , - EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v }; + Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h, + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; + Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h, + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; + Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::Scalar_v, + EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v , + EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v }; // Device vectors Gpu::AsyncVector alpha_eff_d; @@ -678,7 +678,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, // The surface heat flux hfx_z(i,j,-1) is updated in MOSTStress at // each RK stage if using the ERF_EXPLICIT_MOST_STRESS path, but that // does not change the buoyancy production term here. - if (l_use_deardorff && start_comp <= RhoKE_comp && end_comp >=RhoKE_comp) { + if (l_use_KE && (start_comp <= RhoKE_comp) && (end_comp >=RhoKE_comp)) { int qty_index = RhoKE_comp; ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -705,34 +705,4 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, cell_rhs(i,j,k,qty_index) -= diss(i,j,k); }); } - - // Using PBL - if (l_use_QKE && start_comp <= RhoQKE_comp && end_comp >=RhoQKE_comp) { - int qty_index = RhoQKE_comp; - auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1; - - const int rhoqv_comp = solverChoice.RhoQv_comp; - const int rhoqr_comp = solverChoice.RhoQr_comp; - - ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) ); - bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) ); - bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ); - bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) ); - bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ); - bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) ); - - // This computes shear production, buoyancy production, and dissipation terms only. - 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, - 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, - use_most); - }); - } - } diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp index d210b512c..eb8ea8743 100644 --- a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp @@ -91,11 +91,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, const auto& dom_hi = ubound(domain); - bool l_use_QKE = turbChoice.use_QKE; - bool l_use_deardorff = (turbChoice.les_type == LESType::Deardorff); Real l_inv_theta0 = 1.0 / turbChoice.theta_ref; Real l_abs_g = std::abs(grav_gpu[2]); + bool l_use_KE = ( (turbChoice.les_type == LESType::Deardorff) || + (turbChoice.pbl_type == PBLType::MYNN) ); bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha); bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) || (turbChoice.les_type == LESType::Deardorff ) || @@ -109,7 +109,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, const int end_comp = start_comp + num_comp - 1; - // Theta, KE, QKE, Scalar + // Theta, KE, Scalar Vector alpha_eff(NPRIMVAR_max, 0.0); if (l_consA) { for (int i = 0; i < NPRIMVAR_max; ++i) { @@ -177,15 +177,15 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, } } - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , - EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, - EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v , - EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v }; + Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h, + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; + Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::Scalar_h, + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h , + EddyDiff::Q_h , EddyDiff::Q_h, EddyDiff::Q_h }; + Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::Scalar_v, + EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v , + EddyDiff::Q_v , EddyDiff::Q_v, EddyDiff::Q_v }; // Device vectors Gpu::AsyncVector alpha_eff_d; @@ -789,7 +789,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, // The surface heat flux hfx_z(i,j,-1) is updated in MOSTStress at // each RK stage if using the ERF_EXPLICIT_MOST_STRESS path, but that // does not change the buoyancy production term here. - if (l_use_deardorff && start_comp <= RhoKE_comp && end_comp >=RhoKE_comp) { + if (l_use_KE && (start_comp <= RhoKE_comp) && (end_comp >= RhoKE_comp)) { int qty_index = RhoKE_comp; ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -816,34 +816,4 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, cell_rhs(i,j,k,qty_index) -= diss(i,j,k); }); } - - // Using PBL - if (l_use_QKE && start_comp <= RhoQKE_comp && end_comp >=RhoQKE_comp) { - int qty_index = RhoQKE_comp; - auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1; - - const int rhoqv_comp = solverChoice.RhoQv_comp; - const int rhoqr_comp = solverChoice.RhoQr_comp; - - ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) ); - bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) ); - bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ); - bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) ); - bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ); - bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) ); - - const Real met_h_zeta = detJ(i,j,k); - // This computes shear production, buoyancy production, and dissipation terms only. - 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, - 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, - use_most, met_h_zeta); - }); - } } diff --git a/Source/Diffusion/ERF_EddyViscosity.H b/Source/Diffusion/ERF_EddyViscosity.H index cc47191dc..008a39d20 100644 --- a/Source/Diffusion/ERF_EddyViscosity.H +++ b/Source/Diffusion/ERF_EddyViscosity.H @@ -37,14 +37,9 @@ ComputeSmnSmn (int& i, int& j, int& k, const amrex::Array4& tau23, const int& klo, const bool& use_most, - const bool& exp_most) + const bool& exp_most, + const bool& only_pbl) { - amrex::Real s11bar = tau11(i,j,k); - amrex::Real s22bar = tau22(i,j,k); - amrex::Real s33bar = tau33(i,j,k); - amrex::Real s12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k ) - + tau12(i+1, j , k ) + tau12(i+1, j+1, k ) ); - // NOTES: // - If ERF_EXPLICIT_MOST_STRESS is not used, then we do not use the // strains lying on the bottom boundary with MOST. These values are @@ -76,6 +71,15 @@ ComputeSmnSmn (int& i, int& j, int& k, + tau23(i , j+1, k ) + tau23(i , j+1, k+1) ); } + amrex::Real s11bar(0.), s22bar(0.), s33bar(0.), s12bar(0.); + if (!use_pbl) { + s11bar = tau11(i,j,k); + s22bar = tau22(i,j,k); + s33bar = tau33(i,j,k); + s12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k ) + + tau12(i+1, j , k ) + tau12(i+1, j+1, k ) ); + } + amrex::Real SmnSmn = s11bar*s11bar + s22bar*s22bar + s33bar*s33bar + 2.0*s12bar*s12bar + 2.0*s13bar*s13bar + 2.0*s23bar*s23bar; diff --git a/Source/ERF.H b/Source/ERF.H index a977fa35a..e10ba5edf 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -887,12 +887,12 @@ private: amrex::Vector plot_var_names_1; amrex::Vector plot_var_names_2; - const amrex::Vector cons_names {"density", "rhotheta", "rhoKE", "rhoQKE", "rhoadv_0", + const amrex::Vector cons_names {"density", "rhotheta", "rhoKE", "rhoadv_0", "rhoQ1", "rhoQ2", "rhoQ3", "rhoQ4", "rhoQ5", "rhoQ6"}; // Note that the order of variable names here must match the order in ERF_Derive.cpp - const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", + const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "scalar", "vorticity_x","vorticity_y","vorticity_z", "magvel", "divU", "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", diff --git a/Source/ERF.cpp b/Source/ERF.cpp index ef6ec5bfc..d69f9cf27 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -675,24 +675,6 @@ ERF::InitData_post () } // lev } - Real QKE_0; - if (pp.query("QKE_0", QKE_0)) { - Print() << "Initializing uniform QKE=" << QKE_0 << std::endl; - for (int lev = 0; lev <= finest_level; lev++) { - auto& lev_new = vars_new[lev]; - for (MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box &bx = mfi.tilebox(); - const auto &cons_arr = lev_new[Vars::cons].array(mfi); - // We want to set the lateral BC values, too - Box gbx = bx; // Copy constructor - gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cons_arr(i,j,k,RhoQKE_comp) = cons_arr(i,j,k,Rho_comp) * QKE_0; - }); - } // mfi - } - } - if (solverChoice.coupling_type == CouplingType::TwoWay) { AverageDown(); } diff --git a/Source/ERF_Derive.H b/Source/ERF_Derive.H index 0090ef07b..7eb56aa5c 100644 --- a/Source/ERF_Derive.H +++ b/Source/ERF_Derive.H @@ -89,17 +89,6 @@ void erf_derKE ( const int* bcrec, const int level); -void erf_derQKE ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int dcomp, - int ncomp, - const amrex::FArrayBox& datfab, - const amrex::Geometry& geomdata, - amrex::Real time, - const int* bcrec, - const int level); - void erf_dervortx ( const amrex::Box& bx, amrex::FArrayBox& derfab, diff --git a/Source/ERF_Derive.cpp b/Source/ERF_Derive.cpp index 75ff631ba..33fa8f118 100644 --- a/Source/ERF_Derive.cpp +++ b/Source/ERF_Derive.cpp @@ -196,27 +196,6 @@ erf_derKE (const Box& bx, erf_derrhodivide(bx, derfab, datfab, RhoKE_comp); } -/** - * Function to define QKE by dividing (rho QKE) by rho - * - * @params[in] bx box on which to divide by density - * @params[out] derfab array of derived quantity -- here it holds QKE - * @params[in] datfab array of data used to construct derived quantity -*/ -void -erf_derQKE (const Box& bx, - FArrayBox& derfab, - int /*dcomp*/, - int /*ncomp*/, - const FArrayBox& datfab, - const Geometry& /*geomdata*/, - Real /*time*/, - const int* /*bcrec*/, - const int /*level*/) -{ - erf_derrhodivide(bx, derfab, datfab, RhoQKE_comp); -} - void erf_dervortx ( const amrex::Box& bx, diff --git a/Source/ERF_IndexDefines.H b/Source/ERF_IndexDefines.H index 811946464..499380520 100644 --- a/Source/ERF_IndexDefines.H +++ b/Source/ERF_IndexDefines.H @@ -9,8 +9,8 @@ */ // This defines the ACTUAL number of non-moisture vars = -// rho, rhotheta, rhoKE, rhoQKE -#define NDRY 4 +// rho, rhotheta, rhoKE +#define NDRY 3 // This defines the ACTUAL number of non-moisture scalar vars #define NSCALARS 1 @@ -35,10 +35,9 @@ // Cell-centered state variables #define Rho_comp 0 #define RhoTheta_comp (Rho_comp+1) -#define RhoKE_comp (Rho_comp+2) // for Deardorff LES Model -#define RhoQKE_comp (Rho_comp+3) // for MYNN or YSU PBL Model +#define RhoKE_comp (Rho_comp+2) // for Deardorff LES Model or MYNN PBL Model -#define RhoScalar_comp (RhoQKE_comp+1) +#define RhoScalar_comp (RhoKE_comp+1) #define RhoQ1_comp (RhoScalar_comp+NSCALARS) #define RhoQ2_comp (RhoQ1_comp+1) @@ -50,7 +49,6 @@ // Cell-centered primitive variables #define PrimTheta_comp (RhoTheta_comp -1) #define PrimKE_comp (RhoKE_comp -1) -#define PrimQKE_comp (RhoQKE_comp -1) #define PrimScalar_comp (RhoScalar_comp-1) #define PrimQ1_comp (RhoQ1_comp-1) #define PrimQ2_comp (RhoQ2_comp-1) @@ -68,7 +66,6 @@ namespace BCVars { Rho_bc_comp = 0, RhoTheta_bc_comp, RhoKE_bc_comp, - RhoQKE_bc_comp, RhoScalar_bc_comp, RhoQ1_bc_comp, RhoQ2_bc_comp, @@ -143,16 +140,14 @@ namespace EddyDiff { Mom_h = 0, Theta_h, KE_h, - QKE_h, Scalar_h, Q_h, Mom_v, Theta_v, KE_v, - QKE_v, Scalar_v, Q_v, - PBL_lengthscale, + Turb_lengthscale, NumDiffs }; } diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 75557acd4..402e9fa47 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -331,7 +331,6 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } calculate_derived("theta", vars_new[lev][Vars::cons], derived::erf_dertheta); calculate_derived("KE", vars_new[lev][Vars::cons], derived::erf_derKE); - calculate_derived("QKE", vars_new[lev][Vars::cons], derived::erf_derQKE); calculate_derived("scalar", vars_new[lev][Vars::cons], derived::erf_derscalar); calculate_derived("vorticity_x", mf_cc_vel[lev] , derived::erf_dervortx); calculate_derived("vorticity_y", mf_cc_vel[lev] , derived::erf_dervorty); diff --git a/Source/IO/ERF_ReadBndryPlanes.H b/Source/IO/ERF_ReadBndryPlanes.H index 8435daedc..e6dd27e5e 100644 --- a/Source/IO/ERF_ReadBndryPlanes.H +++ b/Source/IO/ERF_ReadBndryPlanes.H @@ -45,7 +45,6 @@ public: [[nodiscard]] int ingested_q1() const {return is_q1_read;} [[nodiscard]] int ingested_q2() const {return is_q2_read;} [[nodiscard]] int ingested_KE() const {return is_KE_read;} - [[nodiscard]] int ingested_QKE() const {return is_QKE_read;} private: @@ -101,7 +100,6 @@ private: int is_q1_read; int is_q2_read; int is_KE_read; - int is_QKE_read; int last_file_read; }; diff --git a/Source/IO/ERF_ReadBndryPlanes.cpp b/Source/IO/ERF_ReadBndryPlanes.cpp index b95710f38..fff5d425c 100644 --- a/Source/IO/ERF_ReadBndryPlanes.cpp +++ b/Source/IO/ERF_ReadBndryPlanes.cpp @@ -159,7 +159,6 @@ ReadBndryPlanes::ReadBndryPlanes (const Geometry& geom, const Real& rdOcp_in) is_q1_read = 0; is_q2_read = 0; is_KE_read = 0; - is_QKE_read = 0; if (pp.contains("bndry_input_var_names")) { @@ -175,7 +174,6 @@ ReadBndryPlanes::ReadBndryPlanes (const Geometry& geom, const Real& rdOcp_in) if (m_var_names[i] == "qv") is_q1_read = 1; if (m_var_names[i] == "qc") is_q2_read = 1; if (m_var_names[i] == "ke") is_KE_read = 1; - if (m_var_names[i] == "qke") is_QKE_read = 1; } } @@ -415,7 +413,6 @@ void ReadBndryPlanes::read_file (const int idx, if (var_name == "theta") n_offset = BCVars::RhoTheta_bc_comp; if (var_name == "temperature") n_offset = BCVars::RhoTheta_bc_comp; if (var_name == "ke") n_offset = BCVars::RhoKE_bc_comp; - if (var_name == "qke") n_offset = BCVars::RhoQKE_bc_comp; if (var_name == "scalar") n_offset = BCVars::RhoScalar_bc_comp; if (var_name == "qv") n_offset = BCVars::RhoQ1_bc_comp; if (var_name == "qc") n_offset = BCVars::RhoQ2_bc_comp; @@ -480,7 +477,7 @@ void ReadBndryPlanes::read_file (const int idx, bndry_mf_arr(i, j, k, 0) = 0.5 * (R1*Th1 + R2*Th2); }); } else if (var_name == "scalar" || var_name == "qv" || var_name == "qc" || - var_name == "ke" || var_name == "qke") { + var_name == "ke") { ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Real R1 = bndry_read_r_arr(i, j, k, 0); @@ -510,7 +507,7 @@ void ReadBndryPlanes::read_file (const int idx, bndry_mf_arr(i, j, k, 0) = 0.5 * (R1*Th1 + R2*Th2); }); } else if (var_name == "scalar" || var_name == "qv" || var_name == "qc" || - var_name == "ke" || var_name == "qke") { + var_name == "ke") { ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Real R1 = l_bc_extdir_vals_d[BCVars::Rho_bc_comp][ori]; diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index 450515a7a..e7a33ff6f 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -210,8 +210,8 @@ ERF::derive_diag_profiles(Real /*time*/, bool l_use_kturb = ((solverChoice.turbChoice[lev].les_type != LESType::None) || (solverChoice.turbChoice[lev].pbl_type != PBLType::None)); - bool l_use_KE = (solverChoice.turbChoice[lev].les_type == LESType::Deardorff); - bool l_use_QKE = solverChoice.turbChoice[lev].use_QKE; + bool l_use_KE = ( (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) || + (solverChoice.turbChoice[lev].pbl_type == PBLType::MYNN) ); // This will hold rho, theta, ksgs, Kmh, Kmv, uu, uv, uw, vv, vw, ww, uth, vth, wth, // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 @@ -287,8 +287,6 @@ ERF::derive_diag_profiles(Real /*time*/, Real ksgs = 0.0; if (l_use_KE) { ksgs = cons_arr(i,j,k,RhoKE_comp) / cons_arr(i,j,k,Rho_comp); - } else if (l_use_QKE) { - ksgs = cons_arr(i,j,k,RhoQKE_comp) / cons_arr(i,j,k,Rho_comp); } fab_arr(i, j, k, 2) = ksgs; #if 1 diff --git a/Source/IO/ERF_Write1DProfiles_stag.cpp b/Source/IO/ERF_Write1DProfiles_stag.cpp index 81fa117f5..3f06b2199 100644 --- a/Source/IO/ERF_Write1DProfiles_stag.cpp +++ b/Source/IO/ERF_Write1DProfiles_stag.cpp @@ -316,8 +316,8 @@ ERF::derive_diag_profiles_stag (Real /*time*/, bool l_use_kturb = ((solverChoice.turbChoice[lev].les_type != LESType::None) || (solverChoice.turbChoice[lev].pbl_type != PBLType::None)); - bool l_use_KE = (solverChoice.turbChoice[lev].les_type == LESType::Deardorff); - bool l_use_QKE = solverChoice.turbChoice[lev].use_QKE; + bool l_use_KE = ( (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) || + (solverChoice.turbChoice[lev].pbl_type == PBLType::MYNN) ); // Note: "uiui" == u_i*u_i = u*u + v*v + w*w // This will hold rho, theta, ksgs, Kmh, Kmv, uu, uv, vv, uth, vth, @@ -375,8 +375,6 @@ ERF::derive_diag_profiles_stag (Real /*time*/, Real ksgs = 0.0; if (l_use_KE) { ksgs = cons_arr(i,j,k,RhoKE_comp) / cons_arr(i,j,k,Rho_comp); - } else if (l_use_QKE) { - ksgs = cons_arr(i,j,k,RhoQKE_comp) / cons_arr(i,j,k,Rho_comp); } fab_arr(i, j, k, 2) = ksgs; if (l_use_kturb) { diff --git a/Source/IO/ERF_WriteBndryPlanes.cpp b/Source/IO/ERF_WriteBndryPlanes.cpp index 5ddc4c458..def6541e5 100644 --- a/Source/IO/ERF_WriteBndryPlanes.cpp +++ b/Source/IO/ERF_WriteBndryPlanes.cpp @@ -196,17 +196,6 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time, derived::erf_derrhodivide(bx, Temp[mfi], S[mfi], RhoKE_comp); } bndry.copyFrom(Temp, nghost, 0, 0, ncomp, m_geom[bndry_lev].periodicity()); - - } else if (var_name == "qke") { - - MultiFab Temp(S.boxArray(),S.DistributionMap(),ncomp,0); - for (MFIter mfi(Temp, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - derived::erf_derrhodivide(bx, Temp[mfi], S[mfi], RhoQKE_comp); - } - bndry.copyFrom(Temp, nghost, 0, 0, ncomp, m_geom[bndry_lev].periodicity()); - } else if (var_name == "qv") { if (S.nComp() > RhoQ2_comp) { MultiFab Temp(S.boxArray(),S.DistributionMap(),ncomp,0); diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp index f4e05cfdb..4f99b069b 100644 --- a/Source/Initialization/ERF_init_bcs.cpp +++ b/Source/Initialization/ERF_init_bcs.cpp @@ -30,7 +30,6 @@ void ERF::init_bcs () m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] = -1.0; // It is important to set this negative // because the sign is tested on below m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori] = 0.0; - m_bc_extdir_vals[BCVars::RhoQKE_bc_comp][ori] = 0.0; m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = 0.0; m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = 0.0; m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = 0.0; @@ -48,7 +47,6 @@ void ERF::init_bcs () m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::RhoKE_bc_comp][ori] = 0.0; - m_bc_neumann_vals[BCVars::RhoQKE_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::RhoScalar_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::RhoQ1_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::RhoQ2_bc_comp][ori] = 0.0; @@ -172,14 +170,6 @@ void ERF::init_bcs () if (pp.query("KE", KE_in)) m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori] = rho_in*KE_in; } - Real QKE_in = 0.; - if (input_bndry_planes && m_r2d->ingested_QKE()) { - m_bc_extdir_vals[BCVars::RhoQKE_bc_comp][ori] = 0.; - } else { - if (pp.query("QKE", QKE_in)) - m_bc_extdir_vals[BCVars::RhoQKE_bc_comp][ori] = rho_in*QKE_in; - } - } else if (bc_type == "noslipwall") { @@ -550,7 +540,6 @@ void ERF::init_bcs () ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) || ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) || ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) || - ( (BCVars::cons_bc+i == BCVars::RhoQKE_bc_comp) && m_r2d->ingested_QKE() ) || ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) || ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) || ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() )) ) @@ -568,7 +557,6 @@ void ERF::init_bcs () ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) || ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) || ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) || - ( (BCVars::cons_bc+i == BCVars::RhoQKE_bc_comp) && m_r2d->ingested_QKE() ) || ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) || ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) || ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() ) diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 984c57004..7e168ec0b 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -89,18 +89,12 @@ ERF::init_custom (int lev) MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoTheta_comp, RhoTheta_comp, 1, cons_pert.nGrow()); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoScalar_comp,RhoScalar_comp,NSCALARS, cons_pert.nGrow()); - // RhoKE is only relevant if using Deardorff with LES - if (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) { + // RhoKE is only relevant if using Deardorff with LES or MYNN with PBL + if ((solverChoice.turbChoice[lev].les_type == LESType::Deardorff) || + (solverChoice.turbChoice[lev].pbl_type == PBLType::MYNN25)) { MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoKE_comp, RhoKE_comp, 1, cons_pert.nGrow()); } - // RhoQKE is only relevant if using MYNN2.5 - if (solverChoice.turbChoice[lev].pbl_type != PBLType::MYNN25) { - lev_new[Vars::cons].setVal(0.0,RhoQKE_comp,1); - } else { - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQKE_comp, RhoQKE_comp, 1, cons_pert.nGrow()); - } - if (solverChoice.moisture_type != MoistureType::None) { int qstate_size = micro->Get_Qstate_Size(); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ1_comp, RhoQ1_comp, 1, cons_pert.nGrow()); diff --git a/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp b/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp index 3faf0dc0f..9458e7ed4 100644 --- a/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp +++ b/Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp @@ -12,6 +12,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, const MultiFab& yvel, const MultiFab& cons_in, MultiFab& eddyViscosity, + MultiFab& diss, const Geometry& geom, const TurbChoice& turbChoice, std::unique_ptr& most, @@ -48,10 +49,11 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box &bx = mfi.growntilebox(1); - const Array4 &cell_data = cons_in.array(mfi); - const Array4 &K_turb = eddyViscosity.array(mfi); - const Array4 &uvel = xvel.array(mfi); - const Array4 &vvel = yvel.array(mfi); + const Array4& cell_data = cons_in.array(mfi); + const Array4& K_turb = eddyViscosity.array(mfi); + const Array4& diss_arr = diss.array(mfi); + const Array4& uvel = xvel.array(mfi); + const Array4& vvel = yvel.array(mfi); // Compute some quantities that are constant in each column // Sbox is shrunk to only include the interior of the domain in the vertical direction to compute integrals @@ -76,20 +78,20 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, const auto invCellSize = geom.InvCellSizeArray(); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp)); - AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value"); + qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp)); + AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "KE must have a positive value"); Real fac = (sbx.contains(i,j,k)) ? 1.0 : 0.0; const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr); - const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr); + const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr); Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz*fac); Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*dz*fac); }); } else { ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp)); - AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value"); + qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp)); + AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "KE must have a positive value"); // Not multiplying by dz: its constant and would fall out when we divide qint0/qint1 anyway @@ -237,7 +239,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, const Real rho = cell_data(i,j,k,Rho_comp); K_turb(i,j,k,EddyDiff::Mom_v) = rho * Lm * qvel(i,j,k) * SM; K_turb(i,j,k,EddyDiff::Theta_v) = rho * Lm * qvel(i,j,k) * SH; - K_turb(i,j,k,EddyDiff::QKE_v) = rho * Lm * qvel(i,j,k) * SQ; + K_turb(i,j,k,EddyDiff::KE_v) = rho * Lm * qvel(i,j,k) * SQ; // TODO: implement partial-condensation scheme? // Currently, implementation matches NN09 without rain (i.e., @@ -250,7 +252,11 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel, K_turb(i,j,k,EddyDiff::Q_v) = rho * Lm * qvel(i,j,k) * SH; } - K_turb(i,j,k,EddyDiff::PBL_lengthscale) = Lm; + K_turb(i,j,k,EddyDiff::Turb_lengthscale) = Lm; + + // NOTE: qvel is sqrt(2 * TKE) + diss_arr(i,j,k) = rho * std::pow(qvel(i,j,k),3.0) / + (mynn.B1 * K_turb(i,j,k,EddyDiff::Turb_lengthscale)); }); } } diff --git a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp index 92503e5f7..f197fdf33 100644 --- a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp +++ b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp @@ -238,7 +238,7 @@ ComputeDiffusivityYSU (const MultiFab& xvel, const Real rhoKmax = rho * Kmax; K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Mom_v) ,rhoKmax), rhoKmin); K_turb(i,j,k,EddyDiff::Theta_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Theta_v) ,rhoKmax), rhoKmin); - K_turb(i,j,k,EddyDiff::PBL_lengthscale) = pblh_arr(i,j,0); + K_turb(i,j,k,EddyDiff::Turb_lengthscale) = pblh_arr(i,j,0); }); // HACK set bottom ghost cell to 1st cell diff --git a/Source/PBL/ERF_MYNNStruct.H b/Source/PBL/ERF_MYNNStruct.H index 15fbde688..8fce2c81c 100644 --- a/Source/PBL/ERF_MYNNStruct.H +++ b/Source/PBL/ERF_MYNNStruct.H @@ -10,7 +10,7 @@ struct MYNNLevel25 { /* * Calculate the stability functions that determine the eddy diffusivities - * of momentum, heat, QKE, and (optionally) moisture. + * of momentum, heat, KE, and (optionally) moisture. */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE diff --git a/Source/PBL/ERF_PBLHeight.H b/Source/PBL/ERF_PBLHeight.H index 772ba6b05..1ca84631f 100644 --- a/Source/PBL/ERF_PBLHeight.H +++ b/Source/PBL/ERF_PBLHeight.H @@ -130,9 +130,9 @@ struct MYNNPBLH { // // Find PBL height based on TKE (for SBLs only) // - amrex::Real tke = 0.5 * cons_arr(i,j,k ,RhoQKE_comp) / cons_arr(i,j,k ,Rho_comp); - amrex::Real tke1 = 0.5 * cons_arr(i,j,k+1,RhoQKE_comp) / cons_arr(i,j,k+1,Rho_comp); - amrex::Real maxtke = 0.5 * cons_arr(i,j,0 ,RhoQKE_comp) / cons_arr(i,j,0 ,Rho_comp); + amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp); + amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp); + amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp); // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS) amrex::Real TKEeps = 0.05 * maxtke; TKEeps = amrex::max(TKEeps, 0.02); // min val from WRF @@ -142,7 +142,7 @@ struct MYNNPBLH { // Interpolate to get lowest height where TKE -> 0 pblh_tke_arr(i,j,0) = zphys_arr(ii,jj,k) + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(tke1-tke) - * (TKEeps - tke); + * (TKEeps - tke); } } }); @@ -202,9 +202,9 @@ struct MYNNPBLH { // // Find PBL height based on TKE (for SBLs only) // - amrex::Real tke = 0.5 * cons_arr(i,j,k ,RhoQKE_comp) / cons_arr(i,j,k ,Rho_comp); - amrex::Real tke1 = 0.5 * cons_arr(i,j,k+1,RhoQKE_comp) / cons_arr(i,j,k+1,Rho_comp); - amrex::Real maxtke = 0.5 * cons_arr(i,j,0 ,RhoQKE_comp) / cons_arr(i,j,0 ,Rho_comp); + amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp); + amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp); + amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp); // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS) amrex::Real TKEeps = 0.05 * maxtke; TKEeps = amrex::max(TKEeps, 0.02); // min val from WRF @@ -249,7 +249,7 @@ struct MYNNPBLH { // // Finally, blend between the two PBLH estimates // - amrex::Real maxqke = cons_arr(i,j,0 ,RhoQKE_comp) / cons_arr(i,j,0 ,Rho_comp); + amrex::Real maxqke = 2.0 * cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp); if (maxqke > 0.05) { amrex::Real wt = 0.5*std::tanh((zi - sbl_lim)/sbl_damp) + 0.5; pblh_arr(i,j,0) = (1.-wt)*pblh_tke_arr(i,j,0) + wt*pblh_arr(i,j,0); diff --git a/Source/PBL/ERF_PBLModels.H b/Source/PBL/ERF_PBLModels.H index d2fdd9d3c..faadc0fb2 100644 --- a/Source/PBL/ERF_PBLModels.H +++ b/Source/PBL/ERF_PBLModels.H @@ -25,6 +25,7 @@ ComputeDiffusivityMYNN25 (const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& cons_in, amrex::MultiFab& eddyViscosity, + amrex::MultiFab& diss, const amrex::Geometry& geom, const TurbChoice& turbChoice, std::unique_ptr& most, @@ -137,82 +138,4 @@ ComputeVerticalDerivativesPBL (int i, int j, int k, dvdz = 0.25*( vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1) )*dz_inv; } } - -/** - * Function for computing the QKE source terms (NN09, Eqn. 5). - * - * @param[in] u velocity in x-dir - * @param[in] v velocity in y-dir - * @param[in] cell_data conserved cell center vars - * @param[in] cell_prim primitive cell center vars - * @param[in] K_turb turbulent viscosity - * @param[in] cellSizeInv inverse cell size array - * @param[in] domain box of the whole domain - * @param[in] pbl_mynn_B1_l a parameter - * @param[in] theta_mean average theta - */ -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -amrex::Real -ComputeQKESourceTerms (int i, int j, int k, - const amrex::Array4& uvel, - const amrex::Array4& vvel, - const amrex::Array4& cell_data, - const amrex::Array4& cell_prim, - const amrex::Array4& K_turb, - const amrex::GpuArray& cellSizeInv, - const amrex::Box& domain, - amrex::Real pbl_mynn_B1_l, - const amrex::Real theta_mean, - const int RhoQv_comp, - const int RhoQr_comp, - bool c_ext_dir_on_zlo, - bool c_ext_dir_on_zhi, - bool u_ext_dir_on_zlo, - bool u_ext_dir_on_zhi, - bool v_ext_dir_on_zlo, - bool v_ext_dir_on_zhi, - const bool use_most=false, - const amrex::Real met_h_zeta=1.0) -{ - // Compute some relevant derivatives - amrex::Real dthetadz, dudz, dvdz; - amrex::Real source_term = 0.0; - - amrex::Real dz_inv = cellSizeInv[2]; - int izmin = domain.smallEnd(2); - int izmax = domain.bigEnd(2); - - // NOTE: With MOST, the ghost cells are filled AFTER k_turb is computed - // so that the non-explicit pathway works. Therefore, at this - // point we DO have valid ghost cells from MOST. We are passing - // the MOST flag to use one-sided diffs here to be consistent with - // the explicit pathway. - - ComputeVerticalDerivativesPBL(i, j, k, - uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta, - 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, use_most); - - // Note: Transport terms due to turbulence and pressure are included when - // DiffusionSrcForState_* is called from ERF_slow_rhs_post. - - // Shear Production - source_term += 2.0*K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz); - - // Buoyancy Production - source_term -= 2.0*(CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz; - - // Dissipation - amrex::Real qke = cell_prim(i,j,k,PrimQKE_comp); - if (std::abs(qke) > 0.0) { - source_term -= 2.0 * cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) / - (pbl_mynn_B1_l * K_turb(i,j,k,EddyDiff::PBL_lengthscale)); - } - - return source_term; -} #endif diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index a2f7000e1..959b9c24d 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -60,8 +60,8 @@ void make_sources (int level, const bool use_terrain = solverChoice.use_terrain; TurbChoice tc = solverChoice.turbChoice[level]; - const bool l_use_deardorff = (tc.les_type == LESType::Deardorff); - const bool l_use_QKE = tc.use_QKE && tc.diffuse_QKE_3D; + const bool l_use_KE = ( (tc.les_type == LESType::Deardorff) || + (tc.pbl_type == PBLType::MYNN) ); const Box& domain = geom.Domain(); @@ -338,18 +338,12 @@ void make_sources (int level, NumericalDiffusion(bx, start_comp, num_comp, dt, solverChoice.NumDiffCoeff, cell_data, cell_src, mf_u, mf_v, false, false); - if (l_use_deardorff) { + if (l_use_KE) { int sc = RhoKE_comp; int nc = 1; NumericalDiffusion(bx, sc, nc, dt, solverChoice.NumDiffCoeff, cell_data, cell_src, mf_u, mf_v, false, false); } - if (l_use_QKE) { - int sc = RhoQKE_comp; - int nc = 1; - NumericalDiffusion(bx, sc, nc, dt, solverChoice.NumDiffCoeff, - cell_data, cell_src, mf_u, mf_v, false, false); - } { int sc = RhoScalar_comp; int nc = NSCALARS; diff --git a/Source/TimeIntegration/ERF_make_tau_terms.cpp b/Source/TimeIntegration/ERF_make_tau_terms.cpp index 30c8251b0..03901782b 100644 --- a/Source/TimeIntegration/ERF_make_tau_terms.cpp +++ b/Source/TimeIntegration/ERF_make_tau_terms.cpp @@ -50,6 +50,10 @@ void erf_make_tau_terms (int level, int nrk, tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); + const bool only_pbl = ( (tc.les_type == LESType::None) && + (tc.pbl_type != PBLType::None) ); + const bool use_KE = ( (tc.les_type == LESType::Deardorff) || + (tc.pbl_type == PBLType::MYNN25) ); const bool use_most = (most != nullptr); const bool exp_most = (solverChoice.use_explicit_most); @@ -232,11 +236,15 @@ void erf_make_tau_terms (int level, int nrk, // Populate SmnSmn if using Deardorff (used as diff src in post) // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) - if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { + if ((nrk==0) && (use_KE)) { SmnSmn_a = SmnSmn->array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); + SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k, + s11,s22,s33, + s12,s13,s23, + domlo_z,use_most, + exp_most,only_pbl); }); } @@ -344,11 +352,15 @@ void erf_make_tau_terms (int level, int nrk, // Populate SmnSmn if using Deardorff (used as diff src in post) // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) - if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { + if ((nrk==0) && (use_KE)) { SmnSmn_a = SmnSmn->array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); + SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k, + s11,s22,s33, + s12,s13,s23, + domlo_z,use_most, + exp_most,only_pbl); }); } diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index a91e7bdde..6ca0bd9df 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -120,9 +120,8 @@ void erf_slow_rhs_post (int level, int finest_level, if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain); const bool l_use_mono_adv = solverChoice.use_mono_adv; - const bool l_use_QKE = tc.use_QKE; - const bool l_advect_QKE = tc.use_QKE && tc.advect_QKE; - const bool l_use_deardorff = (tc.les_type == LESType::Deardorff); + const bool l_use_KE = ( (tc.les_type == LESType::Deardorff) || + (tc.pbl_type == PBLType::MYNN) ); const bool l_use_diff = ((dc.molec_diff_type != MolecDiffType::None) || (tc.les_type != LESType::None) || (tc.pbl_type != PBLType::None) ); @@ -167,9 +166,8 @@ 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); - if (l_use_deardorff) {is_valid_slow_var[ RhoKE_comp] = 1;} - if (l_use_QKE) {is_valid_slow_var[ RhoQKE_comp] = 1;} - is_valid_slow_var[RhoScalar_comp] = 1; + if (l_use_KE) {is_valid_slow_var[ RhoKE_comp] = 1;} + is_valid_slow_var[RhoScalar_comp] = 1; if (solverChoice.moisture_type != MoistureType::None) { is_valid_slow_var[RhoQ1_comp] = 1; } @@ -282,7 +280,7 @@ void erf_slow_rhs_post (int level, int finest_level, const Array4& mf_v = mapfac_v->const_array(mfi); // SmnSmn for KE src with Deardorff - const Array4& SmnSmn_a = l_use_deardorff ? SmnSmn->const_array(mfi) : Array4{}; + const Array4& SmnSmn_a = l_use_KE ? SmnSmn->const_array(mfi) : Array4{}; // ************************************************************************** // Here we fill the "current" data with "new" data because that is the result of the previous RK stage @@ -395,18 +393,6 @@ void erf_slow_rhs_post (int level, int finest_level, num_comp = 1; } - if (( ivar != RhoQKE_comp ) || - ((ivar == RhoQKE_comp) && l_advect_QKE)) - { - AdvectionSrcForScalars(dt, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_cons, cur_prim, cell_rhs, - l_use_mono_adv, max_s_ptr, min_s_ptr, - detJ_arr, dxInv, mf_m, - horiz_adv_type, vert_adv_type, - horiz_upw_frac, vert_upw_frac, - flx_arr, flx_tmp_arr, domain, bc_ptr_h); - } - if (l_use_diff) { const Array4 tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4{}; @@ -475,8 +461,6 @@ void erf_slow_rhs_post (int level, int finest_level, cur_cons(i,j,k,n) = temp_val / detJ_new_arr(i,j,k); if (ivar == RhoKE_comp) { cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps); - } else if (ivar == RhoQKE_comp) { - cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 1e-12); } }); @@ -489,8 +473,6 @@ void erf_slow_rhs_post (int level, int finest_level, cur_cons(i,j,k,n) = old_cons(i,j,k,n) + dt * cell_rhs(i,j,k,n); if (ivar == RhoKE_comp) { cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), eps); - } else if (ivar == RhoQKE_comp) { - cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 1e-12); } else if (ivar >= RhoQ1_comp) { cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 0.0); } diff --git a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp index 0b7915381..0522183f5 100644 --- a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp +++ b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp @@ -51,7 +51,7 @@ EWP::update (const Real& dt_advance, }, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - cons_array(i,j,k,RhoQKE_comp) = cons_array(i,j,k,RhoQKE_comp) + ewp_array(i,j,k,2)*dt_advance; + cons_array(i,j,k,RhoKE_comp) = cons_array(i,j,k,RhoKE_comp) + ewp_array(i,j,k,2)*dt_advance; }); } } diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp index 5fbf55232..7947c847a 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp @@ -91,7 +91,7 @@ Fitch::update (const Real& dt_advance, }, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - cons_array(i,j,k,RhoQKE_comp) = cons_array(i,j,k,RhoQKE_comp) + fitch_array(i,j,k,4)*dt_advance; + cons_array(i,j,k,RhoKE_comp) = cons_array(i,j,k,RhoKE_comp) + fitch_array(i,j,k,4)*dt_advance; }); } }