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.

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

wall model refactor #335

Merged
merged 4 commits into from
Feb 8, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
10 changes: 10 additions & 0 deletions amr-wind/wind_energy/ABLWallFunction.H
Original file line number Diff line number Diff line change
Expand Up @@ -81,8 +81,13 @@ public:

void operator()(Field& velocity, const FieldState rho_state) override;

template <typename ShearStress>
void wall_model(
Field& velocity, const FieldState rho_state, const ShearStress& tau);

private:
const ABLWallFunction& m_wall_func;
std::string m_wall_shear_stress_type{"moeng"};
};

class ABLTempWallFunc : public FieldBCIface
Expand All @@ -92,8 +97,13 @@ public:

void operator()(Field& temperature, const FieldState rho_state) override;

template <typename HeatFlux>
void wall_model(
Field& temperature, const FieldState rho_state, const HeatFlux& tau);

private:
const ABLWallFunction& m_wall_func;
std::string m_wall_shear_stress_type{"moeng"};
};

} // namespace amr_wind
Expand Down
124 changes: 94 additions & 30 deletions amr-wind/wind_energy/ABLWallFunction.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "amr-wind/utilities/tensor_ops.H"
#include "amr-wind/utilities/trig_ops.H"
#include "amr-wind/diffusion/diffusion.H"
#include "amr-wind/wind_energy/ShearStress.H"

#include <cmath>

Expand Down Expand Up @@ -266,12 +267,28 @@ void ABLWallFunction::computeusingheatflux()

ABLVelWallFunc::ABLVelWallFunc(Field&, const ABLWallFunction& wall_func)
: m_wall_func(wall_func)
{}
{
amrex::ParmParse pp("ABL");
pp.query("wall_shear_stress_type", m_wall_shear_stress_type);
m_wall_shear_stress_type = amrex::toLower(m_wall_shear_stress_type);

if (m_wall_shear_stress_type == "constant" ||
m_wall_shear_stress_type == "local" ||
m_wall_shear_stress_type == "schumann" ||
m_wall_shear_stress_type == "moeng")
amrex::Print() << "Shear Stress model: " << m_wall_shear_stress_type
<< std::endl;
else {
amrex::Abort("Shear Stress wall model input mistake");
}
}

void ABLVelWallFunc::operator()(Field& velocity, const FieldState rho_state)
template <typename ShearStress>
void ABLVelWallFunc::wall_model(
Field& velocity, const FieldState rho_state, const ShearStress& tau)
{
#if 1
BL_PROFILE("amr-wind::ABLVelWallFunc");

constexpr bool extrapolate = false;
constexpr int idim = 2;
auto& repo = velocity.repo();
Expand All @@ -287,11 +304,6 @@ void ABLVelWallFunc::operator()(Field& velocity, const FieldState rho_state)

const amrex::Real c0 = (!extrapolate) ? 1.0 : 1.5;
const amrex::Real c1 = (!extrapolate) ? 0.0 : -0.5;
const amrex::Real utau2 = m_wall_func.utau() * m_wall_func.utau();
const auto& mo = m_wall_func.mo();
const amrex::Real umeanx = mo.vel_mean[0];
const amrex::Real umeany = mo.vel_mean[1];
const amrex::Real wspd_mean = mo.vmag_mean;

for (int lev = 0; lev < nlevels; ++lev) {
const auto& geom = repo.mesh().Geom(lev);
Expand Down Expand Up @@ -329,17 +341,41 @@ void ABLVelWallFunc::operator()(Field& velocity, const FieldState rho_state)
varr(i, j, k - 1, 2) = 0.0;

// Shear stress BC
amrex::Real taux =
((uu - umeanx) * wspd_mean + wspd * umeanx) /
(wspd_mean * wspd_mean);
amrex::Real tauy =
((vv - umeany) * wspd_mean + wspd * umeany) /
(wspd_mean * wspd_mean);
varr(i, j, k - 1, 0) = taux * den(i, j, k) * utau2 / mu;
varr(i, j, k - 1, 1) = tauy * den(i, j, k) * utau2 / mu;
varr(i, j, k - 1, 0) =
tau.calc_vel_x(uu, wspd) * den(i, j, k) / mu;
varr(i, j, k - 1, 1) =
tau.calc_vel_y(vv, wspd) * den(i, j, k) / mu;
});
}
}
}

