diff --git a/include/teqp/models/saft/correlation_integrals.hpp b/include/teqp/models/saft/correlation_integrals.hpp
index f988f982cf46f1a7480d14775c6cc82fb10ca0f7..2a2e0416cd8b6cd718161cc8f652fa0f4848b905 100644
--- a/include/teqp/models/saft/correlation_integrals.hpp
+++ b/include/teqp/models/saft/correlation_integrals.hpp
@@ -7,6 +7,12 @@
namespace teqp {
namespace SAFTpolar{
+// |x| = sqrt(x^2), the latter is well-suited to differentiation
+template<typename X>
+inline auto differentiable_abs(const X& x){
+ return forceeval(sqrt(x*x));
+};
+
static const std::map<int, std::array<double, 12>> Luckas_J_coeffs = {
{4, {-1.38410152e00, -7.05792933e-01, 2.60947023e00, 1.96828333e01, 1.13619510e01, -2.98510490e01, -3.15686398e01, -2.00943290e01, 5.11029320e01, 1.44194150e01, 9.40061069e00, -2.36844608e01}},
{5, {-6.89702637e-01, -1.62382602e-01, 1.16302441e00, 1.42067443e01, 4.59642681e00, -1.81421003e01, -2.45012804e01, -8.42839734e00, 3.25579587e01, 1.16339969e01, 4.00080085e00, -1.54419815e01}},
@@ -40,12 +46,9 @@ public:
auto get_J(const TType& Tstar, const RhoType& rhostar) const{
double Z_1 = 0.3 + 0.05*n;
double Z_2 = 1.0/n;
- auto A_0 = forceeval(a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar);
- auto A_1 = forceeval(a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar);
- auto A_2 = forceeval(a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar);
- // |x| = sqrt(x^2), the latter is well-suited to differentiation
- auto differentiable_abs = [](const auto& x){ return forceeval(sqrt(x*x)); };
-// return forceeval(A_2*differentiable_abs(Tstar));
+ RhoType A_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
+ RhoType A_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
+ RhoType A_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
std::common_type_t<TType, RhoType> out = (A_0 + A_1*pow(Tstar, Z_1) + A_2*pow(Tstar, Z_2))*exp(1.0/(Tstar + 4.0/pow(differentiable_abs(log(forceeval(rhostar/sqrt(2.0)))), 3.0)));
return out;
}
@@ -77,10 +80,10 @@ public:
double Z_1 = 2.0;
double Z_2 = 3.0;
double Z_3 = 4.0;
- auto b_0 = forceeval(a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar);
- auto b_1 = forceeval(a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar);
- auto b_2 = forceeval(a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar);
- auto b_3 = forceeval(a03 + a13*rhostar + a23*rhostar*rhostar + a33*rhostar*rhostar*rhostar);
+ RhoType b_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
+ RhoType b_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
+ RhoType b_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
+ RhoType b_3 = a03 + a13*rhostar + a23*rhostar*rhostar + a33*rhostar*rhostar*rhostar;
std::common_type_t<TType, RhoType> out = b_0 + b_1*Tstar + b_2*pow(exp(pow(1.0-rhostar/sqrt(2.0),Z_3)), Z_1) + b_3*pow(exp(pow(1.0-rhostar/sqrt(2.0),Z_3)), Z_2);
return out;
}
diff --git a/include/teqp/models/saft/polar_terms.hpp b/include/teqp/models/saft/polar_terms.hpp
index e45faa51908131b50c2656a541767f0ad5a0b5ea..25ddbcf25a50db5743eb601cb9bc8973b752f177 100644
--- a/include/teqp/models/saft/polar_terms.hpp
+++ b/include/teqp/models/saft/polar_terms.hpp
@@ -406,17 +406,18 @@ public:
const auto& x = mole_fractions; // concision
const std::size_t N = mole_fractions.size();
+ using XTtype = std::common_type_t<TTYPE, decltype(mole_fractions[0])>;
std::common_type_t<TTYPE, RhoType, RhoStarType, decltype(mole_fractions[0])> alpha2_112 = 0.0, alpha2_123 = 0.0, alpha2_224 = 0.0;
- const auto factor_112 = forceeval(-2.0*PI_*rhoN/3.0);
- const auto factor_123 = forceeval(-PI_*rhoN/3.0);
- const auto factor_224 = forceeval(-14.0*PI_*rhoN/5.0);
+ const RhoType factor_112 = -2.0*PI_*rhoN/3.0;
+ const RhoType factor_123 = -PI_*rhoN/3.0;
+ const RhoType factor_224 = -14.0*PI_*rhoN/5.0;
for (std::size_t i = 0; i < N; ++i){
for (std::size_t j = 0; j < N; ++j){
TTYPE Tstari = forceeval(T/EPSKIJ(i, i)), Tstarj = forceeval(T/EPSKIJ(j, j));
- auto leading = forceeval(x[i]*x[j]/(Tstari*Tstarj)); // common for all alpha_2 terms
+ XTtype leading = forceeval(x[i]*x[j]/(Tstari*Tstarj)); // common for all alpha_2 terms
TTYPE Tstarij = forceeval(T/EPSKIJ(i, j));
double sigmaij = SIGMAIJ(i,j);
{
@@ -442,6 +443,7 @@ public:
const VecType& x = mole_fractions; // concision
const std::size_t N = mole_fractions.size();
using type = std::common_type_t<TTYPE, RhoType, RhoStarType, decltype(mole_fractions[0])>;
+ using XTtype = std::common_type_t<TTYPE, decltype(mole_fractions[0])>;
type summerA_112_112_224 = 0.0, summerA_112_123_213 = 0.0, summerA_123_123_224 = 0.0, summerA_224_224_224 = 0.0;
type summerB_112_112_112 = 0.0, summerB_112_123_123 = 0.0, summerB_123_123_224 = 0.0, summerB_224_224_224 = 0.0;
@@ -451,7 +453,7 @@ public:
TTYPE Tstari = forceeval(T/EPSKIJ(i,i)), Tstarj = forceeval(T/EPSKIJ(j,j));
TTYPE Tstarij = forceeval(T/EPSKIJ(i,j));
- auto leading = forceeval(x[i]*x[j]/pow(forceeval(Tstari*Tstarj), 3.0/2.0)); // common for all alpha_3A terms
+ XTtype leading = forceeval(x[i]*x[j]/pow(forceeval(Tstari*Tstarj), 3.0/2.0)); // common for all alpha_3A terms
double sigmaij = SIGMAIJ(i,j);
double POW4sigmaij = powi(sigmaij, 4);
double POW8sigmaij = POW4sigmaij*POW4sigmaij;
@@ -482,7 +484,7 @@ public:
double sigmaik = SIGMAIJ(i,k), sigmajk = SIGMAIJ(j,k);
// Lorentz-Berthelot mixing rules for sigma
- auto leadingijk = forceeval(x[i]*x[j]*x[k]/(Tstari*Tstarj*Tstark));
+ XTtype leadingijk = forceeval(x[i]*x[j]*x[k]/(Tstari*Tstarj*Tstark));
if (std::abs(mubar2[i]*mubar2[j]*mubar2[k]) > 0){
auto K222333 = get_Kijk(K222_333, rhostar, Tstarij, Tstarik, Tstarjk);
@@ -508,21 +510,21 @@ public:
}
}
- type alpha3A_112_112_224 = forceeval(8.0*PI_*rhoN/25.0*summerA_112_112_224);
- type alpha3A_112_123_213 = forceeval(8.0*PI_*rhoN/75.0*summerA_112_123_213);
- type alpha3A_123_123_224 = forceeval(8.0*PI_*rhoN/35.0*summerA_123_123_224);
- type alpha3A_224_224_224 = forceeval(144.0*PI_*rhoN/245.0*summerA_224_224_224);
+ type alpha3A_112_112_224 = 8.0*PI_*rhoN/25.0*summerA_112_112_224;
+ type alpha3A_112_123_213 = 8.0*PI_*rhoN/75.0*summerA_112_123_213;
+ type alpha3A_123_123_224 = 8.0*PI_*rhoN/35.0*summerA_123_123_224;
+ type alpha3A_224_224_224 = 144.0*PI_*rhoN/245.0*summerA_224_224_224;
- type alpha3A = forceeval(3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224);
+ type alpha3A = 3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224;
RhoType rhoN2 = rhoN*rhoN;
- type alpha3B_112_112_112 = forceeval(32.0*POW3(PI_)*rhoN2/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112);
- type alpha3B_112_123_123 = forceeval(64.0*POW3(PI_)*rhoN2/315.0*sqrt(3.0*PI_)*summerB_112_123_123);
- type alpha3B_123_123_224 = forceeval(-32.0*POW3(PI_)*rhoN2/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224);
- type alpha3B_224_224_224 = forceeval(32.0*POW3(PI_)*rhoN2/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224);
+ type alpha3B_112_112_112 = 32.0*POW3(PI_)*rhoN2/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112;
+ type alpha3B_112_123_123 = 64.0*POW3(PI_)*rhoN2/315.0*sqrt(3.0*PI_)*summerB_112_123_123;
+ type alpha3B_123_123_224 = -32.0*POW3(PI_)*rhoN2/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224;
+ type alpha3B_224_224_224 = 32.0*POW3(PI_)*rhoN2/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224;
- type alpha3B = forceeval(alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224);
+ type alpha3B = alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224;
return forceeval(alpha3A + alpha3B);
}
diff --git a/include/teqp/models/saftvrmie.hpp b/include/teqp/models/saftvrmie.hpp
index 2c38cb891fd8d56a2198f19c75b4bc3b83404d2a..8f98a182960f398b3436036829fcbbb0ed65e252 100644
--- a/include/teqp/models/saftvrmie.hpp
+++ b/include/teqp/models/saftvrmie.hpp
@@ -835,19 +835,22 @@ public:
// First values for the Mie chain with dispersion (always included)
error_if_expr(T); error_if_expr(rhomolar);
auto vals = terms.get_core_calcs(T, rhomolar, mole_fractions);
- auto alphar = forceeval(vals.alphar_mono + vals.alphar_chain);
+ using type = std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])>;
+ type alphar = vals.alphar_mono + vals.alphar_chain;
+ type packing_fraction = vals.zeta[3];
if (polar){ // polar term is present
using mas = SAFTpolar::multipolar_argument_spec;
- auto visitor = [&](const auto& contrib){
- if constexpr(std::decay_t<decltype(contrib)>::arg_spec == mas::TK_rhoNA3_packingfraction_molefractions){
- auto rho_A3 = forceeval(rhomolar*N_A*1e-30);
- auto packing_fraction = forceeval(vals.zeta[3]);
- auto alpha = contrib.eval(T, rho_A3, packing_fraction, mole_fractions).alpha;
+ auto visitor = [&T, &rhomolar, &mole_fractions, &packing_fraction](const auto& contrib) -> type {
+
+ constexpr mas arg_spec = std::decay_t<decltype(contrib)>::arg_spec;
+ if constexpr(arg_spec == mas::TK_rhoNA3_packingfraction_molefractions){
+ RhoType rho_A3 = rhomolar*N_A*1e-30;
+ type alpha = contrib.eval(T, rho_A3, packing_fraction, mole_fractions).alpha;
return alpha;
}
- else if constexpr(std::decay_t<decltype(contrib)>::arg_spec == mas::TK_rhoNm3_molefractions){
- auto rhoN_m3 = forceeval(rhomolar*N_A);
+ else if constexpr(arg_spec == mas::TK_rhoNm3_molefractions){
+ RhoType rhoN_m3 = rhomolar*N_A;
return contrib.eval(T, rhoN_m3, mole_fractions).alpha;
}
else{