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

Merge pull request #606 from CoolProp/tab_mix

Browse Source 
Tabular mixtures work with BICUBIC
This commit is contained in:
Ian Bell
2015-04-19 22:06:26 -06:00
parent 545e8761aa 2acfa38744
commit 195707d521
13 changed files with 263 additions and 62 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
Web/coolprop/Tabular.rst
View File

@@ -214,28 +214,13 @@ The primary motivation for the use of tabular interpolation is the improvement i
More Information
----------------
The tables are stored in a zipped format using the msgpack package and miniz. To save space, the uncompressed binary tables are not stored, but if you want them to be stored as well as the compressed tables, you can do something like:
The tables are stored in a zipped format using the msgpack package and miniz. If you want to see what data is serialized in the tabular data, you can unzip and unpack into python (or other high-level languages) using something roughly like::
.. ipython::
In [0]: import CoolProp.CoolProp as CP, json
import msgpack, zlib, StringIO
In [1]: jj = json.loads(CP.get_config_as_json_string())
In [2]: jj['SAVE_RAW_TABLES'] = True
In [3]: CP.set_config_as_json_string(json.dumps(jj))
before you run your code to debug the tables. This can be useful for debugging of the data stored in the tables.
.. warning::
Make sure you delete the tables that are already there, otherwise it will entirely just load the zipped tables
The uncompressed tabled can be read back into python (or other high-level languages) using something roughly like::
with open(r'/path/to/home/.CoolProp/Tables/HelmholtzEOSBackend(R245fa)/single_phase_logph.bin','rb') as fp:
values = msgpack.load(fp)
with open(r'/path/to/home/.CoolProp/Tables/HelmholtzEOSBackend(R245fa)/single_phase_logph.bin.z','rb') as fp:
ph = zlib.decompress(fp.read())
values = msgpack.load(StringIO.StringIO(ph))
revision, matrices = values[0:2]
T,h,p,rho = np.array(matrices['T']), np.array(matrices['hmolar']), np.array(matrices['p']), np.array(matrices['rhomolar'])

12
include/CoolPropTools.h
View File

@@ -302,9 +302,9 @@
{
return (y1-y0)/(x1-x0)*(x-x0)+y0;
};
template<class T> T LinearInterp(const std::vector<T> &x, const std::vector<T> &y, std::size_t i0, std::size_t i1, T val)
template<class T1, class T2> T2 LinearInterp(const std::vector<T1> &x, const std::vector<T1> &y, std::size_t i0, std::size_t i1, T2 val)
{
return LinearInterp(x[i0],x[i1],y[i0],y[i1],val);
return LinearInterp(x[i0],x[i1],y[i0],y[i1], static_cast<T1>(val));
};
template<class T> T QuadInterp(T x0, T x1, T x2, T f0, T f1, T f2, T x)
@@ -319,9 +319,9 @@
L2=((x-x0)*(x-x1))/((x2-x0)*(x2-x1));
return L0*f0+L1*f1+L2*f2;
};
template<class T> T QuadInterp(const std::vector<T> &x, const std::vector<T> &y, std::size_t i0, std::size_t i1, std::size_t i2, T val)
template<class T1, class T2> T2 QuadInterp(const std::vector<T1> &x, const std::vector<T1> &y, std::size_t i0, std::size_t i1, std::size_t i2, T2 val)
{
return QuadInterp(x[i0],x[i1],x[i2],y[i0],y[i1],y[i2],val);
return QuadInterp(x[i0],x[i1],x[i2],y[i0],y[i1],y[i2],static_cast<T1>(val));
};
template<class T> T CubicInterp( T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x)
@@ -337,9 +337,9 @@
L3=((x-x0)*(x-x1)*(x-x2))/((x3-x0)*(x3-x1)*(x3-x2));
return L0*f0+L1*f1+L2*f2+L3*f3;
};
template<class T> T CubicInterp(const std::vector<T> &x, const std::vector<T> &y, std::size_t i0, std::size_t i1, std::size_t i2, std::size_t i3, T val)
template<class T1, class T2> T2 CubicInterp(const std::vector<T1> &x, const std::vector<T1> &y, std::size_t i0, std::size_t i1, std::size_t i2, std::size_t i3, T2 val)
{
return CubicInterp(x[i0],x[i1],x[i2],x[i3],y[i0],y[i1],y[i2],y[i3],val);
return CubicInterp(x[i0],x[i1],x[i2],x[i3],y[i0],y[i1],y[i2],y[i3],static_cast<T1>(val));
};
template<class T> T is_in_closed_range( T x1, T x2, T x)