void ABLVelWallFunc::operator()(Field& velocity, const FieldState rho_state)
{
#if 1

const auto& mo = m_wall_func.mo();

if (m_wall_shear_stress_type == "moeng") {

auto tau = ShearStressMoeng(mo);
wall_model(velocity, rho_state, tau);

} else if (m_wall_shear_stress_type == "constant") {

auto tau = ShearStressConstant(mo);
wall_model(velocity, rho_state, tau);

} else if (m_wall_shear_stress_type == "local") {

auto tau = ShearStressLocal(mo);
wall_model(velocity, rho_state, tau);

} else if (m_wall_shear_stress_type == "schumann") {

auto tau = ShearStressSchumann(mo);
wall_model(velocity, rho_state, tau);
}

#else
diffusion::wall_model_bc_moeng(
Expand All @@ -349,11 +385,18 @@ void ABLVelWallFunc::operator()(Field& velocity, const FieldState rho_state)

ABLTempWallFunc::ABLTempWallFunc(Field&, const ABLWallFunction& wall_fuc)
: m_wall_func(wall_fuc)
{}
{
amrex::ParmParse pp("ABL");
pp.query("wall_shear_stress_type", m_wall_shear_stress_type);
m_wall_shear_stress_type = amrex::toLower(m_wall_shear_stress_type);
amrex::Print() << "Heat Flux model: " << m_wall_shear_stress_type
<< std::endl;
}

void ABLTempWallFunc::operator()(Field& temperature, const FieldState rho_state)
template <typename HeatFlux>
void ABLTempWallFunc::wall_model(
Field& temperature, const FieldState rho_state, const HeatFlux& tau)
{
#if 1
constexpr bool extrapolate = false;
constexpr int idim = 2;
auto& repo = temperature.repo();
Expand All @@ -364,19 +407,14 @@ void ABLTempWallFunc::operator()(Field& temperature, const FieldState rho_state)
repo.mesh().Geom(0).isPeriodic(idim))
return;

BL_PROFILE("amr-wind::ABLVelWallFunc");
BL_PROFILE("amr-wind::ABLTempWallFunc");
auto& velocity = repo.get_field("velocity");
auto& density = repo.get_field("density", rho_state);
auto& alpha = repo.get_field("temperature_mueff");
const int nlevels = repo.num_active_levels();

const amrex::Real c0 = (!extrapolate) ? 1.0 : 1.5;
const amrex::Real c1 = (!extrapolate) ? 0.0 : -0.5;
const auto& mo = m_wall_func.mo();
const amrex::Real wspd_mean = mo.vmag_mean;
const amrex::Real theta_mean = mo.theta_mean;
const amrex::Real theta_surf = mo.surf_temp;
const amrex::Real term1 = (mo.utau * mo.kappa) / (wspd_mean * mo.phi_h());

for (int lev = 0; lev < nlevels; ++lev) {
const auto& geom = repo.mesh().Geom(lev);
Expand Down Expand Up @@ -411,15 +449,41 @@ void ABLTempWallFunc::operator()(Field& temperature, const FieldState rho_state)
const amrex::Real uu = vold_arr(i, j, k, 0);
const amrex::Real vv = vold_arr(i, j, k, 1);
const amrex::Real wspd = std::sqrt(uu * uu + vv * vv);

const amrex::Real theta2 = told_arr(i, j, k);
const amrex::Real num1 = (theta2 - theta_mean) * wspd_mean;
const amrex::Real num2 = (theta_mean - theta_surf) * wspd;
const amrex::Real tauT = term1 * (num1 + num2);
tarr(i, j, k - 1) = den(i, j, k) * tauT / alphaT;
tarr(i, j, k - 1) =
den(i, j, k) * tau.calc_theta(wspd, theta2) / alphaT;
});
}
}
}

