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

Added handling for first partial derivative as a string - added tests - they pass

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-19 00:34:26 +02:00
parent e32dfcf255
commit 881a6c4616
7 changed files with 230 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
include/AbstractState.h
View File

@@ -300,7 +300,11 @@ public:
/** \brief The second partial derivative in homogeneous phases
*
* \sa \ref CoolProp::AbstractState::first_partial_deriv
* The first partial derivative (\ref CoolProp::AbstractState::first_partial_deriv) can be expressed as
*
* \f[ \left(\frac{\partial A}{\partial B}\right)_C = \frac{\left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta}{\left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta} = \frac{N}{D}\f]
*
* and the second derivative can be expressed as
*
* \f[
* \frac{\partial}{\partial D}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_E = \frac{\frac{\partial}{\partial \tau}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\delta\left(\frac{\partial E}{\partial \delta}\right)_\tau-\frac{\partial}{\partial \delta}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_\tau\left(\frac{\partial E}{\partial \tau}\right)_\delta}{\left(\frac{\partial D}{\partial \tau}\right)_\delta\left(\frac{\partial E}{\partial \delta}\right)_\tau-\left(\frac{\partial D}{\partial \delta}\right)_\tau\left(\frac{\partial E}{\partial \tau}\right)_\delta}
@@ -314,6 +318,8 @@ public:
* \f[\left(\frac{\partial D}{\partial \tau}\right)_{\delta} = \left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta\partial\tau}\right)+\left(\frac{\partial^2 B}{\partial \tau^2}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau^2}\right)_\delta-\left(\frac{\partial^2 B}{\partial \delta\partial\tau}\right)\left(\frac{\partial C}{\partial \tau}\right)_\delta\f]
* \f[\frac{\partial}{\partial \tau}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\delta = \frac{D\left(\frac{\partial N}{\partial \tau}\right)_{\delta}-N\left(\frac{\partial D}{\partial \tau}\right)_{\delta}}{D^2}\f]
* \f[\frac{\partial}{\partial \delta}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\tau = \frac{D\left(\frac{\partial N}{\partial \delta}\right)_{\tau}-N\left(\frac{\partial D}{\partial \delta}\right)_{\tau}}{D^2}\f]
*
* The terms \f$ N \f$ and \f$ D \f$ are the numerator and denominator from \ref CoolProp::AbstractState::first_partial_deriv respectively
*/
long double second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Of2, parameters Constant2){return calc_second_partial_deriv(Of1,Wrt1,Constant1,Of2,Constant2);};

5
include/DataStructures.h
View File

@@ -84,6 +84,11 @@ int get_parameter_index(const std::string &param_name);
/// Returns true if the input is trivial (constants, critical parameters, etc.)
bool is_trivial_parameter(int key);
/// Returns true if a valid parameter, and sets value in the variable iOutput
bool is_valid_parameter(const std::string & name, parameters & iOutput);
bool is_valid_first_derivative(const std::string & name, parameters &iOf, parameters &iWrt, parameters &iConstant);
std::string get_csv_parameter_list();
/// These are constants for the compositions

63
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View File

