github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-23 03:00:17 -04:00
Code Issues Packages Projects Releases Wiki Activity

Refactor UNIFAQ code; first steps toward #1310 (#1312)

Browse Source 
* Add derivatives of gE_R

* Finally found this const issue

* Update TEST_CASE

* Move ln_gamma_C into activity_coefficients
This speeds up the computation of VTPR who do not need that part to run.
This commit is contained in:
JonWel
2016-11-05 03:43:30 +01:00
committed by Ian Bell
parent 4ca3865775
commit 1a618a5d8f
4 changed files with 117 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
src/Backends/Cubics/UNIFAQ.cpp
View File

@@ -25,34 +25,9 @@ void UNIFAQ::UNIFAQMixture::set_mole_fractions(const std::vector<double> &z) {
pure_data.clear();
this->mole_fractions = z;
std::size_t N = z.size();
std::vector<double> &r = m_r, &q = m_q, &l = m_l, &phi = m_phi, &theta = m_theta, &ln_Gamma_C = m_ln_Gamma_C;
r.resize(N); q.resize(N); l.resize(N); phi.resize(N); theta.resize(N); ln_Gamma_C.resize(N);
double summerzr = 0, summerzq = 0, summerzl = 0;
for (std::size_t i = 0; i < z.size(); ++i) {
double summerr = 0, summerq = 0;
const UNIFAQLibrary::Component &c = components[i];
for (std::size_t j = 0; j < c.groups.size(); ++j) {
const UNIFAQLibrary::ComponentGroup &cg = c.groups[j];
summerr += cg.count*cg.group.R_k;
summerq += cg.count*cg.group.Q_k;
}
r[i] = summerr;
q[i] = summerq;
summerzr += z[i] * r[i];
summerzq += z[i] * q[i];
}
for (std::size_t i = 0; i < z.size(); ++i) {
phi[i] = z[i] * r[i] / summerzr;
theta[i] = z[i] * q[i] / summerzq;
l[i] = 10.0 / 2.0*(r[i] - q[i]) - (r[i] - 1);
summerzl += z[i] * l[i];
}
for (std::size_t i = 0; i < z.size(); ++i) {
ln_Gamma_C[i] = log(phi[i]/z[i]) + 10.0/2.0*q[i] * log(theta[i]/phi[i]) + l[i] - phi[i]/z[i]*summerzl;
}
/// Calculate the parameters X and theta for the pure components, which does not depend on temperature
for (std::size_t i = 0; i < z.size(); ++i){
for (std::size_t i = 0; i < N; ++i){
int totalgroups = 0;
const UNIFAQLibrary::Component &c = components[i];
ComponentData cd;
@@ -80,7 +55,7 @@ void UNIFAQ::UNIFAQMixture::set_mole_fractions(const std::vector<double> &z) {
}
pure_data.push_back(cd);
}
for (std::size_t i = 0; i < z.size(); ++i) {
for (std::size_t i = 0; i < N; ++i) {
//printf("%g %g %g %g %g %g\n", l[i], phi[i], q[i], r[i], theta[i], ln_Gamma_C[i]);
}
@@ -153,9 +128,9 @@ double UNIFAQ::UNIFAQMixture::theta_pure(std::size_t i, std::size_t sgi) const {
return pure_data[i].theta.find(sgi)->second;
}
void UNIFAQ::UNIFAQMixture::set_temperature(const double T, const std::vector<double> &z){
// // Check whether you are using exactly the same temperature and mole fractions as last time
// if (static_cast<bool>(_T) && std::abs(static_cast<double>(_T) - T) < 1e-15 && maxvectordiff(z, mole_fractions) < 1e-15){
void UNIFAQ::UNIFAQMixture::set_temperature(const double T){
// // Check whether you are using exactly the same temperature as last time
// if (static_cast<bool>(_T) && std::abs(static_cast<double>(_T) - T) < 1e-15 {
// //
// return;
// }
@@ -212,20 +187,52 @@ void UNIFAQ::UNIFAQMixture::set_temperature(const double T, const std::vector<do
}
_T = m_T;
}
double UNIFAQ::UNIFAQMixture::ln_gamma_R(std::size_t i) const{
double summer = 0;
for (std::vector<UNIFAQLibrary::Group>::const_iterator it = unique_groups.begin(); it != unique_groups.end(); ++it) {
std::size_t k = it->sgi;
std::size_t count = group_count(i, k);
if (count > 0){
summer += count*(m_lnGammag.find(k)->second - pure_data[i].lnGamma.find(k)->second);
double UNIFAQ::UNIFAQMixture::ln_gamma_R(const double tau, std::size_t i, std::size_t itau){
if (itau == 0) {
set_temperature(T_r / tau);
double summer = 0;
for (std::vector<UNIFAQLibrary::Group>::const_iterator it = unique_groups.begin(); it != unique_groups.end(); ++it) {
std::size_t k = it->sgi;
std::size_t count = group_count(i, k);
if (count > 0){
summer += count*(m_lnGammag.find(k)->second - pure_data[i].lnGamma.find(k)->second);
}
}
//printf("log(gamma)_{%d}: %g\n", i+1, summer);
return summer;
}
else {
double dtau = 0.01*tau;
return (ln_gamma_R(tau + dtau, i, itau - 1) - ln_gamma_R(tau - dtau, i, itau - 1)) / (2 * dtau);
}
//printf("log(gamma)_{%d}: %g\n", i+1, summer);
return summer;
}
double UNIFAQ::UNIFAQMixture::activity_coefficient(std::size_t i) const {
return exp(ln_gamma_R(i) + m_ln_Gamma_C[i]);
void UNIFAQ::UNIFAQMixture::activity_coefficients(double tau, const std::vector<double> &z, std::vector<double> &gamma){
std::size_t N = z.size();
std::vector<double> r(N), q(N), l(N), phi(N), theta(N), ln_Gamma_C(N);
double summerzr = 0, summerzq = 0, summerzl = 0;
for (std::size_t i = 0; i < N; ++i) {
double summerr = 0, summerq = 0;
const UNIFAQLibrary::Component &c = components[i];
for (std::size_t j = 0; j < c.groups.size(); ++j) {
const UNIFAQLibrary::ComponentGroup &cg = c.groups[j];
summerr += cg.count*cg.group.R_k;
summerq += cg.count*cg.group.Q_k;
}
r[i] = summerr;
q[i] = summerq;
summerzr += z[i] * r[i];
summerzq += z[i] * q[i];
}
for (std::size_t i = 0; i < N; ++i) {
phi[i] = z[i] * r[i] / summerzr;
theta[i] = z[i] * q[i] / summerzq;
l[i] = 10.0 / 2.0*(r[i] - q[i]) - (r[i] - 1);
summerzl += z[i] * l[i];
}
for (std::size_t i = 0; i < N; ++i) {
ln_Gamma_C[i] = log(phi[i] / z[i]) + 10.0 / 2.0*q[i] * log(theta[i] / phi[i]) + l[i] - phi[i] / z[i] * summerzl;
gamma[i] = exp(ln_gamma_R(tau, i, 0) + ln_Gamma_C[i]);
}
}
/// Add a component with the defined groups defined by (count, sgi) pairs

20
src/Backends/Cubics/UNIFAQ.h
View File

@@ -20,13 +20,7 @@ namespace UNIFAQ
CoolProp::CachedElement _T; ///< The cached temperature
double m_T; ///< The temperature in K
std::vector<double> m_r,
m_q,
m_l,
m_phi,
m_theta,
m_ln_Gamma_C;
double T_r; ///< Reduce temperature
std::map<std::size_t, double> m_Xg, ///< Map from sgi to mole fraction of group in the mixture
m_thetag, ///< Map from sgi to theta for the group in the mixture
@@ -52,7 +46,7 @@ namespace UNIFAQ
public:
UNIFAQMixture(const UNIFAQLibrary::UNIFAQParameterLibrary &library) : library(library) {};
UNIFAQMixture(const UNIFAQLibrary::UNIFAQParameterLibrary &library, const double T_r) : library(library), T_r(T_r) {};
/**
* \brief Set all the interaction parameters between groups
@@ -68,8 +62,8 @@ namespace UNIFAQ
/// Get the mole fractions of the components in the mixtures (not the groups)
const std::vector<double> & get_mole_fractions() { return mole_fractions; }
/// Set the mole fractions of the components in the mixtures (not the groups) AND the mole fractions you want to use
void set_temperature(const double T, const std::vector<double> &z);
/// Set the temperature of the components in the mixtures (not the groups)
void set_temperature(const double T);
/// Get the temperature
double get_temperature() const { return m_T; }
@@ -78,9 +72,9 @@ namespace UNIFAQ
double theta_pure(std::size_t i, std::size_t sgi) const;
double activity_coefficient(std::size_t i) const;
void activity_coefficients(double tau, const std::vector<double> &z, std::vector<double> &gamma);
double ln_gamma_R(std::size_t i) const;
double ln_gamma_R(const double tau, std::size_t i, std::size_t itau);
std::size_t group_count(std::size_t i, std::size_t sgi) const;
@@ -94,4 +88,4 @@ namespace UNIFAQ
} /* namespace UNIFAQ */
#endif
#endif

20
src/Backends/Cubics/UNIFAQLibrary.cpp
View File

@@ -144,28 +144,28 @@ TEST_CASE("Check Poling example for UNIFAQ", "[UNIFAQ]")
SECTION("Validate AC for acetone + n-pentane") {
UNIFAQLibrary::UNIFAQParameterLibrary lib;
CHECK_NOTHROW(lib.populate(groups, interactions, acetone_pentane_groups););
UNIFAQ::UNIFAQMixture mix(lib);
UNIFAQ::UNIFAQMixture mix(lib,1.0);
std::vector<std::string> names; names.push_back("Acetone"); names.push_back("n-Pentane");
mix.set_components("name",names);
mix.set_interaction_parameters();
std::vector<double> z(2,0.047); z[1] = 1-z[0];
mix.set_mole_fractions(z);
CHECK_NOTHROW(mix.set_temperature(307, z););
CHECK_NOTHROW(mix.set_temperature(307););
double lngammaR0 = mix.ln_gamma_R(0);
double lngammaR1 = mix.ln_gamma_R(1);
double lngammaR0 = mix.ln_gamma_R(1.0/307,0,0);
double lngammaR1 = mix.ln_gamma_R(1.0/307,1,0);
CAPTURE(lngammaR0);
CAPTURE(lngammaR1);
CHECK(std::abs(lngammaR0 - 1.66) < 1e-2);
CHECK(std::abs(lngammaR1 - 5.68e-3) < 1e-3);
double gamma0 = mix.activity_coefficient(0);
double gamma1 = mix.activity_coefficient(1);
CAPTURE(gamma0);
CAPTURE(gamma1);
CHECK(std::abs(gamma0 - 4.99) < 1e-2);
CHECK(std::abs(gamma1 - 1.005) < 1e-3);
std::vector<double> gamma(2);
mix.activity_coefficients(1.0/307,z,gamma);
CAPTURE(gamma[0]);
CAPTURE(gamma[1]);
CHECK(std::abs(gamma[0] - 4.99) < 1e-2);
CHECK(std::abs(gamma[1] - 1.005) < 1e-3);
};
};

75
src/Backends/Cubics/VTPRCubic.h
View File

@@ -7,6 +7,7 @@
//
#include "GeneralizedCubic.h"
#include "Exceptions.h"
#ifndef VTPRCubic_h
#define VTPRCubic_h
@@ -22,33 +23,33 @@ public:
double R_u,
const UNIFAQLibrary::UNIFAQParameterLibrary & lib
)
: PengRobinson(Tc, pc, acentric, R_u), unifaq(lib) {};
: PengRobinson(Tc, pc, acentric, R_u), unifaq(lib,T_r) {};
VTPRCubic(double Tc,
double pc,
double acentric,
double R_u,
const UNIFAQLibrary::UNIFAQParameterLibrary & lib)
: PengRobinson(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u), unifaq(lib) {};
: PengRobinson(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u), unifaq(lib,T_r) {};
/// Get a reference to the managed UNIFAQ instance
UNIFAQ::UNIFAQMixture &get_unifaq() { return unifaq; }
/// Calculate the non-dimensionalized gE/RT term
double gE_R_RT(const std::vector<double> &x) {
double gE_R_RT(double tau, const std::vector<double> &x, std::size_t itau) {
double summer = 0;
for (std::size_t i = 0; i < x.size(); ++i) {
summer += x[i] * unifaq.ln_gamma_R(i);
summer += x[i] * unifaq.ln_gamma_R(tau, i, itau);
}
return summer;
}
double d_gE_R_RT_dxi(const std::vector<double> &x, std::size_t i, bool xN_independent) {
double d_gE_R_RT_dxi(double tau, const std::vector<double> &x, std::size_t itau, std::size_t i, bool xN_independent) {
if (xN_independent)
{
return unifaq.ln_gamma_R(i);
return unifaq.ln_gamma_R(tau, i, itau);
}
else {
return unifaq.ln_gamma_R(i) - unifaq.ln_gamma_R(N - 1);
return unifaq.ln_gamma_R(tau, i, itau) - unifaq.ln_gamma_R(tau, N-1, itau);
}
}
double gE_R(double tau, const std::vector<double> &x, std::size_t itau) {
@@ -56,13 +57,28 @@ public:
return 0.;
}
else {
if (itau == 0) {
set_temperature(T_r / tau, x);
return R_u*T_r / tau*gE_R_RT(x);
}
else {
double dtau = 0.01*tau;
return (gE_R(tau + dtau, x, itau - 1) - gE_R(tau - dtau, x, itau - 1)) / (2 * dtau);
switch (itau){
case 0:
{
return R_u*T_r / tau*gE_R_RT(tau,x,0);
}
case 1:
{
return R_u*T_r / tau*(-gE_R_RT(tau,x,0)/tau + gE_R_RT(tau,x,1));
}
case 2:
{
return R_u*T_r / tau*( 2*gE_R_RT(tau,x,0)/powInt(tau, 2) - 2*gE_R_RT(tau,x,1)/tau + gE_R_RT(tau,x,2));
}
case 3:
{
return R_u*T_r / tau*(-6*gE_R_RT(tau,x,0)/powInt(tau, 3) + 6*gE_R_RT(tau,x,1)/powInt(tau, 2) - 3*gE_R_RT(tau,x,2)/tau + gE_R_RT(tau,x,3));
}
case 4:
{
return R_u*T_r / tau*(24*gE_R_RT(tau,x,0)/powInt(tau, 4) - 24*gE_R_RT(tau,x,1)/powInt(tau, 3) + 12*gE_R_RT(tau,x,2)/powInt(tau, 2) - 4*gE_R_RT(tau,x,3)/tau + gE_R_RT(tau,x,4));
}
default: throw CoolProp::ValueError(format("itau (%d) is invalid",itau));
}
}
}
@@ -71,13 +87,28 @@ public:
return 0.;
}
else {
if (itau == 0) {
set_temperature(T_r / tau, x);
return R_u*T_r / tau*d_gE_R_RT_dxi(x, i, xN_independent);
}
else {
double dtau = 0.01*tau;
return (d_gE_R_dxi(tau + dtau, x, itau - 1, i, xN_independent) - d_gE_R_dxi(tau - dtau, x, itau - 1, i, xN_independent)) / (2 * dtau);
switch (itau){
case 0:
{
return R_u*T_r / tau*d_gE_R_RT_dxi(tau,x,0,i,xN_independent);
}
case 1:
{
return R_u*T_r / tau*(-d_gE_R_RT_dxi(tau,x,0,i,xN_independent)/tau + d_gE_R_RT_dxi(tau,x,1,i,xN_independent));
}
case 2:
{
return R_u*T_r / tau*( 2*d_gE_R_RT_dxi(tau,x,0,i,xN_independent)/powInt(tau, 2) - 2*d_gE_R_RT_dxi(tau,x,1,i,xN_independent)/tau + d_gE_R_RT_dxi(tau,x,2,i,xN_independent));
}
case 3:
{
return R_u*T_r / tau*(-6*d_gE_R_RT_dxi(tau,x,0,i,xN_independent)/powInt(tau, 3) + 6*d_gE_R_RT_dxi(tau,x,1,i,xN_independent)/powInt(tau, 2) - 3*d_gE_R_RT_dxi(tau,x,2,i,xN_independent)/tau + d_gE_R_RT_dxi(tau,x,3,i,xN_independent));
}
case 4:
{
return R_u*T_r / tau*(24*d_gE_R_RT_dxi(tau,x,0,i,xN_independent)/powInt(tau, 4) - 24*d_gE_R_RT_dxi(tau,x,1,i,xN_independent)/powInt(tau, 3) + 12*d_gE_R_RT_dxi(tau,x,2,i,xN_independent)/powInt(tau, 2) - 4*d_gE_R_RT_dxi(tau,x,3,i,xN_independent)/tau + d_gE_R_RT_dxi(tau,x,4,i,xN_independent));
}
default: throw CoolProp::ValueError(format("itau (%d) is invalid",itau));
}
}
}
@@ -164,8 +195,6 @@ public:
{
return 0;
}
void set_temperature(const double T, const std::vector<double> &z) { unifaq.set_temperature(T, z); }
};
#endif /* VTPRCubic_h */