Program Listing for File velocity.h¶
↰ Return to documentation for file (source/relations/velocity.h)
/*
* velocity.h Atmospheric Boundary Layer velocity profile relations
*
* Author: Tom Clark (thclark @ github)
*
* Copyright (c) 2013-9 Octue Ltd. All Rights Reserved.
*
*/
#ifndef ES_FLOW_VELOCITY_H_
#define ES_FLOW_VELOCITY_H_
#include <Eigen/Dense>
#include <Eigen/Core>
#include "profile.h"
#include <iostream>
namespace es {
template <typename T, typename Talf>
T power_law_speed(T const & z, const double u_ref, const double z_ref, Talf const & alpha){
T z_norm = z / z_ref;
T speed = pow(z_norm, alpha) * u_ref;
return speed;
};
// Remove template specialisations from doc (causes duplicate) @cond
Eigen::VectorXd power_law_speed(Eigen::VectorXd const & z, const double u_ref, const double z_ref, const double alpha) {
Eigen::VectorXd z_norm = z.array() / z_ref;
Eigen::VectorXd speed = pow(z_norm.array(), alpha) * u_ref;
return speed;
};
// @endcond
template <typename T>
T most_law_speed(T const & z, const double kappa, const double d, const double z0, const double L){
std::cout << "MOST Law not implemented yet" << std::endl;
T speed;
return speed;
}
// Remove template specialisation from doc (causes duplicate) @cond
Eigen::VectorXd most_law_speed(Eigen::VectorXd const & z, const double kappa, const double d, const double z0, const double L){
std::cout << "MOST Law not implemented yet" << std::endl;
Eigen::VectorXd speed;
return speed;
};
// @endcond
template <typename T_z, typename T_pi_j>
T_z marusic_jones_speed(T_z const & z, T_pi_j const pi_j, const double kappa, const double z_0,
const double delta, const double u_inf, const double u_tau){
T_z eta = (z + z_0) / (delta + z_0);
T_z eta_cubed = pow(eta, 3.0);
T_z term1 = log(eta) / kappa;
T_z term2 = (eta_cubed - 1.0) / (3.0 * kappa);
T_z term3 = 2.0 * pi_j * (1.0 - pow(eta, 2.0) * 3.0 + eta_cubed * 2.0) / kappa;
T_z u_deficit = term2 - term1 + term3;
T_z speed = u_inf - u_deficit * u_tau;
return speed;
};
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_pi_j>
Eigen::VectorXd marusic_jones_speed(Eigen::VectorXd const & z, T_pi_j const pi_j, const double kappa, const double z_0,
const double delta, const double u_inf, const double u_tau){
Eigen::VectorXd eta = (z.array() + z_0) / (delta + z_0);
Eigen::VectorXd eta_cubed = eta.array().cube();
Eigen::VectorXd term1 = eta.array().log() / kappa;
Eigen::VectorXd term2 = (eta_cubed.array() - 1.0) / (3.0 * kappa);
Eigen::VectorXd term3 = 2.0 * pi_j * (1.0 - eta.array().square() * 3.0 + eta_cubed.array() * 2.0) / kappa;
Eigen::VectorXd u_deficit = term2 - term1 + term3;
Eigen::VectorXd speed = u_inf - u_deficit.array() * u_tau;
return speed;
};
//@endcond
template <typename T_z, typename T_param>
T_z coles_wake(T_z const &eta, T_param const &pi_coles, const bool corrected=true){
T_z wake_param, eta_sqd;
eta_sqd = pow(eta, 2.0);
if (corrected) {
wake_param = 2.0 * eta_sqd * (3.0 - 2.0 * eta)
- eta_sqd * (1.0 - eta) * (1.0 - 2.0 * eta) / pi_coles;
} else {
wake_param = 2.0 * eta_sqd * (3.0 - 2.0 * eta);
}
return wake_param;
};
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_param>
Eigen::VectorXd coles_wake(Eigen::VectorXd const &eta, T_param const &pi_coles, const bool corrected=true){
Eigen::VectorXd wake_param, eta_sqd;
eta_sqd = eta.array().pow(2.0);
if (corrected) {
wake_param = 2.0 * eta_sqd.array() * (3.0 - 2.0 * eta.array())
- eta_sqd.array() * (1.0 - eta.array()) * (1.0 - 2.0*eta.array()) / pi_coles;
} else {
wake_param = 2.0 * eta_sqd.array() * (3.0 - 2.0 * eta.array());
}
return wake_param;
};
//@endcond
// Do not document @cond
/* Template wrapper for std::isinf to kill off problems where ceres::Jet is used (autodifferentiation) instead of
* a double. Call it is_dbl_inf rather than isinf to avoid accidental use.
*/
bool is_double_inf(double in) {
return std::isinf(in);
};
template <typename T>
bool is_double_inf(T in) {
return false;
};
// @endcond
template <typename T_eta, typename T_param>
T_eta deficit(T_eta const &eta, const double kappa, T_param const & pi_coles, const double shear_ratio, const bool lewkowicz=false) {
T_eta f, ones;
ones = 1.0;
f = -1.0 * log(eta)/ kappa;
if (lewkowicz) {
f = f + (pi_coles/kappa) * coles_wake(ones, pi_coles)
- (pi_coles/kappa) * coles_wake(eta, pi_coles);
} else {
f = f + (pow(eta, 3.0) - 1.0)/(3.0*kappa)
+ 2.0*pi_coles*(1.0 - 3.0*pow(eta, 2.0) + 2.0*pow(eta, 3.0))/kappa;
}
if (is_double_inf(f)) {
f = shear_ratio;
};
return f;
}
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_param>
Eigen::ArrayXd deficit(const Eigen::ArrayXd &eta, const double kappa, T_param const & pi_coles, const double shear_ratio, const bool lewkowicz=false) {
Eigen::ArrayXd f;
Eigen::ArrayXd ones;
ones.setOnes(eta.size());
f = -1.0 * log(eta)/ kappa;
if (lewkowicz) {
f = f + (pi_coles/kappa) * coles_wake(ones, pi_coles)
- (pi_coles/kappa) * coles_wake(eta, pi_coles);
} else {
f = f + (eta.pow(3.0) - 1.0)/(3.0*kappa)
+ 2.0*pi_coles*(1.0 - 3.0*eta.pow(2.) + 2.0*eta.pow(3.0))/kappa;
}
for (int k = 0; k < f.size(); k++) {
if (std::isinf(f[k])) {
f(k) = shear_ratio;
}
}
return f;
}
// @endcond
template <typename T_z, typename T_param>
T_z lewkowicz_speed(T_z const & z, T_param const & pi_coles, T_param const & kappa, T_param const & u_inf, T_param const & shear_ratio, T_param const &delta_c) {
T_z f, speed, eta;
eta = z / delta_c;
T_param u_tau = u_inf / shear_ratio;
T_z term1 = log(eta) / (-1.0*kappa);
T_z term2 = pi_coles * coles_wake(T_z(1.0), pi_coles, true) / kappa;
T_z term3 = pi_coles * coles_wake(eta, pi_coles, true) / kappa;
f = term1 + term2 - term3;
// TODO why isn't this set directly to shear_ratio?
if (is_double_inf(f)) {
f = u_inf / u_tau;
};
speed = u_inf - f * u_tau;
return speed;
};
// Remove template specialisation from doc (causes duplicate) @cond
template <typename T_param>
Eigen::VectorXd lewkowicz_speed(Eigen::VectorXd const & z, T_param const &pi_coles, T_param const &kappa, T_param const &u_inf, T_param const &shear_ratio, T_param const &delta_c=1.0){
Eigen::VectorXd f, speed, eta;
eta = z.array() / delta_c;
T_param u_tau = u_inf/shear_ratio;
Eigen::VectorXd term1 = eta.array().log() / (-1.0*kappa);
double term2 = pi_coles * coles_wake(1.0, pi_coles) / kappa;
Eigen::VectorXd term3 = pi_coles * coles_wake(eta, pi_coles).array() / kappa;
f = term1.array() + term2 - term3.array();
for (int k = 0; k < f.size(); k++) {
if (std::isinf(f[k])) {
f(k) = u_inf / u_tau;
}
}
speed = u_inf - f.array() * u_tau;
return speed;
};
template <typename T_param>
Eigen::ArrayXd lewkowicz_speed(Eigen::ArrayXd const & z, T_param const &pi_coles, T_param const &kappa, T_param const &u_inf, T_param const &shear_ratio, T_param const &delta_c=1.0){
Eigen::ArrayXd f, speed, eta;
eta = z / delta_c;
T_param u_tau = u_inf/shear_ratio;
Eigen::ArrayXd term1 = eta.log() / (-1.0*kappa);
double term2 = pi_coles * coles_wake(1.0, pi_coles) / kappa;
Eigen::ArrayXd term3 = pi_coles * coles_wake(eta, pi_coles) / kappa;
f = term1 + term2 - term3;
for (int k = 0; k < f.size(); k++) {
if (std::isinf(f[k])) {
f(k) = u_inf / u_tau;
}
}
speed = u_inf - f * u_tau;
return speed;
};
// @endcond
} /* namespace es */
#endif /* ES_FLOW_VELOCITY_H_ */