Skip to content

Commit

Permalink
Merge branch 'development' into update_build_docs_intel
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi authored Feb 29, 2024
2 parents 0fc113a + 16ab087 commit 5810a87
Show file tree
Hide file tree
Showing 9 changed files with 162 additions and 88 deletions.
4 changes: 4 additions & 0 deletions Docs/sphinx_doc/MOST.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,10 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in
(\rho \theta)_{i,j,-n} = \rho_{i,j,-n} \left[ \frac{(\rho\theta)_{i,j,0}}{\rho_{i,j,0}} - \left. \frac{\tau_{\theta z}}{\rho} \right|_{i,j,0} \frac{\rho_{i,j,0}}{K_{\theta,v,(i,j,0)}} n \Delta z \right]
Finally, it must be noted that complex terrain will modify the surface normal and tangent vectors. Consequently, the MOST implentation with terrain will require local vector rotations. While the ERF dycore accounts for
terrain metrics when computing fluxes (e.g. for advection, diffusion, etc.), the impact of terrain metrics on MOST is still a work in progress. Therefore, running with terrain (``erf.use_terrain = true``) and with MOST
(``zlo.type = "Most"``) should be cautioned.

MOST Inputs
~~~~~~~~~~~~~~~~~~~
To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 <https://doi.org/10.1175/JAMC-D-17-0137.1>`_. The rougness calculation method may be specified with
Expand Down
6 changes: 5 additions & 1 deletion Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ zhi.type = "SlipWall"
# MOST BOUNDARY
#=================================================================
zlo.type = "Most"
erf.most.average_policy = 1 # POLICY FOR AVERAGING
erf.most.average_policy = 0 # POLICY FOR AVERAGING
erf.most.use_normal_vector = false # USE NORMAL VECTOR W/ TERRAIN
erf.most.use_interpolation = true # USE INTERPOLATION ON DESTINATION
#erf.most.time_average = true # USE TIME AVERAGING
Expand Down Expand Up @@ -81,3 +81,7 @@ erf.alpha_T = 0.0 # [m^2/s]
prob.T_0 = 300.0
prob.U_0 = 1.0
prob.rho_0 = 1.16
prob.U_0_Pert_Mag = 0.01
prob.V_0_Pert_Mag = 0.01
prob.W_0_Pert_Mag = 0.0
prob.pert_ref_height = 1.0
20 changes: 19 additions & 1 deletion Exec/RegTests/Terrain3d_Hemisphere/prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,12 +12,30 @@ struct ProbParm : ProbParmDefaults {
amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature
amrex::Real U_0 = 0.0;
amrex::Real V_0 = 0.0;

// random initial perturbations (legacy code)
amrex::Real U_0_Pert_Mag = 0.0;
amrex::Real V_0_Pert_Mag = 0.0;
amrex::Real W_0_Pert_Mag = 0.0;

// divergence-free initial perturbations
amrex::Real pert_deltaU = 0.0;
amrex::Real pert_deltaV = 0.0;
amrex::Real pert_periods_U = 5.0;
amrex::Real pert_periods_V = 5.0;
amrex::Real pert_ref_height = 1.0;

// helper vars
amrex::Real aval;
amrex::Real bval;
amrex::Real ufac;
amrex::Real vfac;
}; // namespace ProbParm

class Problem : public ProblemBase
{
public:
Problem();
Problem (const amrex::Real* problo, const amrex::Real* probhi);

#include "Prob/init_density_hse_dry_terrain.H"

Expand Down
152 changes: 100 additions & 52 deletions Exec/RegTests/Terrain3d_Hemisphere/prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,25 +4,33 @@
using namespace amrex;

std::unique_ptr<ProblemBase>
amrex_probinit(
const amrex_real* /*problo*/,
const amrex_real* /*probhi*/)
amrex_probinit (const amrex_real* problo,
const amrex_real* probhi)
{
return std::make_unique<Problem>();
return std::make_unique<Problem>(problo, probhi);
}

Problem::Problem()
Problem::Problem (const amrex::Real* problo,
const amrex::Real* probhi)
{
// Parse params
amrex::ParmParse pp("prob");
pp.query("rho_0", parms.rho_0);
pp.query("U_0", parms.U_0);
pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag);
pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag);
pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag);
pp.query("pert_ref_height", parms.pert_ref_height);
parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]);
parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]);
parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;

init_base_parms(parms.rho_0, parms.T_0);
}

