github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-08 21:05:14 -05:00
Code Issues Packages Projects Releases Wiki Activity

Added the first incompressible test cases and it seems to work...

Browse Source
This commit is contained in:
jowr
2014-06-23 17:34:37 +02:00
parent 864df6d590
commit de03057f48
7 changed files with 920 additions and 516 deletions
Download Patch File Download Diff File Expand all files Collapse all files

120
src/Backends/Helmholtz/Fluids/Incompressible.cpp
View File

@@ -1,120 +0,0 @@
///*
// * Incompressible.cpp
// *
// * Created on: 20 Dec 2013
// * Author: jowr
// */
//
//#include "Incompressible.h"
//
//
//namespace CoolProp {
//
//
//
//// /* All functions need T and p as input. Might not
//// * be necessary, but gives a clearer structure.
//// */
//// /// Density as a function of temperature, pressure and composition.
//// double Incompressible::rho(double T_K, double p) {
//// return poly.polyval(cRho, getxInput(x), getTInput(T_K));
//// }
//// /// Heat capacities as a function of temperature, pressure and composition.
//// double Incompressible::c(double T_K, double p) {
//// return poly.polyval(cHeat, getxInput(x), getTInput(T_K));
//// }
//// /// Enthalpy as a function of temperature, pressure and composition.
//// double Incompressible::h(double T_K, double p) {
//// return h_u(T_K, p);
//// }
//// /// Entropy as a function of temperature, pressure and composition.
//// double Incompressible::s(double T_K, double p) {
//// return poly.polyfracintcentral(cHeat, getxInput(x), T_K, Tbase)
//// - poly.polyfracintcentral(cHeat, getxInput(x), Tref, Tbase);
//// }
//// /// Viscosity as a function of temperature, pressure and composition.
//// double Incompressible::visc(double T_K, double p) {
//// return expo.expval(cVisc, getxInput(x), getTInput(T_K), 2) / 1e3;
//// }
//// /// Thermal conductivity as a function of temperature, pressure and composition.
//// double Incompressible::cond(double T_K, double p) {
//// return poly.polyval(cCond, getxInput(x), getTInput(T_K));
//// }
//// /// Internal energy as a function of temperature, pressure and composition.
//// double Incompressible::u(double T_K, double p) {
//// return poly.polyint(cHeat, getxInput(x), getTInput(T_K))
//// - poly.polyint(cHeat, getxInput(x), getTInput(Tref));
//// }
//
///// Saturation pressure as a function of temperature and composition.
//double Incompressible::psat(double T_K ){throw NotImplementedError("Psat is not available");};
///// Freezing temperature as a function of pressure and composition.
//double Incompressible::Tfreeze( double p){throw NotImplementedError("Tfreeze is not available");};
//
//
///*
// * Some more functions to provide a single implementation
// * of important routines.
// * We start with the check functions that can validate input
// * in terms of pressure p, temperature T and composition x.
// */
//
///// Check validity of temperature input.
///** Compares the given temperature T to the result of a
// * freezing point calculation. This is not necessarily
// * defined for all fluids, default values do not
// * cause errors. */
//bool Incompressible::checkT(double T_K, double p){
// if( Tmin < 0. ) {
// throw ValueError("Please specify the minimum temperature.");
// } else if( Tmax < 0.) {
// throw ValueError("Please specify the maximum temperature.");
// } else if ( (Tmin>T_K) || (T_K>Tmax) ) {
// throw ValueError(format("Your temperature %f is not between %f and %f.",T_K,Tmin,Tmax));
// } else if (T_K < Tfreeze(p)) {
// throw ValueError(format("Your temperature %f is below the freezing point of %f.",T_K,Tfreeze(p)));
// } else {
// return true;
// }
// return false;
//}
//
///// Check validity of pressure input.
///** Compares the given pressure p to the saturation pressure at
// * temperature T and throws and exception if p is lower than
// * the saturation conditions.
// * The default value for psat is -1 yielding true if psat
// * is not redefined in the subclass.
// * */
//bool Incompressible::checkP(double T_K, double p) {
// double ps = psat(T_K);
// if (p<ps) {
// throw ValueError(format("Equations are valid for solution phase only: %f < %f. ",p,ps));
// } else {
// return true;
// }
//}
//
///// Check validity of composition input.
///** Compares the given composition x to a stored minimum and
// * maximum value. Enforces the redefinition of xmin and
// * xmax since the default values cause an error. */
//bool Incompressible::checkX(double x){
// if( xmin < 0. ) {
// throw ValueError("Please specify the minimum concentration.");
// } else if( xmax < 0.) {
// throw ValueError("Please specify the maximum concentration.");
// } else if ( (xmin>x) || (x>xmax) ) {
// throw ValueError(format("Your composition %f is not between %f and %f.",x,xmin,xmax));
// } else {
// return true;
// }
// return false;
//}
//
///// Check validity of temperature, pressure and composition input.
//bool Incompressible::checkTPX(double T, double p, double x) {
// return (checkT(T,p) && checkP(T,p) && checkX(x));
//}
//
//} /* namespace CoolProp */

169
src/Backends/Helmholtz/Fluids/Incompressible.h
View File

