Skip to content

Commit

Permalink
Limiting of MYNN2.5 for QKE heating. (#1268)
Browse files Browse the repository at this point in the history 
* Limiting of MYNN2.5 for QKE heating.

* Fix cons name in forward declaration.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Oct 16, 2023
1 parent d6e6fe7 commit c1466aa
Show file tree
Hide file tree
Showing 6 changed files with 48 additions and 28 deletions.
3 changes: 2 additions & 1 deletion Exec/DevTests/MiguelDev/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ zlo.type = "Most" # "Most" # "NoSlipWall"
#erf.most.time_average = false
erf.most.z0 = 0.001
erf.most.zref = 15.0 # >=dz/2
erf.most.surf_temp_flux = 0.05
#erf.most.surf_temp = 303.0
zhi.type = "SlipWall"

Expand All @@ -38,7 +39,7 @@ erf.check_file = chk # root name of checkpoint file
erf.check_int =-100 # number of timesteps between checkpoints

# PLOTFILES
erf.plotfile_type = netcdf
#erf.plotfile_type = netcdf
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 100 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhotheta rhoQKE x_velocity y_velocity z_velocity pressure temp theta QKE pres_hse
Expand Down
9 changes: 4 additions & 5 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ void
ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const amrex::MultiFab& yvel,
const amrex::MultiFab& cons_in,
const amrex::MultiFab& cons_old,
amrex::MultiFab& eddyViscosity,
const amrex::Geometry& geom,
const TurbChoice& turbChoice,
Expand Down Expand Up @@ -368,6 +369,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
const amrex::MultiFab& Tau11, const amrex::MultiFab& Tau22, const amrex::MultiFab& Tau33,
const amrex::MultiFab& Tau12, const amrex::MultiFab& Tau13, const amrex::MultiFab& Tau23,
const amrex::MultiFab& cons_in,
const amrex::MultiFab& cons_old,
amrex::MultiFab& eddyViscosity,
amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss,
const amrex::Geometry& geom,
Expand All @@ -388,10 +390,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
// ComputeTurbulentViscosityPBL computes the PBL viscosity just for the vertical component.
//

// We must initialize all the components to 0 because we may not set all of them below
// (which ones depends on which LES / PBL scheme we are using)
eddyViscosity.setVal(0.0);

if (most) {
bool l_use_turb = ( turbChoice.les_type == LESType::Smagorinsky ||
turbChoice.les_type == LESType::Deardorff ||
Expand All @@ -413,7 +411,8 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
}

if (turbChoice.pbl_type != PBLType::None) {
ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity,
// NOTE: state_new is passed in for Cons_old (due to ptr swap in advance)
ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, cons_old, eddyViscosity,
geom, turbChoice, most, bc_ptr, vert_only);
}
}
1 change: 1 addition & 0 deletions Source/Diffusion/EddyViscosity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
const amrex::MultiFab& Tau11, const amrex::MultiFab& Tau22, const amrex::MultiFab& Tau33,
const amrex::MultiFab& Tau12, const amrex::MultiFab& Tau13, const amrex::MultiFab& Tau23,
const amrex::MultiFab& cons_in,
const amrex::MultiFab& cons_old,
amrex::MultiFab& eddyViscosity,
amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss,
const amrex::Geometry& geom,
Expand Down
59 changes: 38 additions & 21 deletions Source/Diffusion/PBLModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ void
ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const amrex::MultiFab& yvel,
const amrex::MultiFab& cons_in,
const amrex::MultiFab& cons_old,
amrex::MultiFab& eddyViscosity,
const amrex::Geometry& geom,
const TurbChoice& turbChoice,
Expand All @@ -31,7 +32,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,

const amrex::Real A1 = turbChoice.pbl_A1;
const amrex::Real A2 = turbChoice.pbl_A2;
//const amrex::Real B1 = turbChoice.pbl_B1;
const amrex::Real B1 = turbChoice.pbl_B1;
const amrex::Real B2 = turbChoice.pbl_B2;
const amrex::Real C1 = turbChoice.pbl_C1;
const amrex::Real C2 = turbChoice.pbl_C2;
Expand All @@ -47,13 +48,17 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
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) );

// Epsilon
amrex::Real eps = std::numeric_limits<amrex::Real>::epsilon();

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( amrex::MFIter mfi(eddyViscosity,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {

const amrex::Box &bx = mfi.growntilebox(1);
const amrex::Array4<amrex::Real const > &cell_data = cons_in.array(mfi);
const amrex::Array4<amrex::Real const > &cell_data = cons_in.array(mfi);
const amrex::Array4<amrex::Real const > &cell_data_old = cons_old.array(mfi);
const amrex::Array4<amrex::Real> &K_turb = eddyViscosity.array(mfi);
const amrex::Array4<amrex::Real const> &uvel = xvel.array(mfi);
const amrex::Array4<amrex::Real const> &vvel = yvel.array(mfi);
Expand All @@ -71,18 +76,20 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const amrex::Box xybx = PerpendicularBox<ZDir>(bx, amrex::IntVect{0,0,0});
amrex::FArrayBox qintegral(xybx,2);
qintegral.setVal<amrex::RunOn::Device>(0.0);
amrex::FArrayBox qturb(bx,1);
amrex::FArrayBox qturb(bx,1); amrex::FArrayBox qturb_old(bx,1);
const amrex::Array4<amrex::Real> qint = qintegral.array();
const amrex::Array4<amrex::Real> qvel= qturb.array();
const amrex::Array4<amrex::Real> qvel = qturb.array();
const amrex::Array4<amrex::Real> qvel_old = qturb_old.array();

amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept
{
const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
const amrex::Real rho = cell_data(i,j,k,Rho_comp);
qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / rho);
// We will divide by qvel later
qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps);
AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value");
AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");

const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
if (sbx.contains(i,j,k)) {
amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k), handler);
amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k), handler);
Expand All @@ -94,7 +101,6 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
int izmax = geom.Domain().bigEnd(2);

