Skip to content

Commit

Permalink
2nd order derivs and fix to shear production. (#1266)
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi authored Oct 13, 2023
1 parent 0fd0e07 commit c79c231
Show file tree
Hide file tree
Showing 7 changed files with 109 additions and 33 deletions.
55 changes: 40 additions & 15 deletions Source/Diffusion/ComputeQKESourceTerm.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,33 +26,58 @@ ComputeQKESourceTerms (int i, int j, int k,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
const amrex::Box& domain,
amrex::Real pbl_B1_l,
const amrex::Real theta_mean)
const amrex::Real theta_mean,
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)
{
// Compute some relevant derivatives
amrex::Real dz_inv = cellSizeInv[2];
int izmin = domain.smallEnd(2);
int izmax = domain.bigEnd(2);
amrex::Real dthetadz, dudz, dvdz;
amrex::Real source_term = 0.0;
if (k == izmax) {
dthetadz = (cell_prim(i,j,k,PrimTheta_comp) - cell_prim(i,j,k-1,PrimTheta_comp))*dz_inv;
dudz = 0.5*(uvel(i,j,k) - uvel(i,j,k-1) + uvel(i+1,j,k) - uvel(i+1,j,k-1))*dz_inv;
dvdz = 0.5*(vvel(i,j,k) - vvel(i,j,k-1) + vvel(i,j+1,k) - vvel(i,j+1,k-1))*dz_inv;
} else if (k == izmin){
dthetadz = (cell_prim(i,j,k+1,PrimTheta_comp) - cell_prim(i,j,k,PrimTheta_comp))*dz_inv;
dudz = 0.5*(uvel(i,j,k+1) - uvel(i,j,k) + uvel(i+1,j,k+1) - uvel(i+1,j,k))*dz_inv;
dvdz = 0.5*(vvel(i,j,k+1) - vvel(i,j,k) + vvel(i,j+1,k+1) - vvel(i,j+1,k))*dz_inv;
if ( k==izmax && c_ext_dir_on_zhi ) {
dthetadz = (1.0/3.0)*(-cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp)
- 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp)
+ 4.0 * cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) )*dz_inv;
} else if ( k==izmin && c_ext_dir_on_zlo ) {
dthetadz = (1.0/3.0)*( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
+ 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp)
- 4.0 * cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dz_inv;
} else {
dthetadz = 0.5*(cell_prim(i,j,k+1,PrimTheta_comp) - cell_prim(i,j,k-1,PrimTheta_comp))*dz_inv;
dudz = 0.25*(uvel(i,j,k+1) - uvel(i,j,k-1) + uvel(i+1,j,k+1) - uvel(i+1,j,k-1))*dz_inv;
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;
dthetadz = 0.5*( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
- cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dz_inv;
}

// Production
source_term += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);
if ( k==izmax && u_ext_dir_on_zhi ) {
dudz = (1.0/6.0)*( (-uvel(i ,j,k-1) - 3.0 * uvel(i ,j,k ) + 4.0 * uvel(i ,j,k+1))
+ (-uvel(i+1,j,k-1) - 3.0 * uvel(i+1,j,k ) + 4.0 * uvel(i+1,j,k+1)) )*dz_inv;
} else if ( k==izmin && u_ext_dir_on_zlo ) {
dudz = (1.0/6.0)*( (uvel(i ,j,k+1) + 3.0 * uvel(i ,j,k ) - 4.0 * uvel(i ,j,k-1))
+ (uvel(i+1,j,k+1) + 3.0 * uvel(i+1,j,k ) - 4.0 * uvel(i+1,j,k-1)) )*dz_inv;
} else {
dudz = 0.25*( uvel(i,j,k+1) - uvel(i,j,k-1) + uvel(i+1,j,k+1) - uvel(i+1,j,k-1) )*dz_inv;
}

