From cb497b443cbbf60be06b3485369837a785fb34c9 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 3 Jun 2014 20:58:24 +0200 Subject: [PATCH] Added crossplatform_shared_ptr to allow shared_ptr cleanly on all architectures without needing variable namespace definition --- include/crossplatform_shared_ptr.h | 20 ++ src/Backends/Helmholtz/ExcessHEFunction.h | 1 + .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 2 +- .../REFPROP/REFPROPMixtureBackend.cpp | 104 ++++----- src/CoolProp.cpp | 48 +++-- src/Helmholtz.cpp | 113 +++++----- src/HumidAirProp.cpp | 200 +++++++++--------- src/PolyMath.cpp | 1 - src/SpeedTest.cpp | 1 + src/Tests/CoolProp-Tests.cpp | 32 +-- src/Tests/Tests.cpp | 6 +- src/main.cxx | 6 +- 12 files changed, 280 insertions(+), 254 deletions(-) create mode 100644 include/crossplatform_shared_ptr.h diff --git a/include/crossplatform_shared_ptr.h b/include/crossplatform_shared_ptr.h new file mode 100644 index 00000000..b375c26b --- /dev/null +++ b/include/crossplatform_shared_ptr.h @@ -0,0 +1,20 @@ +#ifndef CROSSPLATFORM_SHARED_PTR +#define CROSSPLATFORM_SHARED_PTR + +#include "PlatformDetermination.h" + +// Based on the platform and compiler, include the necessary header to give access to std::tr1::shared_ptr + +#if defined(__ISLINUX__) +#include +#elif defined(__ISAPPLE__) +#include +#elif defined(__ISWINDOWS__) && defined(__MINGW32__) +#include +#else +#pragma error +#endif + +using namespace std::tr1; + +#endif diff --git a/src/Backends/Helmholtz/ExcessHEFunction.h b/src/Backends/Helmholtz/ExcessHEFunction.h index f53bb1bf..80271a12 100644 --- a/src/Backends/Helmholtz/ExcessHEFunction.h +++ b/src/Backends/Helmholtz/ExcessHEFunction.h @@ -4,6 +4,7 @@ #include #include #include "CoolPropFluid.h" +#include "crossplatform_shared_ptr.h" namespace CoolProp{ diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index a3caf2fe..e28acdf4 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -164,7 +164,7 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_background(long double et // Residual part long double B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this); long double rho = rhomolar(); - long double initial_part = eta_dilute*B_eta_initial*rhomolar(); + long double initial_part = eta_dilute*B_eta_initial*rho; // Higher order terms long double delta_eta_h; diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp index 01e5d206..d5671183 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp @@ -51,7 +51,7 @@ surface tension N/m #include #include #include - +#include "crossplatform_shared_ptr.h" #if defined(_MSC_VER) @@ -63,24 +63,24 @@ surface tension N/m #include #endif -// Some constants for REFPROP... defined by macros for ease of use +// Some constants for REFPROP... defined by macros for ease of use #define refpropcharlength 255 #define filepathlength 255 #define lengthofreference 3 #define errormessagelength 255 #define ncmax 20 // Note: ncmax is the max number of components -#define numparams 72 +#define numparams 72 #define maxcoefs 50 std::string LoadedREFPROPRef; -// Some constants for REFPROP... defined by macros for ease of use +// Some constants for REFPROP... defined by macros for ease of use #define refpropcharlength 255 #define filepathlength 255 #define lengthofreference 3 #define errormessagelength 255 #define ncmax 20 // Note: ncmax is the max number of components -#define numparams 72 +#define numparams 72 #define maxcoefs 50 // Check windows @@ -371,7 +371,7 @@ bool load_REFPROP() // INCREASE ROBUSTNESS. ALWAYS THROW AN ERROR ON THE ELSE. #error "Must define either ENV32BIT or ENV64BIT" #endif - + #elif defined(__ISLINUX__) RefpropdllInstance = dlopen ("librefprop.so", RTLD_LAZY); #elif defined(__ISAPPLE__) @@ -401,8 +401,8 @@ bool load_REFPROP() } #if defined(__ISWINDOWS__) - - // Get data associated with path using the windows libraries, + + // Get data associated with path using the windows libraries, // and if you can (result == 0), the path exists #ifdef __MINGW32__ struct stat buf; @@ -434,7 +434,7 @@ REFPROPMixtureBackend::REFPROPMixtureBackend(const std::vector& flu // Do the REFPROP instantiation for this fluid _mole_fractions_set = false; - // Try to add this fluid to REFPROP - might want to think about making array of + // Try to add this fluid to REFPROP - might want to think about making array of // components and setting mole fractions if they change a lot. this->set_REFPROP_fluids(fluid_names); @@ -521,8 +521,8 @@ void REFPROPMixtureBackend::set_REFPROP_fluids(const std::vector &f // Load REFPROP if it isn't loaded yet load_REFPROP(); - - // If the name of the refrigerant doesn't match + + // If the name of the refrigerant doesn't match // that of the currently loaded refrigerant if (LoadedREFPROPRef.compare(components_joined)) { @@ -593,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 (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} return Tmax_melt; } @@ -611,7 +611,7 @@ long double REFPROPMixtureBackend::calc_melt_p_T(long double T) MELTTdll(&_T, &(mole_fractions[0]), &p_kPa, &ierr,herr,errormessagelength); // Error message - if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } + if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} return p_kPa*1000; } @@ -632,7 +632,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 (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; @@ -652,7 +652,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 (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(_surface_tension); @@ -667,11 +667,11 @@ 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 (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} return static_cast(fug_cof[i]); } - + void REFPROPMixtureBackend::update(long input_pair, double value1, double value2) { double rho_mol_L=_HUGE, rhoLmol_L=_HUGE, rhoVmol_L=_HUGE, @@ -681,7 +681,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 char herr[255]; clear(); - + // Check that mole fractions have been set, etc. check_status(); @@ -689,8 +689,8 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 WMOLdll(&(mole_fractions[0]), &mm); // returns mole mass in kg/kmol _molar_mass = 0.001*mm; // [kg/mol] - - + + switch(input_pair) { case PT_INPUTS: @@ -724,9 +724,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 TDFLSHdll(&_T,&rho_mol_L,&(mole_fractions[0]),&p_kPa, &rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms &q,&emol,&hmol,&smol,&cvmol,&cpmol,&w, - &ierr,herr,errormessagelength); + &ierr,herr,errormessagelength); if (ierr > 0) { throw ValueError(format("%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; if (0) @@ -757,7 +757,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 if (ierr > 0) { throw ValueError(format("%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; + _rhomolar = value1; _p = value2; if (0) { @@ -783,9 +783,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 DHFLSHdll(&rho_mol_L,&hmol,&(mole_fractions[0]),&_T,&p_kPa, &rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms &q,&emol,&smol,&cvmol,&cpmol,&w, - &ierr,herr,errormessagelength); + &ierr,herr,errormessagelength); if (ierr > 0) { throw ValueError(format("%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; if (0) @@ -813,9 +813,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 DSFLSHdll(&rho_mol_L,&smol,&(mole_fractions[0]),&_T,&p_kPa, &rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms &q,&emol,&hmol,&cvmol,&cpmol,&w, - &ierr,herr,errormessagelength); + &ierr,herr,errormessagelength); if (ierr > 0) { throw ValueError(format("%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; if (0) @@ -843,9 +843,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 DEFLSHdll(&rho_mol_L,&emol,&(mole_fractions[0]),&_T,&p_kPa, &rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms &q,&hmol,&hmol,&cvmol,&cpmol,&w, - &ierr,herr,errormessagelength); + &ierr,herr,errormessagelength); if (ierr > 0) { throw ValueError(format("%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; if (0) @@ -889,7 +889,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 { // Call again, but this time with molar units // H: [J/kg] * [kg/mol] -> [J/mol] - update(HmolarP_INPUTS, value1 * (double)_molar_mass, value2); + update(HmolarP_INPUTS, value1 * (double)_molar_mass, value2); return; } case PSmolar_INPUTS: @@ -919,7 +919,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 { // Call again, but this time with molar units // S: [J/kg/K] * [kg/mol] -> [J/mol/K] - update(PSmolar_INPUTS, value1, value2*(double)_molar_mass); + update(PSmolar_INPUTS, value1, value2*(double)_molar_mass); return; } case PUmolar_INPUTS: @@ -948,9 +948,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 } case PUmass_INPUTS: { - // Call again, but this time with molar units + // Call again, but this time with molar units // U: [J/kg] * [kg/mol] -> [J/mol] - update(PUmolar_INPUTS, value1, value2*(double)_molar_mass); + update(PUmolar_INPUTS, value1, value2*(double)_molar_mass); return; } case HmolarSmolar_INPUTS: @@ -977,16 +977,16 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 } case HmassSmass_INPUTS: { - // Call again, but this time with molar units + // Call again, but this time with molar units // H: [J/kg] * [kg/mol] -> [J/mol/K] // S: [J/kg/K] * [kg/mol] -> [J/mol/K] - update(HmolarSmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass); + update(HmolarSmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass); return; } case SmolarUmolar_INPUTS: { // Unit conversion for REFPROP - smol = value1; emol = value2; + smol = value1; emol = value2; // from REFPROP: subroutine ESFLSH (e,s,z,t,p,D,Dl,Dv,x,y,q,h,cv,cp,w,ierr,herr) ESFLSHdll(&emol,&smol,&(mole_fractions[0]),&_T,&p_kPa,&rho_mol_L, @@ -1008,10 +1008,10 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 } case SmassUmass_INPUTS: { - // Call again, but this time with molar units - // S: [J/kg/K] * [kg/mol] -> [J/mol/K], + // Call again, but this time with molar units + // S: [J/kg/K] * [kg/mol] -> [J/mol/K], // U: [J/kg] * [kg/mol] -> [J/mol] - update(SmolarUmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass); + update(SmolarUmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass); return; } case SmolarT_INPUTS: @@ -1047,9 +1047,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 } case SmassT_INPUTS: { - // Call again, but this time with molar units + // Call again, but this time with molar units // S: [J/kg/K] * [kg/mol] -> [J/mol/K] - update(SmolarT_INPUTS, value1 * (double)_molar_mass, value2 ); + update(SmolarT_INPUTS, value1 * (double)_molar_mass, value2 ); return; } case HmolarT_INPUTS: @@ -1085,9 +1085,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 } case HmassT_INPUTS: { - // Call again, but this time with molar units + // Call again, but this time with molar units // H: [J/kg] * [kg/mol] -> [J/mol] - update(HmolarT_INPUTS, value1 * (double)_molar_mass, value2 ); + update(HmolarT_INPUTS, value1 * (double)_molar_mass, value2 ); return; } case TUmolar_INPUTS: @@ -1123,9 +1123,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 } case TUmass_INPUTS: { - // Call again, but this time with molar units + // Call again, but this time with molar units // U: [J/kg] * [kg/mol] -> [J/mol] - update(TUmolar_INPUTS, value1, value2 * (double)_molar_mass); + update(TUmolar_INPUTS, value1, value2 * (double)_molar_mass); return; } case PQ_INPUTS: @@ -1194,7 +1194,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2 { throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str())); } - + }; // Set these common variables that are used in every flash calculation _hmolar = hmol; @@ -1222,7 +1222,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]") { std::string Name = (*it); std::string RPName = CoolProp::get_fluid_param_string((*it),"REFPROP_name"); - + // Skip fluids not in REFPROP if (RPName.find("N/A") == 0){continue;} @@ -1239,7 +1239,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]") CAPTURE(RPName); CAPTURE(rho_CP); CAPTURE(rho_RP); - + double DH = (rho_RP-rho_CP)/rho_RP; CHECK(fabs(DH) < 0.005); } @@ -1252,7 +1252,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]") { std::string Name = (*it); std::string RPName = CoolProp::get_fluid_param_string((*it),"REFPROP_name"); - + // Skip fluids not in REFPROP if (RPName.find("N/A") == 0){continue;} @@ -1270,7 +1270,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]") CAPTURE(cp_CP); CAPTURE(cp_RP); CAPTURE(0.9*Tr); - + double Dcp = (cp_RP-cp_CP)/cp_RP; CHECK(fabs(Dcp) < 0.005); } @@ -1283,7 +1283,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]") { std::string Name = (*it); std::string RPName = CoolProp::get_fluid_param_string((*it),"REFPROP_name"); - + // Skip fluids not in REFPROP if (RPName.find("N/A") == 0){continue;} @@ -1312,4 +1312,4 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]") } } -#endif \ No newline at end of file +#endif diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp index bdf4a6d4..98e830ce 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -45,8 +45,8 @@ int get_debug_level(void){return debug_level;} //int global_Phase = -1; -void set_warning_string(std::string warning){ - warning_string = warning; +void set_warning_string(std::string warning){ + warning_string = warning; } void set_error_string(std::string error){ error_string = error; @@ -77,7 +77,7 @@ void set_error_string(std::string error){ // // First check whether it is one of the Brines that does // // not have a pure-fluid equivalent in CoolProp // if ( -// strcmp(Ref,"HC-10")==0 || +// strcmp(Ref,"HC-10")==0 || // strncmp(Ref,"PG-",3)==0 || // strncmp(Ref,"EG-",3)==0 || // strncmp(Ref,"EA-",3)==0 || @@ -237,7 +237,7 @@ std::string extract_fractions(const std::string &fluid_string, std::vectoris_trivial_output(iOutput)) - { + { double val = State->trivial_keyed_output(iOutput); return val; };*/ @@ -339,7 +339,7 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop // Return the value return val; } - catch(...){ + catch(...){ throw; } } @@ -355,7 +355,7 @@ double PropsSI(const char *Output, const char *Name1, double Prop1, const char * double PropsSI(const std::string &Output, const std::string &Name1, double Prop1, const std::string &Name2, double Prop2, const std::string &Ref) { // In this function the error catching happens; - try{ + try{ // Extract fractions if they are encoded in the fluid string if (has_fractions_in_string(Ref)) { @@ -392,7 +392,7 @@ double Props1SI(std::string FluidName,std::string Output) return PropsSI(Output,"",0,"",0,FluidName); } - + //EXPORT_CODE double CONVENTION IProps(long iOutput, long iName1, double Prop1, long iName2, double Prop2, long iFluid) //{ // Prop1 = convert_from_unit_system_to_SI(iName1, Prop1, get_standard_unit_system()); @@ -408,10 +408,10 @@ double Props1SI(std::string FluidName,std::string Output) // if (get_debug_level()>3){ // std::cout << format("%s:%d: _CoolProp_Fluid_PropsSI(%d,%d,%g,%d,%g,%s)\n",__FILE__,__LINE__,iOutput,iName1, Prop1, iName2, Prop2, pFluid->get_name().c_str()).c_str(); // } -// if (iName1 == iT){ -// T = Prop1;} +// if (iName1 == iT){ +// T = Prop1;} // else if (iName2 == iT){ -// T = Prop2;} +// T = Prop2;} // // // Generate a State instance wrapped around the Fluid instance // CoolPropStateClassSI CPS(pFluid); @@ -567,7 +567,7 @@ double Props1SI(std::string FluidName,std::string Output) // // // Update the class // CPS.update(iT,T,iD,rho); -// +// // // Get the output value // val = CPS.keyed_output(iTerm); // break; @@ -590,7 +590,7 @@ double Props1SI(std::string FluidName,std::string Output) // val = CPS.keyed_output(iTerm); // break; // } -// +// // default: // throw ValueError(format("Sorry DerivTerms is a work in progress, your derivative term [%d] is not available!",iTerm)); // } @@ -642,7 +642,8 @@ double Props1SI(std::string FluidName,std::string Output) void set_reference_stateS(std::string Ref, std::string reference_state) { std::tr1::shared_ptr HEOS; - HEOS.reset(new CoolProp::HelmholtzEOSMixtureBackend(std::vector(1, Ref))); + std::vector _comps(1, Ref); + HEOS.reset(new CoolProp::HelmholtzEOSMixtureBackend(_comps)); if (!reference_state.compare("IIR")) { @@ -687,7 +688,7 @@ void set_reference_stateS(std::string Ref, std::string reference_state) HEOS->get_components()[0]->pEOS->alpha0.EnthalpyEntropyOffset.set(0, 0, ""); } else - { + { throw ValueError(format("reference state string is invalid: [%s]",reference_state.c_str())); } } @@ -714,10 +715,10 @@ void set_reference_stateS(std::string Ref, std::string reference_state) //} std::string get_BibTeXKey(std::string Ref, std::string key) -{ +{ std::vector names(1, Ref); HelmholtzEOSMixtureBackend HEOS(names); - + if (!key.compare("EOS")){ return HEOS.get_components()[0]->pEOS->BibTeX_EOS; } else if (!key.compare("CP0")){ return HEOS.get_components()[0]->pEOS->BibTeX_CP0; } else if (!key.compare("VISCOSITY")){ return HEOS.get_components()[0]->transport.BibTeX_viscosity; } @@ -756,8 +757,9 @@ std::string get_global_param_string(std::string ParamName) std::string get_fluid_param_string(std::string FluidName, std::string ParamName) { try{ - std::tr1::shared_ptr HEOS(new CoolProp::HelmholtzEOSMixtureBackend(std::vector(1,FluidName))); - + std::vector comps(1,FluidName); + std::tr1::shared_ptr HEOS(new CoolProp::HelmholtzEOSMixtureBackend(comps)); + CoolProp::CoolPropFluid *fluid = HEOS->get_components()[0]; if (!ParamName.compare("aliases")) @@ -790,4 +792,4 @@ std::string get_fluid_param_string(std::string FluidName, std::string ParamName) } } -} /* namespace CoolProp */ \ No newline at end of file +} /* namespace CoolProp */ diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index 03a8f373..2bcef397 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -375,7 +375,7 @@ long double ResidualHelmholtzExponential::dTau3(const long double &tau, const lo void ResidualHelmholtzGaussian::to_json(rapidjson::Value &el, rapidjson::Document &doc) { el.AddMember("type","ResidualHelmholtzGaussian",doc.GetAllocator()); - rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType), + rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType), _eta(rapidjson::kArrayType), _epsilon(rapidjson::kArrayType), _beta(rapidjson::kArrayType), _gamma(rapidjson::kArrayType); for (unsigned int i=0; i<=N; ++i) { @@ -515,8 +515,8 @@ long double ResidualHelmholtzGaussian::dTau3(const long double &tau, const long void ResidualHelmholtzGERG2008Gaussian::to_json(rapidjson::Value &el, rapidjson::Document &doc) { el.AddMember("type","ResidualHelmholtzGERG2008Gaussian",doc.GetAllocator()); - rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType), - _eta(rapidjson::kArrayType), _epsilon(rapidjson::kArrayType), + rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType), + _eta(rapidjson::kArrayType), _epsilon(rapidjson::kArrayType), _beta(rapidjson::kArrayType), _gamma(rapidjson::kArrayType); for (unsigned int i=0; i<=N; ++i) { @@ -581,7 +581,7 @@ long double ResidualHelmholtzGERG2008Gaussian::dDelta3(const long double &tau, c { /** Term derived in sympy using - + from sympy import * n_i, d_i, t_i, tau, delta, eta, gamma, beta, epsilon = symbols('n_i, d_i, t_i, tau, delta, eta, gamma, beta, epsilon') psi = exp(-eta*(delta-epsilon)**2-beta*(delta-gamma)) @@ -653,8 +653,8 @@ long double ResidualHelmholtzGERG2008Gaussian::dTau3(const long double &tau, con void ResidualHelmholtzLemmon2005::to_json(rapidjson::Value &el, rapidjson::Document &doc) { el.AddMember("type","ResidualHelmholtzLemmon2005",doc.GetAllocator()); - rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), - _t(rapidjson::kArrayType), _l(rapidjson::kArrayType), + rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), + _t(rapidjson::kArrayType), _l(rapidjson::kArrayType), _m(rapidjson::kArrayType); for (unsigned int i=0;i<=N; ++i) { @@ -743,7 +743,7 @@ long double ResidualHelmholtzLemmon2005::dDelta2(const long double &tau, const l ResidualHelmholtzLemmon2005Element &el = elements[i]; long double ni = el.n, ti = el.t, di = el.d, md = el.md; int li = el.l, mi = el.m; - if (li != 0 && mi != 0){ + if (li != 0 && mi != 0){ long double pow_delta_li = pow(delta, li); long double pow_tau_mi = pow(tau, md); s[i] = ni*((di-li*pow_delta_li)*(di-1.0-li*pow_delta_li) - li*li*pow_delta_li)*exp(ti*log_tau+(di-2)*log_delta-pow_delta_li-pow_tau_mi); @@ -933,8 +933,8 @@ void ResidualHelmholtzNonAnalytic::to_json(rapidjson::Value &el, rapidjson::Docu { el.AddMember("type","ResidualHelmholtzNonAnalytic",doc.GetAllocator()); - rapidjson::Value _n(rapidjson::kArrayType), _a(rapidjson::kArrayType), _b(rapidjson::kArrayType), - _beta(rapidjson::kArrayType), _A(rapidjson::kArrayType), _B(rapidjson::kArrayType), + rapidjson::Value _n(rapidjson::kArrayType), _a(rapidjson::kArrayType), _b(rapidjson::kArrayType), + _beta(rapidjson::kArrayType), _A(rapidjson::kArrayType), _B(rapidjson::kArrayType), _C(rapidjson::kArrayType), _D(rapidjson::kArrayType); for (unsigned int i=0; i<=N; ++i) { @@ -969,7 +969,7 @@ long double ResidualHelmholtzNonAnalytic::base(const long double &tau, const lon long double theta=(1.0-tau)+Ai*pow(pow(delta-1.0,2),1.0/(2.0*betai)); long double DELTA=pow(theta,2)+Bi*pow(pow(delta-1.0,2),ai); long double PSI=exp(-Ci*pow(delta-1.0,2)-Di*pow(tau-1.0,2)); - + s[i] = ni*pow(DELTA, bi)*delta*PSI; } return std::accumulate(s.begin(), s.end(), 0.0); @@ -988,7 +988,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta(const long double &tau, const l long double PSI=exp(-Ci*pow(delta-1.0,2)-Di*pow(tau-1.0,2)); long double dPSI_dDelta=-2.0*Ci*(delta-1.0)*PSI; long double dDELTA_dDelta=(delta-1.0)*(Ai*theta*2.0/betai*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)+2.0*Bi*ai*pow(pow(delta-1.0,2),ai-1.0)); - + // At critical point, DELTA is 0, and 1/0^n is undefined if (fabs(DELTA) < 10*DBL_EPSILON) { @@ -1016,7 +1016,7 @@ long double ResidualHelmholtzNonAnalytic::dTau(const long double &tau, const lon long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI; long double dDELTA_dDelta=(delta-1.0)*(Ai*theta*2.0/betai*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)+2.0*Bi*ai*pow(pow(delta-1.0,2),ai-1.0)); long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0); - + // At critical point, DELTA is 0, and 1/0^n is undefined if (fabs(DELTA) < 10*DBL_EPSILON) { @@ -1048,7 +1048,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta2(const long double &tau, const long double dPSI2_dDelta2=(2.0*Ci*pow(delta-1.0,2)-1.0)*2.0*Ci*PSI; long double dDELTA2_dDelta2=1.0/(delta-1.0)*dDELTA_dDelta+pow(delta-1.0,2)*(4.0*Bi*ai*(ai-1.0)*pow(pow(delta-1.0,2),ai-2.0)+2.0*pow(Ai/betai,2)*pow(pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0),2)+Ai*theta*4.0/betai*(1.0/(2.0*betai)-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-2.0)); long double dDELTAbi2_dDelta2=bi*(pow(DELTA,bi-1.0)*dDELTA2_dDelta2+(bi-1.0)*pow(DELTA,bi-2.0)*pow(dDELTA_dDelta,2)); - + // At critical point, DELTA is 0, and 1/0^n is undefined if (fabs(DELTA) < 10*DBL_EPSILON) { @@ -1079,7 +1079,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta_dTau(const long double &tau, co long double dPSI2_dDelta_dTau=4.0*Ci*Di*(delta-1.0)*(tau-1.0)*PSI; long double dDELTAbi2_dDelta_dTau=-Ai*bi*2.0/betai*pow(DELTA,bi-1.0)*(delta-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)-2.0*theta*bi*(bi-1.0)*pow(DELTA,bi-2.0)*dDELTA_dDelta; - + long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI; long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0); @@ -1135,12 +1135,12 @@ long double ResidualHelmholtzNonAnalytic::dDelta3(const long double &tau, const long double dPSI2_dDelta2=(2.0*Ci*pow(delta-1.0,2)-1.0)*2.0*Ci*PSI; long double dDELTA2_dDelta2=1.0/(delta-1.0)*dDELTA_dDelta+pow(delta-1.0,2)*(4.0*Bi*ai*(ai-1.0)*pow(pow(delta-1.0,2),ai-2.0)+2.0*pow(Ai/betai,2)*pow(pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0),2)+Ai*theta*4.0/betai*(1.0/(2.0*betai)-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-2.0)); long double dDELTAbi2_dDelta2=bi*(pow(DELTA,bi-1.0)*dDELTA2_dDelta2+(bi-1.0)*pow(DELTA,bi-2.0)*pow(dDELTA_dDelta,2)); - + long double dtheta_dDelta = Ai/(2*betai)*pow(pow(delta-1,2),1/(2*betai)-1)*2*(delta-1); long double dPSI3_dDelta3=2.0*Ci*PSI*(-4*Ci*Ci*pow(delta-1.0,3)+6*Ci*(delta-1)); long double PI = 4*Bi*ai*(ai-1)*pow(pow(delta-1,2),ai-2)+2*pow(Ai/betai,2)*pow(pow(delta-1,2),1/betai-2)+4*Ai*theta/betai*(1/(2*betai)-1)*pow(pow(delta-1,2),1/(2*betai)-2); long double dPI_dDelta = -8*Bi*ai*(ai-1)*(ai-2)*pow(pow(delta-1,2),ai-5.0/2.0)-8*pow(Ai/betai,2)*(1/(2*betai)-1)*pow(pow(delta-1,2),1/betai-5.0/2.0)-(8*Ai*theta)/betai*(1/(2*betai)-1)*(1/(2*betai)-2)*pow(pow(delta-1,2),1/(2*betai)-5.0/2.0)+4*Ai/betai*(1/(2*betai)-1)*pow(pow(delta-1,2),1/(2*betai)-2)*dtheta_dDelta; - long double dDELTA3_dDelta3 = 1/(delta-1)*dDELTA2_dDelta2-1/pow(delta-1,2)*dDELTA_dDelta+pow(delta-1,2)*dPI_dDelta+2*(delta-1)*PI; + long double dDELTA3_dDelta3 = 1/(delta-1)*dDELTA2_dDelta2-1/pow(delta-1,2)*dDELTA_dDelta+pow(delta-1,2)*dPI_dDelta+2*(delta-1)*PI; long double dDELTAbi3_dDelta3 = bi*(pow(DELTA,bi-1)*dDELTA3_dDelta3+dDELTA2_dDelta2*(bi-1)*pow(DELTA,bi-2)*dDELTA_dDelta+(bi-1)*(pow(DELTA,bi-2)*2*dDELTA_dDelta*dDELTA2_dDelta2+pow(dDELTA_dDelta,2)*(bi-2)*pow(DELTA,bi-3)*dDELTA_dDelta)); // At critical point, DELTA is 0, and 1/0^n is undefined @@ -1174,7 +1174,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta_dTau2(const long double &tau, c long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0); long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI; - + long double dtheta_dDelta = Ai/(2*betai)*pow(pow(delta-1,2),1/(2*betai)-1)*2*(delta-1); long double dPSI2_dDelta_dTau=4.0*Ci*Di*(delta-1.0)*(tau-1.0)*PSI; @@ -1239,7 +1239,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta2_dTau(const long double &tau, c long double dDELTA_dTau = -2*theta; long double dDELTA2_dDelta_dTau = 2.0*Ai/(betai)*pow(pow(delta-1,2),1.0/(2.0*betai)-0.5); long double dDELTA3_dDelta2_dTau = 2.0*Ai*(betai-1)/(betai*betai)*pow(pow(delta-1,2),1/(2*betai)-1.0); - + long double dDELTAbim1_dTau = (bi-1)*pow(DELTA,bi-2)*dDELTA_dTau; long double dDELTAbim2_dTau = (bi-2)*pow(DELTA,bi-3)*dDELTA_dTau; long double Line11 = dDELTAbim1_dTau*dDELTA2_dDelta2 + pow(DELTA,bi-1)*dDELTA3_dDelta2_dTau; @@ -1288,7 +1288,7 @@ void ResidualHelmholtzSAFTAssociating::to_json(rapidjson::Value &el, rapidjson:: long double ResidualHelmholtzSAFTAssociating::Deltabar(const long double &tau, const long double &delta) { return this->g(this->eta(delta))*(exp(this->epsilonbar*tau)-1)*this->kappabar; -} +} long double ResidualHelmholtzSAFTAssociating::dDeltabar_ddelta__consttau(const long double &tau, const long double &delta) { return this->dg_deta(this->eta(delta))*(exp(this->epsilonbar*tau)-1)*this->kappabar*this->vbarn; @@ -1382,12 +1382,12 @@ long double ResidualHelmholtzSAFTAssociating::d2X_ddelta2(const long double &tau long double X = this->X(delta, Deltabar); long double alpha = this->dDeltabar_ddelta__consttau(tau, delta); long double dalpha_ddelta = this->d2Deltabar_ddelta2__consttau(tau, delta); - + long double dX_ddelta_constall = X*X*(2*Deltabar*Deltabar*X-alpha)/pow(2*Deltabar*delta*X+1,2); long double d_dXddelta_dX = -(Deltabar+delta*alpha)*2*(Deltabar*delta*X*X+X)/pow(2*Deltabar*delta*X+1,2); long double d_dXddelta_dDeltabar = X*X*(2*delta*delta*X*alpha-1)/pow(2*Deltabar*delta*X+1,2); long double d_dXddelta_dalpha = -delta*X*X/(2*Deltabar*delta*X+1); - + long double dX_dDeltabar = this->dX_dDeltabar__constdelta(delta, Deltabar); long double dX_ddelta = this->dX_ddelta__constDeltabar(delta, Deltabar); @@ -1397,7 +1397,7 @@ long double ResidualHelmholtzSAFTAssociating::d2X_ddelta2(const long double &tau + d_dXddelta_dDeltabar*alpha + d_dXddelta_dalpha*dalpha_ddelta); return val; -} +} long double ResidualHelmholtzSAFTAssociating::d3X_dtau3(const long double &tau, const long double &delta) { long double Delta = this->Deltabar(tau, delta); @@ -1475,7 +1475,7 @@ long double ResidualHelmholtzSAFTAssociating::d3X_ddelta3(const long double &tau long double ResidualHelmholtzSAFTAssociating::g(const long double &eta) { return 0.5*(2-eta)/pow(1-eta,(int)3); -} +} long double ResidualHelmholtzSAFTAssociating::dg_deta(const long double &eta) { return 0.5*(5-2*eta)/pow(1-eta,(int)4); @@ -1483,11 +1483,11 @@ long double ResidualHelmholtzSAFTAssociating::dg_deta(const long double &eta) long double ResidualHelmholtzSAFTAssociating::d2g_deta2(const long double &eta) { return 3*(3-eta)/pow(1-eta,(int)5); -} +} long double ResidualHelmholtzSAFTAssociating::d3g_deta3(const long double &eta) { return 6*(7-2*eta)/pow(1-eta,(int)6); -} +} long double ResidualHelmholtzSAFTAssociating::eta(const long double &delta){ return this->vbarn*delta; } @@ -1592,7 +1592,7 @@ void IdealHelmholtzCP0PolyT::to_json(rapidjson::Value &el, rapidjson::Document & el.AddMember("T0", static_cast(T0), doc.GetAllocator()); } long double IdealHelmholtzCP0PolyT::base(const long double &tau, const long double &delta) throw() -{ +{ double sum=0; for (std::size_t i = 0; i < N; ++i){ if (fabs(t[i])<10*DBL_EPSILON) @@ -1680,7 +1680,7 @@ void IdealHelmholtzCP0AlyLee::to_json(rapidjson::Value &el, rapidjson::Document el.AddMember("T0",static_cast(T0),doc.GetAllocator()); } long double IdealHelmholtzCP0AlyLee::base(const long double &tau, const long double &delta) throw() -{ +{ if (!enabled){ return 0.0;} return -tau*(anti_deriv_cp0_tau2(tau)-anti_deriv_cp0_tau2(tau0))+(anti_deriv_cp0_tau(tau)-anti_deriv_cp0_tau(tau0)); } @@ -1722,13 +1722,14 @@ IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffset; #ifdef ENABLE_CATCH #include #include "catch.hpp" +#include "crossplatform_shared_ptr.h" class HelmholtzConsistencyFixture { public: long double numerical, analytic; - - std::tr1::shared_ptr Lead, LogTau, IGPower, PlanckEinstein, + + std::tr1::shared_ptr Lead, LogTau, IGPower, PlanckEinstein, CP0Constant, CP0PolyT, Gaussian, Lemmon2005, Power, SAFT, NonAnalytic, Exponential, GERG2008; HelmholtzConsistencyFixture(){ @@ -1750,13 +1751,13 @@ public: } CP0Constant.reset(new CoolProp::IdealHelmholtzCP0Constant(4/8.314472,300,250)); { - long double beta[] = {1.24, 0.821, 15.45, 2.21, 437, 0.743}, - d[] = {1, 1, 2, 2, 3, 3}, - epsilon[] = {0.6734, 0.9239, 0.8636, 1.0507, 0.8482, 0.7522}, - eta[] = {0.9667, 1.5154, 1.0591, 1.6642, 12.4856, 0.9662}, + long double beta[] = {1.24, 0.821, 15.45, 2.21, 437, 0.743}, + d[] = {1, 1, 2, 2, 3, 3}, + epsilon[] = {0.6734, 0.9239, 0.8636, 1.0507, 0.8482, 0.7522}, + eta[] = {0.9667, 1.5154, 1.0591, 1.6642, 12.4856, 0.9662}, gamma[] = {1.2827, 0.4317, 1.1217, 1.1871, 1.1243, 0.4203}, n[] = {1.2198, -0.4883, -0.0033293, -0.0035387, -0.51172, -0.16882}, - t[] = {1, 2.124, 0.4, 3.5, 0.5, 2.7}; + t[] = {1, 2.124, 0.4, 3.5, 0.5, 2.7}; Gaussian.reset(new CoolProp::ResidualHelmholtzGaussian(std::vector(n,n+sizeof(n)/sizeof(n[0])), std::vector(d,d+sizeof(d)/sizeof(d[0])), std::vector(t,t+sizeof(t)/sizeof(t[0])), @@ -1767,8 +1768,8 @@ public: )); } { - long double d[] = {1, 1, 1, 2, 4, 1, 1, 2, 2, 3, 4, 5, 1, 5, 1, 2, 3, 5}, - l[] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 2, 3, 3}, + long double d[] = {1, 1, 1, 2, 4, 1, 1, 2, 2, 3, 4, 5, 1, 5, 1, 2, 3, 5}, + l[] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 2, 3, 3}, m[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.7, 7, 6}, n[] = {5.28076, -8.67658, 0.7501127, 0.7590023, 0.01451899, 4.777189, -3.330988, 3.775673, -2.290919, 0.8888268, -0.6234864, -0.04127263, -0.08455389, -0.1308752, 0.008344962, -1.532005, -0.05883649, 0.02296658}, t[]= {0.669, 1.05,2.75, 0.956, 1, 2, 2.75, 2.38, 3.37, 3.47, 2.63, 3.45, 0.72, 4.23, 0.2, 4.5, 29, 24}; @@ -1780,9 +1781,9 @@ public: )); } { - long double d[] = {1, 1, 1, 3, 7, 1, 2, 5, 1, 1, 4, 2}, - l[] = {0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3}, - n[] = {1.0038, -2.7662, 0.42921, 0.081363, 0.00024174, 0.48246, 0.75542, -0.00743, -0.4146, -0.016558, -0.10644, -0.021704}, + long double d[] = {1, 1, 1, 3, 7, 1, 2, 5, 1, 1, 4, 2}, + l[] = {0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3}, + n[] = {1.0038, -2.7662, 0.42921, 0.081363, 0.00024174, 0.48246, 0.75542, -0.00743, -0.4146, -0.016558, -0.10644, -0.021704}, t[] = {0.25, 1.25, 1.5, 0.25, 0.875, 2.375, 2, 2.125, 3.5, 6.5, 4.75, 12.5}; Power.reset(new CoolProp::ResidualHelmholtzPower(std::vector(n, n+sizeof(n)/sizeof(n[0])), std::vector(d, d+sizeof(d)/sizeof(d[0])), @@ -1791,13 +1792,13 @@ public: )); } { - + long double a = 1, epsilonbar = 12.2735737, kappabar = 1.09117041e-05, m = 1.01871348, vbarn = 0.0444215309; SAFT.reset(new CoolProp::ResidualHelmholtzSAFTAssociating(a,m,epsilonbar,vbarn,kappabar)); } { long double n[] = {-0.666422765408, 0.726086323499, 0.0550686686128}, - A[] = {0.7, 0.7, 0.7}, B[] = {0.3, 0.3, 1}, C[] = {10, 10, 12.5}, D[] = {275, 275, 275}, + A[] = {0.7, 0.7, 0.7}, B[] = {0.3, 0.3, 1}, C[] = {10, 10, 12.5}, D[] = {275, 275, 275}, a[] = {3.5, 3.5, 3}, b[] = {0.875, 0.925, 0.875}, beta[] = {0.3, 0.3, 0.3}; NonAnalytic.reset(new CoolProp::ResidualHelmholtzNonAnalytic(std::vector(n, n+sizeof(n)/sizeof(n[0])), std::vector(a, a+sizeof(a)/sizeof(a[0])), @@ -1810,10 +1811,10 @@ public: )); } { - long double d[] = {2, 2, 2, 0, 0, 0}, - g[] = {1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788}, - l[] = {2, 2, 2, 2, 2, 2}, - n[] = {-3.821884669859, 8.30345065618981, -4.4832307260286, -1.02590136933231, 2.20786016506394, -1.07889905203761}, + long double d[] = {2, 2, 2, 0, 0, 0}, + g[] = {1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788}, + l[] = {2, 2, 2, 2, 2, 2}, + n[] = {-3.821884669859, 8.30345065618981, -4.4832307260286, -1.02590136933231, 2.20786016506394, -1.07889905203761}, t[] = {3, 4, 5, 3, 4, 5}; Exponential.reset(new CoolProp::ResidualHelmholtzExponential(std::vector(n, n+sizeof(n)/sizeof(n[0])), std::vector(d, d+sizeof(d)/sizeof(n[0])), @@ -1841,7 +1842,7 @@ public: } } - void call(std::string d, std::shared_ptr term, long double tau, long double delta, long double ddelta) + void call(std::string d, shared_ptr term, long double tau, long double delta, long double ddelta) { if (!d.compare("dTau")){return dTau(term,tau,delta,ddelta);} else if (!d.compare("dTau2")){return dTau2(term,tau,delta,ddelta);} @@ -1856,7 +1857,7 @@ public: throw CoolProp::ValueError("don't understand deriv type"); } } - std::shared_ptr get(std::string t) + shared_ptr get(std::string t) { if (!t.compare("Lead")){return Lead;} else if (!t.compare("LogTau")){return LogTau;} @@ -1876,55 +1877,55 @@ public: throw CoolProp::ValueError(format("don't understand helmholtz type: %s",t.c_str())); } } - void dTau(std::shared_ptr term, long double tau, long double delta, long double dtau){ + void dTau(shared_ptr term, long double tau, long double delta, long double dtau){ long double term_plus = term->base(tau + dtau, delta); long double term_minus = term->base(tau - dtau, delta); numerical = (term_plus - term_minus)/(2*dtau); analytic = term->dTau(tau, delta); }; - void dTau2(std::shared_ptr term, long double tau, long double delta, long double dtau){ + void dTau2(shared_ptr term, long double tau, long double delta, long double dtau){ long double term_plus = term->dTau(tau + dtau, delta); long double term_minus = term->dTau(tau - dtau, delta); numerical = (term_plus - term_minus)/(2*dtau); analytic = term->dTau2(tau, delta); }; - void dTau3(std::shared_ptr term, long double tau, long double delta, long double dtau){ + void dTau3(shared_ptr term, long double tau, long double delta, long double dtau){ long double term_plus = term->dTau2(tau + dtau, delta); long double term_minus = term->dTau2(tau - dtau, delta); numerical = (term_plus - term_minus)/(2*dtau); analytic = term->dTau3(tau, delta); }; - void dDelta(std::shared_ptr term, long double tau, long double delta, long double ddelta){ + void dDelta(shared_ptr term, long double tau, long double delta, long double ddelta){ long double term_plus = term->base(tau, delta + ddelta); long double term_minus = term->base(tau, delta - ddelta); numerical = (term_plus - term_minus)/(2*ddelta); analytic = term->dDelta(tau, delta); }; - void dDelta2(std::shared_ptr term, long double tau, long double delta, long double ddelta){ + void dDelta2(shared_ptr term, long double tau, long double delta, long double ddelta){ long double term_plus = term->dDelta(tau, delta + ddelta); long double term_minus = term->dDelta(tau, delta - ddelta); numerical = (term_plus - term_minus)/(2*ddelta); analytic = term->dDelta2(tau, delta); }; - void dDelta3(std::shared_ptr term, long double tau, long double delta, long double ddelta){ + void dDelta3(shared_ptr term, long double tau, long double delta, long double ddelta){ long double term_plus = term->dDelta2(tau, delta + ddelta); long double term_minus = term->dDelta2(tau, delta - ddelta); numerical = (term_plus - term_minus)/(2*ddelta); analytic = term->dDelta3(tau, delta); }; - void dDelta_dTau(std::shared_ptr term, long double tau, long double delta, long double ddelta){ + void dDelta_dTau(shared_ptr term, long double tau, long double delta, long double ddelta){ long double term_plus = term->dTau(tau, delta + ddelta); long double term_minus = term->dTau(tau, delta - ddelta); numerical = (term_plus - term_minus)/(2*ddelta); analytic = term->dDelta_dTau(tau, delta); }; - void dDelta_dTau2(std::shared_ptr term, long double tau, long double delta, long double ddelta){ + void dDelta_dTau2(shared_ptr term, long double tau, long double delta, long double ddelta){ long double term_plus = term->dTau2(tau, delta + ddelta); long double term_minus = term->dTau2(tau, delta - ddelta); numerical = (term_plus - term_minus)/(2*ddelta); analytic = term->dDelta_dTau2(tau, delta); }; - void dDelta2_dTau(std::shared_ptr term, long double tau, long double delta, long double ddelta){ + void dDelta2_dTau(shared_ptr term, long double tau, long double delta, long double ddelta){ long double term_plus = term->dDelta_dTau(tau, delta + ddelta); long double term_minus = term->dDelta_dTau(tau, delta - ddelta); numerical = (term_plus - term_minus)/(2*ddelta); diff --git a/src/HumidAirProp.cpp b/src/HumidAirProp.cpp index 735497a8..388bf5b5 100644 --- a/src/HumidAirProp.cpp +++ b/src/HumidAirProp.cpp @@ -10,8 +10,7 @@ #include "CoolPropTools.h" #include "Ice.h" #include "CoolProp.h" - - +#include "crossplatform_shared_ptr.h" #include #include "math.h" @@ -20,6 +19,7 @@ #include #include + std::tr1::shared_ptr Water, Air; namespace HumidAir @@ -105,19 +105,19 @@ void UseVirialCorrelations(int flag) if (flag==0 || flag==1) { FlagUseVirialCorrelations=flag; - } + } else { printf("UseVirialCorrelations takes an integer, either 0 (no) or 1 (yes)\n"); } - + } void UseIsothermCompressCorrelation(int flag) { if (flag==0 || flag==1) { FlagUseIsothermCompressCorrelation=flag; - } + } else { printf("UseIsothermCompressCorrelation takes an integer, either 0 (no) or 1 (yes)\n"); @@ -128,7 +128,7 @@ void UseIdealGasEnthalpyCorrelations(int flag) if (flag==0 || flag==1) { FlagUseIdealGasEnthalpyCorrelations=flag; - } + } else { printf("UseIdealGasEnthalpyCorrelations takes an integer, either 0 (no) or 1 (yes)\n"); @@ -154,7 +154,7 @@ static double Brent_HAProps_T(char *OutputName, char *Input1Name, double Input1, this->TargetVal = TargetVal; }; ~BrentSolverResids(){}; - + double call(double T){ return HAPropsSI(OutputName,(char *)"T",T,Input1Name,Input1,Input2Name,Input2)-TargetVal; } @@ -231,7 +231,7 @@ static double _B_aw(double T) static double _dB_aw_dT(double T) { check_fluid_instantiation(); - // Returns value in m^3/mol + // Returns value in m^3/mol double a[]={0,0.665687e2,-0.238834e3,-0.176755e3}; double b[]={0,-0.237,-1.048,-3.183}; double rhobarstar=1000,Tstar=100; @@ -301,11 +301,11 @@ static double B_m(double T, double psi_w) double B_aa,B_ww,B_aw; if (FlagUseVirialCorrelations==1) { - B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2) - +4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5) + B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2) + +4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5) -2.157430412913e-18*pow(T,6) +9.453830907795e-22*pow(T,7); - B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2) - +1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5) + B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2) + +1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5) -8.567808759123e-14*pow(T,6) +3.958203548563e-17*pow(T,7); } else @@ -313,7 +313,7 @@ static double B_m(double T, double psi_w) B_aa = B_Air(T); // [m^3/mol] B_ww = B_Water(T); // [m^3/mol] } - + B_aw=_B_aw(T); // [m^3/mol] return pow(1-psi_w,2)*B_aa+2*(1-psi_w)*psi_w*B_aw+psi_w*psi_w*B_ww; } @@ -342,11 +342,11 @@ static double C_m(double T, double psi_w) double C_aaa,C_www,C_aww,C_aaw; if (FlagUseVirialCorrelations) { - C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2) - -6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5) + C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2) + -6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5) +3.184306136120e-23*pow(T,6) -1.386043640106e-26*pow(T,7); - C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2) - +7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5) + C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2) + +7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5) -5.708156494894e-15*pow(T,6) +2.704605721778e-18*pow(T,7); } else @@ -362,7 +362,7 @@ static double C_m(double T, double psi_w) static double dC_m_dT(double T, double psi_w) { // dCm_dT has units of m^6/mol^2/K - + double Tj,tau_Air,tau_Water,dC_dT_aaa,dC_dT_www,dC_dT_aww,dC_dT_aaw; // NDG for fluid EOS for virial terms Tj=132.6312; @@ -391,7 +391,7 @@ static double HenryConstant(double T) { // Result has units of 1/Pa double p_ws,beta_N2,beta_O2,beta_Ar,beta_a,tau,Tr,Tc=647.096; - Tr=T/Tc; + Tr=T/Tc; tau=1-Tr; Water->update(CoolProp::QT_INPUTS, 1.0, T); p_ws = Water->keyed_output(CoolProp::iP); //[Pa] @@ -404,7 +404,7 @@ static double HenryConstant(double T) double isothermal_compressibility(double T, double p) { double k_T; - + if (T> 273.16) { if (FlagUseIsothermCompressCorrelation) @@ -421,7 +421,7 @@ double isothermal_compressibility(double T, double p) else { k_T = IsothermCompress_Ice(T,p); //[1/Pa] - } + } return k_T; } double f_factor(double T, double p) @@ -432,7 +432,7 @@ double f_factor(double T, double p) double p_ws,tau_Air,tau_Water,B_aa,B_aw,B_ww,C_aaa,C_aaw,C_aww,C_www, line1,line2,line3,line4,line5,line6,line7,line8,k_T,beta_H,LHS,RHS,psi_ws, vbar_ws; - + // Saturation pressure [Pa] if (T>273.16) { @@ -449,33 +449,33 @@ double f_factor(double T, double p) beta_H = 0; vbar_ws = dg_dp_Ice(T,p)*MM_Water(); //[m^3/mol] } - + k_T = isothermal_compressibility(T,p); //[1/Pa] - + // Hermann: In the iteration process of the enhancement factor in Eq. (3.25), k_T is set to zero for pw,s (T) > p. if (p_ws>p) { k_T=0; beta_H=0; } - + // NDG for fluid EOS for virial terms Tj=132.6312; tau_Air=Tj/T; tau_Water=Water->keyed_output(CoolProp::iT_reducing)/T; if (FlagUseVirialCorrelations) { - B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2) - +4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5) + B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2) + +4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5) -2.157430412913e-18*pow(T,6) +9.453830907795e-22*pow(T,7); - B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2) - +1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5) + B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2) + +1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5) -8.567808759123e-14*pow(T,6) +3.958203548563e-17*pow(T,7); - C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2) - -6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5) + C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2) + -6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5) +3.184306136120e-23*pow(T,6) -1.386043640106e-26*pow(T,7); - C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2) - +7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5) + C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2) + +7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5) -5.708156494894e-15*pow(T,6) +2.704605721778e-18*pow(T,7); } else @@ -488,7 +488,7 @@ double f_factor(double T, double p) B_aw = _B_aw(T); //[m^3/mol] C_aaw = _C_aaw(T); //[m^6/mol^2] C_aww = _C_aww(T); //[m^6/mol^2] - + // Use a little secant loop to find f iteratively // Start out with a guess value of 1 for f while ((iter<=3 || change>eps) && iter<100) @@ -496,12 +496,12 @@ double f_factor(double T, double p) if (iter==1){x1=1.00; f=x1;} if (iter==2){x2=1.00+0.000001; f=x2;} if (iter>2) {f=x2;} - + // Left-hand-side of Equation 3.25 LHS=log(f); // Eqn 3.24 psi_ws=f*p_ws/p; - + // All the terms forming the RHS of Eqn 3.25 line1=((1+k_T*p_ws)*(p-p_ws)-k_T*(p*p-p_ws*p_ws)/2.0)/(Rbar*T)*vbar_ws+log(1-beta_H*(1-psi_ws)*p); line2=pow(1-psi_ws,2)*p/(Rbar*T)*B_aa-2*pow(1-psi_ws,2)*p/(Rbar*T)*B_aw-(p-p_ws-pow(1-psi_ws,2)*p)/(Rbar*T)*B_ww; @@ -512,7 +512,7 @@ double f_factor(double T, double p) line7=(6*pow(1-psi_ws,2)*psi_ws*psi_ws*p*p)/pow(Rbar*T,2)*B_ww*B_aw-(3*pow(1-psi_ws,4)*p*p)/(2*pow(Rbar*T,2))*B_aa*B_aa; line8=-(2*pow(1-psi_ws,2)*psi_ws*(-2+3*psi_ws)*p*p)/pow(Rbar*T,2)*B_aw*B_aw-(p_ws*p_ws-(4-3*psi_ws)*pow(psi_ws,3)*p*p)/(2*pow(Rbar*T,2))*B_ww*B_ww; RHS=line1+line2+line3+line4+line5+line6+line7+line8; - + if (iter==1){y1=LHS-RHS;} if (iter>1) { @@ -554,9 +554,9 @@ double Viscosity(double T, double p, double psi_w) { /* Using the method of: - + P.T. Tsilingiris, 2009, Thermophysical and transport properties of humid air at temperature range between 0 and 100 oC, Energy Conversion and Management, 49, 1098-1010 - + but using the detailed measurements for pure fluid from IAPWS formulations */ double mu_a,mu_w,Phi_av,Phi_va,Ma,Mw; @@ -576,9 +576,9 @@ double Conductivity(double T, double p, double psi_w) { /* Using the method of: - + P.T. Tsilingiris, 2009, Thermophysical and transport properties of humid air at temperature range between 0 and 100 oC, Energy Conversion and Management, 49, 1098-1010 - + but using the detailed measurements for pure fluid from IAPWS formulations */ double mu_a,mu_w,k_a,k_w,Phi_av,Phi_va,Ma,Mw; @@ -602,14 +602,14 @@ double MolarVolume(double T, double p, double psi_w) // Output in m^3/mol int iter; double v_bar0, v_bar=0, R_bar=8.314472,x1=0,x2=0,x3,y1=0,y2,resid,eps,Bm,Cm; - + // ----------------------------- // Iteratively find molar volume // ----------------------------- - + // Start by assuming it is an ideal gas to get initial guess v_bar0=R_bar*T/p; - + //Bring outside the loop since not a function of v_bar Bm=B_m(T,psi_w); Cm=C_m(T,psi_w); @@ -620,10 +620,10 @@ double MolarVolume(double T, double p, double psi_w) if (iter==1){x1=v_bar0; v_bar=x1;} if (iter==2){x2=v_bar0+0.000001; v_bar=x2;} if (iter>2) {v_bar=x2;} - + // want v_bar in m^3/mol and R_bar in J/mol-K resid = (p-(R_bar)*T/v_bar*(1+Bm/v_bar+Cm/(v_bar*v_bar)))/p; - + if (iter==1){y1=resid;} if (iter>1) { @@ -653,7 +653,7 @@ double IdealGasMolarEntropy_Water(double T, double p) tau = Water->keyed_output(CoolProp::iT_reducing)/T; Water->update(CoolProp::DmolarT_INPUTS,p/(R_bar*T),T); sbar_w = R_bar*(tau*Water->keyed_output(CoolProp::idalpha0_dtau_constdelta)-Water->keyed_output(CoolProp::ialpha0)); //[kJ/kmol/K] - return sbar_w; + return sbar_w; } double IdealGasMolarEnthalpy_Air(double T, double vmolar) { @@ -672,15 +672,15 @@ double IdealGasMolarEnthalpy_Air(double T, double vmolar) double IdealGasMolarEntropy_Air(double T, double vmolar_a) { double sbar_0_Lem, tau, sbar_a, R_bar_Lemmon = 8.314510, T0=273.15, p0=101325, vmolar_a_0; - + // Ideal-Gas contribution to entropy of air sbar_0_Lem=-196.1375815; //[J/mol/K] - + // Tj and rhoj are given by 132.6312 and 302.5507652 respectively tau = 132.6312/T; //[no units] vmolar_a_0 = R_bar_Lemmon*T0/p0; //[m^3/mol] - + Air->update(CoolProp::DmolarT_INPUTS,1/vmolar_a_0,T); sbar_a=sbar_0_Lem+R_bar_Lemmon*(tau*Air->keyed_output(CoolProp::idalpha0_dtau_constdelta)-Air->keyed_output(CoolProp::ialpha0))+R_bar_Lemmon*log(vmolar_a/vmolar_a_0); //[J/mol/K] @@ -691,9 +691,9 @@ double IdealGasMolarEntropy_Air(double T, double vmolar_a) double MolarEnthalpy(double T, double p, double psi_w, double vmolar) { // In units of kJ/kmol - + // vbar (molar volume) in m^3/kg - + double hbar_0, hbar_a, hbar_w, hbar, R_bar=8.314472; // ---------------------------------------- // Enthalpy @@ -701,7 +701,7 @@ double MolarEnthalpy(double T, double p, double psi_w, double vmolar) // Constant for enthalpy // Not clear why getting rid of this term yields the correct values in the table, but enthalpies are equal to an additive constant, so not a big deal hbar_0=0.0;//2.924425468; //[kJ/kmol] - + if (FlagUseIdealGasEnthalpyCorrelations){ hbar_w = 2.7030251618E-03*T*T + 3.1994361015E+01*T + 3.6123174929E+04; hbar_a = 9.2486716590E-04*T*T + 2.8557221776E+01*T - 7.8616129429E+03; @@ -710,7 +710,7 @@ double MolarEnthalpy(double T, double p, double psi_w, double vmolar) hbar_w = IdealGasMolarEnthalpy_Water(T, vmolar); hbar_a = IdealGasMolarEnthalpy_Air(T, vmolar); } - + hbar = hbar_0+(1-psi_w)*hbar_a+psi_w*hbar_w+R_bar*T*((B_m(T,psi_w)-T*dB_m_dT(T,psi_w))/vmolar+(C_m(T,psi_w)-T/2.0*dC_m_dT(T,psi_w))/(vmolar*vmolar)); return hbar; //[J/mol] } @@ -730,7 +730,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar) // Serious typo in RP-1485 - should use total pressure rather than // reference pressure in density calculation for water vapor molar entropy - + // vbar (molar volume) in m^3/mol double x1=0,x2=0,x3=0,y1=0,y2=0,eps=1e-8,f=999,R_bar_Lem=8.314510; int iter=1; @@ -741,7 +741,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar) // Calculate vbar_a, the molar volume of dry air // B_m, C_m, etc. functions take care of the units - Baa = B_m(T,0); + Baa = B_m(T,0); B = B_m(T,psi_w); dBdT = dB_m_dT(T,psi_w); Caaa = C_m(T,0); @@ -749,7 +749,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar) dCdT = dC_m_dT(T,psi_w); vbar_a_guess = R_bar_Lem*T/p; //[m^3/mol] since p in [Pa] - + while ((iter<=3 || fabs(f)>eps) && iter<100) { if (iter==1){x1=vbar_a_guess; vbar_a=x1;} @@ -766,7 +766,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar) iter=iter+1; if (iter>100){ return _HUGE; } } - + if (FlagUseIdealGasEnthalpyCorrelations){ std::cout << "Not implemented" << std::endl; } @@ -788,7 +788,7 @@ double DewpointTemperature(double T, double p, double psi_w) int iter; double p_w,eps,resid,Tdp=0,x1=0,x2=0,x3,y1=0,y2,T0; double p_ws_dp,f_dp; - + // Make sure it isn't dry air, return an impossible temperature otherwise if ((1-psi_w)<1e-16) { @@ -798,9 +798,9 @@ double DewpointTemperature(double T, double p, double psi_w) // Iteratively find the dewpoint temperature // ------------------------------------------ - // The highest dewpoint temperature possible is the dry-bulb temperature. + // The highest dewpoint temperature possible is the dry-bulb temperature. // When they are equal, the air is saturated (R=1) - + p_w = psi_w*p; // 0.61165... is the triple point pressure of water in kPa @@ -813,14 +813,14 @@ double DewpointTemperature(double T, double p, double psi_w) } // A good guess for Tdp is that enhancement factor is unity, which yields // p_w_s = p_w, and get guess for T from saturation temperature - + iter=1; eps=1e-8; resid=999; while ((iter<=3 || fabs(resid)>eps) && iter<100) { if (iter==1){x1 = T0; Tdp=x1;} if (iter==2){x2 = x1 + 0.1; Tdp=x2;} if (iter>2) {Tdp=x2;} - + if (Tdp >= 273.16) { // Saturation pressure at dewpoint [kPa] @@ -836,7 +836,7 @@ double DewpointTemperature(double T, double p, double psi_w) f_dp=f_factor(Tdp,p); // Error between target and actual pressure [kPa] resid=p_w-p_ws_dp*f_dp; - + if (iter==1){y1=resid;} if (iter>1) { @@ -875,7 +875,7 @@ public: if (Twb > 273.16) { // Saturation pressure at wetbulb temperature [Pa] - Water->update(CoolProp::QT_INPUTS,0,Twb); + Water->update(CoolProp::QT_INPUTS,0,Twb); p_ws_wb= Water->keyed_output(CoolProp::iP); } else @@ -883,7 +883,7 @@ public: // Sublimation pressure at wetbulb temperature [kPa] p_ws_wb=psub_Ice(Twb); } - + // Vapor pressure p_s_wb = f_wb*p_ws_wb; // wetbulb humidity ratio @@ -902,7 +902,7 @@ public: h_w=h_Ice(Twb,_p); } // Mole masses of wetbulb and humid air - + M_ha_wb=MM_Water()*psi_wb+(1-psi_wb)*28.966; v_bar_wb=MolarVolume(Twb,_p,psi_wb); double RHS = (MolarEnthalpy(Twb,_p,psi_wb,v_bar_wb)*(1+W_s_wb)/M_ha_wb+(_W-W_s_wb)*h_w); @@ -934,7 +934,7 @@ double WetbulbTemperature(double T, double p, double psi_w) // ------------------------------------------ // Iteratively find the wetbulb temperature // ------------------------------------------ - // + // // If the temperature is less than the saturation temperature of water // for the given atmospheric pressure, the highest wetbulb temperature that is possible is the dry bulb // temperature @@ -953,17 +953,17 @@ double WetbulbTemperature(double T, double p, double psi_w) WetBulbSolver WBS(T, p, psi_w); std::string errstr; - + double return_val; try{ return_val = Secant(WBS,Tmax,0.0001,1e-8,50,errstr); - + // Solution obtained is out of range (T>Tmax) if (return_val > Tmax) {throw CoolProp::ValueError();} } catch(std::exception &) { - // The lowest wetbulb temperature that is possible for a given dry bulb temperature + // The lowest wetbulb temperature that is possible for a given dry bulb temperature // is the saturated air temperature which yields the enthalpy of dry air at dry bulb temperature try{ @@ -980,7 +980,7 @@ double WetbulbTemperature(double T, double p, double psi_w) return_val = _HUGE; } } - return return_val; + return return_val; } static int Name2Type(const char *Name) { @@ -1023,7 +1023,7 @@ int TypeMatch(int TypeCode, const char *Input1Name, const char *Input2Name, cons double MoleFractionWater(double T, double p, int HumInput, double InVal) { double p_ws,f,W,epsilon=0.621945,Tdp,p_ws_dp,f_dp,p_w_dp,p_s,RH; - + if (HumInput==GIVEN_HUMRAT) //(2) { W=InVal; @@ -1041,7 +1041,7 @@ 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); @@ -1107,21 +1107,21 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, { // 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; 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; char MainInputName[100], SecondaryInputName[100],Name1[100],Name2[100]; - + vals[0]=Input1; vals[1]=Input2; vals[2]=Input3; - + // First figure out what kind of inputs you have, convert names to Macro expansions In1Type=Name2Type(Input1Name); In2Type=Name2Type(Input2Name); In3Type=Name2Type(Input3Name); - + // Pressure must be included ip=TypeMatch(GIVEN_P,Input1Name,Input2Name,Input3Name); if (ip>0) @@ -1132,13 +1132,13 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, // ----------------------------------------------------------------------------------------------------- // Check whether the remaining values give explicit solution for mole fraction of water - nice and fast // ----------------------------------------------------------------------------------------------------- - + // Find the codes if they are there iT= TypeMatch(GIVEN_T,Input1Name,Input2Name,Input3Name); iRH= TypeMatch(GIVEN_RH,Input1Name,Input2Name,Input3Name); iW= TypeMatch(GIVEN_HUMRAT,Input1Name,Input2Name,Input3Name); iTdp= TypeMatch(GIVEN_TDP,Input1Name,Input2Name,Input3Name); - + if (iT>0) // Found T (or alias) as an input { T=vals[iT-1]; @@ -1189,9 +1189,9 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, else { // Need to iterate to find dry bulb temperature since temperature is not provided - + // Pick one input, and alter T to match the other input - + // Get the variables and their values that are NOT pressure for simplicity // because you know you need pressure as an input and you already have // its value in variable p @@ -1219,17 +1219,17 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, else{ return _HUGE; } - + // Get the integer type codes Type1=Name2Type(Name1); Type2=Name2Type(Name2); - - // First, if one of the inputs is something that can potentially yield + + // 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 input variable is a "nice" one - + // MainInput is the one that you are using in the call to HAProps MainInputValue=Value1; strcpy(MainInputName,Name1); @@ -1240,7 +1240,7 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, else if (Type2==GIVEN_RH || Type2==GIVEN_HUMRAT || Type2==GIVEN_TDP) { // 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); @@ -1266,8 +1266,8 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, T = Secant_HAProps_T(SecondaryInputName,(char *)"P",p,MainInputName,MainInputValue,SecondaryInputValue,T_guess); double val = HAPropsSI(SecondaryInputName,(char *)"T",T,(char *)"P",p,MainInputName,MainInputValue); if (!ValidNumber(T) || !ValidNumber(val) || !(T_min < T && T < T_max) || fabs(val-SecondaryInputValue)>1e-6) - { - throw CoolProp::ValueError(); + { + throw CoolProp::ValueError(); } else { @@ -1276,13 +1276,13 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, } catch (std::exception &){}; } - + if (T < 0) // No solution found using secant { // Use the Brent's method solver to find T T = Brent_HAProps_T(SecondaryInputName,(char *)"P",p,MainInputName,MainInputValue,SecondaryInputValue,T_min,T_max); } - + // If you want the temperature, return it if (Name2Type(OutputName)==GIVEN_T) return T; @@ -1294,7 +1294,7 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, } } M_ha=(1-psi_w)*0.028966+MM_Water()*psi_w; //[kg_ha/mol_ha] - + // ----------------------------------------------------------------- // Calculate and return the desired value for known set of T,p,psi_w // ----------------------------------------------------------------- @@ -1395,16 +1395,16 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1, double HAProps_Aux(const char* Name,double T, double p, double W, char *units) { // This function provides some things that are not usually needed, but could be interesting for debug purposes. - + // Requires W since it is nice and fast and always defined. Put a dummy value if you want something that doesn't use humidity - + // Takes temperature, pressure, and humidity ratio W as inputs; double psi_w,Tj,tau_Water,tau_Air,B_aa,C_aaa,B_ww,C_www,B_aw,C_aaw,C_aww,v_bar; - + Tj=132.6312; tau_Air=Tj/T; tau_Water=Water->keyed_output(CoolProp::iT_critical)/T; - + try{ if (!strcmp(Name,"Baa")) { @@ -1589,8 +1589,8 @@ double HAProps_Aux(const char* Name,double T, double p, double W, char *units) } double cair_sat(double T) { - // Air saturation specific heat - // Based on a correlation from EES, good from 250K to 300K. + // Air saturation specific heat + // Based on a correlation from EES, good from 250K to 300K. // No error bound checking is carried out // T: [K] // cair_s: [kJ/kg-K] @@ -1674,7 +1674,7 @@ TEST_CASE("Check HA Virials from Table A.2.1","[RP1485]") CHECK(fabs(HumidAir::_C_aww(0+273.15)/(-224234/1e12)-1) < 1e-3); CHECK(fabs(HumidAir::_C_aww(200+273.15)/(-8436.5/1e12)-1) < 1e-3); CHECK(fabs(HumidAir::_C_aww(350+273.15)/(-2486.9/1e12)-1) < 1e-3); - } + } } TEST_CASE("Enhancement factor from Table A.3","[RP1485]") { @@ -1800,7 +1800,7 @@ TEST_CASE_METHOD(HAPropsConsistencyFixture, "ASHRAE RP1485 Tables", "[RP1485]") CHECK(fabs(actual/expected-1) < 0.01); } } - + } #endif /* CATCH_ENABLED */ diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index dc7003cb..7b109054 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -998,7 +998,6 @@ double BaseExponential::expval(const std::vector< std::vector > &coeffic #include #include "catch.hpp" - class PolynomialConsistencyFixture { public: CoolProp::BasePolynomial poly; diff --git a/src/SpeedTest.cpp b/src/SpeedTest.cpp index ad5f5cd7..92ff3408 100644 --- a/src/SpeedTest.cpp +++ b/src/SpeedTest.cpp @@ -2,6 +2,7 @@ #include "SpeedTest.h" #include "AbstractState.h" #include "DataStructures.h" +#include "crossplatform_shared_ptr.h" #include diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index 71bd065b..8bc6dc53 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -9,6 +9,7 @@ #if defined(ENABLE_CATCH) +#include "crossplatform_shared_ptr.h" #include "catch.hpp" namespace TransportValidation{ @@ -22,7 +23,7 @@ public: vel(std::string fluid, std::string in1, double v1, std::string in2, double v2, std::string out, double expected, double tol) { this->in1 = in1; this->in2 = in2; this->fluid = fluid; - this->v1 = v1; this->v2 = v2; this->expected = expected; + this->v1 = v1; this->v2 = v2; this->expected = expected; this->tol = tol; }; }; @@ -218,7 +219,7 @@ public: void set_backend(std::string backend, std::string fluid_name){ pState.reset(CoolProp::AbstractState::factory(backend, fluid_name)); } - void set_pair(std::string &in1, double v1, std::string &in2, double v2){ + void set_pair(std::string &in1, double v1, std::string &in2, double v2){ double o1, o2; long iin1 = CoolProp::get_parameter_index(in1); long iin2 = CoolProp::get_parameter_index(in2); @@ -470,24 +471,24 @@ TEST_CASE_METHOD(TransportValidationFixture, "Compare thermal conductivities aga static int inputs[] = { CoolProp::DmolarT_INPUTS, //CoolProp::SmolarT_INPUTS, - //CoolProp::HmolarT_INPUTS, + //CoolProp::HmolarT_INPUTS, //CoolProp::TUmolar_INPUTS, - CoolProp::DmolarP_INPUTS, - CoolProp::DmolarHmolar_INPUTS, - CoolProp::DmolarSmolar_INPUTS, + CoolProp::DmolarP_INPUTS, + CoolProp::DmolarHmolar_INPUTS, + CoolProp::DmolarSmolar_INPUTS, CoolProp::DmolarUmolar_INPUTS, - + /* CoolProp::HmolarP_INPUTS, CoolProp::PSmolar_INPUTS, - CoolProp::PUmolar_INPUTS, + CoolProp::PUmolar_INPUTS, */ /* - CoolProp::HmolarSmolar_INPUTS, - CoolProp::HmolarUmolar_INPUTS, - CoolProp::SmolarUmolar_INPUTS + CoolProp::HmolarSmolar_INPUTS, + CoolProp::HmolarUmolar_INPUTS, + CoolProp::SmolarUmolar_INPUTS */ }; @@ -503,7 +504,7 @@ public: void set_backend(std::string backend, std::string fluid_name){ pState.reset(CoolProp::AbstractState::factory(backend, fluid_name)); } - void set_pair(int pair){ + void set_pair(int pair){ this->pair = pair; } void set_TP(long double T, long double p) @@ -513,14 +514,13 @@ public: // Start with T,P as inputs, cycle through all the other pairs that are supported State.update(CoolProp::PT_INPUTS, p, T); - + // Set the other state variables rhomolar = State.rhomolar(); hmolar = State.hmolar(); smolar = State.smolar(); umolar = State.umolar(); } void get_variables() { - CoolProp::AbstractState &State = *pState; - + switch (pair) { /// In this group, T is one of the known inputs, iterate for the other one (easy) @@ -571,7 +571,7 @@ public: TEST_CASE_METHOD(ConsistencyFixture, "Test all input pairs for CO2 using all valid backends", "[]") { CHECK_NOTHROW(set_backend("HEOS", "CO2")); - + int inputsN = sizeof(inputs)/sizeof(inputs[0]); for (double p = 600000; p < 800000000.0; p *= 5) { diff --git a/src/Tests/Tests.cpp b/src/Tests/Tests.cpp index e2b77b9b..d7c46273 100644 --- a/src/Tests/Tests.cpp +++ b/src/Tests/Tests.cpp @@ -1,15 +1,15 @@ /** -This file includes some testing functions that will get built +This file includes some testing functions that will get built into the program. Otherwise CTest can be used by removing this file from the build to avoid double declaration of the main function and Catch clashing */ #include "Tests.h" -#if defined ENABLE_CATCH +#if defined ENABLE_CATCH #define CATCH_CONFIG_RUNNER #include "catch.hpp" - + static Catch::Session session; // There must be exactly one instance #endif // ENABLE_CATCH diff --git a/src/main.cxx b/src/main.cxx index 5f232609..bf940685 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -66,7 +66,8 @@ int main() { try{ //set_reference_stateS(NBP_refs[i],"RESET"); - HelmholtzEOSMixtureBackend HEOS(std::vector(1,NBP_refs[i])); + std::vector comps(1,NBP_refs[i]); + HelmholtzEOSMixtureBackend HEOS(comps); HEOS.update(PQ_INPUTS, 101325, 0); double delta_a1 = HEOS.smass()/(HEOS.gas_constant()/HEOS.molar_mass()); double delta_a2 = -HEOS.hmass()/(HEOS.gas_constant()/HEOS.molar_mass()*HEOS.get_reducing().T); @@ -81,7 +82,8 @@ int main() { try{ //set_reference_stateS(IIR_refs[i],"RESET"); - HelmholtzEOSMixtureBackend HEOS(std::vector(1,IIR_refs[i])); + std::vector comps(1,IIR_refs[i]); + HelmholtzEOSMixtureBackend HEOS(comps); HEOS.update(QT_INPUTS, 0, 273.15); double delta_a1 = (HEOS.smass()-1000)/(HEOS.gas_constant()/HEOS.molar_mass()); double delta_a2 = -(HEOS.hmass()-200000)/(HEOS.gas_constant()/HEOS.molar_mass()*HEOS.get_reducing().T);