Skip to content

Commit

Permalink
Finalize the switch to sidecars for J and K
Browse files Browse the repository at this point in the history 
Seems to make the compilation about 40% faster
  • Loading branch information
@ianhbell
ianhbell committed Jul 31, 2024
1 parent cf71d9c commit 99291b4
Show file tree
Hide file tree
Showing 5 changed files with 148 additions and 65 deletions.
116 changes: 111 additions & 5 deletions include/teqp/models/saft/correlation_integrals.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,21 @@ class LuckasKIntegral{
}
};

class KLuckasSidecar{
public:

const int n1, n2;
const std::array<double, 16> a;
double a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33;

KLuckasSidecar(const int n1, const int n2) : n1(n1), n2(n2), a(Luckas_K_coeffs.at({n1, n2})){
a00 = a[0]; a01 = a[1]; a02 = a[2]; a03 = a[3];
a10 = a[4]; a11 = a[5]; a12 = a[6]; a13 = a[7];
a20 = a[8]; a21 = a[9]; a22 = a[10]; a23 = a[11];
a30 = a[12]; a31 = a[13]; a32 = a[14]; a33 = a[15];
}
};

static const std::map<int, std::array<double, 6>> GubbinsTwu_J_coeffs = {
{4, {-0.257431, 0.439229, 0.414783, -0.457019, -0.145520, 0.299666}},
{5, {-0.396724, 0.690721, 0.628935, -0.652622, -0.201462, -0.231635}},
Expand Down Expand Up @@ -197,7 +212,21 @@ class GubbinsTwuKIntegral{
}
};


class KGubbinsTwuSidecar{
public:

const int n1, n2;
const std::array<double, 6> a;
double A, B, C, D, E, F;
double sign_term;

KGubbinsTwuSidecar(const int n1, const int n2) : n1(n1), n2(n2), a(GubbinsTwu_K_coeffs.at({n1, n2})){
A = a[0]; B = a[1]; C = a[2]; D = a[3]; E = a[4]; F = a[5];

// The {334, 445} term has opposite sign to the others
sign_term = ((n1==334) && (n2==445) ? -1 : 1);
}
};

//Type R^2 a[0, 0] a[0, 1] a[0,2] a[0,3] a[1, 0] a[1, 1] a[1,2] a[1,3] a[2, 0] a[2, 1] a[2,2] a[2,3] a[3, 0] a[3, 1] a[3,2] a[3,3] a[4, 0] a[4, 1] a[4,2] a[4,3] b[0, 0] b[0, 1] b[0,2] b[1, 0] b[1, 1] b[1,2] b[2, 0] b[2, 1] b[2,2] b[3, 0] b[3, 1] b[3,2] b[4, 0] b[4, 1] b[4,2]
static const std::map<int, std::array<double, 35>> Gottschalk_J_coeffs = {
Expand Down Expand Up @@ -421,6 +450,28 @@ class GottschalkKIntegral{
}
};


class KGottschalkSidecar{
private:
std::tuple<int,int,int> int2key(int i){
if (i < 222){
throw teqp::InvalidArgument("");
}
if (i > 999){
throw teqp::InvalidArgument("");
}
return std::make_tuple(i / 100 % 10, i / 10 % 10, i % 10);
}
public:
const std::tuple<int,int,int> k1, k2;
const std::array<double, 40> abc;

/// Constructor taking two tuples of ints
KGottschalkSidecar(std::tuple<int,int,int> k1, std::tuple<int,int,int> k2) : k1(k1), k2(k2), abc(Gottschalk_K_coeffs.at({k1, k2})){}
/// Constructor taking two three digit integers, each of which are split into tuples of ints
KGottschalkSidecar(int k1, int k2) : k1(int2key(k1)), k2(int2key(k2)), abc(Gottschalk_K_coeffs.at({this->k1, this->k2})){}
};

using JSidecar = std::variant<JLuckasSidecar, JGubbinsTwuSidecar, JGottschalkSidecar>;