if ( k==izmax && v_ext_dir_on_zhi ) {
dvdz = (1.0/6.0)*( (-vvel(i,j ,k-1) - 3.0 * vvel(i,j ,k ) + 4.0 * vvel(i,j ,k+1))
+ (-vvel(i,j+1,k-1) - 3.0 * vvel(i,j+1,k ) + 4.0 * vvel(i,j+1,k+1)) )*dz_inv;
} else if ( k==izmin && v_ext_dir_on_zlo ) {
dvdz = (1.0/6.0)*( (vvel(i,j ,k+1) + 3.0 * vvel(i,j ,k ) - 4.0 * vvel(i,j ,k-1))
+ (vvel(i,j+1,k+1) + 3.0 * vvel(i,j+1,k ) - 4.0 * vvel(i,j+1,k-1)) )*dz_inv;
} else {
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;
}

// Production (We store mu_turb, which is 0.5*K_turb)
source_term += 4.0*K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);

// Bouyancy
source_term -= 2*(CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;
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);
Expand Down
4 changes: 3 additions & 1 deletion Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const amrex::Geometry& geom,
const TurbChoice& turbChoice,
std::unique_ptr<ABLMost>& most,
const amrex::BCRec* bc_ptr,
bool /*vert_only*/);

/**
Expand Down Expand Up @@ -373,6 +374,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v,
const TurbChoice& turbChoice, const Real const_grav,
std::unique_ptr<ABLMost>& most,
const amrex::BCRec* bc_ptr,
bool vert_only)
{
BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity);
Expand Down Expand Up @@ -412,6 +414,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi

if (turbChoice.pbl_type != PBLType::None) {
ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity,
geom, turbChoice, most, vert_only);
geom, turbChoice, most, bc_ptr, vert_only);
}
}
12 changes: 11 additions & 1 deletion Source/Diffusion/DiffusionSrcForState_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -404,10 +404,20 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain,
if (l_use_QKE && start_comp <= RhoQKE_comp && end_comp >=RhoQKE_comp) {
int qty_index = RhoQKE_comp;
auto pbl_B1_l = turbChoice.pbl_B1;
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) );
amrex::ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
mu_turb,cellSizeInv,domain,pbl_B1_l,tm_arr(i,j,0));
mu_turb,cellSizeInv,domain,
pbl_B1_l,tm_arr(i,j,0),
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);
});
}

Expand Down
11 changes: 10 additions & 1 deletion Source/Diffusion/DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -583,10 +583,19 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
if (l_use_QKE && start_comp <= RhoQKE_comp && end_comp >=RhoQKE_comp) {
int qty_index = RhoQKE_comp;
auto pbl_B1_l = turbChoice.pbl_B1;
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) );
amrex::ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
mu_turb,dxInv,domain,pbl_B1_l,tm_arr(i,j,0));
mu_turb,dxInv,domain,pbl_B1_l,tm_arr(i,j,0),
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);
});
}
}
2 changes: 2 additions & 0 deletions Source/Diffusion/EddyViscosity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

#include <ABLMost.H>
#include <DataStruct.H>
#include <AMReX_BCRec.H>

void
ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab& yvel ,
Expand All @@ -17,6 +18,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v,
const TurbChoice& turbChoice, const amrex::Real const_grav,
std::unique_ptr<ABLMost>& most,
const amrex::BCRec* bc_ptr,
bool vert_only = false);

AMREX_GPU_DEVICE
Expand Down
55 changes: 41 additions & 14 deletions Source/Diffusion/PBLModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const amrex::Geometry& geom,
const TurbChoice& turbChoice,
std::unique_ptr<ABLMost>& most,
const amrex::BCRec* bc_ptr,
bool /*vert_only*/)
{
// MYNN Level 2.5 PBL Model
Expand All @@ -38,6 +39,14 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
//const amrex::Real C4 = turbChoice.pbl_C4;
const amrex::Real C5 = turbChoice.pbl_C5;

// Dirichlet flags to switch derivative stencil
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) );

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
Expand Down Expand Up @@ -97,25 +106,43 @@ 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 (1st order at domain boundary)
// Compute some partial derivatives that we will need (second order)
// U and V derivatives are interpolated to account for staggered grid
amrex::Real dthetadz, dudz, dvdz;
if (k == izmax) {
dthetadz = (cell_data(i,j,k,RhoTheta_comp)/cell_data(i,j,k,Rho_comp) -
cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp))*dz_inv;
dudz = 0.5*(uvel(i,j,k) - uvel(i,j,k-1) + uvel(i+1,j,k) - uvel(i+1,j,k-1))*dz_inv;
dvdz = 0.5*(vvel(i,j,k) - vvel(i,j,k-1) + vvel(i,j+1,k) - vvel(i,j+1,k-1))*dz_inv;
} else if (k == izmin){
dthetadz = (cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) -
cell_data(i,j,k,RhoTheta_comp)/cell_data(i,j,k,Rho_comp))*dz_inv;
dudz = 0.5*(uvel(i,j,k+1) - uvel(i,j,k) + uvel(i+1,j,k+1) - uvel(i+1,j,k))*dz_inv;
dvdz = 0.5*(vvel(i,j,k+1) - vvel(i,j,k) + vvel(i,j+1,k+1) - vvel(i,j+1,k))*dz_inv;
if ( k==izmax && c_ext_dir_on_zhi ) {
dthetadz = (1.0/3.0)*(-cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp)
- 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp)
+ 4.0 * cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) )*dz_inv;
} else if ( k==izmin && c_ext_dir_on_zlo ) {
dthetadz = (1.0/3.0)*( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
+ 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp)
- 4.0 * cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dz_inv;
} else {
dthetadz = 0.5*(cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
- cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp))*dz_inv;
}

