Merge branch 'master' into two_phase_derivatives_and_splines

This commit is contained in:
Ian Bell
2015-02-02 12:23:12 -05:00
64 changed files with 2249 additions and 540 deletions

View File

@@ -11,6 +11,7 @@ namespace CoolProp{
void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.get_mole_fractions_ref().size() < 2){throw ValueError("Cannot build phase envelope for pure fluid");}
std::size_t failure_count = 0;
// Set some imput options
SaturationSolvers::mixture_VLE_IO io;

View File

@@ -1177,8 +1177,8 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS
}
while(std::abs(f) > 1e-12 && iter < options.Nstep_max);
HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq);
HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);
HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
options.p = HEOS.SatL->p();
options.T = HEOS.SatL->T();

View File

@@ -26,11 +26,17 @@ IncompressibleBackend::IncompressibleBackend() {
IncompressibleBackend::IncompressibleBackend(IncompressibleFluid* fluid) {
//this->_fractions_id = fluid->getxid();
this->fluid = fluid;
if (this->fluid->is_pure()){
this->set_fractions(std::vector<long double>(1,1.0));
}
}
IncompressibleBackend::IncompressibleBackend(const std::string &fluid_name) {
this->fluid = &get_incompressible_fluid(fluid_name);
//this->_fractions_id = this->fluid->getxid();
if (this->fluid->is_pure()){
this->set_fractions(std::vector<long double>(1,1.0));
}
}
IncompressibleBackend::IncompressibleBackend(const std::vector<std::string> &component_names) {
@@ -49,17 +55,22 @@ void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double valu
clear();
if (get_debug_level()>=50) std::cout << format("Incompressible backend: _fractions are %s ",vec_to_string(_fractions).c_str()) << std::endl;
if (get_debug_level()>=50) {
std::cout << format("Incompressible backend: _fractions are %s ",vec_to_string(_fractions).c_str()) << std::endl;
}
if (_fractions.size()!=1){
throw ValueError(format("%s is an incompressible fluid, mass fractions must be set to a vector with ONE entry, not %d.",this->name().c_str(),_fractions.size()));
}
if (fluid->is_pure()){
this->_fluid_type = FLUID_TYPE_INCOMPRESSIBLE_LIQUID;
if (get_debug_level()>=50) std::cout << format("Incompressible backend: Fluid type is %d ",this->_fluid_type) << std::endl;
if ((_fractions.size()!=1) || (_fractions[0]!=0.0)){
throw ValueError(format("%s is a pure fluid. The composition has to be set to a vector with one entry equal to 0.0 or nothing. %s is not valid.",this->name().c_str(),vec_to_string(_fractions).c_str()));
if (_fractions[0]!=1.0){
throw ValueError(format("%s is a pure fluid. The composition has to be set to a vector with one entry equal to 1.0. %s is not valid.",this->name().c_str(),vec_to_string(_fractions).c_str()));
}
} else {
this->_fluid_type = FLUID_TYPE_INCOMPRESSIBLE_SOLUTION;
if (get_debug_level()>=50) std::cout << format("Incompressible backend: Fluid type is %d ",this->_fluid_type) << std::endl;
if (_fractions.size()!=1 || ((_fractions[0]<0.0) || (_fractions[0]>1.0)) ){
if ( (_fractions[0]<0.0) || (_fractions[0]>1.0) ){
throw ValueError(format("%s is a solution or brine. Mass fractions must be set to a vector with one entry between 0 and 1. %s is not valid.",this->name().c_str(),vec_to_string(_fractions).c_str()));
}
}
@@ -136,8 +147,8 @@ void IncompressibleBackend::set_fractions(const std::vector<long double> &fracti
void IncompressibleBackend::set_mole_fractions(const std::vector<long double> &mole_fractions){
if (get_debug_level()>=10) std::cout << format("Incompressible backend: Called set_mole_fractions with %s ",vec_to_string(mole_fractions).c_str()) << std::endl;
if (mole_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the mole fraction vector and not %d.",mole_fractions.size()));
if (fluid->getxid()==IFRAC_PURE) {
this->set_fractions(std::vector<long double>(1,0));
if ((fluid->getxid()==IFRAC_PURE) && true ){//( this->_fractions[0]!=1.0 )){
this->set_fractions(std::vector<long double>(1,1.0));
if (get_debug_level()>=20) std::cout << format("Incompressible backend: Overwriting fractions for pure fluid with %s -> %s",vec_to_string(mole_fractions).c_str(),vec_to_string(this->_fractions).c_str()) << std::endl;
} else if (fluid->getxid()==IFRAC_MOLE) {
this->set_fractions(mole_fractions);
@@ -157,8 +168,8 @@ void IncompressibleBackend::set_mole_fractions(const std::vector<long double> &m
void IncompressibleBackend::set_mass_fractions(const std::vector<long double> &mass_fractions){
if (get_debug_level()>=10) std::cout << format("Incompressible backend: Called set_mass_fractions with %s ",vec_to_string(mass_fractions).c_str()) << std::endl;
if (mass_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the mass fraction vector and not %d.",mass_fractions.size()));
if (fluid->getxid()==IFRAC_PURE) {
this->set_fractions(std::vector<long double>(1,0));
if ((fluid->getxid()==IFRAC_PURE) && true ){// ( this->_fractions[0]!=1.0 )) {
this->set_fractions(std::vector<long double>(1,1.0));
if (get_debug_level()>=20) std::cout << format("Incompressible backend: Overwriting fractions for pure fluid with %s -> %s",vec_to_string(mass_fractions).c_str(),vec_to_string(this->_fractions).c_str()) << std::endl;
} else if (fluid->getxid()==IFRAC_MASS) {
this->set_fractions(mass_fractions);
@@ -178,8 +189,8 @@ void IncompressibleBackend::set_mass_fractions(const std::vector<long double> &m
void IncompressibleBackend::set_volu_fractions(const std::vector<long double> &volu_fractions){
if (get_debug_level()>=10) std::cout << format("Incompressible backend: Called set_volu_fractions with %s ",vec_to_string(volu_fractions).c_str()) << std::endl;
if (volu_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the volume fraction vector and not %d.",volu_fractions.size()));
if (fluid->getxid()==IFRAC_PURE) {
this->set_fractions(std::vector<long double>(1,0));
if ((fluid->getxid()==IFRAC_PURE) && true ){// ( this->_fractions[0]!=1.0 )) {
this->set_fractions(std::vector<long double>(1,1.0));
if (get_debug_level()>=20) std::cout << format("Incompressible backend: Overwriting fractions for pure fluid with %s -> %s",vec_to_string(volu_fractions).c_str(),vec_to_string(this->_fractions).c_str()) << std::endl;
} else if (fluid->getxid()==IFRAC_VOLUME) {
this->set_fractions(volu_fractions);
@@ -425,6 +436,18 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Make sure the ruslt does not change -> reference state...
val = backend.calc_smass();
backend.update( CoolProp::PT_INPUTS, p, T );
res = backend.calc_smass();
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = fluid.u(T, p, x);
res = backend.calc_umass();
@@ -524,6 +547,22 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
CAPTURE(errmsg);
CHECK( check_abs(expected,actual,acc) );
}
// entropy reference state problems
//CoolProp::set_debug_level(51);
expected = CoolProp::PropsSI("S","T",T,"P",p,"INCOMP::"+fluid+format("[%f]",x));
actual = CoolProp::PropsSI("S","T",T,"P",p,"INCOMP::"+fluid+format("[%f]",x));
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(expected);
CAPTURE(actual);
std::string name = "INCOMP::"+fluid+format("[%f]",x);
CAPTURE(name);
std::string errmsg = CoolProp::get_global_param_string("errstring");
CAPTURE(errmsg);
CHECK( check_abs(expected,actual,acc) );
}
}
SECTION("SecCool example")
{

View File

@@ -19,15 +19,23 @@ and transport properties.
//IncompressibleFluid::IncompressibleFluid();
void IncompressibleFluid::set_reference_state(double T0, double p0, double x0, double h0, double s0){
this->Tref = T0;
this->rhoref = rho(T0,p0,x0);
this->pref = p0;
this->href = h0;
// Now we take care of the energy related values
this->uref = 0.0;
this->uref = u(T0,p0,x0) - h0; // (value without ref) - (desired ref)
this->sref = 0.0;
this->sref = s(T0,p0,x0) - s0; // (value without ref) - (desired ref)
this->_href = h0;
this->_pref = p0;
this->_Tref = T0;
this->_xref = x0;
this->_rhoref = rho(T0,p0,x0);
this->_uref = h0; // Reset the old reference value
this->_uref = u(T0,p0,x0); // uses _uref
this->_sref = s0; // Reset the old reference value
this->_sref = s(T0,p0,x0); // uses _sref
// Save the values for later calls at composition changes
this->href = h0;
this->pref = p0;
this->Tref = T0;
this->xref = x0;
this->rhoref = this->_rhoref;
this->uref = this->_uref;
this->sref = s0;
}
void IncompressibleFluid::validate(){
@@ -130,7 +138,7 @@ double IncompressibleFluid::s (double T, double p, double x){
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - sref;
return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - _sref;
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
@@ -147,7 +155,7 @@ double IncompressibleFluid::u (double T, double p, double x){
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - uref;
return poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - _uref;
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
@@ -422,7 +430,7 @@ double IncompressibleFluid::T_c (double Cmass, double p, double x){
}
/// Temperature as a function of entropy as a function of temperature, pressure and composition.
double IncompressibleFluid::T_s (double Smass, double p, double x){
double s_raw = Smass + sref;
double s_raw = Smass + _sref;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limitsInt(specific_heat.coeffs, x, s_raw, Tmin, Tmax, 0, -1, 0, Tbase, xbase, 0);
@@ -438,7 +446,7 @@ double IncompressibleFluid::T_s (double Smass, double p, double x){
}
/// Temperature as a function of internal energy as a function of temperature, pressure and composition.
double IncompressibleFluid::T_u (double Umass, double p, double x){
double u_raw = Umass + uref;
double u_raw = Umass + _uref;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limitsInt(specific_heat.coeffs, x, u_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase, 0);
@@ -664,10 +672,38 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
// Compare reference state
{
res = XLT.h(Tref,pref,xref);
CAPTURE(Tref);
CAPTURE(pref);
CAPTURE(xref);
CAPTURE(href);
CAPTURE(res);
CHECK( check_abs(href,res,acc) );
}
{
res = XLT.s(Tref,pref,xref);
CAPTURE(Tref);
CAPTURE(pref);
CAPTURE(xref);
CAPTURE(sref);
CAPTURE(res);
CHECK( check_abs(sref,res,acc) );
}
double res2 = XLT.s(Tref,pref,xref);
CAPTURE(Tref);
CAPTURE(pref);
CAPTURE(xref);
CAPTURE(res);
CAPTURE(res2);
CHECK( check_abs(res2,res,acc) );
res = XLT.s(Tref,pref,xref);
CAPTURE(Tref);
CAPTURE(pref);
CAPTURE(xref);
CAPTURE(res);
CAPTURE(res2);
CHECK( check_abs(res2,res,acc) );
}
Tref = 25+273.15;
pref = 0.0;
@@ -837,6 +873,15 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
CAPTURE(actual);
CHECK( check_abs(expected,actual,acc) );
}
expected = CH3OH.s(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(expected);
CAPTURE(actual);
CHECK( check_abs(expected,actual,acc) );
}
expected = 0.0;
actual = CH3OH.s(Tref,pref,xref);
@@ -848,6 +893,15 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
CAPTURE(actual);
CHECK( expected==actual );
}
expected = CH3OH.s(Tref,pref,xref);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(expected);
CAPTURE(actual);
CHECK( expected==actual );
}
// Compare u
expected = -60043.78429641827;

View File

@@ -32,8 +32,10 @@ REFPROPBackend::REFPROPBackend(const std::string & fluid_name) {
// Set the mole fraction to 1 in the base class (we can't set the mole fraction in this class,
// otherwise a NotImplementedError will be returned)
std::vector<long double> x(1, 1.0); // (one element with value of 1.0)
REFPROPMixtureBackend::set_mole_fractions(x);
if (get_mole_fractions().empty()){
std::vector<long double> x(1, 1.0); // (one element with value of 1.0)
REFPROPMixtureBackend::set_mole_fractions(x);
}
}
REFPROPBackend::~REFPROPBackend() {

View File

@@ -51,6 +51,7 @@ enum DLLNameManglingStyle{ NO_NAME_MANGLING = 0, LOWERCASE_NAME_MANGLING, LOWERC
#include "REFPROPMixtureBackend.h"
#include "Exceptions.h"
#include "Configuration.h"
#include "CoolProp.h"
#include <stdlib.h>
#include <string>
@@ -250,23 +251,6 @@ bool load_REFPROP()
return false;
}
#if defined(__ISWINDOWS__)
// Get data associated with path using the windows libraries,
// and if you can (result == 0), the path exists
#ifdef __MINGW32__
struct stat buf;
if ( stat( "c:\\Program Files\\REFPROP\\fluids", &buf) != 0){
throw CoolProp::ValueError("REFPROP fluid files must be copied to c:\\Program Files\\REFPROP\\fluids");
}
#else
struct _stat buf;
if ( _stat( "c:\\Program Files\\REFPROP\\fluids", &buf) != 0){
throw CoolProp::ValueError("REFPROP fluid files must be copied to c:\\Program Files\\REFPROP\\fluids");
}
#endif
#endif
if (setFunctionPointers()!=COOLPROP_OK)
{
printf("There was an error setting the REFPROP function pointers, check types and names in header file.\n");
@@ -395,6 +379,40 @@ void REFPROPMixtureBackend::set_REFPROP_fluids(const std::vector<std::string> &f
}
else
{
if (N == 1 && upper(components_joined_raw).find(".MIX") != std::string::npos){
// It's a predefined mixture
ierr = 0;
std::vector<double> x(ncmax);
char mix[255];
strcpy(mix, components_joined_raw.c_str());
SETMIXdll(mix,
"HMX.BNC",
"DEF",
&N,
component_string,
&(x[0]),
&ierr,
herr,
255,
255,
3,
10000,
255);
if (static_cast<int>(ierr) <= 0){
this->Ncomp = N;
mole_fractions.resize(N);
mole_fractions_liq.resize(N);
mole_fractions_vap.resize(N);
LoadedREFPROPRef = components_joined_raw;
if (CoolProp::get_debug_level() > 5){ std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n",__FILE__,__LINE__, components_joined.c_str()); }
if (dbg_refprop) std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n",__FILE__,__LINE__, components_joined.c_str());
set_mole_fractions(std::vector<long double>(x.begin(), x.begin()+N));
return;
}
else{
throw ValueError(format("Unable to load mixture: %s",components_joined_raw.c_str()));
}
}
// Loop over the file names - first we try with nothing, then .fld, then .FLD, then .ppf - means you can't mix and match
for (unsigned int k = 0; k < number_of_endings; k++)
{
@@ -425,22 +443,25 @@ void REFPROPMixtureBackend::set_REFPROP_fluids(const std::vector<std::string> &f
errormessagelength // Length of error message
);
if (ierr <= 0) // Success (or a warning, which is silently squelched for now)
if (static_cast<int>(ierr) <= 0) // Success (or a warning, which is silently squelched for now)
{
this->Ncomp = N;
mole_fractions.resize(N);
mole_fractions_liq.resize(N);
mole_fractions_vap.resize(N);
LoadedREFPROPRef = components_joined_raw;
if (CoolProp::get_debug_level() > 5){ std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n",__FILE__,__LINE__, components_joined.c_str()); }
if (dbg_refprop) std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n",__FILE__,__LINE__, components_joined.c_str());
return;
}
else if (k < number_of_endings-1){ // Keep going
if (CoolProp::get_debug_level() > 5){std::cout << format("REFPROP error/warning [ierr: %d]: %s",ierr, herr) << std::endl;}
continue;
}
else
{
throw ValueError(format("Could not load these fluids: %s; error: %s", components_joined_raw.c_str(), herr));
if (CoolProp::get_debug_level() > 5){std::cout << format("k: %d #endings: %d", k, number_of_endings) << std::endl;}
throw ValueError(format("Could not load these fluids: %s", components_joined_raw.c_str()));
}
}
}
@@ -514,21 +535,21 @@ long double REFPROPMixtureBackend::calc_T_critical(){
char herr[255];
double Tcrit, pcrit_kPa, dcrit_mol_L;
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255);
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(Tcrit);
};
long double REFPROPMixtureBackend::calc_p_critical(){
long ierr = 0;
char herr[255];
double Tcrit, pcrit_kPa, dcrit_mol_L;
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(pcrit_kPa*1000);
};
long double REFPROPMixtureBackend::calc_rhomolar_critical(){
long ierr = 0;
char herr[255];
double Tcrit, pcrit_kPa, dcrit_mol_L;
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(dcrit_mol_L*1000);
};
long double REFPROPMixtureBackend::calc_T_reducing(){
@@ -543,7 +564,7 @@ long double REFPROPMixtureBackend::calc_rhomolar_reducing(){
};
long double REFPROPMixtureBackend::calc_Ttriple(){
if (mole_fractions.size() != 1){throw ValueError("calc_Ttriple cannot be evaluated for mixtures");}
long icomp = 0;
long icomp = 1;
double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
INFOdll(&icomp, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
return static_cast<long double>(ttrp);
@@ -572,7 +593,7 @@ double REFPROPMixtureBackend::calc_melt_Tmax()
MELTPdll(&pmax_kPa, &(mole_fractions[0]),
&Tmax_melt,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return Tmax_melt;
}
@@ -586,7 +607,7 @@ long double REFPROPMixtureBackend::calc_melting_line(int param, int given, long
MELTTdll(&_T, &(mole_fractions[0]),
&p_kPa,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
else if (param == iT && given == iP){
@@ -594,7 +615,7 @@ long double REFPROPMixtureBackend::calc_melting_line(int param, int given, long
MELTPdll(&p_kPa, &(mole_fractions[0]),
&_T,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
else{
@@ -616,7 +637,7 @@ long double REFPROPMixtureBackend::calc_viscosity(void)
TRNPRPdll(&_T,&rhomol_L,&(mole_fractions[0]), // Inputs
&eta,&tcx, // Outputs
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
_viscosity = 1e-6*eta;
_conductivity = tcx;
@@ -636,7 +657,7 @@ long double REFPROPMixtureBackend::calc_surface_tension(void)
SURFTdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
&sigma, // Outputs
&ierr, herr, errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
_surface_tension = sigma;
return static_cast<double>(_surface_tension);
@@ -651,7 +672,7 @@ long double REFPROPMixtureBackend::calc_fugacity_coefficient(int i)
FUGCOFdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
&(fug_cof[0]), // Outputs
&ierr, herr, errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(fug_cof[i]);
}
@@ -662,7 +683,7 @@ void REFPROPMixtureBackend::calc_phase_envelope(const std::string &type)
char herr[255];
SATSPLNdll(&(mole_fractions[0]), // Inputs
&ierr, herr, errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); }
}
long double REFPROPMixtureBackend::calc_cpmolar_idealgas(void)
{
@@ -701,7 +722,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions[0]), // Saturation terms
&q,&emol,&hmol,&smol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength); //
if (ierr > 0) { throw ValueError(format("PT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("PT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = value1;
@@ -723,7 +744,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&hmol,&smol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("DmolarT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("DmolarT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
@@ -752,7 +773,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&hmol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("DmolarP: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("DmolarP: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_rhomolar = value1;
@@ -782,7 +803,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&smol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("DmolarHmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("DmolarHmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
@@ -812,7 +833,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&hmol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("DmolarSmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("DmolarSmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
@@ -842,7 +863,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&hmol,&hmol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("DmolarUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("DmolarUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
@@ -871,7 +892,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("HmolarPmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("HmolarPmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = value2;
@@ -901,7 +922,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&emol,&hmol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("PSmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("PSmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = value1;
@@ -932,7 +953,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&hmol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("PUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("PUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = value1;
@@ -961,7 +982,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&emol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("HmolarSmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("HmolarSmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000; // 1000 for conversion from kPa to Pa
@@ -992,7 +1013,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("SmolarUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("SmolarUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000; // 1000 for conversion from kPa to Pa
@@ -1031,7 +1052,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&emol,&hmol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("SmolarT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("SmolarT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000; // 1000 for conversion from kPa to Pa
@@ -1069,7 +1090,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&emol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("HmolarT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("HmolarT: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000; // 1000 for conversion from kPa to Pa
@@ -1107,7 +1128,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&q,&hmol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("TUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("TUmolar: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000; // 1000 for conversion from kPa to Pa
@@ -1145,7 +1166,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&emol,&hmol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) { throw ValueError(format("PQ: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
if (static_cast<int>(ierr) > 0) { throw ValueError(format("PQ: %s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = value1;
@@ -1176,7 +1197,7 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
&emol,&hmol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) {
if (static_cast<int>(ierr) > 0) {
throw ValueError(format("TQ(%s): %s",LoadedREFPROPRef.c_str(), herr).c_str());
}// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());

View File

@@ -483,7 +483,6 @@ TEST_CASE("Check inputs to PropsSI","[PropsSI]")
SECTION("Single state, trivial output, pure incompressible"){
CHECK(ValidNumber(CoolProp::PropsSI("Tmin","P",0,"T",0,"INCOMP::DowQ")));
};
std::cout << get_global_param_string("errstring");
SECTION("Bad input pair"){
CHECK(!ValidNumber(CoolProp::PropsSI("D","Q",0,"Q",0,"Water")));
};
@@ -499,6 +498,12 @@ TEST_CASE("Check inputs to PropsSI","[PropsSI]")
SECTION("Single state, single output, predefined mixture"){
CHECK(ValidNumber(CoolProp::PropsSI("T","P",101325,"Q",0,"R410A.mix")));
};
SECTION("Single state, single output, predefined mixture from REFPROP"){
CHECK(ValidNumber(CoolProp::PropsSI("T","P",101325,"Q",0,"REFPROP::R410A.mix")));
};
SECTION("Single state, single output, bad predefined mixture from REFPROP"){
CHECK(!ValidNumber(CoolProp::PropsSI("T","P",101325,"Q",0,"REFPROP::RRRRRR.mix")));
};
SECTION("Predefined mixture"){
std::vector<double> p(1, 101325), Q(1, 1.0), z;
std::vector<std::string> outputs(1,"T"); outputs.push_back("Dmolar");

View File

@@ -5,14 +5,6 @@
%ignore CoolProp::set_config_json(rapidjson::Document &);
%ignore CoolProp::get_config_as_json(rapidjson::Document &);
#ifdef SWIGMATLAB
%ignore configuration_keys;
%ignore CoolProp::ConfigurationItem::ConfigurationDataTypes;
%ignore CoolProp::SsatSimpleState::SsatSimpleStateEnum;
%ignore CoolProp::composition_types;
%ignore CoolProp::fluid_types;
#endif
%include "std_string.i" // This %include allows the use of std::string natively
%include "std_vector.i" // This allows for the use of STL vectors natively(ish)
%include "exception.i" //
@@ -46,10 +38,12 @@ namespace std {
#define SWIG
#include "Configuration.h"
#undef SWIG
#include "HumidAirProp.h"
%}
%include "DataStructures.h"
%include "AbstractState.h"
%include "CoolProp.h"
%include "PhaseEnvelope.h"
%include "Configuration.h"
%include "Configuration.h"
%include "HumidAirProp.h"

View File

@@ -228,7 +228,12 @@ EXPORT_CODE void CONVENTION set_debug_level(int level){
CoolProp::set_debug_level(level);
}
EXPORT_CODE long CONVENTION get_param_index(const char * param){
return CoolProp::get_parameter_index(param);
try{
return CoolProp::get_parameter_index(param);
}
catch(std::exception &){
return -1;
}
}
EXPORT_CODE long CONVENTION get_global_param_string(const char *param, char * Output, int n)
{
@@ -249,7 +254,13 @@ EXPORT_CODE long CONVENTION get_global_param_string(const char *param, char * Ou
}
EXPORT_CODE long CONVENTION get_parameter_information_string(const char *param, char * Output, int n)
{
int key = CoolProp::get_parameter_index(param);
int key;
try{
key = CoolProp::get_parameter_index(param);
}
catch(std::exception &){
return 0;
}
if (key >= 0){
std::string s = CoolProp::get_parameter_information(key, Output);
if (s.size() < static_cast<unsigned int>(n)){

View File

@@ -13,12 +13,14 @@
#include "crossplatform_shared_ptr.h"
#include "Exceptions.h"
#include <algorithm> // std::next_permutation
#include <stdlib.h>
#include "math.h"
#include "time.h"
#include "stdio.h"
#include <string.h>
#include <iostream>
#include <list>
/// This is a stub overload to help with all the strcmp calls below and avoid needing to rewrite all of them
std::size_t strcmp(const std::string &s, const std::string e){
@@ -166,14 +168,88 @@ static double Brent_HAProps_T(const std::string &OutputName, const std::string &
~BrentSolverResids(){};
double call(double T){
return HAPropsSI(OutputName,"T",T,Input1Name,Input1,Input2Name,Input2)-TargetVal;
double val = HAPropsSI(OutputName, "T", T, Input1Name, Input1, Input2Name, Input2);
return val - TargetVal;
}
};
BrentSolverResids BSR = BrentSolverResids(OutputName, Input1Name, Input1, Input2Name, Input2, TargetVal);
std::string errstr;
T = CoolProp::Brent(BSR,T_min,T_max,1e-7,1e-4,50,errstr);
// Now we need to check the bounds and make sure that they are ok (don't yield invalid output)
// and actually bound the solution
double r_min = BSR.call(T_min);
bool T_min_valid = ValidNumber(r_min);
double r_max = BSR.call(T_max);
bool T_max_valid = ValidNumber(r_max);
if (!T_min_valid && !T_max_valid){
throw CoolProp::ValueError(format("Both T_min [%g] and T_max [%g] yield invalid output values in Brent_HAProps_T",T_min,T_max).c_str());
}
else if (T_min_valid && !T_max_valid){
while (!T_max_valid){
// Reduce T_max until it works
T_max = 0.95*T_max + 0.05*T_min;
r_max = BSR.call(T_max);
T_max_valid = ValidNumber(r_max);
}
}
else if (!T_min_valid && T_max_valid){
while (!T_min_valid){
// Increase T_min until it works
T_min = 0.95*T_min + 0.05*T_max;
r_min = BSR.call(T_min);
T_min_valid = ValidNumber(r_min);
}
}
// We will do a secant call if the values at T_min and T_max have the same sign
if (r_min*r_max > 0){
if (std::abs(r_min) < std::abs(r_max)){
T = CoolProp::Secant(BSR, T_min, 0.01*T_min, 1e-7, 50, errstr);
}
else{
T = CoolProp::Secant(BSR, T_max, -0.01*T_max, 1e-7, 50, errstr);
}
}
else{
T = CoolProp::Brent(BSR, T_min, T_max, 1e-7, 1e-4, 50, errstr);
}
return T;
}
static double Secant_Tdb_at_saturated_W(double psi_w, double p, double T_guess)
{
double T;
class BrentSolverResids : public CoolProp::FuncWrapper1D
{
private:
double pp_water, psi_w, p, r;
public:
BrentSolverResids(double psi_w, double p) : psi_w(psi_w), p(p) { pp_water = psi_w*p; };
~BrentSolverResids(){};
double call(double T){
double p_ws;
if (T>=273.16){
// Saturation pressure [Pa]
Water->update(CoolProp::QT_INPUTS, 0, T);
p_ws= Water->keyed_output(CoolProp::iP);
}
else{
// Sublimation pressure [Pa]
p_ws=psub_Ice(T);
}
double f = f_factor(T, p);
double pp_water_calc = f*p_ws;
double psi_w_calc = pp_water_calc/p;
r = (psi_w_calc - psi_w)/psi_w;
return r;
}
};
BrentSolverResids Resids(psi_w, p);
std::string errstr;
T = CoolProp::Brent(Resids, 150, 350, 1e-16, 1e-7, 100, errstr);
return T;
}
@@ -207,7 +283,7 @@ static double Secant_HAProps_T(const std::string &OutputName, const std::string
static double Secant_HAProps_W(const std::string &OutputName, const std::string &Input1Name, double Input1, const std::string &Input2Name, double Input2, double TargetVal, double W_guess)
{
// Use a secant solve in order to yield a target output value for HAProps by altering humidity ratio
double x1=0,x2=0,x3=0,y1=0,y2=0,eps=1e-8,f=999,W=0.0001;
double x1=0,x2=0,x3=0,y1=0,y2=0,eps=1e-12,f=999,W=0.0001;
int iter=1;
while ((iter<=3 || std::abs(f)>eps) && iter<100)
@@ -1089,7 +1165,6 @@ double MoleFractionWater(double T, double p, int HumInput, double InVal)
{
// Sublimation pressure [Pa]
p_ws=psub_Ice(T);
}
// Enhancement Factor [-]
f=f_factor(T,p);
@@ -1211,8 +1286,8 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d
{
// Add a check to make sure that Air and Water fluid states have been properly instantiated
check_fluid_instantiation();
int In1Type, In2Type, In3Type,iT,iW,iTdp,iRH,ip,Type1,Type2;
int iT, iW, iTdp, iRH, ip;
givens In1Type, In2Type, In3Type, Type1, Type2, OutputType;
double vals[3],p,T,RH,W,Tdp,psi_w,M_ha,v_bar,h_bar,s_bar,MainInputValue,SecondaryInputValue,T_guess;
double Value1,Value2,W_guess;
std::string MainInputName, SecondaryInputName, Name1, Name2;
@@ -1220,11 +1295,17 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d
vals[0]=Input1;
vals[1]=Input2;
vals[2]=Input3;
OutputType=Name2Type(OutputName.c_str());
// First figure out what kind of inputs you have, convert names to Macro expansions
In1Type=Name2Type(Input1Name.c_str());
In2Type=Name2Type(Input2Name.c_str());
In3Type=Name2Type(Input3Name.c_str());
if (OutputType == In1Type){return Input1;}
if (OutputType == In2Type){return Input2;}
if (OutputType == In3Type){return Input3;}
// Pressure must be included
ip=TypeMatch(GIVEN_P,Input1Name,Input2Name,Input3Name);
@@ -1266,17 +1347,17 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d
// Temperature and pressure are known, figure out which variable holds the other value
if (In1Type!=GIVEN_T && In1Type!=GIVEN_P)
{
strcpy(SecondaryInputName,Input1Name);
strcpy(SecondaryInputName, Input1Name);
SecondaryInputValue=Input1;
}
else if (In2Type!=GIVEN_T && In2Type!=GIVEN_P)
{
strcpy(SecondaryInputName,Input2Name);
strcpy(SecondaryInputName, Input2Name);
SecondaryInputValue=Input2;
}
else if (In3Type!=GIVEN_T && In3Type!=GIVEN_P)
{
strcpy(SecondaryInputName,Input3Name);
strcpy(SecondaryInputName, Input3Name);
SecondaryInputValue=Input3;
}
else{
@@ -1327,43 +1408,70 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d
// Get the integer type codes
Type1=Name2Type(Name1);
Type2=Name2Type(Name2);
// First, if one of the inputs is something that can potentially yield
// an explicit solution at a given iteration of the solver, use it
if (Type1==GIVEN_RH || Type1==GIVEN_HUMRAT || Type1==GIVEN_TDP)
// First see if humidity ratio is provided, will be the fastest option
if (Type1 == GIVEN_HUMRAT)
{
// First input variable is a "nice" one
// MainInput is the one that you are using in the call to HAProps
MainInputValue=Value1;
strcpy(MainInputName,Name1);
MainInputValue=Value1; strcpy(MainInputName,Name1);
// SecondaryInput is the one that you are trying to match
SecondaryInputValue=Value2;
strcpy(SecondaryInputName,Name2);
SecondaryInputValue=Value2; strcpy(SecondaryInputName,Name2);
}
else if (Type2==GIVEN_RH || Type2==GIVEN_HUMRAT || Type2==GIVEN_TDP)
else if (Type2 == GIVEN_HUMRAT)
{
// Second input variable is a "nice" one
// MainInput is the one that you are using in the call to HAProps
MainInputValue=Value2;
strcpy(MainInputName,Name2);
// MainInput is the one that you are using in the call to HAPropsSI
MainInputValue=Value2; strcpy(MainInputName, Name2);
// SecondaryInput is the one that you are trying to match
SecondaryInputValue=Value1;
strcpy(SecondaryInputName,Name1);
SecondaryInputValue=Value1; strcpy(SecondaryInputName, Name1);
}
// Next, if one of the inputs is something that can potentially yield
// an explicit solution at a given iteration of the solver, use it
else if (Type1 == GIVEN_RH || Type1 == GIVEN_TDP)
{
// MainInput is the one that you are using in the call to HAProps
MainInputValue = Value1; strcpy(MainInputName, Name1);
// SecondaryInput is the one that you are trying to match
SecondaryInputValue = Value2; strcpy(SecondaryInputName, Name2);
}
else if (Type2 == GIVEN_RH || Type2 == GIVEN_TDP)
{
// MainInput is the one that you are using in the call to HAProps
MainInputValue = Value2; strcpy(MainInputName, Name2);
// SecondaryInput is the one that you are trying to match
SecondaryInputValue=Value1; strcpy(SecondaryInputName, Name1);
}
else
{
printf("Sorry, but currently at least one of the variables as an input to HAPropsSI() must be temperature, relative humidity, humidity ratio, or dewpoint\n Eventually will add a 2-D NR solver to find T and psi_w simultaneously, but not included now\n");
return -1000;
CoolProp::ValueError("Sorry, but currently at least one of the variables as an input to HAPropsSI() must be temperature, relative humidity, humidity ratio, or dewpoint\n Eventually will add a 2-D NR solver to find T and psi_w simultaneously, but not included now\n");
return _HUGE;
}
double T_min = 210;
double T_min = 200;
double T_max = 450;
if (Name2Type(MainInputName) == GIVEN_RH){
T_max = CoolProp::PropsSI("T","P",p,"Q",0,"Water") - 1;
if (MainInputValue < 1e-10){
T_max = 1000;
}
else{
T_max = CoolProp::PropsSI("T","P",p,"Q",0,"Water") - 1;
}
}
// Minimum drybulb temperature is the drybulb temperature corresponding to saturated air for the humidity ratio
// if the humidity ratio is provided
else if (Name2Type(MainInputName) == GIVEN_HUMRAT){
if (MainInputValue < 1e-10){
T_min = 135; // Around the critical point of dry air
T_max = 1000;
}
else{
// Calculate the saturated humid air water partial pressure;
double psi_w_sat = MoleFractionWater(T_min, p, GIVEN_HUMRAT, MainInputValue);
//double pp_water_sat = psi_w_sat*p; // partial pressure of water, which is equal to f*p_{w_s}
// Iteratively solve for temperature that will give desired pp_water_sat
T_min = Secant_Tdb_at_saturated_W(psi_w_sat, p, T_min);
}
}
try{
// Use the Brent's method solver to find T. Slow but reliable
@@ -1912,7 +2020,94 @@ TEST_CASE_METHOD(HAPropsConsistencyFixture, "ASHRAE RP1485 Tables", "[RP1485]")
CHECK(std::abs(actual/expected-1) < 0.01);
}
}
}
TEST_CASE("Assorted tests","[HAPropsSI]")
{
CHECK(ValidNumber(HumidAir::HAPropsSI("T", "H", 267769, "P", 104300, "W", 0.0)));
CHECK(ValidNumber(HumidAir::HAPropsSI("T", "B", 252.84, "W", 5.097e-4, "P", 101325)));
CHECK(ValidNumber(HumidAir::HAPropsSI("T", "B",290, "R", 1, "P", 101325)));
}
// a predicate implemented as a function:
bool is_not_a_pair (const std::set<std::size_t> &item) { return item.size() != 2; }
const int number_of_inputs = 6;
std::string inputs[number_of_inputs] = {"W","D","B","R","T","V"};//,"H","S"};
class ConsistencyTestData
{
public:
bool is_built;
std::vector<Dictionary> data;
std::list<std::set<std::size_t> > inputs_list;
ConsistencyTestData(){
is_built = false;
};
void build(){
if (is_built){return;}
std::vector<std::size_t> indices(number_of_inputs);
for (std::size_t i = 0; i < number_of_inputs; ++i){ indices[i] = i;}
// Generate a powerset of all the permutations of all lengths of inputs
std::set<std::size_t> indices_set(indices.begin(), indices.end());
std::set<std::set<std::size_t> > inputs_powerset = powerset(indices_set);
inputs_list = std::list<std::set<std::size_t> >(inputs_powerset.begin(), inputs_powerset.end());
inputs_list.remove_if(is_not_a_pair);
const int NT = 10, NW = 5;
double p = 101325;
for (double T = 210; T < 350; T += (350-210)/(NT-1))
{
double Wsat = HumidAir::HAPropsSI("W", "T", T, "P", p, "R", 1.0);
for (double W = 1e-5; W < Wsat; W += (Wsat-1e-5)/(NW-1)){
Dictionary vals;
// Calculate all the values using T, W
for (int i = 0; i < number_of_inputs; ++i){
double v = HumidAir::HAPropsSI(inputs[i], "T", T, "P", p, "W", W);
vals.add_number(inputs[i], v);
}
data.push_back(vals);
std::cout << format("T %g W %g\n",T,W);
}
}
is_built = true;
};
} consistency_data;
/*
* This test is incredibly slow, which is why it is currently commented out. Many of the tests also fail
*
TEST_CASE("HAPropsSI", "[HAPropsSI]")
{
consistency_data.build();
double p = 101325;
for (std::size_t i = 0; i < consistency_data.data.size(); ++i)
{
for (std::list<std::set<std::size_t> >::iterator iter = consistency_data.inputs_list.begin(); iter != consistency_data.inputs_list.end(); ++iter)
{
std::vector<std::size_t> pair(iter->begin(), iter->end());
std::string i0 = inputs[pair[0]], i1 = inputs[pair[1]];
double v0 = consistency_data.data[i].get_double(i0), v1 = consistency_data.data[i].get_double(i1);
if ((i0 == "B" && i1 == "V") || (i1 == "B" && i0 == "V")){continue;}
std::ostringstream ss2;
ss2 << "Inputs: \"" << i0 << "\"," << v0 << ",\"" << i1 << "\"," << v1;
SECTION(ss2.str(), ""){
double T = consistency_data.data[i].get_double("T");
double W = consistency_data.data[i].get_double("W");
double Wcalc = HumidAir::HAPropsSI("W", i0, v0, i1, v1, "P", p);
double Tcalc = HumidAir::HAPropsSI("T", i0, v0, i1, v1, "P", p);
std::string err = CoolProp::get_global_param_string("errstring");
CAPTURE(T);
CAPTURE(W);
CAPTURE(Tcalc);
CAPTURE(Wcalc);
CAPTURE(err);
CHECK(std::abs(Tcalc - T) < 1e-1);
CHECK(std::abs((Wcalc - W)/W) < 1e-3);
}
}
}
}
*/
#endif /* CATCH_ENABLED */