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

Finally, some luck with the Matrix classes... C++ lessons learned.

Browse Source
This commit is contained in:
jowr
2014-06-06 18:39:33 +02:00
parent c2b06e301c
commit f4be09df63
5 changed files with 187 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
dev/TTSE/TTSE_ranges.py
View File

@@ -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

13
include/MatrixMath.h
View File

@@ -34,10 +34,10 @@ namespace CoolProp{
* it is still needed.
*/
/// @param coefficients matrix containing the ordered coefficients
template<class T, int R, int C> void convert(const Eigen::Matrix<T,R,C> &coefficients, std::vector<std::vector<T> > &result){
template<class T, typename Derived> void convert(const Eigen::MatrixBase<Derived> &coefficients, std::vector<std::vector<T> > &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<std::vector<T> > result(R, std::vector<T>(C, 0));
result.resize(R, std::vector<T>(C, 0)); // extends vector if necessary
for (size_t i = 0; i < R; ++i) {
@@ -48,8 +48,9 @@ template<class T, int R, int C> void convert(const Eigen::Matrix<T,R,C> &coeffic
}
}
/// @param coefficients matrix containing the ordered coefficients
template <class T, int R, int C> void convert(const std::vector<std::vector<T> > &coefficients, Eigen::Matrix<T,R,C> &result){
template <class T, typename Derived> void convert(const std::vector<std::vector<T> > &coefficients, Eigen::MatrixBase<Derived> &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 <class T, int R, int C> void convert(const std::vector<std::vector<T> >
}
/// @param coefficients matrix containing the ordered coefficients
template<class T, int R> void convert(const Eigen::Matrix<T,R,1> &coefficients, std::vector<T> &result){
template<class T, typename Derived> void convert(const Eigen::MatrixBase<Derived> &coefficients, std::vector<T> &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 <class T, int R> void convert(const std::vector<T> &coefficients, Eigen::Matrix<T,R,1> &result){
template <class T, typename Derived> void convert(const std::vector<T> &coefficients, Eigen::MatrixBase<Derived> &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];

6
include/PolyMath.h
View File

@@ -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<std::vector<double> > &coefficients);
void setCoefficients(const Eigen::MatrixXd &coefficients);
void setCoefficients(const std::vector<std::vector<double> > &coefficients);
/// Basic checks for coefficient vectors.
/** Starts with only the first coefficient dimension

45
src/MatrixMath.cpp
View File

@@ -20,6 +20,7 @@ namespace CoolProp{
#ifdef ENABLE_CATCH
#include <math.h>
#include <iostream>
#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<std::vector<double> > 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<double,2,1> matrix;
CoolProp::convert(cHeat, matrix);
CHECK_THROWS( CoolProp::convert(cHeat, matrix) );
}{
Eigen::Matrix<double,5,1> matrix;
CHECK_THROWS( CoolProp::convert(cHeat, matrix) );
}{
Eigen::Matrix<double,4,1> 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<double,2,2> matrix;
CHECK_THROWS( CoolProp::convert(cHeat2D, matrix) );
}{
Eigen::Matrix<double,2,5> matrix;
CHECK_THROWS( CoolProp::convert(cHeat2D, matrix) );
}{
Eigen::Matrix<double,2,4> matrix;
CHECK_NOTHROW( CoolProp::convert(cHeat2D, matrix) );
}
}
}
//
//int main() {
//
// std::vector<std::string> tags;
// tags.push_back("[MatrixMath]");
// run_user_defined_tests(tags);
// //run_tests();
//}
#endif /* ENABLE_CATCH */

129
src/PolyMath.cpp
View File

@@ -10,6 +10,7 @@
//#include <sstream>
//#include <numeric>
#include <math.h>
#include <iostream>
#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<std::vector<double> > &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<std::vector<double> > &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<std::vector<double> > Polynomial2D::integrateCoeffs(const std::vector<std::vector<double> > &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<std::vector<double> > Polynomial2D::deriveCoeffs(const std::vector<std::vector<double> > &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 <math.h>
#include <iostream>
#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<double> 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<std::vector<double> > 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<double> > &coefficients){
// return this->setCoefficients(convert(coefficients));
//}
}