@@ -1,169 +0,0 @@
/*
* Incompressible.h
*
* Created on: 20 Dec 2013
* Author: jowr
*/
#ifndef INCOMPRESSIBLE_H_
#define INCOMPRESSIBLE_H_
#include <Eigen/Core>
#include "PolyMath.h"
#include "MatrixMath.h"
namespace CoolProp {
class Incompressible{
protected:
std::string name;
std::string description;
std::string reference;
double Tmin, Tmax;
double xmin, xmax;
double TminPsat;
double xref, Tref;
double xbase, Tbase;
Eigen::MatrixXd cRho;
Eigen::MatrixXd cHeat;
Eigen::MatrixXd cVisc;
Eigen::MatrixXd cCond;
Eigen::MatrixXd cPsat;
Eigen::MatrixXd cTfreeze;
Eigen::MatrixXd cV2M;
Polynomial2DFrac poly;
public:
Incompressible();
virtual ~Incompressible();
std::string getName() const {return name;}
std::string get_name() const {return getName();}// For backwards-compatibility.
std::string getDescription() const {return description;}
std::string getReference() const {return reference;}
double getTmax() const {return Tmax;}
double getTmin() const {return Tmin;}
double getxmax() const {return xmax;}
double getxmin() const {return xmin;}
double getTminPsat() const {return TminPsat;}
double getTref() const {return Tref;}
double getxref() const {return xref;}
double getTbase() const {return Tbase;}
double getxbase() const {return xbase;}
void setName(std::string name) {this->name = name;}
void setDescription(std::string description) {this->description = description;}
void setReference(std::string reference) {this->reference = reference;}
void setTmax(double Tmax) {this->Tmax = Tmax;}
void setTmin(double Tmin) {this->Tmin = Tmin;}
void setxmax(double xmax) {this->xmax = xmax;}
void setxmin(double xmin) {this->xmin = xmin;}
void setTminPsat(double TminPsat) {this->TminPsat = TminPsat;}
void setTref(double Tref) {this->Tref = Tref;}
void setxref(double xref) {this->xref = xref;}
void setTbase(double Tbase) {this->Tbase = Tbase;}
void setxbase(double xbase) {this->xbase = xbase;}
// Setters for the coefficients
void setcRho(Eigen::MatrixXd cRho){this->cRho = cRho;}
void setcHeat(Eigen::MatrixXd cHeat){this->cHeat = cHeat;}
void setcVisc(Eigen::MatrixXd cVisc){this->cVisc = cVisc;}
void setcCond(Eigen::MatrixXd cCond){this->cCond = cCond;}
void setcPsat(Eigen::MatrixXd cPsat){this->cPsat = cPsat;}
void setcTfreeze(Eigen::MatrixXd cTfreeze){this->cTfreeze = cTfreeze;}
void setcV2M(Eigen::MatrixXd cV2M){this->cV2M = cV2M;}
double getTInput(double curTValue){return curTValue-Tbase;}
double getxInput(double curxValue){return (curxValue-xbase)*100.0;}
/* All functions need T and p as input. Might not
* be necessary, but gives a clearer structure.
*/
/// Density as a function of temperature, pressure and composition.
virtual double rho (double T_K, double p, double x);
/// Heat capacities as a function of temperature, pressure and composition.
virtual double c (double T_K, double p, double x);
virtual double cp (double T_K, double p, double x){return c(T_K,p,x);};
virtual double cv (double T_K, double p, double x){return c(T_K,p,x);};
/// Entropy as a function of temperature, pressure and composition.
virtual double s (double T_K, double p, double x);
/// Internal energy as a function of temperature, pressure and composition.
virtual double u (double T_K, double p, double x);
/// Enthalpy as a function of temperature, pressure and composition.
virtual double h (double T_K, double p, double x);
/// Viscosity as a function of temperature, pressure and composition.
virtual double visc(double T_K, double p, double x);
/// Thermal conductivity as a function of temperature, pressure and composition.
virtual double cond(double T_K, double p, double x);
/// Saturation pressure as a function of temperature and composition.
virtual double psat(double T_K, double x );
/// Freezing temperature as a function of pressure and composition.
virtual double Tfreeze( double p, double x);
/// Conversion from volume-based to mass-based composition.
virtual double V2M( double x);
protected:
/* Define internal energy and enthalpy as functions of the
* other properties to provide data in case there are no
* coefficients.
*/
/// Enthalpy from u, p and rho.
/** Calculate enthalpy as a function of temperature and
* pressure employing functions for internal energy and
* density. Provides consistent formulations. */
double h_u(double T_K, double p, double x) {
return u(T_K,p,x)+p/rho(T_K,p,x);
};
/// Internal energy from h, p and rho.
/** Calculate internal energy as a function of temperature
* and pressure employing functions for enthalpy and
* density. Provides consistent formulations. */
double u_h(double T_K, double p, double x) {
return h(T_K,p,x)-p/rho(T_K,p,x);
};
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not cause errors. */
bool checkT(double T_K, double p, double x);
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool checkP(double T_K, double p, double x);
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool checkX(double x);
/// Check validity of temperature, pressure and composition input.
bool checkTPX(double T, double p, double x);
};
} /* namespace CoolProp */
#endif /* INCOMPRESSIBLE_H_ */

174
src/Backends/Incompressible/IncompressibleLibrary.cpp
View File