class JIntegral{
Expand All @@ -430,7 +481,7 @@ class JIntegral{
JIntegral(const JSidecar& sidecar) : sidecar(sidecar) {};

template<typename TType, typename RhoType>
auto call(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t<TType, RhoType>{
auto get_J(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t<TType, RhoType>{

// Runtime type switching
if (std::holds_alternative<JLuckasSidecar>(sidecar)){
Expand Down Expand Up @@ -470,12 +521,67 @@ class JIntegral{
throw teqp::InvalidArgument("don't know what to do with this sidecar");
}
}
};

using KSidecar = std::variant<KLuckasSidecar, KGubbinsTwuSidecar, KGottschalkSidecar>;

class KIntegral{
const KSidecar sidecar;
public:

KIntegral(const KSidecar& sidecar) : sidecar(sidecar) {};

template<typename TType, typename RhoType>
auto get_J(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t<TType, RhoType>{
return call(Tstar, rhostar);
auto get_K(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t<TType, RhoType>{

// Runtime type switching
if (std::holds_alternative<KLuckasSidecar>(sidecar)){
const auto& d = std::get<KLuckasSidecar>(sidecar);
double Z_1 = 2.0;
double Z_2 = 3.0;
double Z_3 = 4.0;
RhoType b_0 = d.a00 + d.a10*rhostar + d.a20*rhostar*rhostar + d.a30*rhostar*rhostar*rhostar;
RhoType b_1 = d.a01 + d.a11*rhostar + d.a21*rhostar*rhostar + d.a31*rhostar*rhostar*rhostar;
RhoType b_2 = d.a02 + d.a12*rhostar + d.a22*rhostar*rhostar + d.a32*rhostar*rhostar*rhostar;
RhoType b_3 = d.a03 + d.a13*rhostar + d.a23*rhostar*rhostar + d.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;
}
else if (std::holds_alternative<KGubbinsTwuSidecar>(sidecar)){
const auto& d = std::get<KGubbinsTwuSidecar>(sidecar);
std::common_type_t<TType, RhoType> out = d.sign_term*exp(d.A*rhostar*rhostar*log(Tstar) + d.B*rhostar*rhostar + d.C*rhostar*log(Tstar) + d.D*rhostar + d.E*log(Tstar) + d.F);
return out;
}
else if (std::holds_alternative<KGottschalkSidecar>(sidecar)){
const auto& d = std::get<KGottschalkSidecar>(sidecar);
const auto& abc = d.abc;
std::common_type_t<TType, RhoType> summer = 0.0;
int N1 = 8, N2 = 8; // N3 = 24

for (auto i = 0; i <= 3; ++i){
for (auto j = 1; j <= 2; ++j){
auto I = 2*i + (j-1); // 2 entries in j for each i
summer += abc[I]*pow(rhostar, i)*pow(exp((1.0-rhostar/3.0)/Tstar), j);
}
}
for (auto i = 0; i <= 3; ++i){
for (auto j = 1; j <= 2; ++j){
auto I = N1 + 2*i + (j-1); // 2 entries in j for each i
summer += abc[I]*pow(rhostar, i)*pow(exp((1.0-rhostar/3.0)*(1.0-rhostar/3.0)/Tstar), j);
}
}
for (auto i = 0; i <= 5; ++i){
for (auto j = 0; j <= 3; ++j){
auto I = N1 + N2 + 4*i + j; // 4 entries in j for each i
summer += abc[I]*pow(rhostar, i)*pow(Tstar, j);
}
}
return summer;
}
else{
throw teqp::InvalidArgument("don't know what to do with this sidecar");
}
}

};

}
Expand Down
64 changes: 23 additions & 41 deletions include/teqp/models/saft/polar_terms.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,36 +38,30 @@ auto get_Kijk_334445(const KType& Kint, const RhoType& rhostar, const Txy &Tstar
};

/**
\tparam KIntegral A type that can be indexed with a two integers a and b to give the K(a,b) integral
The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface)
The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface)
*/
template<class JSidecar, class KIntegral>

class MultipolarContributionGubbinsTwu {
public:
static constexpr multipolar_argument_spec arg_spec = multipolar_argument_spec::TK_rhoNm3_rhostar_molefractions;
private:
const JIntegral J6, J8, J10, J11, J13, J15;
const KIntegral K222_333, K233_344, K334_445, K444_555;

const Eigen::ArrayXd sigma_m, epsilon_over_k, mubar2, Qbar2;
const bool has_a_polar;
const Eigen::ArrayXd sigma_m3, sigma_m5;

const JIntegral J6{JSidecar{6}};
const JIntegral J8{JSidecar{8}};
const JIntegral J10{JSidecar{10}};
const JIntegral J11{JSidecar{11}};
const JIntegral J13{JSidecar{13}};
const JIntegral J15{JSidecar{15}};
const KIntegral K222_333{222, 333};
const KIntegral K233_344{233, 344};
const KIntegral K334_445{334, 445};
const KIntegral K444_555{444, 555};
const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2
const double PI_ = static_cast<double>(EIGEN_PI);
Eigen::MatrixXd SIGMAIJ, EPSKIJ;
multipolar_rhostar_approach approach = multipolar_rhostar_approach::use_packing_fraction;

public:
MultipolarContributionGubbinsTwu(const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mubar2, const Eigen::ArrayX<double> &Qbar2, multipolar_rhostar_approach approach) : sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2), has_a_polar(mubar2.cwiseAbs().sum() > 0 || Qbar2.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(approach) {

MultipolarContributionGubbinsTwu(const auto& sidecarJ, const auto& sidecarK, const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mubar2, const Eigen::ArrayX<double> &Qbar2, multipolar_rhostar_approach approach) :
J6(decltype(sidecarJ){6}), J8(decltype(sidecarJ){8}), J10(decltype(sidecarJ){10}), J11(decltype(sidecarJ){11}), J13(decltype(sidecarJ){13}), J15(decltype(sidecarJ){15}),
K222_333(decltype(sidecarK){222,333}), K233_344(decltype(sidecarK){233,344}), K334_445(decltype(sidecarK){334,445}), K444_555(decltype(sidecarK){444,555}),
sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2), has_a_polar(mubar2.cwiseAbs().sum() > 0 || Qbar2.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(approach) {
// Check lengths match
if (sigma_m.size() != mubar2.size()){
throw teqp::InvalidArgument("bad size of mubar2");
Expand Down Expand Up @@ -111,15 +105,15 @@ class MultipolarContributionGubbinsTwu {
double sigmaij = SIGMAIJ(i,j);
{
double dbl = sigma_m3[i]*sigma_m3[j]/powi(sigmaij,3)*mubar2[i]*mubar2[j];
alpha2_112 += leading*dbl*J6.call(Tstarij, rhostar);
alpha2_112 += leading*dbl*J6.get_J(Tstarij, rhostar);
}
{
double dbl = sigma_m3[i]*sigma_m5[j]/powi(sigmaij,5)*mubar2[i]*Qbar2[j];
alpha2_123 += leading*dbl*J8.call(Tstarij, rhostar);
alpha2_123 += leading*dbl*J8.get_J(Tstarij, rhostar);
}
{
double dbl = sigma_m5[i]*sigma_m5[j]/powi(sigmaij,7)*Qbar2[i]*Qbar2[j];
alpha2_224 += leading*dbl*J10.call(Tstarij, rhostar);
alpha2_224 += leading*dbl*J10.get_J(Tstarij, rhostar);
}
}
}
Expand Down Expand Up @@ -274,31 +268,21 @@ struct PolarizableArrays{
The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface)
*/
template<class JSidecar, class KIntegral>
class MultipolarContributionGrayGubbins {
public:
static constexpr multipolar_argument_spec arg_spec = multipolar_argument_spec::TK_rhoNm3_rhostar_molefractions;
private:
const JIntegral J6, J8, J10, J11, J13, J15;
const KIntegral K222_333, K233_344, K334_445, K444_555;

const Eigen::ArrayXd sigma_m, epsilon_over_k;
Eigen::MatrixXd SIGMAIJ, EPSKIJ;
const Eigen::ArrayXd mu, Q, mu2, Q2, Q3;
const bool has_a_polar;
const Eigen::ArrayXd sigma_m3, sigma_m5;

const JIntegral J6{JSidecar{6}};
const JIntegral J8{JSidecar{8}};
const JIntegral J10{JSidecar{10}};
const JIntegral J11{JSidecar{11}};
const JIntegral J13{JSidecar{13}};
const JIntegral J15{JSidecar{15}};
const KIntegral K222_333{222, 333};
const KIntegral K233_344{233, 344};
const KIntegral K334_445{334, 445};
const KIntegral K444_555{444, 555};

const double PI_ = static_cast<double>(EIGEN_PI);
const double PI3 = PI_*PI_*PI_;
const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2
const double k_e = teqp::constants::k_e; // coulomb constant, with units of N m^2 / C^2
const double k_B = teqp::constants::k_B; // Boltzmann constant, with units of J/K

Expand Down Expand Up @@ -337,9 +321,12 @@ class MultipolarContributionGrayGubbins {
}

public:
MultipolarContributionGrayGubbins(const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::MatrixXd& SIGMAIJ, const Eigen::MatrixXd& EPSKIJ, const Eigen::ArrayX<double> &mu, const Eigen::ArrayX<double> &Q, const std::optional<nlohmann::json>& flags)
MultipolarContributionGrayGubbins(const auto& sidecarJ, const auto& sidecarK, const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::MatrixXd& SIGMAIJ, const Eigen::MatrixXd& EPSKIJ, const Eigen::ArrayX<double> &mu, const Eigen::ArrayX<double> &Q, const std::optional<nlohmann::json>& flags)

: sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), SIGMAIJ(SIGMAIJ), EPSKIJ(EPSKIJ), mu(mu), Q(Q), mu2(mu.pow(2)), Q2(Q.pow(2)), Q3(Q.pow(3)), has_a_polar(Q.cwiseAbs().sum() > 0 || mu.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(get_approach(flags)), C3(get_C3(flags)), C3b(get_C3b(flags)), polarizable(get_polarizable(flags)) {
:
J6(decltype(sidecarJ){6}), J8(decltype(sidecarJ){8}), J10(decltype(sidecarJ){10}), J11(decltype(sidecarJ){11}), J13(decltype(sidecarJ){13}), J15(decltype(sidecarJ){15}),
K222_333(decltype(sidecarK){222,333}), K233_344(decltype(sidecarK){233,344}), K334_445(decltype(sidecarK){334,445}), K444_555(decltype(sidecarK){444,555}),
sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), SIGMAIJ(SIGMAIJ), EPSKIJ(EPSKIJ), mu(mu), Q(Q), mu2(mu.pow(2)), Q2(Q.pow(2)), Q3(Q.pow(3)), has_a_polar(Q.cwiseAbs().sum() > 0 || mu.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(get_approach(flags)), C3(get_C3(flags)), C3b(get_C3b(flags)), polarizable(get_polarizable(flags)) {
// Check lengths match
if (sigma_m.size() != mu.size()){
throw teqp::InvalidArgument("bad size of mu");
Expand Down Expand Up @@ -685,13 +672,8 @@ class MultipolarContributionGrayGubbins {
/// The variant containing the multipolar types that can be provided
using multipolar_contributions_variant = std::variant<
teqp::saft::polar_terms::GrossVrabec::MultipolarContributionGrossVrabec,
MultipolarContributionGrayGubbins<JGubbinsTwuSidecar, GubbinsTwuKIntegral>,
MultipolarContributionGrayGubbins<JGottschalkSidecar, GottschalkKIntegral>,
MultipolarContributionGrayGubbins<JLuckasSidecar, LuckasKIntegral>,

MultipolarContributionGubbinsTwu<JLuckasSidecar, LuckasKIntegral>,
MultipolarContributionGubbinsTwu<JGubbinsTwuSidecar, GubbinsTwuKIntegral>,
MultipolarContributionGubbinsTwu<JGottschalkSidecar, GottschalkKIntegral>
MultipolarContributionGrayGubbins,
MultipolarContributionGubbinsTwu
>;

}
Expand Down
Loading

0 comments on commit 99291b4

Please sign in to comment.