Skip to content

Commit

Permalink
Change default advection scheme (#1358)
Browse files Browse the repository at this point in the history 
* Make 3rd order upwind the default advection scheme

* update docs for advection schemes
  • Loading branch information
@asalmgren
asalmgren authored Jan 4, 2024
1 parent 850aa40 commit 10adbad
Show file tree
Hide file tree
Showing 93 changed files with 3,618 additions and 1,225 deletions.
8 changes: 4 additions & 4 deletions Docs/sphinx_doc/Inputs.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -553,19 +553,19 @@ List of Parameters
| Parameter | Definition | Acceptable | Default |
| | | Values | |
+==================================+====================+=====================+==============+
| **erf.dycore_horiz_adv_type** | Horizontal | see below | Centered_2nd |
| **erf.dycore_horiz_adv_type** | Horizontal | see below | Upwind_3rd |
| | advection type | | |
| | for dycore vars | | |
+----------------------------------+--------------------+---------------------+--------------+
| **erf.dycore_vert_adv_type** | Vertical | see below | Centered_2nd |
| **erf.dycore_vert_adv_type** | Vertical | see below | Upwind_3rd |
| | advection type | | |
| | for dycore vars | | |
+----------------------------------+--------------------+---------------------+--------------+
| **erf.dryscal_horiz_adv_type** | Horizontal | see below | Centered_2nd |
| **erf.dryscal_horiz_adv_type** | Horizontal | see below | Upwind_3rd |
| | advection type | | |
| | for dry scalars | | |
+----------------------------------+--------------------+---------------------+--------------+
| **erf.dryscal_vert_adv_type** | Vertical | see below | Centered_2nd |
| **erf.dryscal_vert_adv_type** | Vertical | see below | Upwind_3rd |
| | advection type | | |
| | for dry scalars | | |
+----------------------------------+--------------------+---------------------+--------------+
Expand Down
10 changes: 5 additions & 5 deletions Source/DataStructs/AdvStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,13 +201,13 @@ struct AdvChoice {
// Default prefix
std::string pp_prefix {"erf"};

// Spatial discretization
// Order and type of spatial discretizations used in advection
// Defaults given below but these can be over-written at run-time
bool use_efficient_advection = false;
AdvType dycore_horiz_adv_type = AdvType::Centered_2nd;
AdvType dycore_vert_adv_type = AdvType::Centered_2nd;
AdvType dryscal_horiz_adv_type = AdvType::Centered_2nd;
AdvType dryscal_vert_adv_type = AdvType::Centered_2nd;
AdvType dycore_horiz_adv_type = AdvType::Upwind_3rd;
AdvType dycore_vert_adv_type = AdvType::Upwind_3rd;
AdvType dryscal_horiz_adv_type = AdvType::Upwind_3rd;
AdvType dryscal_vert_adv_type = AdvType::Upwind_3rd;
AdvType moistscal_horiz_adv_type = AdvType::Weno_3;
AdvType moistscal_vert_adv_type = AdvType::Weno_3;
};
Expand Down
139 changes: 69 additions & 70 deletions Source/Diffusion/PBLModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include "TurbStruct.H"
#include "PBLModels.H"

using namespace amrex;

/**
* Function to compute turbulent viscosity with PBL.
*
Expand Down Expand Up @@ -33,15 +35,15 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
// MYNN Level 2.5 PBL Model
if (turbChoice.pbl_type == PBLType::MYNN25) {

const amrex::Real A1 = turbChoice.pbl_A1;
const amrex::Real A2 = turbChoice.pbl_A2;
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;
const amrex::Real C3 = turbChoice.pbl_C3;
//const amrex::Real C4 = turbChoice.pbl_C4;
const amrex::Real C5 = turbChoice.pbl_C5;
const Real A1 = turbChoice.pbl_A1;
const Real A2 = turbChoice.pbl_A2;
const Real B1 = turbChoice.pbl_B1;
const Real B2 = turbChoice.pbl_B2;
const Real C1 = turbChoice.pbl_C1;
const Real C2 = turbChoice.pbl_C2;
const Real C3 = turbChoice.pbl_C3;
//const Real C4 = turbChoice.pbl_C4;
const 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) );
Expand All @@ -52,18 +54,18 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
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();
Real eps = std::numeric_limits<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 > &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);
const amrex::Array4<Real const> &cell_data = cons_in.array(mfi);
const amrex::Array4<Real > &K_turb = eddyViscosity.array(mfi);
const amrex::Array4<Real const> &uvel = xvel.array(mfi);
const amrex::Array4<Real const> &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
Expand All @@ -79,15 +81,15 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
amrex::FArrayBox qintegral(xybx,2);
qintegral.setVal<amrex::RunOn::Device>(0.0);
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_old = qturb_old.array();
const amrex::Array4<Real> qint = qintegral.array();
const amrex::Array4<Real> qvel = qturb.array();
const amrex::Array4<Real> qvel_old = qturb_old.array();

// vertical integrals to compute lengthscale
if (use_terrain) {
const amrex::Array4<amrex::Real const> &z_nd = z_phys_nd->array(mfi);
const amrex::Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
const auto invCellSize = geom.InvCellSizeArray();
amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept
{
qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
Expand All @@ -96,14 +98,14 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");

if (sbx.contains(i,j,k)) {
const amrex::Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd);
const amrex::Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd);
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);
amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz, handler);
amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k)*dz, handler);
}
});
} else {
amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept
{
qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
Expand All @@ -113,21 +115,21 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,

// Not multiplying by dz: its constant and would fall out when we divide qint0/qint1 anyway
if (sbx.contains(i,j,k)) {
const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
const Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
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);
}
});
}

amrex::Real dz_inv = geom.InvCellSize(2);
Real dz_inv = geom.InvCellSize(2);
const auto& dxInv = geom.InvCellSizeArray();
int izmin = geom.Domain().smallEnd(2);
int izmax = geom.Domain().bigEnd(2);

// Spatially varying MOST
amrex::Real d_kappa = KAPPA;
amrex::Real d_gravity = CONST_GRAV;
Real d_kappa = KAPPA;
Real d_gravity = CONST_GRAV;

const auto& t_mean_mf = most->get_mac_avg(0,2); // TODO: IS THIS ACTUALLY RHOTHETA
const auto& u_star_mf = most->get_u_star(0); // Use coarsest level
Expand All @@ -137,17 +139,14 @@ 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);

const amrex::Array4<amrex::Real const> z_nd;
if (use_terrain) {
const auto& z_nd = z_phys_nd->array(mfi);
}
const amrex::Array4<Real const> z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4<Real>{};

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Compute some partial derivatives that we will need (second order)
// U and V derivatives are interpolated to account for staggered grid
const amrex::Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd) : 1.0;
amrex::Real dthetadz, dudz, dvdz;
const Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : 1.0;
Real dthetadz, dudz, dvdz;
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,
Expand All @@ -156,22 +155,22 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
dthetadz, dudz, dvdz);

// Spatially varying MOST
amrex::Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
amrex::Real theta0 = tm_arr(i,j,0); // TODO: IS THIS ACTUALLY RHOTHETA
amrex::Real l_obukhov;
Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
Real theta0 = tm_arr(i,j,0); // TODO: IS THIS ACTUALLY RHOTHETA
Real l_obukhov;
if (std::abs(surface_heat_flux) > eps) {
l_obukhov = ( theta0 * u_star_arr(i,j,0) * u_star_arr(i,j,0) ) /
( d_kappa * d_gravity * t_star_arr(i,j,0) );
} else {
l_obukhov = std::numeric_limits<amrex::Real>::max();
l_obukhov = std::numeric_limits<Real>::max();
}

// First Length Scale
AMREX_ASSERT(l_obukhov != 0);
int lk = amrex::max(k,0);
const amrex::Real zval = gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
const amrex::Real zeta = zval/l_obukhov;
amrex::Real l_S;
const Real zval = gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
const Real zeta = zval/l_obukhov;
Real l_S;
if (zeta >= 1.0) {
l_S = KAPPA*zval/3.7;
} else if (zeta >= 0) {
Expand All @@ -181,61 +180,61 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
}

// Second Length Scale
amrex::Real l_T;
Real l_T;
if (qint(i,j,0,1) > 0.0) {
l_T = 0.23*qint(i,j,0,0)/qint(i,j,0,1);
} else {
l_T = std::numeric_limits<amrex::Real>::max();
l_T = std::numeric_limits<Real>::max();
}

// Third Length Scale
amrex::Real l_B;
Real l_B;
if (dthetadz > 0) {
amrex::Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
if (zeta < 0) {
amrex::Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T;
Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T;
qc = std::pow(qc,1.0/3.0);
l_B = (1.0 + 5.0*std::sqrt(qc/(N_brunt_vaisala * l_T))) * qvel(i,j,k)/N_brunt_vaisala;
} else {
l_B = qvel(i,j,k) / N_brunt_vaisala;
}
} else {
l_B = std::numeric_limits<amrex::Real>::max();
l_B = std::numeric_limits<Real>::max();
}

// Overall Length Scale
amrex::Real l_comb = 1.0 / (1.0/l_S + 1.0/l_T + 1.0/l_B);
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 l_comb_old = K_turb(i,j,k,EddyDiff::PBL_lengthscale);
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_old * ( lSM * shearProd + lSH * buoyProd );
amrex::Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);
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;
Real l_comb_old = K_turb(i,j,k,EddyDiff::PBL_lengthscale);
Real shearProd = dudz*dudz + dvdz*dvdz;
Real buoyProd = -(CONST_GRAV/theta0) * dthetadz;
Real lSM = K_turb(i,j,k,EddyDiff::Mom_v) / (qvel_old(i,j,k) + eps);
Real lSH = K_turb(i,j,k,EddyDiff::Theta_v) / (qvel_old(i,j,k) + eps);
Real qe2 = B1 * l_comb_old * ( lSM * shearProd + lSH * buoyProd );
Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);
Real one_m_alpha = (qvel(i,j,k) > qe) ? 1.0 : qvel(i,j,k) / (qe + eps);
Real one_m_alpha2 = one_m_alpha * one_m_alpha;

// 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 * 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*A1*A2*GM );
amrex::Real E4 = 1.0 + one_m_alpha2 * ( -12.0*A1*A2*(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 - R2*E1)/(E2*E3 - E1*E4);
amrex::Real SQ = 3.0 * SM; // Nakanishi & Niino 2009
Real l2_over_q2 = l_comb*l_comb/(qvel(i,j,k)*qvel(i,j,k));
Real GM = l2_over_q2 * shearProd;
Real GH = l2_over_q2 * buoyProd;
Real E1 = 1.0 + one_m_alpha2 * ( 6.0*A1*A1*GM - 9.0*A1*A2*(1.0-C2)*GH );
Real E2 = one_m_alpha2 * ( -3.0*A1*(4.0*A1 + 3.0*A2*(1.0-C5))*(1.0-C2)*GH );
Real E3 = one_m_alpha2 * ( 6.0*A1*A2*GM );
Real E4 = 1.0 + one_m_alpha2 * ( -12.0*A1*A2*(1.0-C2)*GH - 3.0*A2*B2*(1.0-C3)*GH );
Real R1 = one_m_alpha * ( A1*(1.0-3.0*C1) );
Real R2 = one_m_alpha * A2;

Real SM = (R2*E2 - R1*E4)/(E2*E3 - E1*E4);
Real SH = (R1*E3 - R2*E1)/(E2*E3 - E1*E4);
Real SQ = 3.0 * SM; // Nakanishi & Niino 2009

// Finally, compute the eddy viscosity/diffusivities
const amrex::Real rho = cell_data(i,j,k,Rho_comp);
const 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; // 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
2 changes: 1 addition & 1 deletion Tests/ERFGoldFiles/DensityCurrent/Header
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@ density
x_velocity
y_velocity
z_velocity
pressure
theta
pres_hse
dens_hse
pressure
3
10
0
Expand Down
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00000
View file
Binary file not shown.
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00001
View file
Binary file not shown.
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00002
View file
Binary file not shown.
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00003
View file
Binary file not shown.
Loading

0 comments on commit 10adbad

Please sign in to comment.