void
Problem::init_custom_pert(
Problem::init_custom_pert (
const Box& bx,
const Box& xbx,
const Box& ybx,
Expand All @@ -34,66 +42,107 @@ Problem::init_custom_pert(
Array4<Real > const&,
Array4<Real > const&,
Array4<Real const> const& z_nd,
Array4<Real const> const& z_cc,
Array4<Real const> const& /*z_cc*/,
GeometryData const& geomdata,
Array4<Real const> const& /*mf_m*/,
Array4<Real const> const& /*mf_u*/,
Array4<Real const> const& /*mf_v*/,
const SolverChoice& sc)
{
const int khi = geomdata.Domain().bigEnd()[2];
const int khi = geomdata.Domain().bigEnd()[2];

const bool use_moisture = (sc.moisture_type != MoistureType::None);
const bool use_terrain = sc.use_terrain;

AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);

// Geometry (note we must include these here to get the data on device)
amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
const auto prob_lo = geomdata.ProbLo();
const auto dx = geomdata.CellSize();
const Real x = prob_lo[0] + (i + 0.5) * dx[0];
const Real z = z_cc(i,j,k);

// Set scalar = 0 everywhere
state(i, j, k, RhoScalar_comp) = 0.0;
// Geometry (note we must include these here to get the data on device)
ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Set scalar = 0 everywhere
state(i, j, k, RhoScalar_comp) = 0.0;

if (use_moisture) {
state(i, j, k, RhoQ1_comp) = 0.0;
state(i, j, k, RhoQ2_comp) = 0.0;
}
});

// Set the x-velocity
amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
x_vel(i, j, k) = 10.0;
});

// Set the y-velocity
amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
y_vel(i, j, k) = 0.0;
});

const auto dx = geomdata.CellSize();
amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv;
dxInv[0] = 1. / dx[0];
dxInv[1] = 1. / dx[1];
dxInv[2] = 1. / dx[2];

// Set the z-velocity from impenetrable condition
amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
z_vel(i, j, k) = 0.0;
});

amrex::Gpu::streamSynchronize();
if (use_moisture) {
state(i, j, k, RhoQ1_comp) = 0.0;
state(i, j, k, RhoQ2_comp) = 0.0;
}
});

// Set the x-velocity
ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
{
const Real* prob_lo = geomdata.ProbLo();
const Real* dx = geomdata.CellSize();
const Real z = use_terrain ? 0.25*( z_nd(i,j ,k) + z_nd(i,j ,k+1)
+ z_nd(i,j+1,k) + z_nd(i,j+1,k+1) )
: prob_lo[2] + (k + 0.5) * dx[2];

x_vel(i, j, k) = 10.0;
if ((z <= parms.pert_ref_height) && (parms.U_0_Pert_Mag != 0.0))
{
Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
Real x_vel_prime = (rand_double*2.0 - 1.0)*parms.U_0_Pert_Mag;
x_vel(i, j, k) += x_vel_prime;
}
});

// Set the y-velocity
ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
{
const Real* prob_lo = geomdata.ProbLo();
const Real* dx = geomdata.CellSize();
const Real x = prob_lo[0] + (i + 0.5) * dx[0];
const Real z = use_terrain ? 0.25*( z_nd(i ,j,k) + z_nd(i ,j,k+1)
+ z_nd(i+1,j,k) + z_nd(i+1,j,k+1) )
: prob_lo[2] + (k + 0.5) * dx[2];

// Set the y-velocity
y_vel(i, j, k) = 0.0;
if ((z <= parms.pert_ref_height) && (parms.V_0_Pert_Mag != 0.0))
{
Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
Real y_vel_prime = (rand_double*2.0 - 1.0)*parms.V_0_Pert_Mag;
y_vel(i, j, k) += y_vel_prime;
}
if (parms.pert_deltaV != 0.0)
{
const amrex::Real xl = x - prob_lo[0];
const amrex::Real zl = z / parms.pert_ref_height;
const amrex::Real damp = std::exp(-0.5 * zl * zl);
y_vel(i, j, k) += parms.vfac * damp * z * std::cos(parms.bval * xl);
}
});

const auto dx = geomdata.CellSize();
amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv;
dxInv[0] = 1. / dx[0];
dxInv[1] = 1. / dx[1];
dxInv[2] = 1. / dx[2];

// Set the z-velocity from impenetrable condition
ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
{
const int dom_lo_z = geomdata.Domain().smallEnd()[2];
const int dom_hi_z = geomdata.Domain().bigEnd()[2];

// Set the z-velocity
if (k == dom_lo_z || k == dom_hi_z+1)
{
z_vel(i, j, k) = 0.0;
}
else if (parms.W_0_Pert_Mag != 0.0)
{
Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
Real z_vel_prime = (rand_double*2.0 - 1.0)*parms.W_0_Pert_Mag;
z_vel(i, j, k) = z_vel_prime;
}
});