73
include/PhaseEnvelope.h
View File

@@ -2,6 +2,10 @@
#define PHASE_ENVELOPE_H
#include "Exceptions.h"
#include "CPmsgpack.h"
#define PHASE_ENVELOPE_MATRICES X(K) X(lnK) X(x) X(y)
#define PHASE_ENVELOPE_VECTORS X(T) X(p) X(lnT) X(lnp) X(rhomolar_liq) X(rhomolar_vap) X(lnrhomolar_liq) X(lnrhomolar_vap) X(hmolar_liq) X(hmolar_vap) X(smolar_liq) X(smolar_vap) X(Q)
namespace CoolProp{
@@ -16,12 +20,25 @@ public:
std::size_t iTsat_max, ///< The index of the point corresponding to the maximum temperature for Type-I mixtures
ipsat_max, ///< The index of the point corresponding to the maximum pressure for Type-I mixtures
icrit; ///< The index of the point corresponding to the critical point
std::vector< std::vector<CoolPropDbl> > K, lnK, x, y;
std::vector<CoolPropDbl> T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap, Q;
// Use X macros to auto-generate the variables;
// each will look something like: std::vector<double> T;
#define X(name) std::vector<double> name;
PHASE_ENVELOPE_VECTORS
#undef X
// Use X macros to auto-generate the variables;
// each will look something like: std::vector<std::vector<double> > K;
#define X(name) std::vector<std::vector<double> > name;
PHASE_ENVELOPE_MATRICES
#undef X
PhaseEnvelopeData() : TypeI(false), built(false), iTsat_max(-1), ipsat_max(-1), icrit(-1) {}
std::map<std::string, std::vector<double> > vectors;
std::map<std::string, std::vector<std::vector<double> > > matrices;
MSGPACK_DEFINE(vectors, matrices); // write the member variables that you want to pack using msgpack
void resize(std::size_t N)
{
K.resize(N);
@@ -34,6 +51,52 @@ public:
lnrhomolar_liq.clear(); lnrhomolar_vap.clear(); hmolar_liq.clear(); hmolar_vap.clear(); smolar_liq.clear(); smolar_vap.clear();
K.clear(); lnK.clear(); x.clear(); y.clear(); Q.clear();
}
/// Take all the vectors that are in the class and pack them into the vectors map for easy unpacking using msgpack
void pack(){
/* Use X macros to auto-generate the packing code; each will look something like: matrices.insert(std::pair<std::string, std::vector<double> >("T", T)); */
#define X(name) vectors.insert(std::pair<std::string, std::vector<double> >(#name, name));
PHASE_ENVELOPE_VECTORS
#undef X
/* Use X macros to auto-generate the packing code; each will look something like: matrices.insert(std::pair<std::string, std::vector<std::vector<CoolPropDbl> > >("T", T)); */
#define X(name) matrices.insert(std::pair<std::string, std::vector<std::vector<double> > >(#name, name));
PHASE_ENVELOPE_MATRICES
#undef X
};
std::map<std::string, std::vector<double> >::iterator get_vector_iterator(const std::string &name){
std::map<std::string, std::vector<double> >::iterator it = vectors.find(name);
if (it == vectors.end()){
throw UnableToLoadError(format("could not find vector %s",name.c_str()));
}
return it;
}
std::map<std::string, std::vector<std::vector<double> > >::iterator get_matrix_iterator(const std::string &name){
std::map<std::string, std::vector<std::vector<double> > >::iterator it = matrices.find(name);
if (it == matrices.end()){
throw UnableToLoadError(format("could not find matrix %s", name.c_str()));
}
return it;
}
/// Take all the vectors that are in the class and unpack them from the vectors map
void unpack(){
/* Use X macros to auto-generate the unpacking code;
* each will look something like: T = get_vector_iterator("T")->second
*/
#define X(name) name = get_vector_iterator(#name)->second;
PHASE_ENVELOPE_VECTORS
#undef X
/* Use X macros to auto-generate the unpacking code;
* each will look something like: T = get_matrix_iterator("T")->second
**/
#define X(name) name = get_matrix_iterator(#name)->second;
PHASE_ENVELOPE_MATRICES
#undef X
};
void deserialize(msgpack::object &deserialized){
PhaseEnvelopeData temp;
deserialized.convert(&temp);
temp.unpack();
std::swap(*this, temp); // Swap if successful
};
void insert_variables(const CoolPropDbl T,
const CoolPropDbl p,
const CoolPropDbl rhomolar_liq,
@@ -68,10 +131,10 @@ public:
this->y[j].insert(this->y[j].begin() + i, y[j]);
}
if (rhomolar_liq > rhomolar_vap){
this->Q.insert(this->Q.begin(), 1);
this->Q.insert(this->Q.begin() + i, 1);
}
else{
this->Q.insert(this->Q.begin(), 0);
this->Q.insert(this->Q.begin() + i, 0);
}
};
void store_variables(const CoolPropDbl T,

6
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -41,7 +41,7 @@ void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS)
// Determine whether you are inside (two-phase) or outside (single-phase)
SimpleState closest_state;
std::size_t i;
bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS, iP, HEOS._p, iT, HEOS._T, i, closest_state);
bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS.PhaseEnvelope, iP, HEOS._p, iT, HEOS._T, i, closest_state);
if (!twophase && HEOS._T > closest_state.T){
// Gas solution - bounded between phase envelope temperature and very high temperature
//
@@ -467,7 +467,7 @@ void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parame
{
// Find the intersections in the phase envelope
std::vector< std::pair<std::size_t, std::size_t> > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS, other, value);
std::vector< std::pair<std::size_t, std::size_t> > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS.get_phase_envelope_data(), other, value);
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
@@ -1156,7 +1156,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
SimpleState closest_state;
std::size_t iclosest;
std::cout << format("pre is inside\n");
bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS, iP, HEOS._p, other, value, iclosest, closest_state);
bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS.PhaseEnvelope, iP, HEOS._p, other, value, iclosest, closest_state);
std::cout << format("post is inside\n");
std::string errstr;

