Skip to content

Commit

Permalink
Merge pull request #68 from andrsd/math-pow
Browse files Browse the repository at this point in the history 
Refactoring math::pow
  • Loading branch information
@andrsd
andrsd authored Aug 8, 2024
2 parents 9172fdf + 1f46b5e commit fc0ad6e
Show file tree
Hide file tree
Showing 6 changed files with 128 additions and 80 deletions.
107 changes: 57 additions & 50 deletions include/fprops/Helmholtz.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); ++i)
sum += this->n[i] * std::pow(tau, this->t[i]);
sum += this->n[i] * math::pow(tau, this->t[i]);
return sum;
}

Expand All @@ -214,7 +214,7 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < n.size(); ++i)
sum += this->n[i] * this->t[i] * std::pow(tau, this->t[i] - 1);
sum += this->n[i] * this->t[i] * math::pow(tau, this->t[i] - 1);
return sum;
}

Expand All @@ -223,7 +223,7 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); ++i)
sum += this->n[i] * this->t[i] * (this->t[i] - 1) * std::pow(tau, this->t[i] - 2);
sum += this->n[i] * this->t[i] * (this->t[i] - 1) * math::pow(tau, this->t[i] - 2);
return sum;
}

Expand Down Expand Up @@ -309,8 +309,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); ++i) {
double exp_theta_tau = std::exp(this->theta[i] * tau);
sum += this->n[i] * sqr(this->theta[i]) * this->c[i] * this->d[i] * exp_theta_tau /
sqr(this->c[i] + this->d[i] * exp_theta_tau);
sum += this->n[i] * math::pow<2>(this->theta[i]) * this->c[i] * this->d[i] *
exp_theta_tau / math::pow<2>(this->c[i] + this->d[i] * exp_theta_tau);
}
return sum;
}
Expand Down Expand Up @@ -429,7 +429,7 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i]);
sum += this->n[i] * math::pow(delta, this->d[i]) * math::pow(tau, this->t[i]);
return sum;
}

Expand All @@ -438,8 +438,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * this->d[i] * std::pow(delta, this->d[i] - 1) *
std::pow(tau, this->t[i]);
sum += this->n[i] * this->d[i] * math::pow(delta, this->d[i] - 1) *
math::pow(tau, this->t[i]);
return sum;
}

Expand All @@ -448,8 +448,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i]) * this->t[i] *
std::pow(tau, this->t[i] - 1);
sum += this->n[i] * math::pow(delta, this->d[i]) * this->t[i] *
math::pow(tau, this->t[i] - 1);
return sum;
}

Expand All @@ -459,7 +459,7 @@ class Helmholtz : public SinglePhaseFluidProperties {
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * this->d[i] * (this->d[i] - 1) *
std::pow(delta, this->d[i] - 2) * std::pow(tau, this->t[i]);
math::pow(delta, this->d[i] - 2) * math::pow(tau, this->t[i]);
return sum;
}

Expand All @@ -468,8 +468,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i]) * this->t[i] * (this->t[i] - 1) *
std::pow(tau, this->t[i] - 2);
sum += this->n[i] * math::pow(delta, this->d[i]) * this->t[i] * (this->t[i] - 1) *
math::pow(tau, this->t[i] - 2);
return sum;
}

Expand All @@ -478,8 +478,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * this->d[i] * std::pow(delta, this->d[i] - 1) * this->t[i] *
std::pow(tau, this->t[i] - 1);
sum += this->n[i] * this->d[i] * math::pow(delta, this->d[i] - 1) * this->t[i] *
math::pow(tau, this->t[i] - 1);
return sum;
}

Expand Down Expand Up @@ -551,7 +551,7 @@ class Helmholtz : public SinglePhaseFluidProperties {
sum +=
this->n[i] * std::pow(delta, this->d[i] - 2) * std::pow(tau, this->t[i]) *
std::exp(-std::pow(delta, this->l[i])) *
(sqr(this->l[i]) * std::pow(delta, 2 * this->l[i]) +
(math::pow<2>(this->l[i]) * std::pow(delta, 2 * this->l[i]) +
(this->d[i] - 1) * this->d[i] -
this->l[i] * (2 * this->d[i] + this->l[i] - 1) * std::pow(delta, this->l[i]));
return sum;
Expand Down Expand Up @@ -623,8 +623,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i]) *
std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) -
this->beta[i] * sqr(tau - this->gamma[i]));
std::exp(-this->eta[i] * math::pow<2>(delta - this->epsilon[i]) -
this->beta[i] * math::pow<2>(tau - this->gamma[i]));
return sum;
}

Expand All @@ -635,8 +635,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i] - 1) * std::pow(tau, this->t[i]) *
(this->d[i] - 2.0 * delta * this->eta[i] * (delta - this->epsilon[i])) *
std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) -
this->beta[i] * sqr(tau - this->gamma[i]));
std::exp(-this->eta[i] * math::pow<2>(delta - this->epsilon[i]) -
this->beta[i] * math::pow<2>(tau - this->gamma[i]));
return sum;
}