// Spatially varying MOST
amrex::Real eps = 1.0e-16;
amrex::Real d_kappa = KAPPA;
amrex::Real d_gravity = CONST_GRAV;

Expand All @@ -106,8 +112,6 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const auto& u_star_arr = u_star_mf->array(mfi);
const auto& t_star_arr = t_star_mf->array(mfi);



amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Compute some partial derivatives that we will need (second order)
Expand Down Expand Up @@ -196,17 +200,30 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
// Overall Length Scale
amrex::Real l_comb = 1.0 / (1.0/l_S + 1.0/l_T + 1.0/l_B);

// NOTE: Level 2 limiting from balance of production and dissipation.
// K_turb has a setval of 0.0 when the MF is created (NOT EACH STEP).
// We do this inline to avoid storing qe^2 at each cell.
amrex::Real shearProd = dudz*dudz + dvdz*dvdz;
amrex::Real buoyProd = -(CONST_GRAV/theta0) * dthetadz;
amrex::Real lSM = K_turb(i,j,k,EddyDiff::Mom_v) / (qvel_old(i,j,k) + eps);
amrex::Real lSH = K_turb(i,j,k,EddyDiff::Theta_v) / (qvel_old(i,j,k) + eps);
amrex::Real qe2 = B1 * l_comb * ( lSM * shearProd + lSH * buoyProd );
amrex::Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);

// Compute non-dimensional parameters
amrex::Real l2_over_q2 = l_comb*l_comb/(qvel(i,j,k)*qvel(i,j,k));
amrex::Real GM = l2_over_q2 * (dudz*dudz + dvdz*dvdz);
amrex::Real GH = -l2_over_q2 * (CONST_GRAV/theta0) * dthetadz;
amrex::Real E1 = 1.0 + 6.0*A1*A1*GM - 9.0*A1*A2*(1.0-C2)*GH;
amrex::Real E2 = -3.0*A1*(4.0*A1 + 3.0*A2*(1.0-C5))*(1.0-C2)*GH;
amrex::Real E3 = 6.0*A2*A1*GM;
amrex::Real E4 = 1.0 - 12.0*A2*A1*(1.0-C2)*GH -3.0*A2*B2*(1.0-C3)*GH;
amrex::Real R1 = A1*(1.0-3.0*C1);

amrex::Real SM = (A2*E2 - R1*E4)/(E2*E3 - E1*E4);
amrex::Real one_m_alpha = (qvel(i,j,k) > qe) ? 1.0 : qvel(i,j,k) / (qe + eps);
amrex::Real one_m_alpha2 = one_m_alpha * one_m_alpha;
amrex::Real l2_over_q2 = l_comb*l_comb/(qvel(i,j,k)*qvel(i,j,k));
amrex::Real GM = l2_over_q2 * shearProd;
amrex::Real GH = l2_over_q2 * buoyProd;
amrex::Real E1 = 1.0 + one_m_alpha2 * ( 6.0*A1*A1*GM - 9.0*A1*A2*(1.0-C2)*GH );
amrex::Real E2 = one_m_alpha2 * ( -3.0*A1*(4.0*A1 + 3.0*A2*(1.0-C5))*(1.0-C2)*GH );
amrex::Real E3 = one_m_alpha2 * ( 6.0*A2*A1*GM );
amrex::Real E4 = 1.0 + one_m_alpha2 * ( -12.0*A2*A1*(1.0-C2)*GH - 3.0*A2*B2*(1.0-C3)*GH );
amrex::Real R1 = one_m_alpha * ( A1*(1.0-3.0*C1) );
amrex::Real R2 = one_m_alpha * A2;

amrex::Real SM = (R2*E2 - R1*E4)/(E2*E3 - E1*E4);
amrex::Real SH = (R1*E3 - A2*E1)/(E2*E3 - E1*E4);
amrex::Real SQ = 3.0 * SM;

Expand Down
1 change: 1 addition & 0 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,6 +392,7 @@ ERF::update_arrays (int lev, const BoxArray& ba, const DistributionMapping& dm)

if (l_use_kturb) {
eddyDiffs_lev[lev] = std::make_unique<MultiFab>( ba, dm, EddyDiff::NumDiffs, 1 );
eddyDiffs_lev[lev]->setVal(0.0);
if(l_use_ddorf) {
SmnSmn_lev[lev] = std::make_unique<MultiFab>( ba, dm, 1, 0 );
} else {
Expand Down
3 changes: 2 additions & 1 deletion Source/TimeIntegration/ERF_advance_dycore.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,11 +216,12 @@ void ERF::advance_dycore(int level,
// *************************************************************************
if (l_use_kturb)
{
// NOTE: state_new transfers to state_old for PBL (due to ptr swap in advance)
const amrex::BCRec* bc_ptr_d = domain_bcs_type_d.data();
ComputeTurbulentViscosity(xvel_old, yvel_old,
*Tau11, *Tau22, *Tau33,
*Tau12, *Tau13, *Tau23,
state_old[IntVar::cons],
state_old[IntVar::cons],state_new[IntVar::cons],
*eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated
fine_geom, *mapfac_u[level], *mapfac_v[level],
tc, solverChoice.gravity, m_most, bc_ptr_d);
Expand Down

0 comments on commit c1466aa

Please sign in to comment.