16
src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp
View File

@@ -421,8 +421,8 @@ void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS)
IO.rhomolar_vap = rho0;
}
else if (Nsoln == 2){
if (is_in_closed_range(env.rhomolar_vap[imax-1], env.rhomolar_vap[imax+1], (CoolPropDbl)rho0)){ IO.rhomolar_vap = rho0; }
if (is_in_closed_range(env.rhomolar_vap[imax-1], env.rhomolar_vap[imax+1], (CoolPropDbl)rho1)){ IO.rhomolar_vap = rho1; }
if (is_in_closed_range(env.rhomolar_vap[imax-1], env.rhomolar_vap[imax+1], rho0)){ IO.rhomolar_vap = rho0; }
if (is_in_closed_range(env.rhomolar_vap[imax-1], env.rhomolar_vap[imax+1], rho1)){ IO.rhomolar_vap = rho1; }
}
else{
throw ValueError("More than 2 solutions found");
@@ -492,11 +492,10 @@ void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS)
env.ipsat_max = std::distance(env.p.begin(), std::max_element(env.p.begin(), env.p.end()));
}
std::vector<std::pair<std::size_t, std::size_t> > PhaseEnvelopeRoutines::find_intersections(HelmholtzEOSMixtureBackend &HEOS, parameters iInput, CoolPropDbl value)
std::vector<std::pair<std::size_t, std::size_t> > PhaseEnvelopeRoutines::find_intersections(const PhaseEnvelopeData &env, parameters iInput, double value)
{
std::vector<std::pair<std::size_t, std::size_t> > intersections;
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
for (std::size_t i = 0; i < env.p.size()-1; ++i){
bool matched = false;
switch(iInput){
@@ -518,11 +517,10 @@ std::vector<std::pair<std::size_t, std::size_t> > PhaseEnvelopeRoutines::find_in
}
return intersections;
}
bool PhaseEnvelopeRoutines::is_inside(HelmholtzEOSMixtureBackend &HEOS, parameters iInput1, CoolPropDbl value1, parameters iInput2, CoolPropDbl value2, std::size_t &iclosest, SimpleState &closest_state)
bool PhaseEnvelopeRoutines::is_inside(const PhaseEnvelopeData &env, parameters iInput1, CoolPropDbl value1, parameters iInput2, CoolPropDbl value2, std::size_t &iclosest, SimpleState &closest_state)
{
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
// Find the indices that bound the solution(s)
std::vector<std::pair<std::size_t, std::size_t> > intersections = find_intersections(HEOS, iInput1, value1);
std::vector<std::pair<std::size_t, std::size_t> > intersections = find_intersections(env, iInput1, value1);
// For now, first input must be p
if (iInput1 != iP){throw ValueError("For now, first input must be p in is_inside");}
@@ -538,8 +536,8 @@ bool PhaseEnvelopeRoutines::is_inside(HelmholtzEOSMixtureBackend &HEOS, paramete
if (intersections.size()%2 == 0){
if (intersections.size() != 2){throw ValueError("for now only even value accepted is 2"); }
std::vector<std::size_t> other_indices(4, 0);
std::vector<CoolPropDbl> *y;
std::vector<CoolPropDbl> other_values(4, 0);
std::vector<double> const *y;
std::vector<double> other_values(4, 0);
other_indices[0] = intersections[0].first; other_indices[1] = intersections[0].second;
other_indices[2] = intersections[1].first; other_indices[3] = intersections[1].second;

8
src/Backends/Helmholtz/PhaseEnvelopeRoutines.h
View File

@@ -29,15 +29,15 @@ class PhaseEnvelopeRoutines{
* can be used to determine whether another input is "inside" or "outside" the phase
* boundary.
*
* @param HEOS The HelmholtzEOSMixtureBackend instance to be used
* @param env The PhaseEnvelopeData instance to be used
* @param iInput The key for the variable type that is to be checked
* @param value The value associated with iInput
*/
static std::vector<std::pair<std::size_t, std::size_t> > find_intersections(HelmholtzEOSMixtureBackend &HEOS, parameters iInput, CoolPropDbl value);
static std::vector<std::pair<std::size_t, std::size_t> > find_intersections(const PhaseEnvelopeData &env, parameters iInput, double value);
/** \brief Determine whether a pair of inputs is inside or outside the phase envelope
*
* @param HEOS The HelmholtzEOSMixtureBackend instance to be used
* @param env The PhaseEnvelopeData instance to be used
* @param iInput1 The key for the first input
* @param value1 The value of the first input
* @param iInput2 The key for the second input
@@ -45,7 +45,7 @@ class PhaseEnvelopeRoutines{
* @param iclosest The index of the phase envelope for the closest point
* @param closest_state A SimpleState corresponding to the closest point found on the phase envelope
*/
static bool is_inside(HelmholtzEOSMixtureBackend &HEOS, parameters iInput1, CoolPropDbl value1, parameters iInput2, CoolPropDbl value2, std::size_t &iclosest, SimpleState &closest_state);
static bool is_inside(const PhaseEnvelopeData &env, parameters iInput1, CoolPropDbl value1, parameters iInput2, CoolPropDbl value2, std::size_t &iclosest, SimpleState &closest_state);
};
} /* namespace CoolProp */

11
src/Backends/REFPROP/REFPROPMixtureBackend.cpp
View File

@@ -750,15 +750,20 @@ void REFPROPMixtureBackend::calc_phase_envelope(const std::string &type)
{
double y; iderv = 0;
PhaseEnvelope.rhomolar_vap.push_back(rho_molL*1000);
PhaseEnvelope.lnrhomolar_vap.push_back(log(rho_molL*1000));
isp = nc + 1;
SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
PhaseEnvelope.T.push_back(y);
PhaseEnvelope.lnT.push_back(log(y));
isp = nc + 2;
SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
PhaseEnvelope.p.push_back(y*1000);
PhaseEnvelope.lnp.push_back(log(y*1000));
isp = nc + 3;
SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
PhaseEnvelope.rhomolar_liq.push_back(y*1000);
PhaseEnvelope.lnrhomolar_liq.push_back(log(y*1000));
PhaseEnvelope.Q.push_back(static_cast<double>(y > rho_molL));
isp = nc + 4;
SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
PhaseEnvelope.hmolar_vap.push_back(y*1000);
@@ -1366,8 +1371,8 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
_cvmolar = cvmol;
_cpmolar = cpmol;
_speed_sound = w;
_tau = calc_T_critical()/_T;
_delta = _rhomolar/calc_rhomolar_critical();
_tau = calc_T_reducing()/_T;
_delta = _rhomolar/calc_rhomolar_reducing();
_Q = q;
}
CoolPropDbl REFPROPMixtureBackend::call_phixdll(long itau, long idel)
@@ -1384,7 +1389,7 @@ CoolPropDbl REFPROPMixtureBackend::call_phi0dll(long itau, long idel)
double val = 0, tau = _tau, delta = _delta, __T = T(), __rho = rhomolar()/1000;
if (PHI0dll == NULL){throw ValueError("PHI0dll function is not available in your version of REFPROP. Please upgrade");}
PHI0dll(&itau, &idel, &__T, &__rho, &(mole_fractions[0]), &val);
return static_cast<CoolPropDbl>(val)/pow(delta,idel)/pow(tau,itau);
return static_cast<CoolPropDbl>(val)/pow(tau,itau); // Not multplied by delta^idel
}
void REFPROP_SETREF(char hrf[3], long ixflag, double x0[1], double &h0, double &s0, double &T0, double &p0, long &ierr, char herr[255], long l1, long l2){

31
src/Backends/Tabular/BicubicBackend.cpp
View File

@@ -2,6 +2,7 @@
#include "BicubicBackend.h"
#include "MatrixMath.h"
#include "PhaseEnvelopeRoutines.h"
/// The inverse of the A matrix for the bicubic interpolation (http://en.wikipedia.org/wiki/Bicubic_interpolation)
/// NOTE: The matrix is transposed below
@@ -26,6 +27,7 @@ static Eigen::Matrix<double, 16, 16> Ainv(Ainv_data);
void CoolProp::BicubicBackend::build_coeffs(SinglePhaseGriddedTableData &table, std::vector<std::vector<CellCoeffs> > &coeffs)
{
if (!coeffs.empty()){ return; }
const bool debug = get_debug_level() > 5 || false;
const int param_count = 6;
parameters param_list[param_count] = {iDmolar, iT, iSmolar, iHmolar, iP, iUmolar};
@@ -144,6 +146,15 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
// Check the tables and build if necessary
check_tables();
bool is_mixture = (this->AS->get_mole_fractions().size() >= 2);
if (is_mixture){
// For mixtures, the construction of the coefficients is delayed until this
// function so that the set_mole_fractions function can be called
build_coeffs(single_phase_logph, coeffs_ph);
build_coeffs(single_phase_logpT, coeffs_pT);
}
// Flush the cached indices (set to large number)
cached_single_phase_i = std::numeric_limits<std::size_t>::max();
cached_single_phase_j = std::numeric_limits<std::size_t>::max();
@@ -161,9 +172,15 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
else{
using_single_phase_table = true; // Use the table !
std::size_t iL, iV;
std::size_t iL, iV, iclosest = 0;
CoolPropDbl hL = 0, hV = 0;
if (pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV)){
SimpleState closest_state;
if (
(is_mixture && PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iHmolar, _hmolar, iclosest, closest_state))
||
(!is_mixture && pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV))
)
{
using_single_phase_table = false;
_Q = (static_cast<double>(_hmolar)-hL)/(hV-hL);
if(!is_in_closed_range(0.0,1.0,static_cast<double>(_Q))){
@@ -205,9 +222,15 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
else{
using_single_phase_table = true; // Use the table !
std::size_t iL = 0, iV = 0;
std::size_t iL = 0, iV = 0, iclosest = 0;
CoolPropDbl TL = 0, TV = 0;
if (pure_saturation.is_inside(iP, _p, iT, _T, iL, iV, TL, TV)){
SimpleState closest_state;
if (
(is_mixture && PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iT, _T, iclosest, closest_state))
||
(!is_mixture && pure_saturation.is_inside(iP, _p, iT, _T, iL, iV, TL, TV))
)
{
using_single_phase_table = false;
throw ValueError(format("P,T with TTSE cannot be two-phase for now"));
}

13
src/Backends/Tabular/TabularBackends.cpp
View File

@@ -55,8 +55,8 @@ template <typename T> void load_table(T &table, const std::string &path_to_table
double toc = clock();
if (get_debug_level() > -1){std::cout << format("Loaded table: %s in %g sec.", path_to_table.c_str(), (toc-tic)/CLOCKS_PER_SEC) << std::endl;}
}
catch(std::exception &){
std::string err = format("Unable to deserialize %s", path_to_table.c_str());
catch(std::exception &e){
std::string err = format("Unable to deserialize %s; err: %s", path_to_table.c_str(), e.what());
if (get_debug_level() > 0){std::cout << err << std::endl;}
throw UnableToLoadError(err);
}
@@ -149,7 +149,7 @@ void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractS
if (debug){
std::cout << format("***********************************************\n");
std::cout << format(" Single-Phase Table (%s) \n", AS->name().c_str());
std::cout << format(" Single-Phase Table (%s) \n", strjoin(AS->fluid_names(), "&").c_str());
std::cout << format("***********************************************\n");
}
// ------------------------
@@ -202,7 +202,7 @@ void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractS
}
// Skip two-phase states - they will remain as _HUGE holes in the table
if (AS->phase() == iphase_twophase){
if (is_in_closed_range(0.0, 1.0, AS->Q())){
if (debug){std::cout << " 2Phase" << std::endl;}
continue;
};
@@ -284,16 +284,21 @@ void CoolProp::TabularBackend::write_tables(){
write_table(single_phase_logph, path_to_tables, "single_phase_logph");
write_table(single_phase_logpT, path_to_tables, "single_phase_logpT");
write_table(pure_saturation, path_to_tables, "pure_saturation");
write_table(phase_envelope, path_to_tables, "phase_envelope");
}
void CoolProp::TabularBackend::load_tables(){
std::string path_to_tables = this->path_to_tables();
single_phase_logph.AS = this->AS;
single_phase_logpT.AS = this->AS;
single_phase_logph.AS->set_mole_fractions(AS->get_mole_fractions());
single_phase_logpT.AS->set_mole_fractions(AS->get_mole_fractions());
single_phase_logph.set_limits();
single_phase_logpT.set_limits();
load_table(single_phase_logph, path_to_tables, "single_phase_logph.bin.z");
load_table(single_phase_logpT, path_to_tables, "single_phase_logpT.bin.z");
load_table(pure_saturation, path_to_tables, "pure_saturation.bin.z");
load_table(phase_envelope, path_to_tables, "phase_envelope.bin.z");
}
#endif // !defined(NO_TABULAR_BACKENDS)

15
src/Backends/Tabular/TabularBackends.h
View File

@@ -512,6 +512,7 @@ class TabularBackend : public AbstractState
LogPHTable single_phase_logph;
LogPTTable single_phase_logpT;
PureFluidSaturationTableData pure_saturation; // This will ultimately be split into pure and mixture backends which derive from this backend
PhaseEnvelopeData phase_envelope;
bool using_mole_fractions(void){return true;}
bool using_mass_fractions(void){return false;}
@@ -534,7 +535,18 @@ class TabularBackend : public AbstractState
void load_tables();
/// Build the tables
void build_tables(){
pure_saturation.build(AS);
// Pure or pseudo-pure fluid
if (AS->get_mole_fractions().size() == 1){
pure_saturation.build(AS);
}
else{
// Call function to actually construct the phase envelope
AS->build_phase_envelope("");
// Copy constructed phase envelope into this class
phase_envelope = AS->get_phase_envelope_data();
// Resize so that it will load properly
pure_saturation.resize(pure_saturation.N);
}
single_phase_logph.build(AS);
single_phase_logpT.build(AS);
}
@@ -542,6 +554,7 @@ class TabularBackend : public AbstractState
single_phase_logph.pack();
single_phase_logpT.pack();
pure_saturation.pack();
phase_envelope.pack();
}
/// Write the tables to file
void write_tables();

24
wrappers/Python/CoolProp/AbstractState.pxd
View File

@@ -94,4 +94,26 @@ cdef class AbstractState:
cpdef PyPhaseEnvelopeData get_phase_envelope_data(self)
cpdef mole_fractions_liquid(self)
cpdef mole_fractions_vapor(self)
cpdef mole_fractions_vapor(self)
cpdef long double alpha0(self) except *
cpdef long double dalpha0_dDelta(self) except *
cpdef long double dalpha0_dTau(v) except *
cpdef long double d2alpha0_dDelta2(self) except *
cpdef long double d2alpha0_dDelta_dTau(self) except *
cpdef long double d2alpha0_dTau2(self) except *
cpdef long double d3alpha0_dTau3(self) except *
cpdef long double d3alpha0_dDelta_dTau2(self) except *
cpdef long double d3alpha0_dDelta2_dTau(self) except *
cpdef long double d3alpha0_dDelta3(self) except *
cpdef long double alphar(self) except *
cpdef long double dalphar_dDelta(self) except *
cpdef long double dalphar_dTau(self) except *
cpdef long double d2alphar_dDelta2(self) except *
cpdef long double d2alphar_dDelta_dTau(self) except *
cpdef long double d2alphar_dTau2(self) except *
cpdef long double d3alphar_dDelta3(self) except *
cpdef long double d3alphar_dDelta2_dTau(self) except *
cpdef long double d3alphar_dDelta_dTau2(self) except *
cpdef long double d3alphar_dTau3(self) except *

67
wrappers/Python/CoolProp/AbstractState.pyx
View File

@@ -244,4 +244,69 @@ cdef class AbstractState:
pe_out.ipsat_max = pe_data.ipsat_max
pe_out.TypeI = pe_data.TypeI
return pe_out
## -----------------------------------------
## Helmholtz energy derivatives
## -----------------------------------------
cpdef long double alpha0(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::alpha0` """
return self.thisptr.alpha0()
cpdef long double dalpha0_dDelta(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::dalpha0_dDelta` """
return self.thisptr.dalpha0_dDelta()
cpdef long double dalpha0_dTau(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::dalpha0_dTau` """
return self.thisptr.dalpha0_dTau()
cpdef long double d2alpha0_dDelta2(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d2alpha0_dDelta2` """
return self.thisptr.d2alpha0_dDelta2()
cpdef long double d2alpha0_dDelta_dTau(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d2alpha0_dDelta_dTau` """
return self.thisptr.d2alpha0_dDelta_dTau()
cpdef long double d2alpha0_dTau2(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d2alpha0_dTau2` """
return self.thisptr.d2alpha0_dTau2()
cpdef long double d3alpha0_dTau3(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alpha0_dTau3` """
return self.thisptr.d3alpha0_dTau3()
cpdef long double d3alpha0_dDelta_dTau2(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alpha0_dDelta_dTau2` """
return self.thisptr.d3alpha0_dDelta_dTau2()
cpdef long double d3alpha0_dDelta2_dTau(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alpha0_dDelta2_dTau` """
return self.thisptr.d3alpha0_dDelta2_dTau()
cpdef long double d3alpha0_dDelta3(self) except *:
""" Get the ideal-gas reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alpha0_dDelta3` """
return self.thisptr.d3alpha0_dDelta3()
cpdef long double alphar(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::alphar` """
return self.thisptr.alphar()
cpdef long double dalphar_dDelta(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::dalphar_dDelta` """
return self.thisptr.dalphar_dDelta()
cpdef long double dalphar_dTau(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::dalphar_dTau` """
return self.thisptr.dalphar_dTau()
cpdef long double d2alphar_dDelta2(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d2alphar_dDelta2` """
return self.thisptr.d2alphar_dDelta2()
cpdef long double d2alphar_dDelta_dTau(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d2alphar_dDelta_dTau` """
return self.thisptr.d2alphar_dDelta_dTau()
cpdef long double d2alphar_dTau2(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d2alphar_dTau2` """
return self.thisptr.d2alphar_dTau2()
cpdef long double d3alphar_dTau3(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alphar_dTau3` """
return self.thisptr.d3alphar_dTau3()
cpdef long double d3alphar_dDelta_dTau2(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alphar_dDelta_dTau2` """
return self.thisptr.d3alphar_dDelta_dTau2()
cpdef long double d3alphar_dDelta2_dTau(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alphar_dDelta2_dTau` """
return self.thisptr.d3alphar_dDelta2_dTau()
cpdef long double d3alphar_dDelta3(self) except *:
""" Get the residual reduced Helmholtz energy - wrapper of c++ function :cpapi:`CoolProp::AbstractState::d3alphar_dDelta3` """
return self.thisptr.d3alphar_dDelta3()

24
wrappers/Python/CoolProp/cAbstractState.pxd
View File

@@ -8,7 +8,7 @@ cdef extern from "PhaseEnvelope.h" namespace "CoolProp":
cdef cppclass PhaseEnvelopeData:
bool TypeI
size_t iTsat_max, ipsat_max, icrit
vector[long double] T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap, Q
vector[double] T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap, Q
cdef extern from "AbstractState.h" namespace "CoolProp":
@@ -105,6 +105,28 @@ cdef extern from "AbstractState.h" namespace "CoolProp":
double build_phase_envelope() except+ValueError
void build_phase_envelope(string) except+ValueError
PhaseEnvelopeData get_phase_envelope_data() except+ValueError
long double alpha0() except+ValueError
long double dalpha0_dDelta() except+ValueError
long double dalpha0_dTau() except+ValueError
long double d2alpha0_dDelta2() except+ValueError
long double d2alpha0_dDelta_dTau() except+ValueError
long double d2alpha0_dTau2() except+ValueError
long double d3alpha0_dTau3() except+ValueError
long double d3alpha0_dDelta_dTau2() except+ValueError
long double d3alpha0_dDelta2_dTau() except+ValueError
long double d3alpha0_dDelta3() except+ValueError
long double alphar() except+ValueError
long double dalphar_dDelta() except+ValueError
long double dalphar_dTau() except+ValueError
long double d2alphar_dDelta2() except+ValueError
long double d2alphar_dDelta_dTau() except+ValueError
long double d2alphar_dTau2() except+ValueError
long double d3alphar_dDelta3() except+ValueError
long double d3alphar_dDelta2_dTau() except+ValueError
long double d3alphar_dDelta_dTau2() except+ValueError
long double d3alphar_dTau3() except+ValueError
# The static factory method for the AbstractState
cdef extern from "AbstractState.h" namespace "CoolProp::AbstractState":