From de03057f48af66560b2a211f8a8b458893e6f09b Mon Sep 17 00:00:00 2001 From: jowr Date: Mon, 23 Jun 2014 17:34:37 +0200 Subject: [PATCH] Added the first incompressible test cases and it seems to work... --- include/IncompressibleFluid.h | 234 +++++++-- .../Helmholtz/Fluids/Incompressible.cpp | 120 ----- .../Helmholtz/Fluids/Incompressible.h | 169 ------ .../Incompressible/IncompressibleLibrary.cpp | 174 +++++++ .../Incompressible/IncompressibleLibrary.h | 162 +----- src/IncompressibleFluid.cpp | 492 ++++++++++++++++++ src/PolyMath.cpp | 85 +-- 7 files changed, 920 insertions(+), 516 deletions(-) delete mode 100644 src/Backends/Helmholtz/Fluids/Incompressible.cpp delete mode 100644 src/Backends/Helmholtz/Fluids/Incompressible.h create mode 100644 src/IncompressibleFluid.cpp diff --git a/include/IncompressibleFluid.h b/include/IncompressibleFluid.h index a6e509a9..2c25730a 100644 --- a/include/IncompressibleFluid.h +++ b/include/IncompressibleFluid.h @@ -19,70 +19,200 @@ #include #include +#include +#include "PolyMath.h" +#include "MatrixMath.h" + namespace CoolProp { -struct IncompressiblePolynomialData{ - std::vector coeffs; -}; - -struct IncompressibleViscosityVariables{ - enum IncompressibleViscosityEnum{ - INCOMPRESSIBLE_VISCOSITY_POLYNOMIAL, - INCOMPRESSIBLE_VISCOSITY_NOT_SET - }; - int type; - IncompressiblePolynomialData poly; - IncompressibleViscosityVariables(){type = INCOMPRESSIBLE_VISCOSITY_NOT_SET;}; -}; - -struct IncompressibleConductivityVariables{ - enum IncompressibleConductivityEnum{ - INCOMPRESSIBLE_CONDUCTIVITY_POLYNOMIAL, - INCOMPRESSIBLE_CONDUCTIVITY_NOT_SET - }; - int type; - IncompressiblePolynomialData poly; - IncompressibleConductivityVariables(){type = INCOMPRESSIBLE_CONDUCTIVITY_NOT_SET;}; -}; - -struct IncompressibleDensityVariables{ - enum IncompressibleDensityEnum{ - INCOMPRESSIBLE_DENSITY_POLYNOMIAL, - INCOMPRESSIBLE_DENSITY_NOT_SET - }; - int type; - IncompressiblePolynomialData poly; - IncompressibleDensityVariables(){type = INCOMPRESSIBLE_DENSITY_NOT_SET;}; -}; - -struct IncompressibleSpecificHeatVariables{ - enum IncompressibleSpecificHeatEnum{ - INCOMPRESSIBLE_SPECIFIC_HEAT_POLYNOMIAL, - INCOMPRESSIBLE_SPECIFIC_HEAT_NOT_SET - }; - int type; - IncompressiblePolynomialData poly; - IncompressibleSpecificHeatVariables(){type = INCOMPRESSIBLE_SPECIFIC_HEAT_NOT_SET;}; +struct IncompressibleData { + int type; + enum IncompressibleTypeEnum { + INCOMPRESSIBLE_POLYNOMIAL, + INCOMPRESSIBLE_EXPONENTIAL, + INCOMPRESSIBLE_EXPPOLYNOMIAL, + INCOMPRESSIBLE_NOT_SET + }; + Eigen::MatrixXd coeffs; //TODO: Can we store the Eigen::Matrix objects more efficiently? + //std::vector > coeffs; + IncompressibleData() { + type = INCOMPRESSIBLE_NOT_SET; + }; }; /// 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 */ -struct IncompressibleFluid { +class IncompressibleFluid{ - std::string name; - std::string reference; +protected: + std::string name; + std::string description; + std::string reference; - double Tmin, Tmax; + double Tmin, Tmax; + double xmin, xmax; - IncompressibleViscosityVariables viscosity; - IncompressibleConductivityVariables conductivity; - IncompressibleSpecificHeatVariables specific_heat; - IncompressibleDensityVariables density; + double TminPsat; + double xref, Tref, pref; + double href, sref; + double uref, rhoref; + double xbase, Tbase; + IncompressibleData density; + IncompressibleData specific_heat; + IncompressibleData viscosity; + IncompressibleData conductivity; + IncompressibleData p_sat; + IncompressibleData T_freeze; + IncompressibleData volToMass; + IncompressibleData massToMole; + + Polynomial2DFrac poly; + + // Forward declaration of the some internal functions + //double h_u(double T, double p, double x); + //double u_h(double T, double p, double x); + +public: + IncompressibleFluid(){}; + virtual ~IncompressibleFluid(){}; + + 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 getpref() const {return pref;} + double getxref() const {return xref;} + double gethref() const {return href;} + double getsref() const {return sref;} + 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 setpref(double pref) {this->pref = pref;} + //void setxref(double xref) {this->xref = xref;} + void set_reference_state(double T0, double p0, double x0, double h0, double s0); + void setTbase(double Tbase) {this->Tbase = Tbase;} + void setxbase(double xbase) {this->xbase = xbase;} + + /// Setters for the coefficients + void setDensity(IncompressibleData density){this->density = density;} + void setSpecificHeat(IncompressibleData specific_heat){this->specific_heat = specific_heat;} + void setViscosity(IncompressibleData viscosity){this->viscosity = viscosity;} + void setConductivity(IncompressibleData conductivity){this->conductivity = conductivity;} + void setPsat(IncompressibleData p_sat){this->p_sat = p_sat;} + void setTfreeze(IncompressibleData T_freeze){this->T_freeze = T_freeze;} + void setVolToMass(IncompressibleData volToMass){this->volToMass = volToMass;} + void setMassToMole(IncompressibleData massToMole){this->massToMole = massToMole;} + + /// A function to check coefficients and equation types. + void validate(); + + +protected: + /// Base function that handles the custom data type, just a place holder to show the structure. + double baseFunction(IncompressibleData data, double x_in, double y_in); + +public: + /* 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 rho (double T, double p, double x=0.0){return baseFunction(density, T, x);}; + /// Heat capacities as a function of temperature, pressure and composition. + double c (double T, double p, double x=0.0){return baseFunction(specific_heat, T, x);}; + double cp (double T, double p, double x=0.0){return c(T,p,x);}; + double cv (double T, double p, double x=0.0){return c(T,p,x);}; + /// Entropy as a function of temperature, pressure and composition. + double s (double T, double p, double x); + /// Internal energy as a function of temperature, pressure and composition. + double u (double T, double p, double x); + /// Enthalpy as a function of temperature, pressure and composition. + double h (double T, double p, double x=0.0){return h_u(T,p,x);}; + /// Viscosity as a function of temperature, pressure and composition. + double visc(double T, double p, double x=0.0){return baseFunction(viscosity, T, x);}; + /// Thermal conductivity as a function of temperature, pressure and composition. + double cond(double T, double p, double x=0.0){return baseFunction(conductivity, T, x);}; + /// Saturation pressure as a function of temperature and composition. + double psat(double T, double x=0.0){return baseFunction(p_sat, T, x);}; + /// Freezing temperature as a function of pressure and composition. + double Tfreeze( double p, double x); + /// Conversion from volume-based to mass-based composition. + double V2M (double T, double y); + /// Conversion from mass-based to mole-based composition. + double M2M (double T, 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, double p, double x) { + return u(T,p,x)+p/rho(T,p,x)-href; + }; + + /// 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, double p, double x) { + return h(T,p,x)-p/rho(T,p,x)+href; + }; + + + /* + * 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, 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, 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 /* COOLPROPFLUID_H_ */ +#endif /* INCOMPRESSIBLEFLUID_H_ */ diff --git a/src/Backends/Helmholtz/Fluids/Incompressible.cpp b/src/Backends/Helmholtz/Fluids/Incompressible.cpp deleted file mode 100644 index d176a7ba..00000000 --- a/src/Backends/Helmholtz/Fluids/Incompressible.cpp +++ /dev/null @@ -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 (px) || (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 */ diff --git a/src/Backends/Helmholtz/Fluids/Incompressible.h b/src/Backends/Helmholtz/Fluids/Incompressible.h deleted file mode 100644 index 675228ef..00000000 --- a/src/Backends/Helmholtz/Fluids/Incompressible.h +++ /dev/null @@ -1,169 +0,0 @@ -/* - * Incompressible.h - * - * Created on: 20 Dec 2013 - * Author: jowr - */ - -#ifndef INCOMPRESSIBLE_H_ -#define INCOMPRESSIBLE_H_ - -#include -#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_ */ diff --git a/src/Backends/Incompressible/IncompressibleLibrary.cpp b/src/Backends/Incompressible/IncompressibleLibrary.cpp index 70343e37..0bc9ebcf 100644 --- a/src/Backends/Incompressible/IncompressibleLibrary.cpp +++ b/src/Backends/Incompressible/IncompressibleLibrary.cpp @@ -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::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::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; diff --git a/src/Backends/Incompressible/IncompressibleLibrary.h b/src/Backends/Incompressible/IncompressibleLibrary.h index fe3fa41f..8d0f2087 100644 --- a/src/Backends/Incompressible/IncompressibleLibrary.h +++ b/src/Backends/Incompressible/IncompressibleLibrary.h @@ -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 fluid_map; std::vector name_vector; std::map 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::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::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 diff --git a/src/IncompressibleFluid.cpp b/src/IncompressibleFluid.cpp new file mode 100644 index 00000000..4f145e31 --- /dev/null +++ b/src/IncompressibleFluid.cpp @@ -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 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 +#include +#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 tmpVector; + std::vector< std::vector > tmpMatrix; + + + SECTION("Test case for \"SylthermXLT\" by Dow Chemicals") { + + std::vector 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 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 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 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 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 > 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 > vec2Dtmp; +// +// SECTION("Coefficient parsing") { +// CoolProp::Polynomial2D poly; +// CHECK_THROWS(poly.checkCoefficients(matrix2D,4,5)); +// CHECK( poly.checkCoefficients(matrix2D,3,4) ); +// } + + +} + + +#endif /* ENABLE_CATCH */ diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 73baeac6..8a967953 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -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 × = 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 × = 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 ×, 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)= -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)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