Expand All @@ -647,8 +647,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i] - 1) *
(2.0 * this->beta[i] * (this->gamma[i] - tau) * tau + this->t[i]) *
std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) -
this->beta[i] * sqr(tau - this->gamma[i]));
std::exp(-this->eta[i] * math::pow<2>(delta - this->epsilon[i]) -
this->beta[i] * math::pow<2>(tau - this->gamma[i]));
return sum;
}

Expand All @@ -658,12 +658,12 @@ class Helmholtz : public SinglePhaseFluidProperties {
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i] - 2) * std::pow(tau, this->t[i]) *
(2 * sqr(delta) * this->eta[i] *
(2 * this->eta[i] * sqr(delta - this->epsilon[i]) - 1) +
sqr(this->d[i]) +
(2 * math::pow<2>(delta) * this->eta[i] *
(2 * this->eta[i] * math::pow<2>(delta - this->epsilon[i]) - 1) +
math::pow<2>(this->d[i]) +
this->d[i] * (4 * delta * this->eta[i] * (this->epsilon[i] - delta) - 1)) *
std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) -
this->beta[i] * sqr(tau - this->gamma[i]));
std::exp(-this->eta[i] * math::pow<2>(delta - this->epsilon[i]) -
this->beta[i] * math::pow<2>(tau - this->gamma[i]));
return sum;
}

Expand All @@ -673,12 +673,12 @@ class Helmholtz : public SinglePhaseFluidProperties {
T sum = 0;
for (std::size_t i = 0; i < this->n.size(); i++)
sum += this->n[i] * std::pow(delta, this->d[i]) * std::pow(tau, this->t[i] - 2) *
(2 * this->beta[i] * sqr(tau) *
(2 * this->beta[i] * sqr(this->gamma[i] - tau) - 1) +
sqr(this->t[i]) +
(2 * this->beta[i] * math::pow<2>(tau) *
(2 * this->beta[i] * math::pow<2>(this->gamma[i] - tau) - 1) +
math::pow<2>(this->t[i]) +
this->t[i] * (4 * this->beta[i] * (this->gamma[i] - tau) * tau - 1)) *
std::exp(-this->eta[i] * sqr(delta - this->epsilon[i]) -
this->beta[i] * sqr(tau - this->gamma[i]));
std::exp(-this->eta[i] * math::pow<2>(delta - this->epsilon[i]) -
this->beta[i] * math::pow<2>(tau - this->gamma[i]));
return sum;
}

Expand All @@ -691,8 +691,8 @@ class Helmholtz : public SinglePhaseFluidProperties {
std::pow(tau, this->t[i] - 1) *
(2. * this->beta[i] * (this->gamma[i] - tau) * tau + this->t[i]) *
(2. * delta * this->eta[i] * (this->epsilon[i] - delta) + this->d[i]) *
std::exp(-this->beta[i] * sqr(this->gamma[i] - tau) -
this->eta[i] * sqr(delta - this->epsilon[i]));
std::exp(-this->beta[i] * math::pow<2>(this->gamma[i] - tau) -
this->eta[i] * math::pow<2>(delta - this->epsilon[i]));
return sum;
}

Expand Down Expand Up @@ -837,13 +837,13 @@ class Helmholtz : public SinglePhaseFluidProperties {
theta(std::size_t i, T delta, T tau) const
{
return (1.0 - tau) +
this->A[i] * std::pow(sqr(delta - 1.0), 1.0 / (2.0 * this->beta[i]));
this->A[i] * std::pow(math::pow<2>(delta - 1.0), 1.0 / (2.0 * this->beta[i]));
}

T
dtheta_ddelta(std::size_t i, T delta, T tau) const
{
return this->A[i] * std::pow(sqr(delta - 1), 1. / (2. * this->beta[i])) /
return this->A[i] * std::pow(math::pow<2>(delta - 1), 1. / (2. * this->beta[i])) /
(this->beta[i] * (delta - 1));
}

Expand All @@ -858,14 +858,15 @@ class Helmholtz : public SinglePhaseFluidProperties {
{
return -1 *
(this->A[i] * (this->beta[i] - 1) *
std::pow(sqr(delta - 1), 1. / (2 * this->beta[i]) - 1.)) /
(sqr(this->beta[i]));
std::pow(math::pow<2>(delta - 1), 1. / (2 * this->beta[i]) - 1.)) /
(math::pow<2>(this->beta[i]));
}

T
DELTA(std::size_t i, T delta, T tau) const
{
return sqr(theta(i, delta, tau)) + this->B[i] * std::pow(sqr(delta - 1.0), this->a[i]);
return math::pow<2>(theta(i, delta, tau)) +
this->B[i] * std::pow(math::pow<2>(delta - 1.0), this->a[i]);
}

T
Expand All @@ -884,16 +885,16 @@ class Helmholtz : public SinglePhaseFluidProperties {
T
d2DELTA_ddelta2(std::size_t i, T delta, T tau) const
{
return 2. * (sqr(dtheta_ddelta(i, delta, tau) +
theta(i, delta, tau) * d2theta_ddelta2(i, delta, tau))) +
return 2. * (math::pow<2>(dtheta_ddelta(i, delta, tau) +
theta(i, delta, tau) * d2theta_ddelta2(i, delta, tau))) +
2 * this->a[i] * (2 * this->a[i] - 1) * this->B[i] *
std::pow(delta - 1, 2 * this->a[i] - 2);
}

T
d2DELTA_dtau2(std::size_t i, T delta, T tau) const
{
return 2. * sqr(dtheta_dtau(i, delta, tau));
return 2. * math::pow<2>(dtheta_dtau(i, delta, tau));
}

T
Expand Down Expand Up @@ -946,42 +947,48 @@ class Helmholtz : public SinglePhaseFluidProperties {
T
PSI(std::size_t i, T delta, T tau) const
{
return std::exp(-this->C[i] * sqr(delta - 1.0) - this->D[i] * sqr(tau - 1.0));
return std::exp(-this->C[i] * math::pow<2>(delta - 1.0) -
this->D[i] * math::pow<2>(tau - 1.0));
}

T
dPSI_ddelta(std::size_t i, T delta, T tau) const
{
return -2 * this->C[i] * (delta - 1) *
std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1));
std::exp(-this->C[i] * math::pow<2>(delta - 1) -
this->D[i] * math::pow<2>(tau - 1));
}

T
dPSI_dtau(std::size_t i, T delta, T tau) const
{
return -2 * this->D[i] * (tau - 1) *
std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1));
std::exp(-this->C[i] * math::pow<2>(delta - 1) -
this->D[i] * math::pow<2>(tau - 1));
}

