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
Files
267d64533a9682377532c98b79ee84f59a387290
CoolProp/src/AbstractState.cpp
Ian Bell 267d64533a Superancillaries for pure fluids (#2511) 
* Expansions are fully wrapped, looking good. Next step is the set of expansions that is the 1D approximation

* Get 1D approx working via cython

* Count solutions

* SuperAncillary class is working

>1000x speedup for water
Time for C++!

* Superancillaries are working!

In C++, speedup is more than 2000x. In Python, more like 150x because of Python <-> C++ overhead

* Add pmax check for PQ superancillary calls

* Update tests

* Allow T limits to be obtained

* Implement get_fluid_parameter_double for getting superanc value

* Add tests for getting parameters from superanc

* Script for testing superancillaries for sphinx

* Microoptimizations; don't help speed

The limiting factor remains the clear function, which takes about 30 ns

* Add R125 superancillary

* Use the release from fastchebpure for the files

* Drop a .gitignore in the unzipped folder

* Update superancillary injection script

* Turn on superancillaries by default

* Missing header

* Many int conversions in superancillary

* Another int cast

* More annoying solution for boost iter max

* Fix warnings

* One more warning

* Clear up the calculation of rho

* Update docs_docker-build.yml

Use arm64 since the containers were built on mac

* Superfluous ;

* Update backend.py

* Get the critical points working for superancillaries

* Fix wrapping changes of xmin&xmax methods

* squelch warnings

* Version 0 of jupyter notebook for docs

* Try to add the notebook to the docs

* Add jupyter notebook for superancillary

* Lots of updates to superancillary notebook

* More updates to docs

* Skip pseudo-pure for superancillary docs

* Fix output of superancillary figures

* Add superancillary plots to docs for the page for each fluid

* Make a placeholder figure for fluids without superancillary

* Add superancillary plots to task list

* Bump to release fixing m-xylene

* Relax the location of the REFPROP stuff

* Change default name for R-1336mzz(E)

* No need for figures to be so large

* Don't need REFPROP setting

* Bump to fastchebpure release with methanol

* Benchmark caching options

* Benchmark more granularly

* Add the fast methods to public API for HEOS class

* Back to memset - can memset with 0 but no other value

* Fix how caching is managed in Helmholtz class

* Close to final implementation

Perhaps a tiny bit more optimization possible?

* Update function name

* Make message more accurate

* Fix init order

* Expose update_QT_pure_superanc to Python

* Fix when _reducing is set for pures

* Fix the post_update

* Indent

* Notebook

* Notebook

* Make ln(p) construction lazy

Only really matters for debug builds

* Also make reference non-const

* Inject superancillary for methanol

* Make the superancillary loading entirely lazy in debug

* Fix PH bug for Nitrogen

 Closes #2470

* Force the clear to be called on SatL and SatV

To invalidate them at start

* Default is non-lazy superancillary loading

* Add CMake option to have lazy-loading superancillaries [skip ci]

Not a good idea unless doing very narrow testing
2025-05-17 20:27:19 -04:00

1145 lines
49 KiB
C++
Raw Blame History

/*
* AbstractState.cpp
*
* Created on: 21 Dec 2013
* Author: jowr
*/
#ifndef _CRT_SECURE_NO_WARNINGS
#define _CRT_SECURE_NO_WARNINGS
#endif
#include <stdlib.h>
#include "math.h"
#include "AbstractState.h"
#include "DataStructures.h"
#include "Backends/IF97/IF97Backend.h"
#include "Backends/Cubics/CubicBackend.h"
#include "Backends/Cubics/VTPRBackend.h"
#include "Backends/Incompressible/IncompressibleBackend.h"
#include "Backends/PCSAFT/PCSAFTBackend.h"
#if !defined(NO_TABULAR_BACKENDS)
#include "Backends/Tabular/TTSEBackend.h"
#include "Backends/Tabular/BicubicBackend.h"
#endif
namespace CoolProp {
/// This tiny class holds pointers to generators for the backends and can be used to look up
/// generators at runtime. This class should be populated through the use of static initialized
class BackendLibrary
{
private:
std::map<backend_families, shared_ptr<AbstractStateGenerator>> backends;
public:
void add_backend(const backend_families& bg, const shared_ptr<AbstractStateGenerator>& asg) {
backends[bg] = asg;
};
void get_generator_iterators(const backend_families& bg,
std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& generator,
std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& end) {
generator = backends.find(bg);
end = backends.end();
};
std::size_t size() {
return backends.size();
};
};
inline BackendLibrary& get_backend_library() {
static BackendLibrary the_library;
return the_library;
}
void register_backend(const backend_families& bf, shared_ptr<AbstractStateGenerator> gen) {
get_backend_library().add_backend(bf, gen);
};
class IF97BackendGenerator : public AbstractStateGenerator
{
public:
AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
if (fluid_names.size() == 1) { // Check that fluid_names[0] has only one component
std::string str = fluid_names[0]; // Check that the fluid name is an alias for "Water"
if ((upper(str) == "WATER") || (upper(str) == "H2O")) {
return new IF97Backend();
} else {
throw ValueError(format("The IF97 backend returns Water props only; fluid name [%s] not allowed", fluid_names[0].c_str()));
}
} else {
throw ValueError(format("The IF97 backend does not support mixtures, only Water"));
};
};
};
// This static initialization will cause the generator to register
static GeneratorInitializer<IF97BackendGenerator> if97_gen(IF97_BACKEND_FAMILY);
class SRKGenerator : public AbstractStateGenerator
{
public:
AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
return new SRKBackend(fluid_names, get_config_double(R_U_CODATA));
};
};
static GeneratorInitializer<SRKGenerator> srk_gen(CoolProp::SRK_BACKEND_FAMILY);
class PRGenerator : public AbstractStateGenerator
{
public:
AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
return new PengRobinsonBackend(fluid_names, get_config_double(R_U_CODATA));
};
};
static GeneratorInitializer<PRGenerator> pr_gen(CoolProp::PR_BACKEND_FAMILY);
class IncompressibleBackendGenerator : public AbstractStateGenerator
{
public:
AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
if (fluid_names.size() != 1) {
throw ValueError(format("For INCOMP backend, name vector must be one element long"));
}
return new IncompressibleBackend(fluid_names[0]);
};
};
// This static initialization will cause the generator to register
static GeneratorInitializer<IncompressibleBackendGenerator> incomp_gen(INCOMP_BACKEND_FAMILY);
class VTPRGenerator : public CoolProp::AbstractStateGenerator
{
public:
CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
return new CoolProp::VTPRBackend(fluid_names, CoolProp::get_config_double(R_U_CODATA));
};
};
// This static initialization will cause the generator to register
static CoolProp::GeneratorInitializer<VTPRGenerator> vtpr_gen(CoolProp::VTPR_BACKEND_FAMILY);
class PCSAFTGenerator : public CoolProp::AbstractStateGenerator
{
public:
CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
return new CoolProp::PCSAFTBackend(fluid_names);
};
};
// This static initialization will cause the generator to register
static CoolProp::GeneratorInitializer<PCSAFTGenerator> pcsaft_gen(CoolProp::PCSAFT_BACKEND_FAMILY);
AbstractState* AbstractState::factory(const std::string& backend, const std::vector<std::string>& fluid_names) {
if (get_debug_level() > 0) {
std::cout << "AbstractState::factory(" << backend << "," << stringvec_to_string(fluid_names) << ")" << std::endl;
}
backend_families f1;
std::string f2;
extract_backend_families_string(backend, f1, f2);
std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator gen, end;
get_backend_library().get_generator_iterators(f1, gen, end);
if (get_debug_level() > 0) {
std::cout << "AbstractState::factory backend_library size: " << get_backend_library().size() << std::endl;
}
if (gen != end) {
// One of the registered backends was able to match the given backend family
return gen->second->get_AbstractState(fluid_names);
}
#if !defined(NO_TABULAR_BACKENDS)
else if (f1 == TTSE_BACKEND_FAMILY) {
// Will throw if there is a problem with this backend
shared_ptr<AbstractState> AS(factory(f2, fluid_names));
return new TTSEBackend(AS);
} else if (f1 == BICUBIC_BACKEND_FAMILY) {
// Will throw if there is a problem with this backend
shared_ptr<AbstractState> AS(factory(f2, fluid_names));
return new BicubicBackend(AS);
}
#endif
else if (!backend.compare("?") || backend.empty()) {
std::size_t idel = fluid_names[0].find("::");
// Backend has not been specified, and we have to figure out what the backend is by parsing the string
if (idel == std::string::npos) // No '::' found, no backend specified, try HEOS, otherwise a failure
{
// Figure out what backend to use
return factory("HEOS", fluid_names);
} else {
// Split string at the '::' into two std::string, call again
return factory(std::string(fluid_names[0].begin(), fluid_names[0].begin() + idel),
std::string(fluid_names[0].begin() + (idel + 2), fluid_names[0].end()));
}
} else {
throw ValueError(format("Invalid backend name [%s] to factory function", backend.c_str()));
}
}
std::vector<std::string> AbstractState::fluid_names(void) {
return calc_fluid_names();
}
bool AbstractState::clear_comp_change() {
// Reset all instances of CachedElement and overwrite
// the internal double values with -_HUGE
this->_molar_mass.clear();
this->_critical.fill(_HUGE);
this->_reducing.fill(_HUGE);
calc_reducing_state();
this->_gas_constant = calc_gas_constant();
return true;
}
bool AbstractState::clear() {
// Reset all instances of CachedElement
cache.clear();
this->_critical.fill(_HUGE);
/// Bulk values
this->_rhomolar = -_HUGE;
this->_T = -_HUGE;
this->_p = -_HUGE;
this->_Q = -_HUGE;
return true;
}
void AbstractState::mass_to_molar_inputs(CoolProp::input_pairs& input_pair, CoolPropDbl& value1, CoolPropDbl& value2) {
// Check if a mass based input, convert it to molar units
switch (input_pair) {
case DmassT_INPUTS: ///< Mass density in kg/m^3, Temperature in K
//case HmassT_INPUTS: ///< Enthalpy in J/kg, Temperature in K (NOT CURRENTLY IMPLEMENTED)
case SmassT_INPUTS: ///< Entropy in J/kg/K, Temperature in K
//case TUmass_INPUTS: ///< Temperature in K, Internal energy in J/kg (NOT CURRENTLY IMPLEMENTED)
case DmassP_INPUTS: ///< Mass density in kg/m^3, Pressure in Pa
case DmassQ_INPUTS: ///< Mass density in kg/m^3, molar quality
case HmassP_INPUTS: ///< Enthalpy in J/kg, Pressure in Pa
case PSmass_INPUTS: ///< Pressure in Pa, Entropy in J/kg/K
case PUmass_INPUTS: ///< Pressure in Pa, Internal energy in J/kg
case HmassSmass_INPUTS: ///< Enthalpy in J/kg, Entropy in J/kg/K
case SmassUmass_INPUTS: ///< Entropy in J/kg/K, Internal energy in J/kg
case DmassHmass_INPUTS: ///< Mass density in kg/m^3, Enthalpy in J/kg
case DmassSmass_INPUTS: ///< Mass density in kg/m^3, Entropy in J/kg/K
case DmassUmass_INPUTS: ///< Mass density in kg/m^3, Internal energy in J/kg
{
// Set the cache value for the molar mass if it hasn't been set yet
molar_mass();
// Molar mass (just for compactness of the following switch)
CoolPropDbl mm = static_cast<CoolPropDbl>(_molar_mass);
switch (input_pair) {
case DmassT_INPUTS:
input_pair = DmolarT_INPUTS;
value1 /= mm;
break;
//case HmassT_INPUTS: input_pair = HmolarT_INPUTS; value1 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
case SmassT_INPUTS:
input_pair = SmolarT_INPUTS;
value1 *= mm;
break;
//case TUmass_INPUTS: input_pair = TUmolar_INPUTS; value2 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
case DmassP_INPUTS:
input_pair = DmolarP_INPUTS;
value1 /= mm;
break;
case DmassQ_INPUTS:
input_pair = DmolarQ_INPUTS;
value1 /= mm;
break;
case HmassP_INPUTS:
input_pair = HmolarP_INPUTS;
value1 *= mm;
break;
case PSmass_INPUTS:
input_pair = PSmolar_INPUTS;
value2 *= mm;
break;
case PUmass_INPUTS:
input_pair = PUmolar_INPUTS;
value2 *= mm;
break;
case HmassSmass_INPUTS:
input_pair = HmolarSmolar_INPUTS;
value1 *= mm;
value2 *= mm;
break;
case SmassUmass_INPUTS:
input_pair = SmolarUmolar_INPUTS;
value1 *= mm;
value2 *= mm;
break;
case DmassHmass_INPUTS:
input_pair = DmolarHmolar_INPUTS;
value1 /= mm;
value2 *= mm;
break;
case DmassSmass_INPUTS:
input_pair = DmolarSmolar_INPUTS;
value1 /= mm;
value2 *= mm;
break;
case DmassUmass_INPUTS:
input_pair = DmolarUmolar_INPUTS;
value1 /= mm;
value2 *= mm;
break;
default:
break;
}
break;
}
default:
return;
}
}
double AbstractState::trivial_keyed_output(parameters key) {
if (get_debug_level() >= 50)
std::cout << format("AbstractState: trivial_keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
switch (key) {
case imolar_mass:
return molar_mass();
case iacentric_factor:
return acentric_factor();
case igas_constant:
return gas_constant();
case iT_min:
return Tmin();
case iT_triple:
return Ttriple();
case iT_max:
return Tmax();
case iP_max:
return pmax();
case iP_min:
case iP_triple:
return this->p_triple();
case iT_reducing:
return calc_T_reducing();
case irhomolar_reducing:
return calc_rhomolar_reducing();
case iP_reducing:
return calc_p_reducing();
case iP_critical:
return this->p_critical();
case iT_critical:
return this->T_critical();
case irhomolar_critical:
return this->rhomolar_critical();
case irhomass_critical:
return this->rhomass_critical();
case iODP:
return this->calc_ODP();
case iGWP100:
return this->calc_GWP100();
case iGWP20:
return this->calc_GWP20();
case iGWP500:
return this->calc_GWP500();
case ifraction_min:
return this->calc_fraction_min();
case ifraction_max:
return this->calc_fraction_max();
case iT_freeze:
return this->calc_T_freeze();
case iFH:
return this->calc_flame_hazard();
case iHH:
return this->calc_health_hazard();
case iPH:
return this->calc_physical_hazard();
case idipole_moment:
return this->calc_dipole_moment();
default:
throw ValueError(
format("This input [%d: \"%s\"] is not valid for trivial_keyed_output", key, get_parameter_information(key, "short").c_str()));
}
}
double AbstractState::keyed_output(parameters key) {
if (get_debug_level() >= 50)
std::cout << format("AbstractState: keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
// Handle trivial inputs
if (is_trivial_parameter(key)) {
return trivial_keyed_output(key);
}
switch (key) {
case iQ:
return Q();
case iT:
return T();
case iP:
return p();
case iDmolar:
return rhomolar();
case iDmass:
return rhomass();
case iHmolar:
return hmolar();
case iHmolar_residual:
return hmolar_residual();
case iHmass:
return hmass();
case iSmolar:
return smolar();
case iSmolar_residual:
return smolar_residual();
case iSmass:
return smass();
case iUmolar:
return umolar();
case iUmass:
return umass();
case iGmolar:
return gibbsmolar();
case iGmolar_residual:
return gibbsmolar_residual();
case iGmass:
return gibbsmass();
case iHelmholtzmolar:
return helmholtzmolar();
case iHelmholtzmass:
return helmholtzmass();
case iCvmolar:
return cvmolar();
case iCvmass:
return cvmass();
case iCpmolar:
return cpmolar();
case iCp0molar:
return cp0molar();
case iCpmass:
return cpmass();
case iCp0mass:
return cp0mass();
case imolar_mass:
return molar_mass();
case iT_reducing:
return get_reducing_state().T;
case irhomolar_reducing:
return get_reducing_state().rhomolar;
case ispeed_sound:
return speed_sound();
case ialphar:
return alphar();
case ialpha0:
return alpha0();
case idalpha0_ddelta_consttau:
return dalpha0_dDelta();
case id2alpha0_ddelta2_consttau:
return d2alpha0_dDelta2();
case id3alpha0_ddelta3_consttau:
return d3alpha0_dDelta3();
case idalpha0_dtau_constdelta:
return dalpha0_dTau();
case idalphar_ddelta_consttau:
return dalphar_dDelta();
case idalphar_dtau_constdelta:
return dalphar_dTau();
case iBvirial:
return Bvirial();
case idBvirial_dT:
return dBvirial_dT();
case iCvirial:
return Cvirial();
case idCvirial_dT:
return dCvirial_dT();
case iisothermal_compressibility:
return isothermal_compressibility();
case iisobaric_expansion_coefficient:
return isobaric_expansion_coefficient();
case iisentropic_expansion_coefficient:
return isentropic_expansion_coefficient();
case iviscosity:
return viscosity();
case iconductivity:
return conductivity();
case iPrandtl:
return Prandtl();
case isurface_tension:
return surface_tension();
case iPhase:
return phase();
case iZ:
return compressibility_factor();
case iPIP:
return PIP();
case ifundamental_derivative_of_gas_dynamics:
return fundamental_derivative_of_gas_dynamics();
default:
throw ValueError(format("This input [%d: \"%s\"] is not valid for keyed_output", key, get_parameter_information(key, "short").c_str()));
}
}
double AbstractState::tau(void) {
if (!_tau) _tau = calc_reciprocal_reduced_temperature();
return _tau;
}
double AbstractState::delta(void) {
if (!_delta) _delta = calc_reduced_density();
return _delta;
}
double AbstractState::Tmin(void) {
return calc_Tmin();
}
double AbstractState::Tmax(void) {
return calc_Tmax();
}
double AbstractState::Ttriple(void) {
return calc_Ttriple();
}
double AbstractState::pmax(void) {
return calc_pmax();
}
double AbstractState::T_critical(void) {
return calc_T_critical();
}
double AbstractState::T_reducing(void) {
if (!ValidNumber(_reducing.T)) {
calc_reducing_state();
}
return _reducing.T;
}
double AbstractState::p_critical(void) {
return calc_p_critical();
}
double AbstractState::p_triple(void) {
return calc_p_triple();
}
double AbstractState::rhomolar_critical(void) {
return calc_rhomolar_critical();
}
double AbstractState::rhomass_critical(void) {
return calc_rhomolar_critical() * molar_mass();
}
double AbstractState::rhomolar_reducing(void) {
if (!ValidNumber(_reducing.rhomolar)) {
calc_reducing_state();
}
return _reducing.rhomolar;
}
double AbstractState::rhomass_reducing(void) {
return rhomolar_reducing() * molar_mass();
}
double AbstractState::hmolar(void) {
if (!_hmolar) _hmolar = calc_hmolar();
return _hmolar;
}
double AbstractState::hmolar_residual(void) {
if (!_hmolar_residual) _hmolar_residual = calc_hmolar_residual();
return _hmolar_residual;
}
double AbstractState::hmolar_excess(void) {
if (!_hmolar_excess) calc_excess_properties();
return _hmolar_excess;
}
double AbstractState::smolar(void) {
if (!_smolar) _smolar = calc_smolar();
return _smolar;
}
double AbstractState::smolar_residual(void) {
if (!_smolar_residual) _smolar_residual = calc_smolar_residual();
return _smolar_residual;
}
double AbstractState::neff(void) {
double tau = calc_T_reducing()/_T;
double delta = _rhomolar/calc_rhomolar_reducing();
double Ar01 = delta*dalphar_dDelta();
double Ar11 = tau*delta*d2alphar_dDelta_dTau();
double Ar20 = tau*tau*d2alphar_dTau2();
return -3.0*(Ar01-Ar11)/Ar20;
}
double AbstractState::smolar_excess(void) {
if (!_smolar_excess) calc_excess_properties();
return _smolar_excess;
}
double AbstractState::umolar(void) {
if (!_umolar) _umolar = calc_umolar();
return _umolar;
}
double AbstractState::umolar_excess(void) {
if (!_umolar_excess) calc_excess_properties();
return _umolar_excess;
}
double AbstractState::gibbsmolar(void) {
if (!_gibbsmolar) _gibbsmolar = calc_gibbsmolar();
return _gibbsmolar;
}
double AbstractState::gibbsmolar_residual(void) {
if (!_gibbsmolar_residual) _gibbsmolar_residual = calc_gibbsmolar_residual();
return _gibbsmolar_residual;
}
double AbstractState::gibbsmolar_excess(void) {
if (!_gibbsmolar_excess) calc_excess_properties();
return _gibbsmolar_excess;
}
double AbstractState::helmholtzmolar(void) {
if (!_helmholtzmolar) _helmholtzmolar = calc_helmholtzmolar();
return _helmholtzmolar;
}
double AbstractState::helmholtzmolar_excess(void) {
if (!_helmholtzmolar_excess) calc_excess_properties();
return _helmholtzmolar_excess;
}
double AbstractState::volumemolar_excess(void) {
if (!_volumemolar_excess) calc_excess_properties();
return _volumemolar_excess;
}
double AbstractState::cpmolar(void) {
if (!_cpmolar) _cpmolar = calc_cpmolar();
return _cpmolar;
}
double AbstractState::cp0molar(void) {
return calc_cpmolar_idealgas();
}
double AbstractState::cvmolar(void) {
if (!_cvmolar) _cvmolar = calc_cvmolar();
return _cvmolar;
}
double AbstractState::speed_sound(void) {
if (!_speed_sound) _speed_sound = calc_speed_sound();
return _speed_sound;
}
double AbstractState::viscosity(void) {
if (!_viscosity) _viscosity = calc_viscosity();
return _viscosity;
}
double AbstractState::conductivity(void) {
if (!_conductivity) _conductivity = calc_conductivity();
return _conductivity;
}
double AbstractState::melting_line(int param, int given, double value) {
return calc_melting_line(param, given, value);
}
double AbstractState::acentric_factor() {
return calc_acentric_factor();
}
double AbstractState::saturation_ancillary(parameters param, int Q, parameters given, double value) {
return calc_saturation_ancillary(param, Q, given, value);
}
double AbstractState::surface_tension(void) {
if (!_surface_tension) _surface_tension = calc_surface_tension();
return _surface_tension;
}
double AbstractState::molar_mass(void) {
if (!_molar_mass) _molar_mass = calc_molar_mass();
return _molar_mass;
}
double AbstractState::gas_constant(void) {
if (!_gas_constant) _gas_constant = calc_gas_constant();
return _gas_constant;
}
double AbstractState::fugacity_coefficient(std::size_t i) {
// TODO: Cache the fug. coeff for each component
return calc_fugacity_coefficient(i);
}
std::vector<double> AbstractState::fugacity_coefficients() {
// TODO: Cache the fug. coeff for each component
return calc_fugacity_coefficients();
}
double AbstractState::fugacity(std::size_t i) {
// TODO: Cache the fug. coeff for each component
return calc_fugacity(i);
}
double AbstractState::chemical_potential(std::size_t i) {
// TODO: Cache the chemical potential for each component
return calc_chemical_potential(i);
}
void AbstractState::build_phase_envelope(const std::string& type) {
calc_phase_envelope(type);
}
double AbstractState::isothermal_compressibility(void) {
return 1.0 / _rhomolar * first_partial_deriv(iDmolar, iP, iT);
}
double AbstractState::isobaric_expansion_coefficient(void) {
return -1.0 / _rhomolar * first_partial_deriv(iDmolar, iT, iP);
}
double AbstractState::isentropic_expansion_coefficient(void) {
return _rhomolar / _p * first_partial_deriv(iP, iDmolar, iSmolar);
}
double AbstractState::Bvirial(void) {
return calc_Bvirial();
}
double AbstractState::Cvirial(void) {
return calc_Cvirial();
}
double AbstractState::dBvirial_dT(void) {
return calc_dBvirial_dT();
}
double AbstractState::dCvirial_dT(void) {
return calc_dCvirial_dT();
}
double AbstractState::compressibility_factor(void) {
return calc_compressibility_factor();
}
double AbstractState::fundamental_derivative_of_gas_dynamics() {
// See Colonna, FPE, 2010, Eq. 1
return 1 + this->second_partial_deriv(iP, iDmass, iSmolar, iDmass, iSmolar) * this->rhomass() / (2 * powInt(speed_sound(), 2));
};
// Get the derivatives of the parameters in the partial derivative with respect to T and rho
void get_dT_drho(AbstractState& AS, parameters index, CoolPropDbl& dT, CoolPropDbl& drho) {
CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), dT_dtau = -pow(T, 2) / Tr,
R = AS.gas_constant(), delta = rho / rhor, tau = Tr / T;
switch (index) {
case iT:
dT = 1;
drho = 0;
break;
case iDmolar:
dT = 0;
drho = 1;
break;
case iDmass:
dT = 0;
drho = AS.molar_mass();
break;
case iP: {
// dp/drho|T
drho = R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
// dp/dT|rho
dT = rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau());
break;
}
case iHmass:
case iHmolar: {
// dh/dT|rho
dT = R
* (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
+ (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
// dh/drhomolar|T
drho = T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
if (index == iHmass) {
// dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass
drho /= AS.molar_mass();
dT /= AS.molar_mass();
}
break;
}
case iSmass:
case iSmolar: {
// ds/dT|rho
dT = R / T * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
// ds/drho|T
drho = R / rho * (-(1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
if (index == iSmass) {
// ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
drho /= AS.molar_mass();
dT /= AS.molar_mass();
}
break;
}
case iUmass:
case iUmolar: {
// du/dT|rho
dT = R * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
// du/drho|T
drho = AS.T() * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau());
if (index == iUmass) {
// du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
drho /= AS.molar_mass();
dT /= AS.molar_mass();
}
break;
}
case iGmass:
case iGmolar: {
// dg/dT|rho
double dTau_dT = 1 / dT_dtau;
dT = R * AS.T() * (AS.dalpha0_dTau() + AS.dalphar_dTau() + AS.delta() * AS.d2alphar_dDelta_dTau()) * dTau_dT
+ R * (1 + AS.alpha0() + AS.alphar() + AS.delta() * AS.dalphar_dDelta());
// dg/drho|T
double dDelta_drho = 1 / rhor;
drho = AS.T() * R * (AS.dalpha0_dDelta() + AS.dalphar_dDelta() + AS.delta() * AS.d2alphar_dDelta2() + AS.dalphar_dDelta()) * dDelta_drho;
if (index == iGmass) {
// dg/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
drho /= AS.molar_mass();
dT /= AS.molar_mass();
}
break;
}
case iTau:
dT = 1 / dT_dtau;
drho = 0;
break;
case iDelta:
dT = 0;
drho = 1 / rhor;
break;
case iCvmolar:
case iCvmass: {
// use the second order derivative of internal energy
// make it cleaner by using the function get_dT_drho_second_derivatives directly?
// dcvdT|rho = d2u/dT2|rho
dT = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
// dcvdrho|T = d2u/dT/drho
drho = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
if (index == iCvmass) {
drho /= AS.molar_mass();
dT /= AS.molar_mass();
}
break;
}
case iCpmolar:
case iCpmass: {
// dcp/dT|rho = d2h/dT2 + dh/drho * dP/dT * d2P/drhodT / ( dp/drho )^2 - ( d2h/dTdrho * dP/dT + dh/drho * d2P/dT2 ) / ( dP/drho )
dT = R / T * pow(tau, 2)
* (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
+ delta * AS.d3alphar_dDelta_dTau2());
dT += (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
* (R
* (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
- tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau()))
/ pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
dT -= ((R / rho * delta
* (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
- tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()))
* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
+ (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
* (rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2())))
/ (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
// dcpdrho|T = d2h/dTdrho + dh/drho * dP/dT * d2P/drho2 / ( dp/drho )^2 - ( d2h/drho2 * dP/dT + dh/drho * d2P/dTdrho ) / ( dP/drho )
drho = R / rho * delta
* (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
- tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()); //d2h/dTdrho
drho +=
(T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
* (T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3()))
/ pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
drho -= ((R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3()))
* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
+ (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
* (R
* (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()
- 2 * delta * tau * AS.d2alphar_dDelta_dTau() - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau())))
/ (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
if (index == iCpmass) {
drho /= AS.molar_mass();
dT /= AS.molar_mass();
}
break;
}
case ispeed_sound: {
//dwdT
double aa = 1.0 + delta * AS.dalphar_dDelta() - delta * tau * AS.d2alphar_dDelta_dTau();
double bb = pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
double daa_dTau = -delta * tau * AS.d3alphar_dDelta_dTau2();
double dbb_dTau = pow(tau, 2) * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2.0 * tau * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
double w = AS.speed_sound();
dT = 1.0 / 2.0 / w / T
* (pow(w, 2)
- R * Tr / AS.molar_mass()
* (2.0 * delta * AS.d2alphar_dDelta_dTau() + pow(delta, 2) * AS.d3alphar_dDelta2_dTau()
- (2 * aa / bb * daa_dTau - pow(aa / bb, 2) * dbb_dTau)));
//dwdrho
double daa_dDelta =
AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2() - tau * (AS.d2alphar_dDelta_dTau() + delta * AS.d3alphar_dDelta2_dTau());
double dbb_dDelta = pow(tau, 2) * (AS.d3alpha0_dDelta_dTau2() + AS.d3alphar_dDelta_dTau2());
drho = R * T / 2.0 / AS.molar_mass() / w / rhor
* (2.0 * (AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2())
+ (2.0 * delta * AS.d2alphar_dDelta2() + pow(delta, 2) * AS.d3alphar_dDelta3())
- (2 * aa / bb * daa_dDelta - pow(aa / bb, 2) * dbb_dDelta));
break;
}
default:
throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));
}
}
void get_dT_drho_second_derivatives(AbstractState& AS, int index, CoolPropDbl& dT2, CoolPropDbl& drho_dT, CoolPropDbl& drho2) {
CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), R = AS.gas_constant(), delta = rho / rhor,
tau = Tr / T;
// Here we use T and rho as independent variables since derivations are already done by Thorade, 2013,
// Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z
switch (index) {
case iT:
case iDmass:
case iDmolar:
dT2 = 0; // d2rhomolar_dtau2
drho2 = 0;
drho_dT = 0;
break;
case iTau:
dT2 = 2 * Tr / pow(T, 3);
drho_dT = 0;
drho2 = 0;
break;
case iDelta:
dT2 = 0;
drho_dT = 0;
drho2 = 0;
break;
case iP: {
drho2 =
T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3());
dT2 = rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
drho_dT = R
* (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
- tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
break;
}
case iHmass:
case iHmolar: {
// d2h/drho2|T
drho2 = R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3());
// d2h/dT2|rho
dT2 = R / T * pow(tau, 2)
* (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
+ delta * AS.d3alphar_dDelta_dTau2());
// d2h/drho/dT
drho_dT = R / rho * delta
* (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
- tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau());
if (index == iHmass) {
drho2 /= AS.molar_mass();
drho_dT /= AS.molar_mass();
dT2 /= AS.molar_mass();
}
break;
}
case iSmass:
case iSmolar: {
// d2s/rho2|T
drho2 = R / pow(rho, 2) * (1 - pow(delta, 2) * AS.d2alphar_dDelta2() + tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
// d2s/dT2|rho
dT2 = R * pow(tau / T, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 3 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
// d2s/drho/dT
drho_dT = R / (T * rho) * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
if (index == iSmass) {
drho2 /= AS.molar_mass();
drho_dT /= AS.molar_mass();
dT2 /= AS.molar_mass();
}
break;
}
case iUmass:
case iUmolar: {
// d2u/rho2|T
drho2 = R * T * tau * pow(delta / rho, 2) * AS.d3alphar_dDelta2_dTau();
// d2u/dT2|rho
dT2 = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
// d2u/drho/dT
drho_dT = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
if (index == iUmass) {
drho2 /= AS.molar_mass();
drho_dT /= AS.molar_mass();
dT2 /= AS.molar_mass();
}
break;
}
default:
throw ValueError(format("input to get_dT_drho_second_derivatives[%s] is invalid", get_parameter_information(index, "short").c_str()));
}
}
CoolPropDbl AbstractState::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) {
CoolPropDbl dOf_dT, dOf_drho, dWrt_dT, dWrt_drho, dConstant_dT, dConstant_drho;
get_dT_drho(*this, Of, dOf_dT, dOf_drho);
get_dT_drho(*this, Wrt, dWrt_dT, dWrt_drho);
get_dT_drho(*this, Constant, dConstant_dT, dConstant_drho);
return (dOf_dT * dConstant_drho - dOf_drho * dConstant_dT) / (dWrt_dT * dConstant_drho - dWrt_drho * dConstant_dT);
}
CoolPropDbl AbstractState::calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2) {
CoolPropDbl dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2, d2Wrt1_dT2,
d2Wrt1_drhodT, d2Wrt1_drho2, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2, dWrt2_dT, dWrt2_drho, dConstant2_dT, dConstant2_drho, N, D,
dNdrho__T, dDdrho__T, dNdT__rho, dDdT__rho, dderiv1_drho, dderiv1_dT, second;
// First and second partials needed for terms involved in first derivative
get_dT_drho(*this, Of1, dOf1_dT, dOf1_drho);
get_dT_drho(*this, Wrt1, dWrt1_dT, dWrt1_drho);
get_dT_drho(*this, Constant1, dConstant1_dT, dConstant1_drho);
get_dT_drho_second_derivatives(*this, Of1, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2);
get_dT_drho_second_derivatives(*this, Wrt1, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2);
get_dT_drho_second_derivatives(*this, Constant1, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2);
// First derivatives of terms involved in the second derivative
get_dT_drho(*this, Wrt2, dWrt2_dT, dWrt2_drho);
get_dT_drho(*this, Constant2, dConstant2_dT, dConstant2_drho);
// Numerator and denominator of first partial derivative term
N = dOf1_dT * dConstant1_drho - dOf1_drho * dConstant1_dT;
D = dWrt1_dT * dConstant1_drho - dWrt1_drho * dConstant1_dT;
// Derivatives of the numerator and denominator of the first partial derivative term with respect to rho, T held constant
// They are of similar form, with Of1 and Wrt1 swapped
dNdrho__T = dOf1_dT * d2Constant1_drho2 + d2Of1_drhodT * dConstant1_drho - dOf1_drho * d2Constant1_drhodT - d2Of1_drho2 * dConstant1_dT;
dDdrho__T = dWrt1_dT * d2Constant1_drho2 + d2Wrt1_drhodT * dConstant1_drho - dWrt1_drho * d2Constant1_drhodT - d2Wrt1_drho2 * dConstant1_dT;
// Derivatives of the numerator and denominator of the first partial derivative term with respect to T, rho held constant
// They are of similar form, with Of1 and Wrt1 swapped
dNdT__rho = dOf1_dT * d2Constant1_drhodT + d2Of1_dT2 * dConstant1_drho - dOf1_drho * d2Constant1_dT2 - d2Of1_drhodT * dConstant1_dT;
dDdT__rho = dWrt1_dT * d2Constant1_drhodT + d2Wrt1_dT2 * dConstant1_drho - dWrt1_drho * d2Constant1_dT2 - d2Wrt1_drhodT * dConstant1_dT;
// First partial of first derivative term with respect to T
dderiv1_drho = (D * dNdrho__T - N * dDdrho__T) / pow(D, 2);
// First partial of first derivative term with respect to rho
dderiv1_dT = (D * dNdT__rho - N * dDdT__rho) / pow(D, 2);
// Complete second derivative
second = (dderiv1_dT * dConstant2_drho - dderiv1_drho * dConstant2_dT) / (dWrt2_dT * dConstant2_drho - dWrt2_drho * dConstant2_dT);
return second;
}
// // ----------------------------------------
// // Smoothing functions for density
// // ----------------------------------------
// /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
// virtual double AbstractState::drhodh_constp_smoothed(double xend);
// /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
// virtual double AbstractState::drhodp_consth_smoothed(double xend);
// /// Density corresponding to the smoothed derivatives in the region of x=0 to x=xend
// virtual void AbstractState::rho_smoothed(double xend, double *rho_spline, double *dsplinedh, double *dsplinedp);
} /* namespace CoolProp */
#ifdef ENABLE_CATCH
#include <catch2/catch_all.hpp>
TEST_CASE("Check AbstractState", "[AbstractState]") {
SECTION("bad backend") {
CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water")));
}
SECTION("good backend - bad fluid") {
CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "DEFINITELY_A_BAD_FLUID")));
}
SECTION("good backend - helmholtz") {
CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Water")));
}
SECTION("good backend - incomp") {
CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("INCOMP", "DEB")));
}
SECTION("good backend - REFPROP") {
CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("REFPROP", "Water")));
}
}
TEST_CASE("Check derivatives in first_partial_deriv", "[derivs_in_first_partial_deriv]") {
shared_ptr<CoolProp::AbstractState> Water(CoolProp::AbstractState::factory("HEOS", "Water"));
shared_ptr<CoolProp::AbstractState> WaterplusT(CoolProp::AbstractState::factory("HEOS", "Water"));
shared_ptr<CoolProp::AbstractState> WaterminusT(CoolProp::AbstractState::factory("HEOS", "Water"));
shared_ptr<CoolProp::AbstractState> Waterplusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
shared_ptr<CoolProp::AbstractState> Waterminusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
double dT = 1e-3, drho = 1;
Water->update(CoolProp::PT_INPUTS, 101325, 300);
WaterplusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 + dT);
WaterminusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 - dT);
Waterplusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() + drho, 300);
Waterminusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() - drho, 300);
// Numerical derivatives
CoolPropDbl dP_dT_num = (WaterplusT->p() - WaterminusT->p()) / (2 * dT);
CoolPropDbl dP_drho_num = (Waterplusrho->p() - Waterminusrho->p()) / (2 * drho);
CoolPropDbl dHmolar_dT_num = (WaterplusT->hmolar() - WaterminusT->hmolar()) / (2 * dT);
CoolPropDbl dHmolar_drho_num = (Waterplusrho->hmolar() - Waterminusrho->hmolar()) / (2 * drho);
CoolPropDbl dHmass_dT_num = (WaterplusT->hmass() - WaterminusT->hmass()) / (2 * dT);
CoolPropDbl dHmass_drho_num = (Waterplusrho->hmass() - Waterminusrho->hmass()) / (2 * drho);
CoolPropDbl dSmolar_dT_num = (WaterplusT->smolar() - WaterminusT->smolar()) / (2 * dT);
CoolPropDbl dSmolar_drho_num = (Waterplusrho->smolar() - Waterminusrho->smolar()) / (2 * drho);
CoolPropDbl dSmass_dT_num = (WaterplusT->smass() - WaterminusT->smass()) / (2 * dT);
CoolPropDbl dSmass_drho_num = (Waterplusrho->smass() - Waterminusrho->smass()) / (2 * drho);
CoolPropDbl dUmolar_dT_num = (WaterplusT->umolar() - WaterminusT->umolar()) / (2 * dT);
CoolPropDbl dUmolar_drho_num = (Waterplusrho->umolar() - Waterminusrho->umolar()) / (2 * drho);
CoolPropDbl dUmass_dT_num = (WaterplusT->umass() - WaterminusT->umass()) / (2 * dT);
CoolPropDbl dUmass_drho_num = (Waterplusrho->umass() - Waterminusrho->umass()) / (2 * drho);
CoolPropDbl dGmolar_dT_num = (WaterplusT->gibbsmolar() - WaterminusT->gibbsmolar()) / (2 * dT);
CoolPropDbl dGmolar_drho_num = (Waterplusrho->gibbsmolar() - Waterminusrho->gibbsmolar()) / (2 * drho);
CoolPropDbl dGmass_dT_num = (WaterplusT->gibbsmass() - WaterminusT->gibbsmass()) / (2 * dT);
CoolPropDbl dGmass_drho_num = (Waterplusrho->gibbsmass() - Waterminusrho->gibbsmass()) / (2 * drho);
CoolPropDbl dCvmolar_dT_num = (WaterplusT->cvmolar() - WaterminusT->cvmolar()) / (2 * dT);
CoolPropDbl dCvmolar_drho_num = (Waterplusrho->cvmolar() - Waterminusrho->cvmolar()) / (2 * drho);
CoolPropDbl dCvmass_dT_num = (WaterplusT->cvmass() - WaterminusT->cvmass()) / (2 * dT);
CoolPropDbl dCvmass_drho_num = (Waterplusrho->cvmass() - Waterminusrho->cvmass()) / (2 * drho);
CoolPropDbl dCpmolar_dT_num = (WaterplusT->cpmolar() - WaterminusT->cpmolar()) / (2 * dT);
CoolPropDbl dCpmolar_drho_num = (Waterplusrho->cpmolar() - Waterminusrho->cpmolar()) / (2 * drho);
CoolPropDbl dCpmass_dT_num = (WaterplusT->cpmass() - WaterminusT->cpmass()) / (2 * dT);
CoolPropDbl dCpmass_drho_num = (Waterplusrho->cpmass() - Waterminusrho->cpmass()) / (2 * drho);
CoolPropDbl dspeed_sound_dT_num = (WaterplusT->speed_sound() - WaterminusT->speed_sound()) / (2 * dT);
CoolPropDbl dspeed_sound_drho_num = (Waterplusrho->speed_sound() - Waterminusrho->speed_sound()) / (2 * drho);
// Pressure
CoolPropDbl dP_dT_analyt, dP_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iP, dP_dT_analyt, dP_drho_analyt);
// Enthalpy
CoolPropDbl dHmolar_dT_analyt, dHmolar_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iHmolar, dHmolar_dT_analyt, dHmolar_drho_analyt);
CoolPropDbl dHmass_dT_analyt, dHmass_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iHmass, dHmass_dT_analyt, dHmass_drho_analyt);
// Entropy
CoolPropDbl dSmolar_dT_analyt, dSmolar_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iSmolar, dSmolar_dT_analyt, dSmolar_drho_analyt);
CoolPropDbl dSmass_dT_analyt, dSmass_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iSmass, dSmass_dT_analyt, dSmass_drho_analyt);
// Internal energy
CoolPropDbl dUmolar_dT_analyt, dUmolar_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iUmolar, dUmolar_dT_analyt, dUmolar_drho_analyt);
CoolPropDbl dUmass_dT_analyt, dUmass_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iUmass, dUmass_dT_analyt, dUmass_drho_analyt);
// Gibbs
CoolPropDbl dGmolar_dT_analyt, dGmolar_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iGmolar, dGmolar_dT_analyt, dGmolar_drho_analyt);
CoolPropDbl dGmass_dT_analyt, dGmass_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iGmass, dGmass_dT_analyt, dGmass_drho_analyt);
// Isochoric heat capacity
CoolPropDbl dCvmolar_dT_analyt, dCvmolar_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iCvmolar, dCvmolar_dT_analyt, dCvmolar_drho_analyt);
CoolPropDbl dCvmass_dT_analyt, dCvmass_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iCvmass, dCvmass_dT_analyt, dCvmass_drho_analyt);
// Isobaric heat capacity
CoolPropDbl dCpmolar_dT_analyt, dCpmolar_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iCpmolar, dCpmolar_dT_analyt, dCpmolar_drho_analyt);
CoolPropDbl dCpmass_dT_analyt, dCpmass_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::iCpmass, dCpmass_dT_analyt, dCpmass_drho_analyt);
// Speed of sound
CoolPropDbl dspeed_sound_dT_analyt, dspeed_sound_drho_analyt;
CoolProp::get_dT_drho(*Water, CoolProp::ispeed_sound, dspeed_sound_dT_analyt, dspeed_sound_drho_analyt);
double eps = 1e-3;
CHECK(std::abs(dP_dT_analyt / dP_dT_num - 1) < eps);
CHECK(std::abs(dP_drho_analyt / dP_drho_num - 1) < eps);
CHECK(std::abs(dHmolar_dT_analyt / dHmolar_dT_num - 1) < eps);
CHECK(std::abs(dHmolar_drho_analyt / dHmolar_drho_num - 1) < eps);
CHECK(std::abs(dHmass_dT_analyt / dHmass_dT_num - 1) < eps);
CHECK(std::abs(dHmass_drho_analyt / dHmass_drho_num - 1) < eps);
CHECK(std::abs(dSmolar_dT_analyt / dSmolar_dT_num - 1) < eps);
CHECK(std::abs(dSmolar_drho_analyt / dSmolar_drho_num - 1) < eps);
CHECK(std::abs(dSmass_dT_analyt / dSmass_dT_num - 1) < eps);
CHECK(std::abs(dSmass_drho_analyt / dSmass_drho_num - 1) < eps);
CHECK(std::abs(dUmolar_dT_analyt / dUmolar_dT_num - 1) < eps);
CHECK(std::abs(dUmolar_drho_analyt / dUmolar_drho_num - 1) < eps);
CHECK(std::abs(dUmass_dT_analyt / dUmass_dT_num - 1) < eps);
CHECK(std::abs(dUmass_drho_analyt / dUmass_drho_num - 1) < eps);
CHECK(std::abs(dGmolar_dT_analyt / dGmolar_dT_num - 1) < eps);
CHECK(std::abs(dGmolar_drho_analyt / dGmolar_drho_num - 1) < eps);
CHECK(std::abs(dGmass_dT_analyt / dGmass_dT_num - 1) < eps);
CHECK(std::abs(dGmass_drho_analyt / dGmass_drho_num - 1) < eps);
CHECK(std::abs(dCvmolar_dT_analyt / dCvmolar_dT_num - 1) < eps);
CHECK(std::abs(dCvmolar_drho_analyt / dCvmolar_drho_num - 1) < eps);
CHECK(std::abs(dCvmass_dT_analyt / dCvmass_dT_num - 1) < eps);
CHECK(std::abs(dCvmass_drho_analyt / dCvmass_drho_num - 1) < eps);
CHECK(std::abs(dCpmolar_dT_analyt / dCpmolar_dT_num - 1) < eps);
CHECK(std::abs(dCpmolar_drho_analyt / dCpmolar_drho_num - 1) < eps);
CHECK(std::abs(dCpmass_dT_analyt / dCpmass_dT_num - 1) < eps);
CHECK(std::abs(dCpmass_drho_analyt / dCpmass_drho_num - 1) < eps);
CHECK(std::abs(dspeed_sound_dT_analyt / dspeed_sound_dT_num - 1) < eps);
CHECK(std::abs(dspeed_sound_drho_analyt / dspeed_sound_drho_num - 1) < eps);
}
#endif
Reference in New Issue View Git Blame Copy Permalink