Skip to content

Commit

Permalink
First pass at removing QKE.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Oct 23, 2024
1 parent be89c40 commit 134cf1d
Show file tree
Hide file tree
Showing 29 changed files with 118 additions and 385 deletions.
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ERF_BoundaryConditions_realbdy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ ERF::fill_from_realbdy (const Vector<MultiFab*>& mfs,
0, 0, 0,
0, 0};
Vector<int> 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};

Expand Down
12 changes: 0 additions & 12 deletions Source/DataStructs/ERF_TurbStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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
42 changes: 10 additions & 32 deletions Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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<Real> 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<Real> Factors = {inv_Pr_t, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr
Gpu::AsyncVector<Real> 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())
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down
56 changes: 13 additions & 43 deletions Source/Diffusion/ERF_DiffusionSrcForState_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ) ||
Expand All @@ -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<Real> alpha_eff(NPRIMVAR_max, 0.0);
if (l_consA) {
for (int i = 0; i < NPRIMVAR_max; ++i) {
Expand Down Expand Up @@ -164,15 +164,15 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
}
}

Vector<int> 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<int> 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<int> 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<int> 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<int> 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<int> 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<Real> alpha_eff_d;
Expand Down Expand Up @@ -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
{
Expand All @@ -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);
});
}

}
56 changes: 13 additions & 43 deletions Source/Diffusion/ERF_DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ) ||
Expand All @@ -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<Real> alpha_eff(NPRIMVAR_max, 0.0);
if (l_consA) {
for (int i = 0; i < NPRIMVAR_max; ++i) {
Expand Down Expand Up @@ -177,15 +177,15 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
}
}

Vector<int> 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<int> 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<int> 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<int> 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<int> 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<int> 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<Real> alpha_eff_d;
Expand Down Expand Up @@ -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
{
Expand All @@ -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);
});
}
}
18 changes: 11 additions & 7 deletions Source/Diffusion/ERF_EddyViscosity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,14 +37,9 @@ ComputeSmnSmn (int& i, int& j, int& k,
const amrex::Array4<amrex::Real const>& 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
Expand Down Expand Up @@ -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;

Expand Down
4 changes: 2 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -887,12 +887,12 @@ private:

amrex::Vector<std::string> plot_var_names_1;
amrex::Vector<std::string> plot_var_names_2;
const amrex::Vector<std::string> cons_names {"density", "rhotheta", "rhoKE", "rhoQKE", "rhoadv_0",
const amrex::Vector<std::string> 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<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar",
const amrex::Vector<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "scalar",
"vorticity_x","vorticity_y","vorticity_z",
"magvel", "divU",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens",
Expand Down
Loading

0 comments on commit 134cf1d

Please sign in to comment.