github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-09 13:25:12 -05:00
Code Issues Packages Projects Releases Wiki Activity

Implement new, more general solvers for the polynomials. Added solvers for rho, c, u and s to the incompressibles. T as function of h is still missing, should got to the backend...

Browse Source
This commit is contained in:
jowr
2014-07-10 12:08:09 +02:00
parent e486f3f78e
commit ed35139d8a
5 changed files with 330 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
src/Backends/Incompressible/IncompressibleBackend.cpp
View File

@@ -27,9 +27,9 @@ IncompressibleBackend::IncompressibleBackend(const std::vector<std::string> &com
}
void IncompressibleBackend::update(long input_pair, double value1, double value2) {
if (mass_fractions.empty()){
throw ValueError("mass fractions have not been set");
}
//if (mass_fractions.empty()){
// throw ValueError("mass fractions have not been set");
//}
switch (input_pair)
{
case PT_INPUTS: {
@@ -65,7 +65,7 @@ void IncompressibleBackend::update(long input_pair, double value1, double value2
@param mole_fractions The vector of mole fractions of the components
*/
void IncompressibleBackend::set_mole_fractions(const std::vector<long double> &mole_fractions) {
throw CoolProp::NotImplementedError("Cannot set mole fractions for incompressible fluid");
throw NotImplementedError("Cannot set mole fractions for incompressible fluid");
}
/// Set the mass fractions
@@ -73,12 +73,54 @@ void IncompressibleBackend::set_mole_fractions(const std::vector<long double> &m
@param mass_fractions The vector of mass fractions of the components
*/
void IncompressibleBackend::set_mass_fractions(const std::vector<long double> &mass_fractions) {
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;
}
/// Check if the mole fractions have been set, etc.
void IncompressibleBackend::check_status() {
throw CoolProp::NotImplementedError("Cannot check status for incompressible fluid");
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);
//
}

83
src/Backends/Incompressible/IncompressibleFluid.cpp
View File

@@ -105,7 +105,7 @@ double IncompressibleFluid::c (double T, double p, double x){
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.",__FILE__,__LINE__,specific_heat.type));
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for specific heat.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
@@ -263,6 +263,87 @@ double IncompressibleFluid::V2M (double T, double y){throw NotImplemen
/// Conversion from mass-based to mole-based composition.
double IncompressibleFluid::M2M (double T, double x){throw NotImplementedError("TODO");}
/* Some functions can be inverted directly, those are listed
* here. It is also possible to solve for other quantities, but
* that involves some more sophisticated processing and is not
* done here, but in the backend, T(h,p) for example.
*/
/// Temperature as a function of density, pressure and composition.
double IncompressibleFluid::T_rho (double Dmass, double p, double x){
double d_raw = Dmass; // No changes needed, no reference values...
switch (density.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limits(density.coeffs, x, d_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse density.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
/// Temperature as a function of heat capacities as a function of temperature, pressure and composition.
double IncompressibleFluid::T_c (double Cmass, double p, double x){
double c_raw = Cmass; // No changes needed, no reference values...
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limits(specific_heat.coeffs, x, c_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse specific heat.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
/// Temperature as a function of entropy as a function of temperature, pressure and composition.
double IncompressibleFluid::T_s (double Smass, double p, double x){
double s_raw = Smass + sref;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limitsInt(specific_heat.coeffs, x, s_raw, Tmin, Tmax, 0, -1, 0, Tbase, xbase, 0);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
/// Temperature as a function of internal energy as a function of temperature, pressure and composition.
double IncompressibleFluid::T_u (double Umass, double p, double x){
double u_raw = Umass + uref;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limitsInt(specific_heat.coeffs, x, u_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase, 0);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
///// Temperature as a function of enthalpy, pressure and composition.
//double IncompressibleFluid::T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));}
///// Viscosity as a function of temperature, pressure and composition.
//double IncompressibleFluid::T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));}
///// Thermal conductivity as a function of temperature, pressure and composition.
//double IncompressibleFluid::T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));}
///// Saturation pressure as a function of temperature and composition.
//double IncompressibleFluid::T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));}
///// Composition as a function of freezing temperature and pressure.
//double IncompressibleFluid::x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));}
/*
* Some more functions to provide a single implementation
* of important routines.

120
src/PolyMath.cpp
View File

@@ -196,6 +196,33 @@ double Polynomial2D::integral(const Eigen::MatrixXd &coefficients, const double
return this->evaluate(this->integrateCoeffs(coefficients, axis, 1), x_in,y_in);
}
/// 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){
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;
}
/// 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){
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;
}
/// @param coefficients vector containing the ordered coefficients
/// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension)
@@ -236,13 +263,7 @@ Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd &coefficients, const d
/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y)
double Polynomial2D::solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis){
Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis);
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;
return solve_limits(res, min, max);
}
/// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension)
@@ -251,13 +272,7 @@ double Polynomial2D::solve_limits(const Eigen::MatrixXd &coefficients, const dou
/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y)
double Polynomial2D::solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis){
Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis);
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;
return solve_guess(res, guess);
}
@@ -650,7 +665,7 @@ double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const dou
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @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
Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base = 0.0, const double &y_base = 0.0){
Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){
Eigen::MatrixXd newCoefficients;
Eigen::VectorXd tmpCoefficients;
@@ -715,16 +730,10 @@ Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, con
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @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 = 0.0, const double &y_base = 0.0){
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){
Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base);
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;
}
return Polynomial2D::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
/// @param coefficients vector containing the ordered coefficients
@@ -736,17 +745,43 @@ double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd &coefficients, const
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @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 = 0.0, const double &y_base = 0.0){
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){
Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base);
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;
}
return Polynomial2D::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
/// @param coefficients vector containing the ordered coefficients
/// @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 min double value that represents the minimum value
/// @param max double value that represents the maximum value
/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y)
/// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @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
/// @param int_axis axis for the integration (0=x, 1=y)
double Polynomial2DFrac::solve_limitsInt(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, const int &int_axis){
Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis);
return Polynomial2D::solve_limits(res, min, max);
} //TODO: Implement tests for this solver
/// Uses the Newton solver to find the roots of Int(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)
/// @param z_in double value that represents the current output in the 3rd dimension
/// @param guess double value that represents the start value
/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y)
/// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @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
/// @param int_axis axis for the integration (0=x, 1=y)
double Polynomial2DFrac::solve_guessInt(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, const int &int_axis){
Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis);
return Polynomial2D::solve_guess(res, guess);
} //TODO: Implement tests for this solver
@@ -822,16 +857,31 @@ Poly2DFracResidual::Poly2DFracResidual(Polynomial2DFrac &poly, const Eigen::Matr
double Poly2DFracResidual::call(double target){
if (axis==iX) return poly.evaluate(coefficients,target,in,x_exp,y_exp,x_base,y_base)-z_in;
if (axis==iY) return poly.evaluate(coefficients,in,target,x_exp,y_exp,x_base,y_base)-z_in;
return -_HUGE;
return _HUGE;
}
double Poly2DFracResidual::deriv(double target){
if (axis==iX) return poly.derivative(coefficients,target,in,axis,x_exp,y_exp,x_base,y_base);
if (axis==iY) return poly.derivative(coefficients,in,target,axis,x_exp,y_exp,x_base,y_base);
return -_HUGE;
return _HUGE;
}
Poly2DFracIntResidual::Poly2DFracIntResidual(Polynomial2DFrac &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, const int &int_axis)
: Poly2DFracResidual(poly, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base){
this->int_axis = int_axis;
}
double Poly2DFracIntResidual::call(double target){
if (axis==iX) return poly.integral(coefficients,target,in,int_axis,x_exp,y_exp,x_base,y_base)-z_in;
if (axis==iY) return poly.integral(coefficients,in,target,int_axis,x_exp,y_exp,x_base,y_base)-z_in;
return _HUGE;
}
double Poly2DFracIntResidual::deriv(double target){
if (axis==iX) return poly.evaluate(coefficients,target,in,x_exp,y_exp,x_base,y_base);
if (axis==iY) return poly.evaluate(coefficients,in,target,x_exp,y_exp,x_base,y_base);
return _HUGE;
}