Skip to content

Commit

Permalink
Revert the src terms for mynn to avoid breaking the ctests for now.
Browse files Browse the repository at this point in the history
  • Loading branch information
Aaron Lattanzi committed Oct 24, 2024
1 parent 7aa4b13 commit 92c1687
Show file tree
Hide file tree
Showing 12 changed files with 190 additions and 32 deletions.
2 changes: 1 addition & 1 deletion Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -487,7 +487,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
}

if (turbChoice.pbl_type == PBLType::MYNN25) {
ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity, Diss,
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);
Expand Down
6 changes: 6 additions & 0 deletions Source/Diffusion/ERF_Diffusion.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@ void DiffusionSrcForMom_T (const amrex::Box& bxx, const amrex::Box& bxy, const a
void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain,
int start_comp, int num_comp,
const bool& exp_most,
const amrex::Array4<const amrex::Real>& u,
const amrex::Array4<const amrex::Real>& v,
const amrex::Array4<const amrex::Real>& cell_data,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<amrex::Real>& cell_rhs,
Expand All @@ -61,6 +63,7 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain,
const amrex::Array4<const amrex::Real>& mu_turb,
const SolverChoice& solverChoice,
const int level,
const amrex::Array4<const amrex::Real>& tm_arr,
const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> grav_gpu,
const amrex::BCRec* bc_ptr,
const bool use_most);
Expand All @@ -69,6 +72,8 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
int start_comp, int num_comp,
const bool& exp_most,
const bool& rot_most,
const amrex::Array4<const amrex::Real>& u,
const amrex::Array4<const amrex::Real>& v,
const amrex::Array4<const amrex::Real>& cell_data,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<amrex::Real>& cell_rhs,
Expand Down Expand Up @@ -96,6 +101,7 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
const amrex::Array4<const amrex::Real>& mu_turb,
const SolverChoice& solverChoice,
const int level,
const amrex::Array4<const amrex::Real>& tm_arr,
const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> grav_gpu,
const amrex::BCRec* bc_ptr,
const bool use_most);
Expand Down
39 changes: 36 additions & 3 deletions Source/Diffusion/ERF_DiffusionSrcForState_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@ void
DiffusionSrcForState_N (const Box& bx, const Box& domain,
int start_comp, int num_comp,
const bool& exp_most,
const Array4<const Real>& u,
const Array4<const Real>& v,
const Array4<const Real>& cell_data,
const Array4<const Real>& cell_prim,
const Array4<Real>& cell_rhs,
Expand All @@ -58,6 +60,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
const Array4<const Real>& mu_turb,
const SolverChoice &solverChoice,
const int level,
const Array4<const Real>& tm_arr,
const GpuArray<Real,AMREX_SPACEDIM> grav_gpu,
const BCRec* bc_ptr,
const bool use_most)
Expand All @@ -79,8 +82,9 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
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::MYNN25) );
bool l_use_ddorf = (turbChoice.les_type == LESType::Deardorff);
bool l_use_mynn = (turbChoice.pbl_type == PBLType::MYNN25);

bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha);
bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) ||
(turbChoice.les_type == LESType::Deardorff ) ||
Expand Down Expand Up @@ -675,7 +679,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_KE && (start_comp <= RhoKE_comp) && (end_comp >=RhoKE_comp)) {
if (l_use_ddorf && (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 @@ -702,4 +706,33 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
cell_rhs(i,j,k,qty_index) -= diss(i,j,k);
});
}

// Using PBL
if (l_use_mynn && start_comp <= RhoKE_comp && end_comp >=RhoKE_comp) {
int qty_index = RhoKE_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);
});
}
}
39 changes: 36 additions & 3 deletions Source/Diffusion/ERF_DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
int start_comp, int num_comp,
const bool& exp_most,
const bool& rot_most,
const Array4<const Real>& u,
const Array4<const Real>& v,
const Array4<const Real>& cell_data,
const Array4<const Real>& cell_prim,
const Array4<Real>& cell_rhs,
Expand Down Expand Up @@ -71,6 +73,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
const Array4<const Real>& mu_turb,
const SolverChoice &solverChoice,
const int level,
const Array4<const Real>& tm_arr,
const GpuArray<Real,AMREX_SPACEDIM> grav_gpu,
const BCRec* bc_ptr,
const bool use_most)
Expand All @@ -91,8 +94,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
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::MYNN25) );
bool l_use_ddorf = (turbChoice.les_type == LESType::Deardorff);
bool l_use_mynn = (turbChoice.pbl_type == PBLType::MYNN25);

bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha);
bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) ||
(turbChoice.les_type == LESType::Deardorff ) ||
Expand Down Expand Up @@ -786,7 +790,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_KE && (start_comp <= RhoKE_comp) && (end_comp >= RhoKE_comp)) {
if (l_use_ddorf && (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 @@ -813,4 +817,33 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
cell_rhs(i,j,k,qty_index) -= diss(i,j,k);
});
}

// Using PBL
if (l_use_mynn && start_comp <= RhoKE_comp && end_comp >=RhoKE_comp) {
int qty_index = RhoKE_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);
});
}
}
5 changes: 2 additions & 3 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -362,8 +362,7 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
(solverChoice.turbChoice[lev].pbl_type != PBLType::None) );
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) ||
(solverChoice.turbChoice[lev].pbl_type == PBLType::MYNN25) );
bool l_use_ddorff = (solverChoice.turbChoice[lev].les_type == LESType::Deardorff);
bool l_use_moist = ( solverChoice.moisture_type != MoistureType::None );

BoxArray ba12 = convert(ba, IntVect(1,1,0));
Expand Down Expand Up @@ -430,7 +429,7 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
if (l_use_kturb) {
eddyDiffs_lev[lev] = std::make_unique<MultiFab>(ba, dm, EddyDiff::NumDiffs, 2);
eddyDiffs_lev[lev]->setVal(0.0);
if(l_use_ke) {
if(l_use_ddorff) {
SmnSmn_lev[lev] = std::make_unique<MultiFab>( ba, dm, 1, 0 );
} else {
SmnSmn_lev[lev] = nullptr;
Expand Down
6 changes: 0 additions & 6 deletions Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
const MultiFab& yvel,
const MultiFab& cons_in,
MultiFab& eddyViscosity,
MultiFab& diss,
const Geometry& geom,
const TurbChoice& turbChoice,
std::unique_ptr<ABLMost>& most,
Expand Down Expand Up @@ -51,7 +50,6 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
const Box &bx = mfi.growntilebox(1);
const Array4<Real const>& cell_data = cons_in.array(mfi);
const Array4<Real >& K_turb = eddyViscosity.array(mfi);
const Array4<Real >& diss_arr = diss.array(mfi);
const Array4<Real const>& uvel = xvel.array(mfi);
const Array4<Real const>& vvel = yvel.array(mfi);

Expand Down Expand Up @@ -253,10 +251,6 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
}

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));
});
}
}
79 changes: 78 additions & 1 deletion Source/PBL/ERF_PBLModels.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,6 @@ 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<ABLMost>& most,
Expand Down Expand Up @@ -138,4 +137,82 @@ 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<const amrex::Real>& uvel,
const amrex::Array4<const amrex::Real>& vvel,
const amrex::Array4<const amrex::Real>& cell_data,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<const amrex::Real>& K_turb,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& 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 += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);

// Buoyancy Production
source_term -= (CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;

// Dissipation
amrex::Real qke = 2.0 * cell_prim(i,j,k,PrimKE_comp);
if (std::abs(qke) > 0.0) {
source_term -= cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) /
(pbl_mynn_B1_l * K_turb(i,j,k,EddyDiff::Turb_lengthscale));
}

return source_term;
}
#endif
2 changes: 2 additions & 0 deletions Source/TimeIntegration/ERF_TI_slow_headers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,8 @@ void erf_slow_rhs_post (int level, int finest_level, int nrk,
amrex::Vector<amrex::MultiFab>& S_data,
const amrex::MultiFab& S_prim,
amrex::Vector<amrex::MultiFab >& S_scratch,
const amrex::MultiFab& xvel,
const amrex::MultiFab& yvel,
const amrex::MultiFab& zvel,
const amrex::MultiFab& source,
const amrex::MultiFab* SmnSmn,
Expand Down
4 changes: 2 additions & 2 deletions Source/TimeIntegration/ERF_TI_slow_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,7 +342,7 @@
if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) {
erf_slow_rhs_post(level, finest_level, nrk, slow_dt, n_qstate,
S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
zvel_new, cc_src, SmnSmn, eddyDiffs,
xvel_new, yvel_new, zvel_new, cc_src, SmnSmn, eddyDiffs,
Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc_new[level],
Expand All @@ -359,7 +359,7 @@
} else {
erf_slow_rhs_post(level, finest_level, nrk, slow_dt, n_qstate,
S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
zvel_new, cc_src, SmnSmn, eddyDiffs,
xvel_new, yvel_new, zvel_new, cc_src, SmnSmn, eddyDiffs,
Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc[level],
Expand Down
Loading

0 comments on commit 92c1687

Please sign in to comment.