if ( k==izmax && u_ext_dir_on_zhi ) {
dudz = (1.0/6.0)*( (-uvel(i ,j,k-1) - 3.0 * uvel(i ,j,k ) + 4.0 * uvel(i ,j,k+1))
+ (-uvel(i+1,j,k-1) - 3.0 * uvel(i+1,j,k ) + 4.0 * uvel(i+1,j,k+1)) )*dz_inv;
} else if ( k==izmin && u_ext_dir_on_zlo ) {
dudz = (1.0/6.0)*( (uvel(i ,j,k+1) + 3.0 * uvel(i ,j,k ) - 4.0 * uvel(i ,j,k-1))
+ (uvel(i+1,j,k+1) + 3.0 * uvel(i+1,j,k ) - 4.0 * uvel(i+1,j,k-1)) )*dz_inv;
} else {
dthetadz = 0.5*(cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) -
cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp))*dz_inv;
dudz = 0.25*(uvel(i,j,k+1) - uvel(i,j,k-1) + uvel(i+1,j,k+1) - uvel(i+1,j,k-1))*dz_inv;
}

if ( k==izmax && v_ext_dir_on_zhi ) {
dvdz = (1.0/6.0)*( (-vvel(i,j ,k-1) - 3.0 * vvel(i,j ,k ) + 4.0 * vvel(i,j ,k+1))
+ (-vvel(i,j+1,k-1) - 3.0 * vvel(i,j+1,k ) + 4.0 * vvel(i,j+1,k+1)) )*dz_inv;
} else if ( k==izmin && v_ext_dir_on_zlo ) {
dvdz = (1.0/6.0)*( (vvel(i,j ,k+1) + 3.0 * vvel(i,j ,k ) - 4.0 * vvel(i,j ,k-1))
+ (vvel(i,j+1,k+1) + 3.0 * vvel(i,j+1,k ) - 4.0 * vvel(i,j+1,k-1)) )*dz_inv;
} else {
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;
}

Expand Down Expand Up @@ -185,7 +212,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,

// Finally, compute the eddy viscosity/diffusivities
const amrex::Real rho = cell_data(i,j,k,Rho_comp);
K_turb(i,j,k,EddyDiff::Mom_v) = rho * l_comb * qvel(i,j,k) * SM * 0.5;
K_turb(i,j,k,EddyDiff::Mom_v) = rho * l_comb * qvel(i,j,k) * SM * 0.5; // 0.5 for mu_turb
K_turb(i,j,k,EddyDiff::Theta_v) = rho * l_comb * qvel(i,j,k) * SH;
K_turb(i,j,k,EddyDiff::QKE_v) = rho * l_comb * qvel(i,j,k) * SQ;

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,13 +216,14 @@ void ERF::advance_dycore(int level,
// *************************************************************************
if (l_use_kturb)
{
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],
*eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated
fine_geom, *mapfac_u[level], *mapfac_v[level],
tc, solverChoice.gravity, m_most);
tc, solverChoice.gravity, m_most, bc_ptr_d);
}

// ***********************************************************************************************
Expand Down

0 comments on commit c79c231

Please sign in to comment.