T
d2PSI_ddelta2(std::size_t i, T delta, T tau) const
{
return 2 * this->C[i] * (2 * this->C[i] * sqr(delta - 1) - 1) *
std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1));
return 2 * this->C[i] * (2 * this->C[i] * math::pow<2>(delta - 1) - 1) *
std::exp(-this->C[i] * math::pow<2>(delta - 1) -
this->D[i] * math::pow<2>(tau - 1));
}

T
d2PSI_dtau2(std::size_t i, T delta, T tau) const
{
return 2 * this->D[i] * (2 * this->D[i] * sqr(tau - 1) - 1) *
std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1));
return 2 * this->D[i] * (2 * this->D[i] * math::pow<2>(tau - 1) - 1) *
std::exp(-this->C[i] * math::pow<2>(delta - 1) -
this->D[i] * math::pow<2>(tau - 1));
}

T
d2PSI_ddeltatau(std::size_t i, T delta, T tau) const
{
return 4. * this->C[i] * (delta - 1) * this->D[i] * (tau - 1) *
std::exp(-this->C[i] * sqr(delta - 1) - this->D[i] * sqr(tau - 1));
std::exp(-this->C[i] * math::pow<2>(delta - 1) -
this->D[i] * math::pow<2>(tau - 1));
}

std::vector<T> n;
Expand Down
63 changes: 52 additions & 11 deletions include/fprops/Numerics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,25 +4,58 @@
#pragma once

#include <functional>
#include <stdexcept>
#include <type_traits>
#include <cmath>

namespace fprops {

inline double
sqr(double x)
{
return x * x;
}
namespace math {

template <int N, typename T>
struct pow_impl {
static inline T
value(T x)
{
if (N % 2)
return x * pow_impl<N - 1, T>::value(x);
T x_n_half = pow_impl<N / 2, T>::value(x);
return x_n_half * x_n_half;
}
};

inline double
cb(double x)
template <typename T>
struct pow_impl<1, T> {
static inline T
value(T x)
{
return x;
}
};

template <typename T>
struct pow_impl<0, T> {
static inline T
value(T)
{
return 1;
}
};

template <int N, typename T>
inline T
pow(T x)
{
return x * x * x;
return pow_impl<N, T>::value(x);
}

namespace math {

/// Compute `e` power of `x`
///
/// @param x Operand
/// @param e Exponent (integer)
/// @return `x` to the power of `e`
template <typename T>
T
inline T
pow(T x, int e)
{
bool neg = false;
Expand All @@ -46,6 +79,14 @@ pow(T x, int e)
return neg ? 1.0 / result : result;
}

/// Alias for covenience, or potentially for better implementation
template <typename T, typename EXP>
inline T
pow(T x, EXP exp)
{
return std::pow(x, exp);
}

} // namespace math

namespace newton {
Expand Down
Loading

0 comments on commit fc0ad6e

Please sign in to comment.