Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change default advection scheme #1358

Merged
merged 3 commits into from
Jan 4, 2024
Merged

Change default advection scheme #1358

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
804a63a
Make 3rd order upwind the default advection scheme
asalmgren Jan 3, 2024
0849ba6
update docs for advection schemes
asalmgren Jan 3, 2024
6068390
Merge branch 'development' into change_default_advection_scheme
asalmgren Jan 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions Docs/sphinx_doc/Inputs.rst
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
Binary file not shown.
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00001
View file
Open in desktop
Binary file not shown.
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00002
View file
Open in desktop
Binary file not shown.
Binary file modified Tests/ERFGoldFiles/DensityCurrent/Level_0/Cell_D_00003
View file
Open in desktop
Binary file not shown.
Loading
Loading