diff --git a/include/MatrixMath.h b/include/MatrixMath.h index ccb05288..2f63b810 100644 --- a/include/MatrixMath.h +++ b/include/MatrixMath.h @@ -10,7 +10,7 @@ #include #include "float.h" -#include +#include "Eigen/Core" /// A wrapper around std::vector /** This wrapper makes the standard vector multi-dimensional. diff --git a/include/PolyMath.h b/include/PolyMath.h index aa248d8e..0a40e0e9 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -9,12 +9,13 @@ #include #include "Solvers.h" #include "MatrixMath.h" -#include +#include "unsupported/Eigen/Polynomials" namespace CoolProp{ // Just a forward declaration class Poly2DResidual; +class Poly2DFracResidual; /// The base class for all Polynomials /** A clear and straight forward implementation of polynomial operations. Still @@ -103,12 +104,14 @@ public: double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis); protected: + // TODO: Why doe these base definitions not work with derived classes? /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in /// @param res Poly2DResidual object to calculate residuals and derivatives /// @param min double value that represents the minimum value /// @param max double value that represents the maximum value double solve_limits(Poly2DResidual res, const double &min, const double &max); + // TODO: Why doe these base definitions not work with derived classes? /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in /// @param res Poly2DResidual object to calculate residuals and derivatives /// @param guess double value that represents the start value @@ -287,6 +290,21 @@ public: /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base); +protected: + // TODO: Why doe these base definitions not work with derived classes? + /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in + /// @param res Poly2DResidual object to calculate residuals and derivatives + /// @param min double value that represents the minimum value + /// @param max double value that represents the maximum value + double solve_limits(Poly2DFracResidual res, const double &min, const double &max); + + // TODO: Why doe these base definitions not work with derived classes? + /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in + /// @param res Poly2DResidual object to calculate residuals and derivatives + /// @param guess double value that represents the start value + double solve_guess(Poly2DFracResidual res, const double &guess); + +public: /// Returns a vector with ALL the real roots of p(x_in,y_in)-z_in /// @param coefficients vector containing the ordered coefficients /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) diff --git a/include/TestObjects.h b/include/TestObjects.h new file mode 100644 index 00000000..f13513a0 --- /dev/null +++ b/include/TestObjects.h @@ -0,0 +1,20 @@ +/** + * This file contains some basic methods to generate + * objects that can be used in the test routines. + * This makes the tests themselves much more readable + * and assures that the objects used for testing are the + * same in all places. + */ +#include "IncompressibleFluid.h" +#include "Eigen/Core" +#include "MatrixMath.h" + +#if defined ENABLE_CATCH +namespace CoolPropTesting { + +Eigen::MatrixXd makeMatrix(const std::vector &coefficients); +CoolProp::IncompressibleFluid incompressibleFluidObject(); +//IncompressibleBackend incompressibleBackendObject(); + +} // namespace CoolPropTesting +#endif // ENABLE_CATCH diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp index 29402e7e..77396caa 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.cpp +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -18,6 +18,10 @@ namespace CoolProp { +IncompressibleBackend::IncompressibleBackend(IncompressibleFluid* fluid) { + this->fluid = fluid; +} + IncompressibleBackend::IncompressibleBackend(const std::string &fluid_name) { fluid = &get_incompressible_fluid(fluid_name); } @@ -76,51 +80,342 @@ void IncompressibleBackend::set_mass_fractions(const std::vector &m if (mass_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the mass fraction vector and not %d.",mass_fractions.size())); this->mass_fractions = mass_fractions; } +/// Set the mass fractions +/** +@param mass_fraction The mass fraction of the component other than water +*/ +void IncompressibleBackend::set_mass_fractions(const long double &mass_fraction) { + this->mass_fractions.clear(); + this->mass_fractions.push_back(mass_fraction); +} /// Check if the mole fractions have been set, etc. void IncompressibleBackend::check_status() { throw NotImplementedError("Cannot check status for incompressible fluid"); } -///// Calculate T given pressure and density -///** -//@param rhomass The mass density in kg/m^3 -//@param p The pressure in Pa -//@returns T The temperature in K -//*/ -//long double IncompressibleBackend::DmassP_flash(long double rhomass, long double p){ -// -//} -///// Calculate T given pressure and enthalpy -///** -//@param hmass The mass enthalpy in J/kg -//@param p The pressure in Pa -//@returns T The temperature in K -//*/ -//long double IncompressibleBackend::HmassP_flash(long double hmass, long double p); -///// Calculate T given pressure and entropy -///** -//@param smass The mass entropy in J/kg/K -//@param p The pressure in Pa -//@returns T The temperature in K -//*/ -//long double IncompressibleBackend::PSmass_flash(long double p, long double smass); -// -///// Calculate T given pressure and internal energy -///** -//@param umass The mass internal energy in J/kg -//@param p The pressure in Pa -//@returns T The temperature in K -//*/ -//long double IncompressibleBackend::PUmass_flash(long double p, long double umass); -// - - +/// Calculate T given pressure and density +/** +@param rhomass The mass density in kg/m^3 +@param p The pressure in Pa +@returns T The temperature in K +*/ +long double IncompressibleBackend::DmassP_flash(long double rhomass, long double p){ + return fluid->T_rho(rhomass, p, mass_fractions[0]); +} +/// Calculate T given pressure and enthalpy +/** +@param hmass The mass enthalpy in J/kg +@param p The pressure in Pa +@returns T The temperature in K +*/ +long double IncompressibleBackend::HmassP_flash(long double hmass, long double p){ + class HmassP_residual : public FuncWrapper1D { + protected: + double p,x,h_in; + IncompressibleFluid* fluid; + protected: + HmassP_residual(); + public: + HmassP_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &h_in){ + this->p = p; + this->x = x; + this->h_in = h_in; + this->fluid = fluid; + } + virtual ~HmassP_residual(){}; + double call(double target){ + return fluid->h(target,p,x) - h_in; //fluid.u(target,p,x)+ p / fluid.rho(target,p,x) - h_in; + } + //double deriv(double target); + }; + //double T_tmp = this->PUmass_flash(p, hmass); // guess value from u=h + HmassP_residual res = HmassP_residual(fluid, p, mass_fractions[0], hmass); + std::string errstring; + double macheps = DBL_EPSILON; + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Brent(res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring); + //if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; + return result; +} +/// Calculate T given pressure and entropy +/** +@param smass The mass entropy in J/kg/K +@param p The pressure in Pa +@returns T The temperature in K +*/ +long double IncompressibleBackend::PSmass_flash(long double p, long double smass){ + return fluid->T_s(smass, p, mass_fractions[0]); +} +/// Calculate T given pressure and internal energy +/** +@param umass The mass internal energy in J/kg +@param p The pressure in Pa +@returns T The temperature in K +*/ +long double IncompressibleBackend::PUmass_flash(long double p, long double umass){ + return fluid->T_u(umass, p, mass_fractions[0]); +} } + + +// Testing routines with fixed parameters and known results +/* These functions try to cover as much as possible, but + * they still need some serious additions. + */ + +#ifdef ENABLE_CATCH +#include +#include +#include "catch.hpp" + +#include "TestObjects.h" + +TEST_CASE("Internal consistency checks and example use cases for the incompressible backend","[IncompressibleBackend]") +{ + CoolProp::IncompressibleFluid fluid = CoolPropTesting::incompressibleFluidObject(); + CoolProp::IncompressibleBackend backend = CoolProp::IncompressibleBackend(&fluid); + + SECTION("Test case for Methanol from SecCool") { + + // Some basic functions + // has to return false + CHECK( backend.using_mole_fractions()==false ); + + //void update(long input_pair, double value1, double value2); + + std::vector fractions; + fractions.push_back(0.4); + CHECK_THROWS( backend.set_mole_fractions(fractions) ); + CHECK_NOTHROW( backend.set_mass_fractions(fractions) ); + fractions.push_back(0.4); + CHECK_THROWS( backend.set_mass_fractions(fractions) ); + CHECK_NOTHROW( backend.set_mass_fractions(0.4) ); + CHECK_THROWS( backend.check_status() ); + + + + // Prepare the results and compare them to the calculated values + double acc = 0.0001; + double T = 273.15+10; + double p = 10e5; + double x = 0.25; + backend.set_mass_fractions(x); + double val = 0; + double res = 0; + + //CoolProp::set_debug_level(100); + + // Compare density flash + val = fluid.rho(T,p,x); + //res = backend.DmassP_flash(val, p); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(val); + CAPTURE(res); + CHECK( check_abs(T,backend.DmassP_flash(val, p),acc) ); + } + +// +// +// +// +// /// Calculate T given pressure and density +// /** +// @param rhomass The mass density in kg/m^3 +// @param p The pressure in Pa +// @returns T The temperature in K +// */ +// long double DmassP_flash(long double rhomass, long double p); +// /// Calculate T given pressure and enthalpy +// /** +// @param hmass The mass enthalpy in J/kg +// @param p The pressure in Pa +// @returns T The temperature in K +// */ +// long double HmassP_flash(long double hmass, long double p); +// /// Calculate T given pressure and entropy +// /** +// @param smass The mass entropy in J/kg/K +// @param p The pressure in Pa +// @returns T The temperature in K +// */ +// long double PSmass_flash(long double p, long double smass); +// +// /// Calculate T given pressure and internal energy +// /** +// @param umass The mass internal energy in J/kg +// @param p The pressure in Pa +// @returns T The temperature in K +// */ +// long double PUmass_flash(long double p, long double umass); +// +// +// +//// /// Get the viscosity [Pa-s] +//// long double calc_viscosity(void){return fluid->visc(_T, _p, mass_fractions[0]);}; +//// /// Get the thermal conductivity [W/m/K] (based on the temperature and pressure in the state class) +//// long double calc_conductivity(void){return fluid->cond(_T, _p, mass_fractions[0]);}; +//// +//// long double calc_rhomass(void){return fluid->rho(_T, _p, mass_fractions[0]);}; +//// long double calc_hmass(void){return fluid->h(_T, _p, mass_fractions[0]);}; +//// long double calc_smass(void){return fluid->s(_T, _p, mass_fractions[0]);}; +//// long double calc_umass(void){return fluid->u(_T, _p, mass_fractions[0]);}; +//// long double calc_cpmass(void){return fluid->cp(_T, _p, mass_fractions[0]);}; +//// long double calc_cvmass(void){return fluid->cv(_T, _p, mass_fractions[0]);}; +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// +// // Compare cp +// val = 3993.9748117022423; +// res = CH3OH.c(T,p,x); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// // Compare s +// val = -206.62646783739274; +// res = CH3OH.s(T,p,x); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// val = 0.0; +// res = CH3OH.s(Tref,pref,xref); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( val==res ); +// } +// +// // Compare u +// val = -60043.78429641827; +// res = CH3OH.u(T,p,x); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// val = href - pref/CH3OH.rho(Tref,pref,xref); +// res = CH3OH.u(Tref,pref,xref); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( val==res ); +// } +// +// // Compare h +// val = -59005.67386390795; +// res = CH3OH.h(T,p,x); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// val = 0.0; +// res = CH3OH.h(Tref,pref,xref); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( val==res ); +// } +// +// // Compare v +// val = 0.0023970245009602097; +// res = CH3OH.visc(T,p,x)/1e3; +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// // Compare l +// val = 0.44791148414693727; +// res = CH3OH.cond(T,p,x); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// // Compare Tfreeze +// val = -20.02+273.15;// 253.1293105454671; +// res = CH3OH.Tfreeze(p,x)+273.15; +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } + + + } + + +} + +#endif /* ENABLE_CATCH */ diff --git a/src/Backends/Incompressible/IncompressibleBackend.h b/src/Backends/Incompressible/IncompressibleBackend.h index 5ca99d84..36308981 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.h +++ b/src/Backends/Incompressible/IncompressibleBackend.h @@ -21,6 +21,9 @@ public: IncompressibleBackend(){}; virtual ~IncompressibleBackend(){}; + /// The instantiator + /// @param fluid object, mostly for testing purposes + IncompressibleBackend(IncompressibleFluid* fluid); /// The instantiator /// @param fluid_name the string with the fluid name IncompressibleBackend(const std::string &fluid_name); @@ -54,6 +57,12 @@ public: */ void set_mass_fractions(const std::vector &mass_fractions); + /// Set the mass fraction + /** + @param mass_fractions The mass fraction of the component other than water + */ + void set_mass_fractions(const long double &mass_fraction); + /// Check if the mole fractions have been set, etc. void check_status(); diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 3c6e1fa2..769b53b2 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -14,7 +14,7 @@ #include "Solvers.h" -#include +#include "unsupported/Eigen/Polynomials" namespace CoolProp{ @@ -161,7 +161,7 @@ double Polynomial2D::evaluate(const Eigen::MatrixXd &coefficients, const double throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",__FILE__,__LINE__,coefficients.rows(),coefficients.cols())); } double result = Eigen::poly_eval( Eigen::RowVectorXd(coefficients), x_in ); - if (this->do_debug()) std::cout << "Running evaluate(" << mat_to_string(coefficients) << ", " << vec_to_string(x_in) << "): " << result << std::endl; + if (this->do_debug()) std::cout << "Running 1D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in) << "): " << result << std::endl; return result; } /// @param coefficients vector containing the ordered coefficients @@ -174,7 +174,7 @@ double Polynomial2D::evaluate(const Eigen::MatrixXd &coefficients, const double result *= x_in; result += evaluate(coefficients.row(i), y_in); } - if (this->do_debug()) std::cout << "Running evaluate(" << mat_to_string(coefficients) << ", " << vec_to_string(x_in) << ", " << vec_to_string(y_in) << "): " << result << std::endl; + if (this->do_debug()) std::cout << "Running 2D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in) << ", y_in:" << vec_to_string(y_in) << "): " << result << std::endl; return result; } @@ -196,11 +196,13 @@ double Polynomial2D::integral(const Eigen::MatrixXd &coefficients, const double return this->evaluate(this->integrateCoeffs(coefficients, axis, 1), x_in,y_in); } +// TODO: Why doe these base definitions not work with derived classes? /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in /// @param res Poly2DResidual object to calculate residuals and derivatives /// @param min double value that represents the minimum value /// @param max double value that represents the maximum value double Polynomial2D::solve_limits(Poly2DResidual res, const double &min, const double &max){ + if (do_debug()) std::cout << format("Called solve_limits with: min=%f and max=%f", min, max) << std::endl; std::string errstring; double macheps = DBL_EPSILON; double tol = DBL_EPSILON*1e3; @@ -210,10 +212,12 @@ double Polynomial2D::solve_limits(Poly2DResidual res, const double &min, const d return result; } +// TODO: Why doe these base definitions not work with derived classes? /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in /// @param res Poly2DResidual object to calculate residuals and derivatives /// @param guess double value that represents the start value double Polynomial2D::solve_guess(Poly2DResidual res, const double &guess){ + if (do_debug()) std::cout << format("Called solve_guess with: guess=%f ", guess) << std::endl; std::string errstring; //set_debug_level(1000); double tol = DBL_EPSILON*1e3; @@ -531,6 +535,8 @@ double Polynomial2DFrac::evaluate(const Eigen::MatrixXd &coefficients, const dou posExp *= x_in-x_base; posExp += evaluate(tmpCoeffs.row(i), y_in, y_exp, y_base); } + if (this->do_debug()) std::cout << "Running 2D evaluate(" << mat_to_string(coefficients) << ", " << std::endl; + if (this->do_debug()) std::cout << "x_in:" << vec_to_string(x_in) << ", y_in:" << vec_to_string(y_in) << ", x_exp:" << vec_to_string(x_exp) << ", y_exp:" << vec_to_string(y_exp) << ", x_base:" << vec_to_string(x_base) << ", y_base:" << vec_to_string(y_base) << "): " << negExp+posExp << std::endl; return negExp+posExp; } @@ -655,6 +661,37 @@ double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const dou } +// TODO: Why doe these base definitions not work with derived classes? +/// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in +/// @param res Poly2DResidual object to calculate residuals and derivatives +/// @param min double value that represents the minimum value +/// @param max double value that represents the maximum value +double Polynomial2DFrac::solve_limits(Poly2DFracResidual res, const double &min, const double &max){ + if (do_debug()) std::cout << format("Called solve_limits with: min=%f and max=%f", min, max) << std::endl; + std::string errstring; + double macheps = DBL_EPSILON; + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Brent(res, min, max, macheps, tol, maxiter, errstring); + if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; + return result; +} + +// TODO: Why doe these base definitions not work with derived classes? +/// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in +/// @param res Poly2DResidual object to calculate residuals and derivatives +/// @param guess double value that represents the start value +double Polynomial2DFrac::solve_guess(Poly2DFracResidual res, const double &guess){ + if (do_debug()) std::cout << format("Called solve_guess with: guess=%f ", guess) << std::endl; + std::string errstring; + //set_debug_level(1000); + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Newton(res, guess, tol, maxiter, errstring); + if (this->do_debug()) std::cout << "Newton solver message: " << errstring << std::endl; + return result; +} + /// Returns a vector with ALL the real roots of p(x_in,y_in)-z_in /// @param coefficients vector containing the ordered coefficients @@ -731,8 +768,9 @@ Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, con /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){ + if (do_debug()) std::cout << format("Called solve_limits with: %f, %f, %f, %f, %d, %d, %d, %f, %f",in, z_in, min, max, axis, x_exp, y_exp, x_base, y_base) << std::endl; Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base); - return Polynomial2D::solve_limits(res, min, max); + return this->solve_limits(res, min, max); } //TODO: Implement tests for this solver /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in @@ -746,8 +784,9 @@ double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd &coefficients, const /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension double Polynomial2DFrac::solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){ + if (do_debug()) std::cout << format("Called solve_guess with: %f, %f, %f, %d, %d, %d, %f, %f",in, z_in, guess, axis, x_exp, y_exp, x_base, y_base) << std::endl; Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base); - return Polynomial2D::solve_guess(res, guess); + return this->solve_guess(res, guess); } //TODO: Implement tests for this solver /// Uses the Brent solver to find the roots of Int(p(x_in,y_in))-z_in diff --git a/src/Tests/TestObjects.cpp b/src/Tests/TestObjects.cpp new file mode 100644 index 00000000..050561e7 --- /dev/null +++ b/src/Tests/TestObjects.cpp @@ -0,0 +1,215 @@ +/** + * This file contains some basic methods to generate + * objects that can be used in the test routines. + * This makes the tests themselves much more readable + * and assures that the objects used for testing are the + * same in all places. + */ +#include "TestObjects.h" +#include "IncompressibleFluid.h" +#include "Eigen/Core" + +#if defined ENABLE_CATCH + +Eigen::MatrixXd CoolPropTesting::makeMatrix(const std::vector &coefficients){ + //IncompressibleClass::checkCoefficients(coefficients,18); + std::vector< std::vector > matrix; + std::vector tmpVector; + + tmpVector.clear(); + tmpVector.push_back(coefficients[0]); + tmpVector.push_back(coefficients[6]); + tmpVector.push_back(coefficients[11]); + tmpVector.push_back(coefficients[15]); + matrix.push_back(tmpVector); + + tmpVector.clear(); + tmpVector.push_back(coefficients[1]*100.0); + tmpVector.push_back(coefficients[7]*100.0); + tmpVector.push_back(coefficients[12]*100.0); + tmpVector.push_back(coefficients[16]*100.0); + matrix.push_back(tmpVector); + + tmpVector.clear(); + tmpVector.push_back(coefficients[2]*100.0*100.0); + tmpVector.push_back(coefficients[8]*100.0*100.0); + tmpVector.push_back(coefficients[13]*100.0*100.0); + tmpVector.push_back(coefficients[17]*100.0*100.0); + matrix.push_back(tmpVector); + + tmpVector.clear(); + tmpVector.push_back(coefficients[3]*100.0*100.0*100.0); + tmpVector.push_back(coefficients[9]*100.0*100.0*100.0); + tmpVector.push_back(coefficients[14]*100.0*100.0*100.0); + tmpVector.push_back(0.0); + matrix.push_back(tmpVector); + + tmpVector.clear(); + tmpVector.push_back(coefficients[4]*100.0*100.0*100.0*100.0); + tmpVector.push_back(coefficients[10]*100.0*100.0*100.0*100.0); + tmpVector.push_back(0.0); + tmpVector.push_back(0.0); + matrix.push_back(tmpVector); + + tmpVector.clear(); + tmpVector.push_back(coefficients[5]*100.0*100.0*100.0*100.0*100.0); + tmpVector.push_back(0.0); + tmpVector.push_back(0.0); + tmpVector.push_back(0.0); + matrix.push_back(tmpVector); + + tmpVector.clear(); + return CoolProp::vec_to_eigen(matrix).transpose(); +} + +CoolProp::IncompressibleFluid CoolPropTesting::incompressibleFluidObject(){ + bool PRINT = false; + std::string tmpStr; + std::vector tmpVector; + std::vector< std::vector > tmpMatrix; + + tmpVector.clear(); + tmpVector.push_back( 960.24665800); + tmpVector.push_back(-1.2903839100); + tmpVector.push_back(-0.0161042520); + tmpVector.push_back(-0.0001969888); + tmpVector.push_back( 1.131559E-05); + tmpVector.push_back( 9.181999E-08); + tmpVector.push_back(-0.4020348270); + tmpVector.push_back(-0.0162463989); + tmpVector.push_back( 0.0001623301); + tmpVector.push_back( 4.367343E-06); + tmpVector.push_back( 1.199000E-08); + tmpVector.push_back(-0.0025204776); + tmpVector.push_back( 0.0001101514); + tmpVector.push_back(-2.320217E-07); + tmpVector.push_back( 7.794999E-08); + tmpVector.push_back( 9.937483E-06); + tmpVector.push_back(-1.346886E-06); + tmpVector.push_back( 4.141999E-08); + CoolProp::IncompressibleData density; + density.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL; + density.coeffs = makeMatrix(tmpVector); + + tmpVector.clear(); + tmpVector.push_back( 3822.9712300); + tmpVector.push_back(-23.122409500); + tmpVector.push_back( 0.0678775826); + tmpVector.push_back( 0.0022413893); + tmpVector.push_back(-0.0003045332); + tmpVector.push_back(-4.758000E-06); + tmpVector.push_back( 2.3501449500); + tmpVector.push_back( 0.1788839410); + tmpVector.push_back( 0.0006828000); + tmpVector.push_back( 0.0002101166); + tmpVector.push_back(-9.812000E-06); + tmpVector.push_back(-0.0004724176); + tmpVector.push_back(-0.0003317949); + tmpVector.push_back( 0.0001002032); + tmpVector.push_back(-5.306000E-06); + tmpVector.push_back( 4.242194E-05); + tmpVector.push_back( 2.347190E-05); + tmpVector.push_back(-1.894000E-06); + CoolProp::IncompressibleData specific_heat; + specific_heat.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL; + specific_heat.coeffs = makeMatrix(tmpVector); + + tmpVector.clear(); + tmpVector.push_back( 0.4082066700); + tmpVector.push_back(-0.0039816870); + tmpVector.push_back( 1.583368E-05); + tmpVector.push_back(-3.552049E-07); + tmpVector.push_back(-9.884176E-10); + tmpVector.push_back( 4.460000E-10); + tmpVector.push_back( 0.0006629321); + tmpVector.push_back(-2.686475E-05); + tmpVector.push_back( 9.039150E-07); + tmpVector.push_back(-2.128257E-08); + tmpVector.push_back(-5.562000E-10); + tmpVector.push_back( 3.685975E-07); + tmpVector.push_back( 7.188416E-08); + tmpVector.push_back(-1.041773E-08); + tmpVector.push_back( 2.278001E-10); + tmpVector.push_back( 4.703395E-08); + tmpVector.push_back( 7.612361E-11); + tmpVector.push_back(-2.734000E-10); + CoolProp::IncompressibleData conductivity; + conductivity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL; + conductivity.coeffs = makeMatrix(tmpVector); + + tmpVector.clear(); + tmpVector.push_back( 1.4725525500); + tmpVector.push_back( 0.0022218998); + tmpVector.push_back(-0.0004406139); + tmpVector.push_back( 6.047984E-06); + tmpVector.push_back(-1.954730E-07); + tmpVector.push_back(-2.372000E-09); + tmpVector.push_back(-0.0411841566); + tmpVector.push_back( 0.0001784479); + tmpVector.push_back(-3.564413E-06); + tmpVector.push_back( 4.064671E-08); + tmpVector.push_back( 1.915000E-08); + tmpVector.push_back( 0.0002572862); + tmpVector.push_back(-9.226343E-07); + tmpVector.push_back(-2.178577E-08); + tmpVector.push_back(-9.529999E-10); + tmpVector.push_back(-1.699844E-06); + tmpVector.push_back(-1.023552E-07); + tmpVector.push_back( 4.482000E-09); + CoolProp::IncompressibleData viscosity; + viscosity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL; + viscosity.coeffs = makeMatrix(tmpVector); + + + tmpVector.clear(); + tmpVector.push_back( 27.755555600/100.0); // reference concentration in per cent + tmpVector.push_back(-22.973221700); + tmpVector.push_back(-1.1040507200*100.0); + tmpVector.push_back(-0.0120762281*100.0*100.0); + tmpVector.push_back(-9.343458E-05*100.0*100.0*100.0); + CoolProp::IncompressibleData T_freeze; + T_freeze.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYOFFSET; + T_freeze.coeffs = CoolProp::vec_to_eigen(tmpVector); + + // After preparing the coefficients, we have to create the objects + CoolProp::IncompressibleFluid CH3OH; + CH3OH.setName("CH3OH"); + CH3OH.setDescription("Methanol solution"); + CH3OH.setReference("SecCool software"); + CH3OH.setTmax( 20 + 273.15); + CH3OH.setTmin(-50 + 273.15); + CH3OH.setxmax(0.5); + CH3OH.setxmin(0.0); + CH3OH.setTminPsat( 20 + 273.15); + + CH3OH.setTbase(-4.48 + 273.15); + CH3OH.setxbase(31.57 / 100.0); + + /// Setters for the coefficients + CH3OH.setDensity(density); + CH3OH.setSpecificHeat(specific_heat); + CH3OH.setViscosity(viscosity); + CH3OH.setConductivity(conductivity); + //CH3OH.setPsat(saturation_pressure); + CH3OH.setTfreeze(T_freeze); + //CH3OH.setVolToMass(volume2mass); + //CH3OH.setMassToMole(mass2mole); + + //XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0); + double Tref = 25+273.15; + double pref = 0.0; + double xref = 0.25; + double href = 0.0; + double sref = 0.0; + CH3OH.set_reference_state(Tref, pref, xref, href, sref); + + /// A function to check coefficients and equation types. + CH3OH.validate(); + + return CH3OH; +} +//CoolProp::IncompressibleBackend CoolProp::Testing::incompressibleBackendObject(){ +// return CoolProp::IncompressibleBackend(CoolProp::Testing::incompressibleFluidObject()); +//} + +#endif // ENABLE_CATCH