amrex::Gpu::streamSynchronize();
}

void
Problem::init_custom_terrain(
Problem::init_custom_terrain (
const Geometry& geom,
MultiFab& z_phys_nd,
const Real& /*time*/)
Expand All @@ -111,7 +160,6 @@ Problem::init_custom_terrain(

// User function parameters
Real a = 0.5;
Real num = 8 * a * a * a;
Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);

Expand Down
6 changes: 3 additions & 3 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,9 +286,9 @@ public:
MODIFIED_CHARNOCK
};

FluxCalcType flux_type;
ThetaCalcType theta_type;
RoughCalcType rough_type;
FluxCalcType flux_type{FluxCalcType::MOENG};
ThetaCalcType theta_type{ThetaCalcType::ADIABATIC};
RoughCalcType rough_type{RoughCalcType::CONSTANT};

private:
bool use_moisture;
Expand Down
4 changes: 2 additions & 2 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,8 @@ ABLMost::compute_most_bcs (const int& lev,
const auto *const u_mean = m_ma.get_average(lev,0);
const auto *const v_mean = m_ma.get_average(lev,1);
const auto *const t_mean = m_ma.get_average(lev,2);
const auto *const u_mag_mean = m_ma.get_average(lev,3);
const auto *const q_mean = m_ma.get_average(lev,3);
const auto *const u_mag_mean = m_ma.get_average(lev,4);

const auto um_arr = u_mean->array(mfi);
const auto vm_arr = v_mean->array(mfi);
Expand Down Expand Up @@ -202,7 +203,6 @@ ABLMost::compute_most_bcs (const int& lev,
Box b2d = bx;
b2d.setBig(2,zlo-1);
int n = RhoTheta_comp;

ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
Expand Down
16 changes: 11 additions & 5 deletions Source/BoundaryConditions/MOSTAverage.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -317,13 +317,15 @@ MOSTAverage::set_k_indices_T()
for (int lev(0); lev < m_maxlev; lev++) {
int kmax = m_geom[lev].Domain().bigEnd(2);
for (MFIter mfi(*m_k_indx[lev], TileNoZ()); mfi.isValid(); ++mfi) {
Box npbx = mfi.tilebox(); npbx.convert({1,1,0});
Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
auto k_arr = m_k_indx[lev]->array(mfi);
ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
k_arr(i,j,k) = 0;
Real z_target = d_zref + z_phys_arr(i,j,k);
Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
+ z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
Real z_target = z_bot_face + d_zref;
for (int lk(0); lk<=kmax; ++lk) {
Real z_lo = 0.25 * ( z_phys_arr(i,j ,lk ) + z_phys_arr(i+1,j ,lk )
+ z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) );
Expand Down Expand Up @@ -392,7 +394,9 @@ MOSTAverage::set_norm_indices_T()
j_arr(i,j,k) = j_new;

// Search for k (grid is stretched in z)
Real z_target = delta_z + z_phys_arr(i,j,k);
Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
+ z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
Real z_target = z_bot_face + delta_z;
for (int lk(0); lk<=kmax; ++lk) {
Real z_lo = 0.25 * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk )
+ z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
Expand Down Expand Up @@ -479,7 +483,7 @@ MOSTAverage::set_norm_positions_T()
const auto dxInv = m_geom[lev].InvCellSizeArray();
IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
for (MFIter mfi(*m_x_pos[lev], TileNoZ()); mfi.isValid(); ++mfi) {
Box npbx = mfi.tilebox(); npbx.convert({1,1,0});
Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
Box gpbx = mfi.growntilebox(ng);
RealBox grb{gpbx,dx.data(),base.dataPtr()};

Expand All @@ -506,7 +510,9 @@ MOSTAverage::set_norm_positions_T()
// Final position at end of vector
x_pos_arr(i,j,k) = x0 + delta_x;
y_pos_arr(i,j,k) = y0 + delta_y;
z_pos_arr(i,j,k) = z_phys_arr(i,j,k) + delta_z;
Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
+ z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
z_pos_arr(i,j,k) = z_bot_face + delta_z;

// Destination position must be contained on the current process!
Real pos[] = {x_pos_arr(i,j,k),y_pos_arr(i,j,k),0.5*dx[2]};
Expand Down
Loading

0 comments on commit 5810a87

Please sign in to comment.