void ABLTempWallFunc::operator()(Field& temperature, const FieldState rho_state)
{
#if 1

const auto& mo = m_wall_func.mo();

if (m_wall_shear_stress_type == "moeng") {

auto tau = ShearStressMoeng(mo);
wall_model(temperature, rho_state, tau);

} else if (m_wall_shear_stress_type == "constant") {

auto tau = ShearStressConstant(mo);
wall_model(temperature, rho_state, tau);

} else if (m_wall_shear_stress_type == "local") {

auto tau = ShearStressLocal(mo);
wall_model(temperature, rho_state, tau);

} else if (m_wall_shear_stress_type == "schumann") {

auto tau = ShearStressSchumann(mo);
wall_model(temperature, rho_state, tau);
}

#else
diffusion::temp_wall_model_bc(
temperature, m_wall_func.instplanar(), rho_state);
Expand Down
164 changes: 164 additions & 0 deletions amr-wind/wind_energy/ShearStress.H
Original file line number Diff line number Diff line change
@@ -0,0 +1,164 @@
#ifndef ShearStress_H
#define ShearStress_H

#include "amr-wind/wind_energy/MOData.H"

/**
* \defgroup shear_stress Shear Stress Wall Models
*
* ShearStress contains functions to compute velocity and temperature shear
* stress wall models the default is the Moeng wall model specifying the wall
* model is done through the input file using ABL.wall_shear_stress_type options
* include "constant", "local", "Schumann", and "Moeng"
*
* \ingroup we_abl
*/

struct ShearStressConstant
{
explicit ShearStressConstant(const amr_wind::MOData& mo)
: utau2(mo.utau * mo.utau)
, u_mean(mo.vel_mean[0])
, v_mean(mo.vel_mean[1])
, wspd_mean(mo.vmag_mean)
, theta_mean(mo.theta_mean)
, theta_surface(mo.surf_temp)
, term1(mo.utau * mo.kappa / mo.phi_h())
{}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_x(amrex::Real /* u */, amrex::Real /* wspd */) const
{
return u_mean / wspd_mean * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_y(amrex::Real /* u */, amrex::Real /* wspd */) const
{
return v_mean / wspd_mean * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_theta(amrex::Real /* wspd */, amrex::Real /* theta */) const
{
return term1 * (theta_mean - theta_surface);
};

amrex::Real utau2;
amrex::Real u_mean;
amrex::Real v_mean;
amrex::Real wspd_mean;
amrex::Real theta_mean;
amrex::Real theta_surface;
amrex::Real term1;
};

struct ShearStressLocal
{
explicit ShearStressLocal(const amr_wind::MOData& mo)
: utau2(mo.utau * mo.utau)
, theta_surface(mo.surf_temp)
, term1(mo.utau * mo.kappa / mo.phi_h())
{}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_x(amrex::Real u, amrex::Real wspd) const
{
return u / amrex::max(wspd, small_vel) * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_y(amrex::Real v, amrex::Real wspd) const
{
return v / amrex::max(wspd, small_vel) * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_theta(amrex::Real /* wspd */, amrex::Real theta) const
{
return term1 * (theta - theta_surface);
};

amrex::Real utau2;
amrex::Real theta_surface;
amrex::Real term1;
amrex::Real small_vel{1.0e-6};
};

struct ShearStressSchumann
{
explicit ShearStressSchumann(const amr_wind::MOData& mo)
: utau2(mo.utau * mo.utau)
, wspd_mean(mo.vmag_mean)
, theta_surface(mo.surf_temp)
, term1(mo.utau * mo.kappa / mo.phi_h())
{}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_x(amrex::Real u, amrex::Real /* wspd */) const
{
return u / wspd_mean * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_y(amrex::Real v, amrex::Real /* wspd */) const
{
return v / wspd_mean * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_theta(amrex::Real /* wspd */, amrex::Real theta) const
{
return term1 * (theta - theta_surface);
};

amrex::Real utau2;
amrex::Real wspd_mean;
amrex::Real theta_surface;
amrex::Real term1;
};

struct ShearStressMoeng
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

May be add doxygen docs and assign them to ABL group?

{
explicit ShearStressMoeng(const amr_wind::MOData& mo)
: utau2(mo.utau * mo.utau)
, u_mean(mo.vel_mean[0])
, v_mean(mo.vel_mean[1])
, wspd_mean(mo.vmag_mean)
, theta_surface(mo.surf_temp)
, theta_mean(mo.theta_mean)
, term1(mo.utau * mo.kappa / (mo.vmag_mean * mo.phi_h()))
{}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_x(amrex::Real u, amrex::Real wspd) const
{
return ((u - u_mean) * wspd_mean + wspd * u_mean) /
(wspd_mean * wspd_mean) * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_vel_y(amrex::Real v, amrex::Real wspd) const
{
return ((v - v_mean) * wspd_mean + wspd * v_mean) /
(wspd_mean * wspd_mean) * utau2;
};

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
calc_theta(amrex::Real wspd, amrex::Real theta) const
{
const amrex::Real num1 = (theta - theta_mean) * wspd_mean;
const amrex::Real num2 = (theta_mean - theta_surface) * wspd;
return term1 * (num1 + num2);
};

amrex::Real utau2;
amrex::Real u_mean;
amrex::Real v_mean;
amrex::Real wspd_mean;
amrex::Real theta_surface;
amrex::Real theta_mean;
amrex::Real term1;
};

#endif /* ShearStress_H */
6 changes: 6 additions & 0 deletions docs/sphinx/user/inputs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -580,7 +580,13 @@ This section is for setting atmospheric boundary layer parameters.
**type:** String, optional, default = ""

Variables for IO for ABL inflow

.. input_param:: ABL.wall_shear_stress_type

**type:** String, optional, default = "Moeng"

Wall shear stress model: options include
"constant", "local", "Schumann", and "Moeng"

Section: Momentum Sources
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand Down