@@ -1298,6 +1298,69 @@ void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long doubl
throw ValueError(format("input to get_dtau_ddelta[%s] is invalid",get_parameter_information(index,"short").c_str()));
}
}
void get_dtau_ddelta_second_derivatives(HelmholtzEOSMixtureBackend *HEOS, long double T, long double rho, int index, long double &dtau2, long double &ddelta_dtau, long double &ddelta2)
{
// long double rhor = HEOS->get_reducing_state().rhomolar,
// Tr = HEOS->get_reducing_state().T,
// dT_dtau = -pow(T, 2)/Tr,
// R = HEOS->gas_constant(),
// delta = rho/rhor,
// tau = Tr/T;
//
// HelmholtzDerivatives derivs;
// components[0]->pEOS->alphar.all(tau, delta, derivs);
//
// switch (index)
// {
// case iT:
// dtau2 = d2T_dtau2; // d2T_dtau2
// ddelta_dtau = 0; // d2T_ddelta_dtau
// ddelta2 = 0; // d2T_ddelta2
// break;
// case iDmolar:
// dtau2 = 0; // d2rhomolar_dtau2
// ddelta2 = rhor;
// break;
// case iTau:
// dtau2 = 0; ddelta_dtau = 0; ddelta2 = 0; break;
// case iDelta:
// dtau2 = 0; ddelta_dtau = 0; ddelta2 = 0; break;
// case iP:
// {
// // dp/ddelta|tau
// d2pdrho2__T = R*T/rhomolar*(1 + 2*delta*dalphar_dDelta + pow(delta, 2)*d2alphar_dDelta2);
// ddelta2 =
// long double dalphar_dDelta = HEOS->calc_alphar_deriv_nocache(0, 1, HEOS->get_mole_fractions(), tau, delta);
// long double d2alphar_dDelta2 = HEOS->calc_alphar_deriv_nocache(0, 2, HEOS->get_mole_fractions(), tau, delta);
// long double d2alphar_dDelta_dTau = HEOS->calc_alphar_deriv_nocache(1, 1, HEOS->get_mole_fractions(), tau, delta);
//
// // dp/dtau|delta
// dtau = dT_dtau*rho*R*(1+delta*dalphar_dDelta-tau*delta*d2alphar_dDelta_dTau);
// break;
// }
// case iHmolar:
// // dh/dtau|delta
// dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau()));
// // dh/ddelta|tau
// ddelta = rhor*T*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()+delta*HEOS->dalphar_dDelta()+pow(delta,2)*HEOS->d2alphar_dDelta2());
// break;
// case iSmolar:
// // ds/dtau|delta
// dtau = dT_dtau*R/T*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()));
// // ds/ddelta|tau
// ddelta = rhor*R/rho*(-(1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau()));
// break;
// case iUmolar:
// // du/dtau|delta
// dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()));
// // du/ddelta|tau
// ddelta = rhor*HEOS->T()*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau());
// break;
//
// default:
// throw ValueError(format("input to get_dtau_ddelta[%s] is invalid",get_parameter_information(index,"short").c_str()));
// }
}
long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv_nocache(long double T, long double rhomolar, int Of, int Wrt, int Constant)
{
long double dOf_dtau, dOf_ddelta, dWrt_dtau, dWrt_ddelta, dConstant_dtau, dConstant_ddelta;

6
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
View File

@@ -220,8 +220,6 @@ public:
void p_phase_determination_pure_or_pseudopure(int other, long double value, bool &saturation_called);
void DmolarP_phase_determination();
// ***************************************************************
// ***************************************************************
@@ -236,10 +234,6 @@ public:
// ***************************************************************
// ***************************************************************
// ***************** MIXTURE DERIVATIVES ***********************

38
src/CoolProp.cpp
View File

@@ -334,15 +334,16 @@ std::string extract_fractions(const std::string &fluid_string, std::vector<doubl
// that error bubbling can be done properly
double _PropsSI(const std::string &Output, const std::string &Name1, double Prop1, const std::string &Name2, double Prop2, const std::string &backend, const std::string &Ref, const std::vector<double> &z)
{
parameters iOutput = iundefined_parameter;
parameters iOf = iundefined_parameter, iWrt = iundefined_parameter, iConstant = iundefined_parameter,
iOf1 = iundefined_parameter, iWrt1 = iundefined_parameter, iConstant1 = iundefined_parameter,
iWrt2 = iundefined_parameter, iConstant2 = iundefined_parameter;
double x1, x2;
if (get_debug_level()>5){
std::cout << format("%s:%d: _PropsSI(%s,%s,%g,%s,%g,%s,%s)\n",__FILE__,__LINE__,Output.c_str(),Name1.c_str(),Prop1,Name2.c_str(),Prop2,backend.c_str(),Ref.c_str(), vec_to_string(z).c_str()).c_str();
}
// Convert all the input and output parameters to integers
long iOutput = get_parameter_index(Output);
// The state we are going to use
shared_ptr<AbstractState> State;
@@ -363,13 +364,13 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop
// We are going to let the factory function load the state
State.reset(AbstractState::factory(backend, fluid_string));
// First check if it is a trivial input (critical/max parameters for instance)
if (is_trivial_parameter(iOutput))
if (is_valid_parameter(Output, iOutput) && is_trivial_parameter(iOutput))
{
double val = State->trivial_keyed_output(iOutput);
return val;
};
}
long iName1 = get_parameter_index(Name1);
long iName2 = get_parameter_index(Name2);
@@ -389,12 +390,25 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop
// Update the state
State->update(pair, x1, x2);
// Return the desired output
double val = State->keyed_output(iOutput);
// Return the value
return val;
if (iOutput != iundefined_parameter){
// Return the desired output
double val = State->keyed_output(iOutput);
// Return the value
return val;
}
else if (is_valid_first_derivative(Output, iOf, iWrt, iConstant)){
// Return the desired output
double val = State->first_partial_deriv(iOf, iWrt, iConstant);
// Return the value
return val;
}
else{// if is_valid_second_derivative(Output, iOf1, iWrt1, iConstant1, iWrt2, iConstant2){
return _HUGE;
};
}
double PropsSI(const std::string &Output, const std::string &Name1, double Prop1, const std::string &Name2, double Prop2, const std::string &Ref, const std::vector<double> &z)
{

56
src/DataStructures.cpp
View File

@@ -3,6 +3,7 @@
#include "DataStructures.h"
#include "Exceptions.h"
#include "CoolPropTools.h"
#include "CoolProp.h"
namespace CoolProp{
@@ -192,19 +193,62 @@ std::string get_csv_parameter_list()
}
return strjoin(strings, ",");
}
int get_parameter_index(const std::string &param_name)
bool is_valid_parameter(const std::string &param_name, parameters &iOutput)
{
std::map<std::string, int>::iterator it;
// Try to find it
it = parameter_information.index_map.find(param_name);
// If equal to end, not found
if (it != parameter_information.index_map.end())
{
if (it != parameter_information.index_map.end()){
// Found, return it
return it->second;
iOutput = static_cast<parameters>(it->second);
return true;
}
else
{
else{
return false;
}
}
bool is_valid_first_derivative(const std::string & name, parameters &iOf, parameters &iWrt, parameters &iConstant)
{
std::size_t iN0, iN1, iD0, iD1;
parameters Of, Wrt, Constant;
if (get_debug_level() > 5){std::cout << format("is_valid_first_derivative(%s)",name.c_str());}
// There should be exactly one /
// There should be exactly one |
// Suppose we start with "d(P)/d(T)|Dmolar"
std::vector<std::string> split_at_bar = strsplit(name, '|'); // "d(P)/d(T)" and "Dmolar"
if (split_at_bar.size() != 2){return false;}
std::vector<std::string> split_at_slash = strsplit(split_at_bar[0], '/'); // "d(P)" and "d(T)"
if (split_at_slash.size() != 2){return false;}
iN0 = split_at_slash[0].find("(");
iN1 = split_at_slash[0].find(")", iN0);
if (!(iN0 > 0 && iN1 > 0 && iN1 > iN0)){return false;}
std::string num = split_at_slash[0].substr(iN0+1, iN1-2);
iD0 = split_at_slash[1].find("(");
iD1 = split_at_slash[1].find(")", iD0);
if (!(iD0 > 0 && iD1 > 0 && iD1 > iD0)){return false;}
std::string den = split_at_slash[1].substr(iD0+1, iD1-2);
if (is_valid_parameter(num, Of) && is_valid_parameter(den, Wrt) && is_valid_parameter(split_at_bar[1], Constant)){
iOf = Of; iWrt = Wrt; iConstant = Constant; return true;
}
else{
return false;
}
}
int get_parameter_index(const std::string &param_name)
{
parameters iOutput;
if (is_valid_parameter(param_name, iOutput)){
return iOutput;
}
else{
throw ValueError(format("Your input name [%s] is not valid in get_parameter_index (names are case sensitive)",param_name.c_str()));
}
}

79
src/Tests/CoolProp-Tests.cpp
View File

@@ -824,6 +824,85 @@ TEST_CASE("Tests for solvers in P,H flash using Propane", "[flashdups],[flash],[
}
}
TEST_CASE("Multiple calls to state class are consistent", "[flashdups],[flash],[PH],[consistency]")
{
double hmolar, hmass;
SECTION("3 times PH with HEOS AbstractState yields same results every time","")
{
shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "n-Propane"));
CHECK_NOTHROW(AS->update(CoolProp::PT_INPUTS, 101325, 300));
hmolar = AS->hmolar();
hmass = AS->hmass();
CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
hmolar = AS->hmolar();
hmass = AS->hmass();
CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
hmolar = AS->hmolar();
hmass = AS->hmass();
CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
}
}
TEST_CASE("Test partial derivatives using PropsSI", "[derivatives]")
{
double hmolar, hmass, T = 300;
SECTION("Check drhodp|T 3 ways","")
{
shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "n-Propane"));
AS->update(CoolProp::PT_INPUTS, 101325, 300);
double drhomolardp__T_AbstractState = AS->first_partial_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iT);
double drhomolardp__T_PropsSI_num = (PropsSI("Dmolar","T",T,"P",101325+1e-3,"n-Propane") - PropsSI("Dmolar","T",T,"P",101325-1e-3,"n-Propane"))/(2*1e-3);
double drhomolardp__T_PropsSI = PropsSI("d(Dmolar)/d(P)|T","T",T,"P",101325,"n-Propane");
CAPTURE(drhomolardp__T_AbstractState);
CAPTURE(drhomolardp__T_PropsSI_num);
CAPTURE(drhomolardp__T_PropsSI);
double rel_err_exact = std::abs(drhomolardp__T_AbstractState-drhomolardp__T_PropsSI)/drhomolardp__T_PropsSI;
CHECK(rel_err_exact < 0.001);
}
SECTION("Invalid first partial derivatives","")
{
CHECK(!ValidNumber(PropsSI("d()/d(P)|T","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(Dmolar)/d()|T","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(Dmolar)/d(P)|","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(XXXX)/d(P)|T","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(Dmolar)d(P)|T","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(Dmolar)/d(P)T","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(Bvirial)/d(P)T","T",300,"P",101325,"n-Propane")));
CHECK(!ValidNumber(PropsSI("d(Tcrit)/d(P)T","T",300,"P",101325,"n-Propane")));
}
// SECTION{
//
// CHECK_NOTHROW(AS->update(CoolProp::PT_INPUTS, 101325, 300));
// hmolar = AS->hmolar();
// hmass = AS->hmass();
// CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
// CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
// hmolar = AS->hmolar();
// hmass = AS->hmass();
// CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
// CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
// hmolar = AS->hmolar();
// hmass = AS->hmass();
// CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
// CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
// hmolar = AS->hmolar();
// hmass = AS->hmass();
// CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
// CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
// hmolar = AS->hmolar();
// hmass = AS->hmass();
// CHECK_NOTHROW(AS->update(CoolProp::HmassP_INPUTS, hmass, 101325));
// CHECK_NOTHROW(AS->update(CoolProp::HmolarP_INPUTS, hmolar, 101325));
// }
}
//TEST_CASE("Test that states agree with CoolProp", "[states]")
//{
// std::vector<std::string> fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),',');