From f4be09df63cc0b3248022524d0d75fb0d7bdb04e Mon Sep 17 00:00:00 2001 From: jowr Date: Fri, 6 Jun 2014 18:39:33 +0200 Subject: [PATCH] Finally, some luck with the Matrix classes... C++ lessons learned. --- dev/TTSE/TTSE_ranges.py | 21 ++++--- include/MatrixMath.h | 13 ++-- include/PolyMath.h | 6 +- src/MatrixMath.cpp | 45 +++++++++++--- src/PolyMath.cpp | 129 +++++++++++++++++++++++++++++++++++++++- 5 files changed, 187 insertions(+), 27 deletions(-) diff --git a/dev/TTSE/TTSE_ranges.py b/dev/TTSE/TTSE_ranges.py index eb4fbcd1..967aca72 100644 --- a/dev/TTSE/TTSE_ranges.py +++ b/dev/TTSE/TTSE_ranges.py @@ -224,7 +224,7 @@ def getlim(key,dicts,fac=1): # Here starts the real script ######################################################################### -fluid = 'water' +fluid = 'R245fa' CP.enable_TTSE_LUT(fluid) PRINT = False @@ -241,7 +241,7 @@ else: Tmin = CP.PropsSI('Tmin' ,'T',0,'P',0,fluid) + 5.0 Tcri = CP.PropsSI('Tcrit','T',0,'P',0,fluid) -Tmax = CP.PropsSI('Tcrit','T',0,'P',0,fluid) * 2.0 +#Tmax = CP.PropsSI('Tcrit','T',0,'P',0,fluid) * 2.0 ## Get phase boundary T_TP = np.linspace(Tmin, Tcri-0.5, 0.5*points) @@ -275,8 +275,12 @@ Hmax = H_L + (H_V-H_L)*2.0 Pmin = np.min([P_L,P_V]) Pmax = CP.PropsSI('pcrit','T',0,'P',0,fluid) * 2.0 # should be p_reduce -Dmin = CP.PropsSI('D','H',Hmin,'P',Pmax,fluid) -Dmax = CP.PropsSI('D','H',Hmax,'P',Pmin,fluid) +Tmin = CP.PropsSI('T','H',Hmin,'P',Pmax,fluid) +Tmax = CP.PropsSI('T','H',Hmax,'P',Pmin,fluid) + +Dmin = CP.PropsSI('D','H',Hmax,'P',Pmin,fluid) +Dmax = CP.PropsSI('D','H',Hmin,'P',Pmax,fluid) + @@ -326,7 +330,6 @@ Tdata = fill_Z(X,Y) HPSdict = {'H': X, 'P': Y, 'S': Z, 'V': logVdata, 'T': Tdata, 'H_TP': X_TP, 'P_TP': Y_TP, 'S_TP': Z_TP} - ######################################################################### # Start with the next diagram, vTp ######################################################################### @@ -338,10 +341,10 @@ Z_TP = logP_TP #Xmin = np.log10(1.0/Dmax) #Xmax = np.log10(1.0/Dmin) -Xmin = np.min(X_TP) -Xmax = np.max(X_TP) -Ymin = np.min(Y_TP) -Ymax = Tmax +Xmin = np.min(logVdata) +Xmax = np.max(logVdata) +Ymin = np.min(Tdata) +Ymax = np.max(Tdata) X = np.linspace(Xmin, Xmax, points) # Volume Y = np.linspace(Ymin, Ymax, points) # Temperature diff --git a/include/MatrixMath.h b/include/MatrixMath.h index a9622d37..a52ce5db 100644 --- a/include/MatrixMath.h +++ b/include/MatrixMath.h @@ -34,10 +34,10 @@ namespace CoolProp{ * it is still needed. */ /// @param coefficients matrix containing the ordered coefficients -template void convert(const Eigen::Matrix &coefficients, std::vector > &result){ +template void convert(const Eigen::MatrixBase &coefficients, std::vector > &result){ // Eigen uses columns as major axis, this might be faster than the row iteration. // However, the 2D vector stores things differently, no idea what is faster... - //size_t nRows = coefficients.rows(), nCols = coefficients.cols(); + size_t R = coefficients.rows(), C = coefficients.cols(); //std::vector > result(R, std::vector(C, 0)); result.resize(R, std::vector(C, 0)); // extends vector if necessary for (size_t i = 0; i < R; ++i) { @@ -48,8 +48,9 @@ template void convert(const Eigen::Matrix &coeffic } } /// @param coefficients matrix containing the ordered coefficients -template void convert(const std::vector > &coefficients, Eigen::Matrix &result){ +template void convert(const std::vector > &coefficients, Eigen::MatrixBase &result){ size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients); + size_t R = result.rows(), C = result.cols(); if (nRows!=R) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",nRows,R)); if (nCols!=C) throw ValueError(format("You have to provide matrices with the same number of columns: %d is not %d. ",nCols,C)); for (size_t i = 0; i < nCols; ++i) { @@ -60,15 +61,17 @@ template void convert(const std::vector > } /// @param coefficients matrix containing the ordered coefficients -template void convert(const Eigen::Matrix &coefficients, std::vector &result){ +template void convert(const Eigen::MatrixBase &coefficients, std::vector &result){ + size_t R = result.rows(); result.resize(R, 0); // extends vector if necessary for (size_t i = 0; i < R; ++i) { result[i] = coefficients(i,0); } } /// @param coefficients matrix containing the ordered coefficients -template void convert(const std::vector &coefficients, Eigen::Matrix &result){ +template void convert(const std::vector &coefficients, Eigen::MatrixBase &result){ size_t nRows = num_rows(coefficients); + size_t R = result.rows(); if (nRows!=R) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",nRows,R)); for (size_t i = 0; i < nRows; ++i) { result(i,0) = coefficients[i]; diff --git a/include/PolyMath.h b/include/PolyMath.h index 06fe4260..1cc21b9f 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -100,7 +100,7 @@ private: public: /// Constructors - Polynomial2D(); + Polynomial2D(){}; Polynomial2D(const Eigen::MatrixXd &coefficients){ this->setCoefficients(coefficients); } @@ -114,8 +114,8 @@ public: public: /// Set the coefficient matrix. /// @param coefficients matrix containing the ordered coefficients - bool setCoefficients(const Eigen::MatrixXd &coefficients); - bool setCoefficients(const std::vector > &coefficients); + void setCoefficients(const Eigen::MatrixXd &coefficients); + void setCoefficients(const std::vector > &coefficients); /// Basic checks for coefficient vectors. /** Starts with only the first coefficient dimension diff --git a/src/MatrixMath.cpp b/src/MatrixMath.cpp index edd87c39..c8d93182 100644 --- a/src/MatrixMath.cpp +++ b/src/MatrixMath.cpp @@ -20,6 +20,7 @@ namespace CoolProp{ #ifdef ENABLE_CATCH #include +#include #include "catch.hpp" TEST_CASE("Internal consistency checks and example use cases for MatrixMath.h","[MatrixMath]") @@ -32,19 +33,45 @@ TEST_CASE("Internal consistency checks and example use cases for MatrixMath.h"," cHeat.push_back(+7.7175381057E-07); cHeat.push_back(-3.7008444051E-20); + std::vector > cHeat2D; + cHeat2D.push_back(cHeat); + cHeat2D.push_back(cHeat); + + SECTION("Pretty printing tests") { + std::cout << std::endl; + CHECK_NOTHROW( CoolProp::vec_to_string(cHeat[0]) ); + CHECK_NOTHROW( CoolProp::vec_to_string(cHeat) ); + CHECK_NOTHROW( CoolProp::vec_to_string(cHeat2D) ); + } + SECTION("Eigen::Vector from std::vector") { + { Eigen::Matrix matrix; - CoolProp::convert(cHeat, matrix); + CHECK_THROWS( CoolProp::convert(cHeat, matrix) ); + }{ + Eigen::Matrix matrix; + CHECK_THROWS( CoolProp::convert(cHeat, matrix) ); + }{ + Eigen::Matrix matrix; + CHECK_NOTHROW( CoolProp::convert(cHeat, matrix) ); + }{ + Eigen::Vector2d matrix; + CHECK_THROWS( CoolProp::convert(cHeat, matrix) ); + }{ + Eigen::Vector4d matrix; + CHECK_NOTHROW( CoolProp::convert(cHeat, matrix) ); + }{ + Eigen::Matrix matrix; + CHECK_THROWS( CoolProp::convert(cHeat2D, matrix) ); + }{ + Eigen::Matrix matrix; + CHECK_THROWS( CoolProp::convert(cHeat2D, matrix) ); + }{ + Eigen::Matrix matrix; + CHECK_NOTHROW( CoolProp::convert(cHeat2D, matrix) ); + } } } -// -//int main() { -// -// std::vector tags; -// tags.push_back("[MatrixMath]"); -// run_user_defined_tests(tags); -// //run_tests(); -//} #endif /* ENABLE_CATCH */ diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 95082e29..926b79b0 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -10,6 +10,7 @@ //#include //#include #include +#include #include "Solvers.h" @@ -17,6 +18,133 @@ namespace CoolProp{ +/// Set the coefficient matrix. +/// @param coefficients matrix containing the ordered coefficients +void Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){ + this->coefficients = coefficients; +} +void Polynomial2D::setCoefficients(const std::vector > &coefficients){ + int c = num_cols(coefficients); + int r = num_rows(coefficients); + Eigen::MatrixXd tmp = Eigen::MatrixXd::Constant(r,c,0.0); + std::cout << "Rows: " << r << " Columns: " << c << std::endl; + std::cout << "Rows: " << tmp.rows() << " Columns: " << tmp.cols() << std::endl; + convert(coefficients,tmp); + this->setCoefficients(tmp); +} + +///// Basic checks for coefficient vectors. +///** Starts with only the first coefficient dimension +// * and checks the matrix size against the parameters rows and columns. */ +///// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension +///// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension +//bool Polynomial2D::checkCoefficients(const unsigned int rows, const unsigned int columns); +///// @param coefficients matrix containing the ordered coefficients +///// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension +///// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension +//bool Polynomial2D::checkCoefficients(const Eigen::MatrixXd &coefficients, const unsigned int rows, const unsigned int columns); +///// @param coefficients matrix containing the ordered coefficients +///// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension +///// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension +//bool Polynomial2D::checkCoefficients(const std::vector > &coefficients, const unsigned int rows, const unsigned int columns); +// +///// Integration functions +///** Integrating coefficients for polynomials is done by dividing the +// * original coefficients by (i+1) and elevating the order by 1 +// * through adding a zero as first coefficient. +// * Some reslicing needs to be applied to integrate along the x-axis. +// * In the brine/solution equations, reordering of the parameters +// * avoids this expensive operation. However, it is included for the +// * sake of completeness. +// */ +///// @param coefficients matrix containing the ordered coefficients +///// @param axis unsigned integer value that represents the desired direction of integration +//Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficients, unsigned int axis = 1); +///// @param coefficients matrix containing the ordered coefficients +///// @param axis unsigned integer value that represents the desired direction of integration +//std::vector > Polynomial2D::integrateCoeffs(const std::vector > &coefficients, unsigned int axis = 1); +// +///// Derivative coefficients calculation +///** Deriving coefficients for polynomials is done by multiplying the +// * original coefficients with i and lowering the order by 1. +// * +// * It is not really deprecated, but untested and therefore a warning +// * is issued. Please check this method before you use it. +// */ +///// @param coefficients matrix containing the ordered coefficients +///// @param axis unsigned integer value that represents the desired direction of derivation +//Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients, unsigned int axis = 1); +///// @param coefficients matrix containing the ordered coefficients +///// @param axis unsigned integer value that represents the desired direction of derivation +//std::vector > Polynomial2D::deriveCoeffs(const std::vector > &coefficients, unsigned int axis = 1); +// +///// The core functions to evaluate the polynomial +///** It is here we implement the different special +// * functions that allow us to specify certain +// * types of polynomials. +// * The derivative might bee needed during the +// * solution process of the solver. It could also +// * be a protected function... +// */ +///// @param x_in double value that represents the current input in the 1st dimension +///// @param y_in double value that represents the current input in the 2nd dimension +//double Polynomial2D::evaluate(const double &x_in, const double &y_in); +///// @param x_in double value that represents the current input in the 1st dimension +///// @param y_in double value that represents the current input in the 2nd dimension +///// @param axis unsigned integer value that represents the axis to solve for (1=x, 2=y) +//double Polynomial2D::derivative(const double &x_in, const double &y_in, unsigned int axis = 1); +///// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) +///// @param z_in double value that represents the current output in the 3rd dimension +///// @param axis unsigned integer value that represents the axis to solve for (1=x, 2=y) +//double Polynomial2D::solve(const double &in, const double &z_in, unsigned int axis = 1); + + +}; + + +#ifdef ENABLE_CATCH +#include +#include +#include "catch.hpp" + +TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","[PolyMath]") +{ + /// Test case for "SylthermXLT" by "Dow Chemicals" + std::vector cHeat; + cHeat.clear(); + cHeat.push_back(+1.1562261074E+03); + cHeat.push_back(+2.0994549103E+00); + cHeat.push_back(+7.7175381057E-07); + cHeat.push_back(-3.7008444051E-20); + + std::vector > cHeat2D; + cHeat2D.push_back(cHeat); + cHeat2D.push_back(cHeat); + + CoolProp::Polynomial2D poly2D; + + SECTION("Coefficient parsing and setting") { + poly2D.setCoefficients(cHeat2D); + } +} + + +#endif /* ENABLE_CATCH */ + + + + + + + + + + + + + + + /// Constructors for the base class for all Polynomials //Polynomial1D::Polynomial1D(); //bool Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){ @@ -26,7 +154,6 @@ namespace CoolProp{ //bool Polynomial2D::setCoefficients(const std::vector< std::vector > &coefficients){ // return this->setCoefficients(convert(coefficients)); //} -}