From 36d601490fcd9e6003fd415222b3f70205251ce0 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 10 Aug 2014 19:28:04 +0200 Subject: [PATCH] Lots more work on melting curves, can now make pretty phase diagram for water including the melting line --- dev/fluids/Water.json | 6 +- include/Ancillaries.h | 12 +- src/AbstractState.cpp | 5 + src/Backends/Helmholtz/Fluids/Ancillaries.cpp | 123 +++++++++++++----- src/Tests/CoolProp-Tests.cpp | 7 +- src/main.cxx | 23 +++- wrappers/Python/CoolProp5/AbstractState.pxd | 4 +- wrappers/Python/CoolProp5/AbstractState.pyx | 6 +- wrappers/Python/CoolProp5/CoolProp.pyx | 6 + wrappers/Python/CoolProp5/cAbstractState.pxd | 2 + wrappers/Python/CoolProp5/constants.pyx | 8 ++ .../Python/CoolProp5/constants_header.pxd | 8 ++ 12 files changed, 160 insertions(+), 50 deletions(-) diff --git a/dev/fluids/Water.json b/dev/fluids/Water.json index 094f6414..75240c53 100644 --- a/dev/fluids/Water.json +++ b/dev/fluids/Water.json @@ -55,9 +55,9 @@ "T_max": 251.165, "T_min": 273.16, "a": [ - -1195393.37e-6, - -80818.3159e-6, - -3338.2686e-6 + -1195393.37, + -80818.3159, + -3338.2686 ], "p_0": 611.657, "t": [ diff --git a/include/Ancillaries.h b/include/Ancillaries.h index 1c5fa259..864a0b10 100644 --- a/include/Ancillaries.h +++ b/include/Ancillaries.h @@ -114,7 +114,7 @@ public: struct MeltingLinePiecewiseSimonSegment { - long double T_0, a, c, p_0, T_max, T_min; + long double T_0, a, c, p_0, T_max, T_min, p_min, p_max; }; struct MeltingLinePiecewiseSimonData { @@ -123,7 +123,15 @@ struct MeltingLinePiecewiseSimonData struct MeltingLinePiecewisePolynomialInTrSegment { std::vector a, t; - long double T_0, p_0, T_max, T_min; + long double T_0, p_0, T_max, T_min, p_min, p_max; + long double evaluate(long double T) + { + long double summer = 0; + for (std::size_t i =0; i < a.size(); ++i){ + summer += a[i]*(pow(T/T_0,t[i])-1); + } + return p_0*(1+summer); + } }; struct MeltingLinePiecewisePolynomialInTrData { diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 878e6205..aa71e645 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -186,6 +186,11 @@ double AbstractState::trivial_keyed_output(int key) double AbstractState::keyed_output(int key) { if (get_debug_level()>=50) std::cout << format("AbstractState: keyed_output called for %s ",get_parameter_information(key,"short").c_str()) << std::endl; + // Handle trivial inputs + if (is_trivial_parameter(key)) + { + return trivial_keyed_output(key); + } switch (key) { case iQ: diff --git a/src/Backends/Helmholtz/Fluids/Ancillaries.cpp b/src/Backends/Helmholtz/Fluids/Ancillaries.cpp index aa2ad30b..038020b0 100644 --- a/src/Backends/Helmholtz/Fluids/Ancillaries.cpp +++ b/src/Backends/Helmholtz/Fluids/Ancillaries.cpp @@ -1,5 +1,13 @@ #include "Ancillaries.h" #include "DataStructures.h" +#include "AbstractState.h" + +#if defined(ENABLE_CATCH) + +#include "crossplatform_shared_ptr.h" +#include "catch.hpp" + +#endif namespace CoolProp{ @@ -103,20 +111,29 @@ double SaturationAncillaryFunction::invert(double value) void MeltingLineVariables::set_limits(void) { if (type == MELTING_LINE_SIMON_TYPE){ - MeltingLinePiecewiseSimonSegment &partmin = simon.parts[0]; - MeltingLinePiecewiseSimonSegment &partmax = simon.parts[simon.parts.size()-1]; - Tmin = partmin.T_0; - Tmax = partmax.T_max; - pmin = partmin.p_0; - pmax = evaluate(iP, iT, Tmax); + + // Fill in the min and max pressures for each part + for (std::size_t i = 0; i < simon.parts.size(); ++i){ + MeltingLinePiecewiseSimonSegment &part = simon.parts[i]; + part.p_min = part.p_0 + part.a*(pow(part.T_min/part.T_0,part.c)-1); + part.p_max = part.p_0 + part.a*(pow(part.T_max/part.T_0,part.c)-1); + } + pmin = simon.parts.front().p_min; + pmax = simon.parts.back().p_max; + Tmin = simon.parts.front().T_min; + Tmax = simon.parts.back().T_max; } else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE){ - MeltingLinePiecewisePolynomialInTrSegment &partmin = polynomial_in_Tr.parts[0]; - MeltingLinePiecewisePolynomialInTrSegment &partmax = polynomial_in_Tr.parts[polynomial_in_Tr.parts.size() - 1]; - Tmin = partmin.T_0; - Tmax = partmax.T_max; - pmin = partmin.p_0; - pmax = evaluate(iP, iT, Tmax); + // Fill in the min and max pressures for each part + for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){ + MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i]; + part.p_min = part.evaluate(part.T_min); + part.p_max = part.evaluate(part.T_max); + } + Tmin = polynomial_in_Tr.parts.front().T_min; + pmin = polynomial_in_Tr.parts.front().p_min; + Tmax = polynomial_in_Tr.parts.back().T_max; + pmax = polynomial_in_Tr.parts.back().p_max; } else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE){ MeltingLinePiecewisePolynomialInThetaSegment &partmin = polynomial_in_Theta.parts[0]; @@ -124,7 +141,7 @@ void MeltingLineVariables::set_limits(void) Tmin = partmin.T_0; Tmax = partmax.T_max; pmin = partmin.p_0; - pmax = evaluate(iP, iT, Tmax); + //pmax = evaluate(iP, iT, Tmax); } else{ throw ValueError("only Simon supported now"); @@ -151,11 +168,7 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value) for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){ MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i]; if (is_in_closed_range(part.T_min, part.T_max, T)){ - long double summer = 0; - for (std::size_t i =0; i < part.a.size(); ++i){ - summer += part.a[i]*(pow(T/part.T_0,part.t[i])-1); - } - return part.p_0*(1+summer); + return part.evaluate(T); } } throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Tr curve"); @@ -189,23 +202,21 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value) return T; } } - throw ValueError("unable to calculate melting line (p,T) for Simon curve"); + throw ValueError("unable to calculate melting line p(T) for Simon curve"); } - else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE || type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE) + else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE) { - class solver_resid : public FuncWrapper1D + class solver_resid : public FuncWrapper1D { public: - - MeltingLineVariables *line; + MeltingLinePiecewisePolynomialInTrSegment *part; long double r, given_p, calc_p, T; - solver_resid(MeltingLineVariables *line, long double p) : line(line), given_p(p){}; + solver_resid(MeltingLinePiecewisePolynomialInTrSegment *part, long double p) : part(part), given_p(p){}; double call(double T){ this->T = T; - // Calculate p using melting line - calc_p = line->evaluate(iP, iT, T); + calc_p = part->evaluate(T); // Difference between the two is to be driven to zero r = given_p - calc_p; @@ -213,17 +224,63 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value) return r; }; }; - solver_resid resid(this, value); - - double pmin = evaluate(iP, iT, Tmin); - double pmax = evaluate(iP, iT, Tmax); - std::string errstr; - return Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-12, 100, errstr); + + // Need to find the right segment + for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){ + MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i]; + if (is_in_closed_range(part.p_min, part.p_max, value)){ + std::string errstr; + solver_resid resid(&part, value); + double T = Brent(resid, part.T_min, part.T_max, DBL_EPSILON, 1e-12, 100, errstr); + return T; + } + } + throw ValueError("unable to calculate melting line T(p) for polynomial_in_Tr curve"); } else{ - throw ValueError(format("Invalid melting line type (T,p) [%d]",type)); + throw ValueError(format("Invalid melting line type T(p) [%d]",type)); } } } }; /* namespace CoolProp */ + +#if defined(ENABLE_CATCH) +TEST_CASE("Water melting line", "") +{ + shared_ptr AS(CoolProp::AbstractState::factory("HEOS","water")); + int iT = CoolProp::iT, iP = CoolProp::iP; + SECTION("Ice Ih-liquid") + { + double actual = AS->melting_line(iT, iP, 138.268e6); + double expected = 260.0; + CAPTURE(actual); + CAPTURE(expected); + CHECK(std::abs(actual-expected) < 0.01); + } + SECTION("Ice III-liquid") + { + double actual = AS->melting_line(iT, iP, 268.685e6); + double expected = 254; + CAPTURE(actual); + CAPTURE(expected); + CHECK(std::abs(actual-expected) < 0.01); + } + SECTION("Ice V-liquid") + { + double actual = AS->melting_line(iT, iP, 479.640e6); + double expected = 265; + CAPTURE(actual); + CAPTURE(expected); + CHECK(std::abs(actual-expected) < 0.01); + } + SECTION("Ice VI-liquid") + { + double actual = AS->melting_line(iT, iP, 1356.76e6); + double expected = 320; + CAPTURE(actual); + CAPTURE(expected); + CHECK(std::abs(actual-expected) < 1); + } +} +#endif \ No newline at end of file diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index 3aed36a2..69c46acd 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -778,8 +778,8 @@ TEST_CASE("Tests for solvers in P,Y flash using Water", "[flash],[PH],[PS],[PU]" { double Tc = Props1SI("Water","Tcrit"); double pc = Props1SI("Water","pcrit"); - double p = pc*1.3; - double T = Tc*0.7; + double p = pc*2; + double T = Tc*0.5; CAPTURE(T); CAPTURE(p); CHECK(ValidNumber(T)); @@ -793,9 +793,6 @@ TEST_CASE("Tests for solvers in P,Y flash using Water", "[flash],[PH],[PS],[PU]" CHECK(ValidNumber(T2)); } } - - - } #endif diff --git a/src/main.cxx b/src/main.cxx index 72c17a61..329307a8 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -18,7 +18,7 @@ using namespace CoolProp; #endif #include "SpeedTest.h" #include "HumidAirProp.h" -//#include "CoolPropLib.h" +#include "CoolPropLib.h" #include "crossplatform_shared_ptr.h" @@ -400,23 +400,36 @@ int main() } #endif #if 0 + { + double Tc = Props1SI("Water","Tcrit"); + double pc = Props1SI("Water","pcrit"); + double p = pc*2; + double T = Tc*0.5; + char ykey[] = "H"; + double y = PropsSI(ykey,"P",p,"T",T,"Water"); + double TT = PropsSI("T","P",p,ykey,y,"Water"); + int rr = 0; + } + #endif + #if 1 { run_tests(); char c; std::cin >> c; } #endif - #if 1 + #if 0 { char ykey[] = "H"; double Ts, y, T2, dT = -1; - double dd = PropsSI("T","P",101325,"T",114.357,"n-Propane"); shared_ptr AS(AbstractState::factory("HEOS","water")); - double ptt = AS->melting_line(iP, iT, 273.159); + double ptt = AS->melting_line(iT, iP, 138.268e6); - Ts = PropsSI("T","P",101325,"Q",0,"n-Propane"); + + + Ts = PropsSI("H","T",841.225,"P",2.86832e+007,"Water"); std::cout << get_global_param_string("errstring"); y = PropsSI(ykey,"T",Ts+dT,"P",101325,"n-Propane"); T2 = PropsSI("T",ykey,y,"P",101325,"n-Propane"); diff --git a/wrappers/Python/CoolProp5/AbstractState.pxd b/wrappers/Python/CoolProp5/AbstractState.pxd index b94e7288..56758969 100644 --- a/wrappers/Python/CoolProp5/AbstractState.pxd +++ b/wrappers/Python/CoolProp5/AbstractState.pxd @@ -27,4 +27,6 @@ cdef class AbstractState: cpdef double speed_sound(self) except * cpdef double molar_mass(self) except * - cpdef double keyed_output(self, long) except * \ No newline at end of file + cpdef double keyed_output(self, long) except * + + cpdef double melting_line(self, int, int, double) except * \ No newline at end of file diff --git a/wrappers/Python/CoolProp5/AbstractState.pyx b/wrappers/Python/CoolProp5/AbstractState.pyx index d2ab63e0..0d31d2c5 100644 --- a/wrappers/Python/CoolProp5/AbstractState.pyx +++ b/wrappers/Python/CoolProp5/AbstractState.pyx @@ -65,4 +65,8 @@ cdef class AbstractState: cpdef double molar_mass(self) except *: """ Get the molar mass kg/mol - wrapper of c++ function :cpapi:`AbstractState::molar_mass` """ - return self.thisptr.molar_mass() \ No newline at end of file + return self.thisptr.molar_mass() + + cpdef double melting_line(self, int param, int given, double value) except *: + """ Get values from the melting line - wrapper of c++ function :cpapi:`AbstractState::melting_line` """ + return self.thisptr.melting_line(param, given, value) \ No newline at end of file diff --git a/wrappers/Python/CoolProp5/CoolProp.pyx b/wrappers/Python/CoolProp5/CoolProp.pyx index d53f1bda..96d3d42b 100644 --- a/wrappers/Python/CoolProp5/CoolProp.pyx +++ b/wrappers/Python/CoolProp5/CoolProp.pyx @@ -12,6 +12,12 @@ cimport cython import math import warnings +try: + import numpy as np + _numpy_supported = True +except ImportError: + _numpy_supported = False + from libcpp.string cimport string from libcpp.vector cimport vector diff --git a/wrappers/Python/CoolProp5/cAbstractState.pxd b/wrappers/Python/CoolProp5/cAbstractState.pxd index 9a08fa88..b64da527 100644 --- a/wrappers/Python/CoolProp5/cAbstractState.pxd +++ b/wrappers/Python/CoolProp5/cAbstractState.pxd @@ -38,6 +38,8 @@ cdef extern from "AbstractState.h" namespace "CoolProp": double viscosity() except+ double conductivity() except+ double surface_tension() except+ + + double melting_line(int,int,double) except+ # The static factory method for the AbstractState cdef extern from "AbstractState.h" namespace "CoolProp::AbstractState": diff --git a/wrappers/Python/CoolProp5/constants.pyx b/wrappers/Python/CoolProp5/constants.pyx index 5d1c8aee..22bd88be 100644 --- a/wrappers/Python/CoolProp5/constants.pyx +++ b/wrappers/Python/CoolProp5/constants.pyx @@ -7,6 +7,14 @@ irhomolar_reducing = constants_header.irhomolar_reducing irhomolar_critical = constants_header.irhomolar_critical iT_reducing = constants_header.iT_reducing iT_critical = constants_header.iT_critical +irhomass_reducing = constants_header.irhomass_reducing +irhomass_critical = constants_header.irhomass_critical +iP_critical = constants_header.iP_critical +iT_triple = constants_header.iT_triple +iP_triple = constants_header.iP_triple +iT_min = constants_header.iT_min +iT_max = constants_header.iT_max +iP_max = constants_header.iP_max iT = constants_header.iT iP = constants_header.iP iQ = constants_header.iQ diff --git a/wrappers/Python/CoolProp5/constants_header.pxd b/wrappers/Python/CoolProp5/constants_header.pxd index f095295e..9242ba8b 100644 --- a/wrappers/Python/CoolProp5/constants_header.pxd +++ b/wrappers/Python/CoolProp5/constants_header.pxd @@ -8,6 +8,14 @@ cdef extern from "DataStructures.h" namespace "CoolProp": irhomolar_critical iT_reducing iT_critical + irhomass_reducing + irhomass_critical + iP_critical + iT_triple + iP_triple + iT_min + iT_max + iP_max iT iP iQ