@@ -1,7 +1,181 @@
#include "IncompressibleLibrary.h"
#include "MatrixMath.h"
#include "rapidjson/rapidjson_include.h"
#include "all_incompressibles_JSON.h" // Makes a std::string variable called all_fluids_JSON
namespace CoolProp{
/// A general function to parse the json files that hold the coefficient matrices
IncompressibleData JSONIncompressibleLibrary::parse_coefficients(rapidjson::Value &obj, std::string id, bool vital){
IncompressibleData fluidData;
if (obj.HasMember(id.c_str())) {
//rapidjson::Value value = obj[id.c_str()];
if (obj[id.c_str()].HasMember("type")){
std::string type = cpjson::get_string(obj[id.c_str()], "type");
if (!type.compare("polynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exponential")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exppolynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else{
throw ValueError(format("The type [%s] is not understood for [%s] of incompressible fluids. Please check your JSON file.", type.c_str(), id.c_str()));
}
}
else{
throw ValueError(format("Your file does not have an entry for \"type\" of [%s], which is vital for this function.", id.c_str()));
}
}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
}
return fluidData;
}
/// Get a double from the JSON storage if it is defined, otherwise return def
double JSONIncompressibleLibrary::parse_value(rapidjson::Value &obj, std::string id, bool vital, double def=0.0){
if (obj.HasMember(id.c_str())) {return cpjson::get_double(obj, id);}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
else{
return def;
}
}
}
/// Add all the fluid entries in the rapidjson::Value instance passed in
void JSONIncompressibleLibrary::add_many(rapidjson::Value &listing) {
for (rapidjson::Value::ValueIterator itr = listing.Begin();
itr != listing.End(); ++itr) {
add_one(*itr);
}
};
void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
_is_empty = false;
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Add index->fluid mapping
fluid_map[index] = IncompressibleFluid();
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
fluid.setName(cpjson::get_string(fluid_json, "name"));
fluid.setDescription(cpjson::get_string(fluid_json, "description"));
fluid.setReference(cpjson::get_string(fluid_json, "reference"));
fluid.setTmax(parse_value(fluid_json, "Tmax", true, 0.0));
fluid.setTmin(parse_value(fluid_json, "Tmin", true, 0.0));
fluid.setxmax(parse_value(fluid_json, "xmax", false, 1.0));
fluid.setxmin(parse_value(fluid_json, "xmin", false, 0.0));
fluid.setTminPsat(parse_value(fluid_json, "TminPsat", false, 0.0));
fluid.setTbase(parse_value(fluid_json, "Tbase", false, 0.0));
fluid.setxbase(parse_value(fluid_json, "xbase", false, 0.0));
/// Setters for the coefficients
fluid.setDensity(parse_coefficients(fluid_json, "density", true));
fluid.setSpecificHeat(parse_coefficients(fluid_json, "specific_heat", true));
fluid.setViscosity(parse_coefficients(fluid_json, "viscosity", false));
fluid.setConductivity(parse_coefficients(fluid_json, "conductivity", false));
fluid.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
fluid.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
fluid.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false));
fluid.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false));
fluid.set_reference_state(
parse_value(fluid_json, "Tref", false, 25+273.15) ,
parse_value(fluid_json, "pref", false, 1.01325e5) ,
parse_value(fluid_json, "xref", false, 0.0) ,
parse_value(fluid_json, "href", false, 0.0) ,
parse_value(fluid_json, "sref", false, 0.0)
);
/// A function to check coefficients and equation types.
fluid.validate();
// Add name->index mapping
string_to_index_map[fluid.getName()] = index;
};
/// Get an IncompressibleFluid instance stored in this library
/**
@param name Name of the fluid
*/
IncompressibleFluid& JSONIncompressibleLibrary::get(std::string key) {
std::map<std::string, std::size_t>::iterator it;
// Try to find it
it = string_to_index_map.find(key);
// If it is found
if (it != string_to_index_map.end()) {
return get(it->second);
} else {
throw ValueError(
format(
"key [%s] was not found in string_to_index_map in JSONIncompressibleLibrary",
key.c_str()
)
);
}
};
/// Get a CoolPropFluid instance stored in this library
/**
@param key The index of the fluid in the map
*/
IncompressibleFluid& JSONIncompressibleLibrary::get(std::size_t key) {
std::map<std::size_t, IncompressibleFluid>::iterator it;
// Try to find it
it = fluid_map.find(key);
// If it is found
if (it != fluid_map.end()) {
return it->second;
} else {
throw ValueError(
format("key [%d] was not found in JSONIncompressibleLibrary",key));
}
};
static JSONIncompressibleLibrary library;

162
src/Backends/Incompressible/IncompressibleLibrary.h
View File

@@ -22,172 +22,44 @@ a rapidjson array of fluids to the add_many function.
*/
class JSONIncompressibleLibrary
{
/// Map from CAS code to JSON instance. For pseudo-pure fluids, use name in place of CAS code since no CASE number is defined for mixtures
/// Map from CAS code to JSON instance.
/** This is not practical for the incomressibles, the CAS may not be
* defined for blends of heat transfer fluids and solutions.
*/
std::map<std::size_t, IncompressibleFluid> fluid_map;
std::vector<std::string> name_vector;
std::map<std::string, std::size_t> string_to_index_map;
bool _is_empty;
protected:
/// A general function to parse the json files that hold the coefficient matrices
IncompressibleData parse_coefficients(rapidjson::Value &obj, std::string id, bool vital);
double parse_value(rapidjson::Value &obj, std::string id, bool vital, double def);
/// Parse the viscosity
void parse_viscosity(rapidjson::Value &viscosity, IncompressibleFluid & fluid)
{
if (viscosity.HasMember("type")){
std::string type = cpjson::get_string(viscosity, "type");
if (!type.compare("polynomial")){
fluid.viscosity.type = CoolProp::IncompressibleViscosityVariables::INCOMPRESSIBLE_VISCOSITY_POLYNOMIAL;
fluid.viscosity.poly.coeffs = cpjson::get_double_array(viscosity["coeffs"]);
return;
}
else{
throw ValueError(format("viscosity type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("viscosity does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the conductivity
void parse_conductivity(rapidjson::Value &conductivity, IncompressibleFluid & fluid)
{
if (conductivity.HasMember("type")){
std::string type = cpjson::get_string(conductivity, "type");
if (!type.compare("polynomial")){
fluid.conductivity.type = CoolProp::IncompressibleConductivityVariables::INCOMPRESSIBLE_CONDUCTIVITY_POLYNOMIAL;
fluid.conductivity.poly.coeffs = cpjson::get_double_array(conductivity["coeffs"]);
return;
}
else{
throw ValueError(format("conductivity type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("conductivity does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the specific_heat
void parse_specific_heat(rapidjson::Value &specific_heat, IncompressibleFluid & fluid)
{
if (specific_heat.HasMember("type")){
std::string type = cpjson::get_string(specific_heat, "type");
if (!type.compare("polynomial")){
fluid.specific_heat.type = CoolProp::IncompressibleSpecificHeatVariables::INCOMPRESSIBLE_SPECIFIC_HEAT_POLYNOMIAL; return;
fluid.specific_heat.poly.coeffs = cpjson::get_double_array(specific_heat["coeffs"]);
}
else{
throw ValueError(format("specific_heat type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("specific_heat does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the density
void parse_density(rapidjson::Value &density, IncompressibleFluid & fluid)
{
if (density.HasMember("type")){
std::string type = cpjson::get_string(density, "type");
if (!type.compare("polynomial")){
fluid.density.type = CoolProp::IncompressibleDensityVariables::INCOMPRESSIBLE_DENSITY_POLYNOMIAL; return;
fluid.density.poly.coeffs = cpjson::get_double_array(density["coeffs"]);
}
else{
throw ValueError(format("density type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("density does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Validate the fluid file that was just constructed
void validate(IncompressibleFluid & fluid)
{
}
public:
// Default constructor;
JSONIncompressibleLibrary(){
_is_empty = true;
};
JSONIncompressibleLibrary(){ _is_empty = true;};
bool is_empty(void){ return _is_empty;};
/// Add all the fluid entries in the rapidjson::Value instance passed in
void add_many(rapidjson::Value &listing)
{
for (rapidjson::Value::ValueIterator itr = listing.Begin(); itr != listing.End(); ++itr)
{
add_one(*itr);
}
};
void add_one(rapidjson::Value &fluid_json)
{
_is_empty = false;
void add_many(rapidjson::Value &listing);
void add_one(rapidjson::Value &fluid_json);
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Add index->fluid mapping
fluid_map[index] = IncompressibleFluid();
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
fluid.name = cpjson::get_string(fluid_json, "name");
fluid.Tmin = cpjson::get_double(fluid_json, "Tmin");
fluid.Tmax = cpjson::get_double(fluid_json, "Tmax");
parse_conductivity(fluid_json["conductivity"], fluid);
parse_density(fluid_json["density"], fluid);
parse_viscosity(fluid_json["viscosity"], fluid);
parse_specific_heat(fluid_json["specific_heat"], fluid);
// Add name->index mapping
string_to_index_map[fluid.name] = index;
};
/// Get an IncompressibleFluid instance stored in this library
/**
@param name Name of the fluid
*/
IncompressibleFluid& get(std::string key)
{
std::map<std::string, std::size_t>::iterator it;
// Try to find it
it = string_to_index_map.find(key);
// If it is found
if (it != string_to_index_map.end()){
return get(it->second);
}
else{
throw ValueError(format("key [%s] was not found in string_to_index_map in JSONIncompressibleLibrary",key.c_str()));
}
};
IncompressibleFluid& get(std::string key);
/// Get a CoolPropFluid instance stored in this library
/**
@param key The index of the fluid in the map
*/
IncompressibleFluid& get(std::size_t key)
{
std::map<std::size_t, IncompressibleFluid>::iterator it;
// Try to find it
it = fluid_map.find(key);
// If it is found
if (it != fluid_map.end()){
return it->second;
}
else{
throw ValueError(format("key [%d] was not found in JSONIncompressibleLibrary",key));
}
};
IncompressibleFluid& get(std::size_t key);
/// Return a comma-separated list of fluid names
std::string get_fluid_list(void)
{
return strjoin(name_vector, ",");
};
std::string get_fluid_list(void){ return strjoin(name_vector, ",");};
};
/// Get a reference to the library instance

492
src/IncompressibleFluid.cpp Normal file
View File

@@ -0,0 +1,492 @@
#include "IncompressibleFluid.h"
#include "math.h"
#include "MatrixMath.h"
#include "PolyMath.h"
namespace CoolProp {
/// A thermophysical property provider for critical and reducing values as well as derivatives of Helmholtz energy
/**
This fluid instance is populated using an entry from a JSON file
*/
//IncompressibleFluid::IncompressibleFluid();
void IncompressibleFluid::set_reference_state(double T0, double p0, double x0=0.0, double h0=0.0, double s0=0.0){
this->rhoref = rho(T0,p0,x0);
this->pref = p0;
this->uref = 0.0;
this->uref = u(T0,p0,x0);//(href - pref/rhoref);
this->href = h0; // set new reference value
this->sref = s0; // set new reference value
this->href = h(T0,p0,x0); // adjust offset to fit to equations
this->sref = s(T0,p0,x0); // adjust offset to fit to equations
}
void IncompressibleFluid::validate(){throw NotImplementedError("TODO");}
/// Base function that handles the custom data type
double IncompressibleFluid::baseFunction(IncompressibleData data, double x_in, double y_in=0.0){
switch (data.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.evaluate(data.coeffs, x_in, y_in, 0, 0, Tbase, xbase);
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
//throw NotImplementedError("Here you should implement the exponential.");
poly.checkCoefficients(data.coeffs, 3,1);
return exp( data.coeffs(0,0) / ( (x_in-Tbase)+data.coeffs(1,0) ) - data.coeffs(2,0) );
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
//throw NotImplementedError("Here you should implement the exponential polynomial.");
return exp(poly.evaluate(data.coeffs, x_in, y_in, 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__,data.type));
break;
default:
throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,data.type));
break;
}
return -_HUGE;
}
/// Entropy as a function of temperature, pressure and composition.
double IncompressibleFluid::s (double T, double p, double x=0.0){
IncompressibleData data = specific_heat;
switch (data.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(data.coeffs, T, x, 0, -1, 0, Tbase, xbase) - sref;
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential function.",__FILE__,__LINE__));
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential polynomial function.",__FILE__,__LINE__));
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__,data.type));
break;
default:
throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,data.type));
break;
}
return -_HUGE;
}
/// Internal energy as a function of temperature, pressure and composition.
double IncompressibleFluid::u (double T, double p, double x=0.0){
IncompressibleData data = specific_heat;
switch (data.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(data.coeffs, T, x, 0, 0, 0, Tbase, xbase) - uref;
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential function.",__FILE__,__LINE__));
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential polynomial function.",__FILE__,__LINE__));
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__,data.type));
break;
default:
throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,data.type));
break;
}
return -_HUGE;
}
/// Freezing temperature as a function of pressure and composition.
double Tfreeze( double p, double x){throw NotImplementedError("TODO");}
/// Define freezing point calculations
//double Tfreeze(double p, double x){
// IncompressibleClass::checkCoefficients(cTfreeze,5);
// std::vector<double> tmpVector(cTfreeze.begin()+1,cTfreeze.end());
// return polyval(tmpVector, x*100.0-cTfreeze[0])+273.15;
//}
/// Conversion from volume-based to mass-based composition.
double V2M (double T, double y){throw NotImplementedError("TODO");}
/// Conversion from mass-based to mole-based composition.
double M2M (double T, double x){throw NotImplementedError("TODO");}
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not cause errors. */
bool IncompressibleFluid::checkT(double T, double p, double x=0.0){throw NotImplementedError("TODO");}
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool IncompressibleFluid::checkP(double T, double p, double x=0.0){throw NotImplementedError("TODO");}
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool IncompressibleFluid::checkX(double x){throw NotImplementedError("TODO");}
/// Check validity of temperature, pressure and composition input.
bool IncompressibleFluid::checkTPX(double T, double p, double x=0.0){throw NotImplementedError("TODO");}
} /* namespace CoolProp */
// Testing still needs to be enhanced.
/* Below, I try to carry out some basic tests for both 2D and 1D
* polynomials as well as the exponential functions for vapour
* pressure etc.
*/
#ifdef ENABLE_CATCH
#include <math.h>
#include <iostream>
#include "catch.hpp"
TEST_CASE("Internal consistency checks and example use cases for the incompressible fluids","[IncompressibleFluids]")
{
bool PRINT = false;
std::string tmpStr;
std::vector<double> tmpVector;
std::vector< std::vector<double> > tmpMatrix;
SECTION("Test case for \"SylthermXLT\" by Dow Chemicals") {
std::vector<double> cRho;
cRho.push_back(+1.1563685145E+03);
cRho.push_back(-1.0269048032E+00);
cRho.push_back(-9.3506079577E-07);
cRho.push_back(+1.0368116627E-09);
CoolProp::IncompressibleData density;
density.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
density.coeffs = CoolProp::vec_to_eigen(cRho);
std::vector<double> cHeat;
cHeat.push_back(+1.1562261074E+03);
cHeat.push_back(+2.0994549103E+00);
cHeat.push_back(+7.7175381057E-07);
cHeat.push_back(-3.7008444051E-20);
CoolProp::IncompressibleData specific_heat;
specific_heat.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
specific_heat.coeffs = CoolProp::vec_to_eigen(cHeat);
std::vector<double> cCond;
cCond.push_back(+1.6121957379E-01);
cCond.push_back(-1.3023781944E-04);
cCond.push_back(-1.4395238766E-07);
CoolProp::IncompressibleData conductivity;
conductivity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
conductivity.coeffs = CoolProp::vec_to_eigen(cCond);
std::vector<double> cVisc;
cVisc.push_back(+1.0337654989E+03);
cVisc.push_back(-4.3322764383E+01);
cVisc.push_back(+1.0715062356E+01);
CoolProp::IncompressibleData viscosity;
viscosity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
viscosity.coeffs = CoolProp::vec_to_eigen(cVisc);
CoolProp::IncompressibleFluid XLT;
XLT.setName("XLT");
XLT.setDescription("SylthermXLT");
XLT.setReference("Dow Chemicals data sheet");
XLT.setTmax(533.15);
XLT.setTmin(173.15);
XLT.setxmax(0.0);
XLT.setxmin(0.0);
XLT.setTminPsat(533.15);
XLT.setTbase(0.0);
XLT.setxbase(0.0);
/// Setters for the coefficients
XLT.setDensity(density);
XLT.setSpecificHeat(specific_heat);
XLT.setViscosity(viscosity);
XLT.setConductivity(conductivity);
//XLT.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
//XLT.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
//XLT.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false));
//XLT.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false));
//XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0);
double Tref = 25+273.15;
double pref = 0.0;
XLT.set_reference_state(Tref, pref, 0.0, 0.0, 0.0);
/// A function to check coefficients and equation types.
//XLT.validate();
// Prepare the results and compare them to the calculated values
double acc = 0.0001;
double T = 273.15+50;
double p = 10e5;
double val = 0;
double res = 0;
// Compare density
val = 824.4615702148608;
res = XLT.rho(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare cp
val = 1834.7455527670554;
res = XLT.c(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare s
val = 145.59157247249246;
res = XLT.s(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = XLT.s(Tref,pref);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare u
val = 45212.407309106304;
res = XLT.u(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = XLT.u(Tref,pref);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare h
val = 46425.32011926845;
res = XLT.h(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = XLT.h(Tref,pref);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare v
val = 0.0008931435169681835;
res = XLT.visc(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare l
val = 0.10410086156049088;
res = XLT.cond(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
}
/// Test case for Methanol from SecCool
// name = std::string("SecCoolSolution");
// description = std::string("Test Methanol SecCool");
// reference = std::string("Test");
//
// Tmin = -50 + 273.15;
// Tmax = 20 + 273.15;
// TminPsat = Tmax;
//
// xmin = 0.0;
// xmax = 0.5;
//
// Tbase = -4.48 + 273.15;
// xbase = 31.57 / 100.0;
//
// 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);
// cRho.clear();
// cRho = makeMatrix(tmpVector);
//
//
//
//
// std::vector<double> 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]);
// tmpVector.push_back(coefficients[7]);
// tmpVector.push_back(coefficients[12]);
// tmpVector.push_back(coefficients[16]);
// matrix.push_back(tmpVector);
//
// tmpVector.clear();
// tmpVector.push_back(coefficients[2]);
// tmpVector.push_back(coefficients[8]);
// tmpVector.push_back(coefficients[13]);
// tmpVector.push_back(coefficients[17]);
// matrix.push_back(tmpVector);
//
// tmpVector.clear();
// tmpVector.push_back(coefficients[3]);
// tmpVector.push_back(coefficients[9]);
// tmpVector.push_back(coefficients[14]);
// tmpVector.push_back(0.0);
// matrix.push_back(tmpVector);
//
// tmpVector.clear();
// tmpVector.push_back(coefficients[4]);
// tmpVector.push_back(coefficients[10]);
// tmpVector.push_back(0.0);
// tmpVector.push_back(0.0);
// matrix.push_back(tmpVector);
//
// tmpVector.clear();
// tmpVector.push_back(coefficients[5]);
// tmpVector.push_back(0.0);
// tmpVector.push_back(0.0);
// tmpVector.push_back(0.0);
// matrix.push_back(tmpVector);
//
//
//
//
//
//
//
//
//
//
//
//
//
//
//
//
//
//
// cTfreeze.clear();
// cTfreeze.push_back( 27.755555600); // reference concentration in per cent
// cTfreeze.push_back(-22.973221700);
// cTfreeze.push_back(-1.1040507200);
// cTfreeze.push_back(-0.0120762281);
// cTfreeze.push_back(-9.343458E-05);
//
//
//
//
//
//
//
//
// double deltaT = 0.1;
// double Tmin = 273.15- 50;
// double Tmax = 273.15+250;
// double Tinc = 200;
//
// std::vector<std::vector<double> > cHeat2D;
// cHeat2D.push_back(cHeat);
// cHeat2D.push_back(cHeat);
// cHeat2D.push_back(cHeat);
//
// Eigen::MatrixXd matrix2D = CoolProp::vec_to_eigen(cHeat2D);
//
// Eigen::MatrixXd matrix2Dtmp;
// std::vector<std::vector<double> > vec2Dtmp;
//
// SECTION("Coefficient parsing") {
// CoolProp::Polynomial2D poly;
// CHECK_THROWS(poly.checkCoefficients(matrix2D,4,5));
// CHECK( poly.checkCoefficients(matrix2D,3,4) );
// }
}
#endif /* ENABLE_CATCH */

85
src/PolyMath.cpp
View File

@@ -30,10 +30,10 @@ bool Polynomial2D::checkCoefficients(const Eigen::MatrixXd &coefficients, const
if (coefficients.cols() == columns) {
return true;
} else {
throw ValueError(format("The number of columns %d does not match with %d. ",coefficients.cols(),columns));
throw ValueError(format("%s (%d): The number of columns %d does not match with %d. ",__FILE__,__LINE__,coefficients.cols(),columns));
}
} else {
throw ValueError(format("The number of rows %d does not match with %d. ",coefficients.rows(),rows));
throw ValueError(format("%s (%d): The number of rows %d does not match with %d. ",__FILE__,__LINE__,coefficients.rows(),rows));
}
return false;
}
@@ -52,7 +52,7 @@ bool Polynomial2D::checkCoefficients(const Eigen::MatrixXd &coefficients, const
/// @param axis integer value that represents the desired direction of integration
/// @param times integer value that represents the desired order of integration
Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis = -1, const int &times = 1){
if (times < 0) throw ValueError(format("You have to provide a positive order for integration, %d is not valid. ",times));
if (times < 0) throw ValueError(format("%s (%d): You have to provide a positive order for integration, %d is not valid. ",__FILE__,__LINE__,times));
if (times == 0) return Eigen::MatrixXd(coefficients);
Eigen::MatrixXd oldCoefficients;
Eigen::MatrixXd newCoefficients(coefficients);
@@ -65,7 +65,7 @@ Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficient
newCoefficients = Eigen::MatrixXd(coefficients.transpose());
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -87,7 +87,7 @@ Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficient
newCoefficients.transposeInPlace();
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -103,7 +103,7 @@ Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficient
/// @param axis integer value that represents the desired direction of derivation
/// @param times integer value that represents the desired order of integration
Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients, const int &axis = -1, const int &times = 1){
if (times < 0) throw ValueError(format("You have to provide a positive order for derivation, %d is not valid. ",times));
if (times < 0) throw ValueError(format("%s (%d): You have to provide a positive order for derivation, %d is not valid. ",__FILE__,__LINE__,times));
if (times == 0) return Eigen::MatrixXd(coefficients);
// Recursion is also possible, but not recommended
//Eigen::MatrixXd newCoefficients;
@@ -119,7 +119,7 @@ Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients,
newCoefficients = Eigen::MatrixXd(coefficients.transpose());
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -139,7 +139,7 @@ Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients,
newCoefficients.transposeInPlace();
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -158,7 +158,7 @@ Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients,
/// @param x_in double value that represents the current input
double Polynomial2D::evaluate(const Eigen::MatrixXd &coefficients, const double &x_in){
if (coefficients.rows() != 1) {
throw ValueError(format("You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",coefficients.rows(),coefficients.cols()));
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;
@@ -218,7 +218,7 @@ Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd &coefficients, const d
}
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
tmpCoefficients(0,0) -= z_in;
@@ -317,7 +317,7 @@ Poly2DResidual::Poly2DResidual(Polynomial2D &poly, const Eigen::MatrixXd &coeffi
this->axis = axis;
break;
default:
throw ValueError(format("You have to provide a dimension to the solver, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension to the solver, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -376,7 +376,7 @@ double Poly2DResidual::deriv(double target){
/// @param times integer value that represents the desired order of derivation
/// @param firstExponent integer value that represents the lowest exponent of the polynomial in axis direction
Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int &times, const int &firstExponent){
if (times < 0) throw ValueError(format("You have to provide a positive order for derivation, %d is not valid. ",times)); if (times < 0) throw ValueError(format("You have to provide a positive order for derivation, %d is not valid. ",times));
if (times < 0) throw ValueError(format("%s (%d): You have to provide a positive order for derivation, %d is not valid. ",__FILE__,__LINE__,times));
if (times == 0) return Eigen::MatrixXd(coefficients);
// Recursion is also possible, but not recommended
//Eigen::MatrixXd newCoefficients;
@@ -392,7 +392,7 @@ Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd &coefficien
newCoefficients = Eigen::MatrixXd(coefficients.transpose());
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -413,7 +413,7 @@ Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd &coefficien
newCoefficients.transposeInPlace();
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -438,10 +438,10 @@ Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd &coefficien
/// @param x_base double value that represents the base value for a centred fit in the 1st dimension
double Polynomial2DFrac::evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const int &firstExponent = 0, const double &x_base = 0.0){
if (coefficients.rows() != 1) {
throw ValueError(format("You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",coefficients.rows(),coefficients.cols()));
throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",__FILE__,__LINE__,coefficients.rows(),coefficients.cols()));
}
if ( (firstExponent<0) && (fabs(x_in-x_base)<DBL_EPSILON)) {
throw ValueError(format("A fraction cannot be evaluated with zero as denominator, x_in-x_base=%f ",x_in-x_base));
throw ValueError(format("%s (%d): A fraction cannot be evaluated with zero as denominator, x_in-x_base=%f ",__FILE__,__LINE__,x_in-x_base));
}
Eigen::MatrixXd tmpCoeffs(coefficients);
Eigen::MatrixXd newCoeffs;
@@ -480,10 +480,10 @@ double Polynomial2DFrac::evaluate(const Eigen::MatrixXd &coefficients, const dou
/// @param y_base double value that represents the base value for a centred fit in the 2nd dimension
double Polynomial2DFrac::evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &x_exp, const int &y_exp, const double &x_base = 0.0, const double &y_base = 0.0){
if ( (x_exp<0) && (fabs(x_in-x_base)<DBL_EPSILON)) {
throw ValueError(format("A fraction cannot be evaluated with zero as denominator, x_in-x_base=%f ",x_in-x_base));
throw ValueError(format("%s (%d): A fraction cannot be evaluated with zero as denominator, x_in-x_base=%f ",__FILE__,__LINE__,x_in-x_base));
}
if ( (y_exp<0) && (fabs(y_in-y_base)<DBL_EPSILON)) {
throw ValueError(format("A fraction cannot be evaluated with zero as denominator, y_in-y_base=%f ",y_in-y_base));
throw ValueError(format("%s (%d): A fraction cannot be evaluated with zero as denominator, y_in-y_base=%f ",__FILE__,__LINE__,y_in-y_base));
}
Eigen::MatrixXd tmpCoeffs(coefficients);
@@ -554,7 +554,7 @@ double Polynomial2DFrac::derivative(const Eigen::MatrixXd &coefficients, const d
other_base = x_base;
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -600,21 +600,32 @@ double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const dou
other_base = x_base;
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
if (int_exp<-1) throw NotImplementedError(format("This function is only implemented for lowest exponents >= -1, int_exp=%d ",int_exp));
if (int_exp<-1) throw NotImplementedError(format("%s (%d): This function is only implemented for lowest exponents >= -1, int_exp=%d ",__FILE__,__LINE__,int_exp));
double result = 0;
size_t r = newCoefficients.rows();
size_t c = newCoefficients.cols();
if (int_exp==-1) {
Eigen::MatrixXd tmpCoefficients = newCoefficients.row(0) * log(int_val-int_base);
newCoefficients = integrateCoeffs(newCoefficients.block(1,0,r-1,c), 0, 1);
newCoefficients.row(0) = tmpCoefficients;
return evaluate(newCoefficients,int_val,other_val,0,other_exp,int_base,other_base);
if (fabs(int_base)<DBL_EPSILON){
Eigen::MatrixXd tmpCoefficients = newCoefficients.row(0) * log(int_val-int_base);
newCoefficients = integrateCoeffs(newCoefficients.block(1,0,r-1,c), 0, 1);
newCoefficients.row(0) = tmpCoefficients;
return evaluate(newCoefficients,int_val,other_val,0,other_exp,int_base,other_base);
}
else {
// Reduce the coefficients to the integration dimension:
newCoefficients = Eigen::MatrixXd(r,1);
for (int i=0; i<r; i++){
newCoefficients(i,0) = evaluate(coefficients.row(i).transpose(), other_val, other_exp, other_base);
}
return fracIntCentral(newCoefficients.transpose(),int_val,int_base);
}
}
Eigen::MatrixXd tmpCoeffs;
@@ -661,7 +672,7 @@ Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, con
input = in - x_base;
break;
default:
throw ValueError(format("You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ",axis));
throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ",__FILE__,__LINE__,axis));
break;
}
@@ -682,7 +693,7 @@ Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, con
tmpCoefficients = Eigen::VectorXd::Zero(r+std::max(diff,0));
tmpCoefficients.block(0,0,r,1) = newCoefficients.block(0,0,r,1);
tmpCoefficients(r+diff-1,0) -= z_in;
throw NotImplementedError(format("Currently, there is no solver for the generalised polynomial, an exponent of %d is not valid. ",solve_exp));
throw NotImplementedError(format("%s (%d): Currently, there is no solver for the generalised polynomial, an exponent of %d is not valid. ",__FILE__,__LINE__,solve_exp));
}
if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string( Eigen::MatrixXd(tmpCoefficients) ) << std::endl;
@@ -768,11 +779,11 @@ double Polynomial2DFrac::binom(const int &nValue, const int &nValue2){
/// @param x_in double value that represents the current input
/// @param x_base double value that represents the basis for the fit
Eigen::MatrixXd Polynomial2DFrac::fracIntCentralDvector(const int &m, const double &x_in, const double &x_base){
if (m<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",m));
if (m<1) throw ValueError(format("%s (%d): You have to provide coefficients, a vector length of %d is not a valid. ",__FILE__,__LINE__,m));
Eigen::MatrixXd D = Eigen::MatrixXd::Zero(1,m);
double tmp;
// TODO: This can be optimized using the Horner scheme!
for (int j=0; j<m; j++){ // loop through row
tmp = pow(-1.0,j) * log(x_in) * pow(x_base,j);
for(int k=0; k<j; k++) { // internal loop for every entry
@@ -788,7 +799,7 @@ Eigen::MatrixXd Polynomial2DFrac::fracIntCentralDvector(const int &m, const doub
/// @param x_base double value that represents the basis for the fit
double Polynomial2DFrac::fracIntCentral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &x_base){
if (coefficients.rows() != 1) {
throw ValueError(format("You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",coefficients.rows(),coefficients.cols()));
throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",__FILE__,__LINE__,coefficients.rows(),coefficients.cols()));
}
int m = coefficients.cols();
Eigen::MatrixXd D = fracIntCentralDvector(m, x_in, x_base);
@@ -1163,6 +1174,20 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","
}
T = 423.15;
c = 3460.895272;
d = frac.integral(matrix, T, 0.0, 0, -1, 0, 348.15, 0.0);
{
CAPTURE(T);
CAPTURE(c);
CAPTURE(d);
tmpStr = CoolProp::mat_to_string(matrix);
CAPTURE(tmpStr);
CHECK( check_abs(c,d,acc) );
}
deltaT = 0.01;
for (T = Tmin; T<Tmax; T+=Tinc) {
a = poly.evaluate(matrix, T-deltaT, y);