diff --git a/.gitignore b/.gitignore index 5350b05c..70df2081 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,4 @@ /include/mixture_reducing_parameters_JSON.h /include/version.h /include/catch.hpp +/build/ diff --git a/dev/fluids/CycloHexane.json b/dev/fluids/CycloHexane.json index e6b8f72d..909c82e5 100644 --- a/dev/fluids/CycloHexane.json +++ b/dev/fluids/CycloHexane.json @@ -107,17 +107,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 553.64, - "a": [ - 0.06485 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.263 - ] } }, "CAS": "110-82-7", diff --git a/dev/fluids/IsoButane.json b/dev/fluids/IsoButane.json index 4d5a36b9..41e7ce52 100644 --- a/dev/fluids/IsoButane.json +++ b/dev/fluids/IsoButane.json @@ -110,19 +110,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 407.81, - "a": [ - -0.01639, - 0.06121 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 2.102, - 1.304 - ] } }, "CAS": "75-28-5", diff --git a/dev/fluids/IsoButene.json b/dev/fluids/IsoButene.json index 9f3ca344..9e3143ad 100644 --- a/dev/fluids/IsoButene.json +++ b/dev/fluids/IsoButene.json @@ -107,17 +107,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 418.09, - "a": [ - 0.0545 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.23 - ] } }, "CAS": "115-11-7", diff --git a/dev/fluids/ParaHydrogen.json b/dev/fluids/ParaHydrogen.json index 677524bb..71754c4b 100644 --- a/dev/fluids/ParaHydrogen.json +++ b/dev/fluids/ParaHydrogen.json @@ -108,17 +108,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 32.938, - "a": [ - 0.005314 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.06 - ] } }, "CAS": "1333-74-0p", diff --git a/dev/fluids/R152A.json b/dev/fluids/R152A.json index 179be08c..5f81bcf3 100644 --- a/dev/fluids/R152A.json +++ b/dev/fluids/R152A.json @@ -106,17 +106,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 386.411, - "a": [ - 0.05808 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.2115 - ] } }, "CAS": "75-37-6", diff --git a/dev/fluids/R227EA.json b/dev/fluids/R227EA.json index ed82915f..0a5f7605 100644 --- a/dev/fluids/R227EA.json +++ b/dev/fluids/R227EA.json @@ -106,21 +106,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 374.9, - "a": [ - 0.06127, - -0.009516, - -0.00192 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.192, - 0.9795, - 1.421 - ] } }, "CAS": "431-89-0", diff --git a/dev/fluids/R236EA.json b/dev/fluids/R236EA.json index 8bb1075f..52dc8df0 100644 --- a/dev/fluids/R236EA.json +++ b/dev/fluids/R236EA.json @@ -106,19 +106,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 412.44, - "a": [ - 0.306974, - -0.247277 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.12614, - 1.09899 - ] } }, "CAS": "431-63-0", diff --git a/dev/fluids/R236FA.json b/dev/fluids/R236FA.json index 41f01557..c6ed7fe2 100644 --- a/dev/fluids/R236FA.json +++ b/dev/fluids/R236FA.json @@ -106,17 +106,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 398.07, - "a": [ - 0.05389 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.249 - ] } }, "CAS": "690-39-1", diff --git a/dev/fluids/R365MFC.json b/dev/fluids/R365MFC.json index 6a8e4922..6dd25025 100644 --- a/dev/fluids/R365MFC.json +++ b/dev/fluids/R365MFC.json @@ -106,17 +106,6 @@ ], "type": "rhoV", "using_tau_r": true - }, - "surface_tension": { - "BibTeX": "Mulero-JPCRD-2012", - "Tc": 460, - "a": [ - 0.0534 - ], - "description": "sigma = sum(a_i*(1-T/Tc)^n_i)", - "n": [ - 1.21 - ] } }, "CAS": "406-58-6", diff --git a/include/PolyMath.h b/include/PolyMath.h index ae85fb5e..e02acb3d 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -209,7 +209,7 @@ public: /// @param coefficients vector containing the ordered coefficients /// @param x double value that represents the current position virtual inline double polyfracval(const std::vector &coefficients, double x){ - return baseHornerFracInt(coefficients,x); + return baseHorner(coefficients,x)/x; } /// Evaluates the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable @@ -217,7 +217,7 @@ public: /// @param x double value that represents the current input in the 1st dimension /// @param y double value that represents the current input in the 2nd dimension virtual inline double polyfracval(const std::vector< std::vector > &coefficients, double x, double y){ - return baseHornerFracInt(coefficients,x,y); + return baseHorner(coefficients,x,y)/y; } @@ -333,6 +333,15 @@ public: virtual double call(double x); virtual double deriv(double x); }; +class PolyFracIntCentralResidual : public PolyResidual { +protected: + double baseVal; +public: + PolyFracIntCentralResidual(const std::vector &coefficients, double y, double xBase):PolyResidual(coefficients, y){this->baseVal = xBase;}; + PolyFracIntCentralResidual(const std::vector< std::vector > &coefficients, double x, double z, double yBase): PolyResidual(coefficients, x, z){this->baseVal = yBase;}; + virtual double call(double x); + virtual double deriv(double x); +}; class PolyDerResidual : public PolyResidual { public: PolyDerResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; diff --git a/src/Backends/Helmholtz/Fluids/Incompressible.cpp b/src/Backends/Helmholtz/Fluids/Incompressible.cpp index f2382558..edd50552 100644 --- a/src/Backends/Helmholtz/Fluids/Incompressible.cpp +++ b/src/Backends/Helmholtz/Fluids/Incompressible.cpp @@ -7,15 +7,124 @@ #include "Incompressible.h" + namespace CoolProp { -Incompressible::Incompressible() { - // TODO Auto-generated constructor stub - +std::vector > Incompressible::changeAxis(const std::vector &input){ + std::vector > output; + std::vector tmp; + size_t sizeX = input.size(); + for (size_t i = 0; i < sizeX; ++i){ + tmp.clear(); + tmp.push_back(input[i]); + output.push_back(tmp); + } + return output; } -Incompressible::~Incompressible() { - // TODO Auto-generated destructor stub +/* 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 index 24e6b759..79713bf8 100644 --- a/src/Backends/Helmholtz/Fluids/Incompressible.h +++ b/src/Backends/Helmholtz/Fluids/Incompressible.h @@ -8,14 +8,169 @@ #ifndef INCOMPRESSIBLE_H_ #define INCOMPRESSIBLE_H_ -//#include "AbstractFluid.h" +#include "PolyMath.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; + double x; + + std::vector< std::vector > cRho; + std::vector< std::vector > cHeat; + std::vector< std::vector > cVisc; + std::vector< std::vector > cCond; + std::vector< std::vector > cPsat; + std::vector cTfreeze; + + std::vector > changeAxis(const std::vector &input); + + BasePolynomial poly; + PolynomialSolver solver; + BaseExponential expo; + 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;} + double getx() const {return x;} + + 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;} + void setx(double x) {this->x = x;} + + // Setters for the coefficients + void setcRho(std::vector > cRho){this->cRho = cRho;} + void setcHeat(std::vector > cHeat){this->cHeat = cHeat;} + void setcVisc(std::vector > cVisc){this->cVisc = cVisc;} + void setcCond(std::vector > cCond){this->cCond = cCond;} + void setcPsat(std::vector > cPsat){this->cPsat = cPsat;} + void setcTfreeze(std::vector cTfreeze){this->cTfreeze = cTfreeze;} + + void setcRho(std::vector cRho){this->cRho = changeAxis(cRho);} + void setcHeat(std::vector cHeat){this->cHeat = changeAxis(cHeat);} + void setcVisc(std::vector cVisc){this->cVisc = changeAxis(cVisc);} + void setcCond(std::vector cCond){this->cCond = changeAxis(cCond);} + void setcPsat(std::vector cPsat){this->cPsat = changeAxis(cPsat);} + + 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); + /// Heat capacities as a function of temperature, pressure and composition. + virtual double c (double T_K, double p); + virtual double cp (double T_K, double p){return c(T_K,p);}; + virtual double cv (double T_K, double p){return c(T_K,p);}; + /// Entropy as a function of temperature, pressure and composition. + virtual double s (double T_K, double p); + /// Internal energy as a function of temperature, pressure and composition. + virtual double u (double T_K, double p); + /// Enthalpy as a function of temperature, pressure and composition. + virtual double h (double T_K, double p); + /// Viscosity as a function of temperature, pressure and composition. + virtual double visc(double T_K, double p); + /// Thermal conductivity as a function of temperature, pressure and composition. + virtual double cond(double T_K, double p); + /// Saturation pressure as a function of temperature and composition. + virtual double psat(double T_K ); + /// Freezing temperature as a function of pressure and composition. + virtual double Tfreeze( double p); + + +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) { + return u(T_K,p)+p/rho(T_K,p); + }; + + /// 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) { + return h(T_K,p)-p/rho(T_K,p); + }; + + + /* + * 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); + + /// 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); + + /// 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 */ diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp new file mode 100644 index 00000000..21b59c1c --- /dev/null +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -0,0 +1,98 @@ +#if defined(_MSC_VER) +#define _CRTDBG_MAP_ALLOC +#define _CRT_SECURE_NO_WARNINGS +#include +#include +#else +#include +#endif + +#include +//#include "CoolProp.h" + +#include "IncompressibleBackend.h" +#include "../Fluids/FluidLibrary.h" +#include "Solvers.h" +#include "MatrixMath.h" + +namespace CoolProp { + +IncompressibleBackend::IncompressibleBackend(const std::string &fluid_name) { + throw CoolProp::NotImplementedError("Mixture-style constructor is not implemented yet for incompressible fluids"); +} + +IncompressibleBackend::IncompressibleBackend(const std::vector &component_names) { + throw CoolProp::NotImplementedError("Mixture-style constructor is not implemented yet for incompressible fluids"); +} + +bool IncompressibleBackend::using_mole_fractions(){return false;} + + +void IncompressibleBackend::update(long input_pair, double value1, double value2) { + switch (input_pair) { + case PT_INPUTS: { + + break; + } +// case DmassP_INPUTS: { +// +// } +// break; +// } +// case HmassP_INPUTS: { +// // Call again, but this time with molar units +// // H: [J/kg] * [kg/mol] -> [J/mol] +// update(HmolarP_INPUTS, value1 * (double) _molar_mass, value2); +// return; +// } +// case PUmass_INPUTS: { +// // Call again, but this time with molar units +// // U: [J/kg] * [kg/mol] -> [J/mol] +// update(PUmolar_INPUTS, value1, value2 * (double) _molar_mass); +// return; +// } +// case PSmass_INPUTS: { +// // Call again, but this time with molar units +// // U: [J/kg] * [kg/mol] -> [J/mol] +// update(PUmolar_INPUTS, value1, value2 * (double) _molar_mass); +// return; +// } + default: { + throw ValueError( + format("This pair of inputs [%s] is not yet supported", + get_input_pair_short_desc(input_pair).c_str())); + } + } +} + +/// Set the mole fractions +/** +@param mole_fractions The vector of mole fractions of the components +*/ +void IncompressibleBackend::set_mole_fractions(const std::vector &mole_fractions) { + throw CoolProp::NotImplementedError("Cannot set mole fractions for incompressible fluid"); +} + +/// Set the mass fractions +/** +@param mass_fractions The vector of mass fractions of the components +*/ +void IncompressibleBackend::set_mass_fractions(const std::vector &mass_fractions) { + this->mass_fractions = mass_fractions; +} + +/// Check if the mole fractions have been set, etc. +void IncompressibleBackend::check_status() { + throw CoolProp::NotImplementedError("Cannot check status for incompressible fluid"); +} + +/// Get the viscosity [Pa-s] +long double IncompressibleBackend::calc_viscosity(void){ + throw NotImplementedError(); +} +/// Get the thermal conductivity [W/m/K] (based on the temperature and density in the state class) +long double IncompressibleBackend::calc_conductivity(void){ + throw NotImplementedError(); +} + +} diff --git a/src/Backends/Incompressible/IncompressibleBackend.h b/src/Backends/Incompressible/IncompressibleBackend.h index 4bd10263..c4f788fa 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.h +++ b/src/Backends/Incompressible/IncompressibleBackend.h @@ -20,44 +20,45 @@ public: virtual ~IncompressibleBackend(){}; /// The instantiator - /// @param fluid_names The vector of strings of the fluid components, without file ending - IncompressibleBackend(const std::string & fluid_names){ }; + /// @param fluid_name the string with the fluid name + IncompressibleBackend(const std::string &fluid_name); + /// The instantiator + /// @param fluid_names The vector of strings of the fluid components, without file ending + IncompressibleBackend(const std::vector &component_names); // Incompressible backend uses mole fractions - bool using_mole_fractions(){return false;} + bool using_mole_fractions(); /// Updating function for incompressible fluid /** In this function we take a pair of thermodynamic states, those defined in the input_pairs - enumeration and update all the internal variables that we can. REFPROP calculates - a lot of other state variables each time you use a flash routine so we cache all the - outputs that we can, which saves on computational time. + enumeration and update all the internal variables that we can. @param input_pair Integer key from CoolProp::input_pairs to the two inputs that will be passed to the function @param value1 First input value @param value2 Second input value */ - void update(long input_pair, double value1, double value2){throw CoolProp::NotImplementedError("Update not implemented yet for incompressible fluids");}; + void update(long input_pair, double value1, double value2); /// Set the mole fractions /** @param mole_fractions The vector of mole fractions of the components */ - void set_mole_fractions(const std::vector &mole_fractions){throw CoolProp::NotImplementedError("Cannot set mole fractions for incompressible fluid");}; + void set_mole_fractions(const std::vector &mole_fractions); /// Set the mass fractions /** @param mass_fractions The vector of mass fractions of the components */ - void set_mass_fractions(const std::vector &mass_fractions){this->mass_fractions = mass_fractions;}; + void set_mass_fractions(const std::vector &mass_fractions); /// Check if the mole fractions have been set, etc. void check_status(); /// Get the viscosity [Pa-s] - long double calc_viscosity(void){throw NotImplementedError();}; + long double calc_viscosity(void); /// Get the thermal conductivity [W/m/K] (based on the temperature and density in the state class) - long double calc_conductivity(void){throw NotImplementedError();}; + long double calc_conductivity(void); }; } /* namespace CoolProp */ diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index d9300d5f..e8015939 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -656,6 +656,55 @@ double PolyIntResidual::deriv(double x){ return polyRes; } +double PolyFracIntResidual::call(double x){ + double polyRes = -1; + switch (this->dim) { + case i1D: + polyRes = this->poly.polyfracint(this->coefficients[0], x); + break; + case i2D: + polyRes = this->poly.polyfracint(this->coefficients, this->firstDim, x); + break; + default: + throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); + } + return polyRes - this->output; +} + +double PolyFracIntResidual::deriv(double x){ + double polyRes = -1; + switch (this->dim) { + case i1D: + polyRes = this->poly.polyfracval(this->coefficients[0], x); + break; + case i2D: + polyRes = this->poly.polyfracval(this->coefficients, this->firstDim, x); + break; + default: + throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); + } + return polyRes; +} + +double PolyFracIntCentralResidual::call(double x){ + double polyRes = -1; + switch (this->dim) { + case i1D: + polyRes = this->poly.polyfracintcentral(this->coefficients[0], x, this->baseVal); + break; + case i2D: + polyRes = this->poly.polyfracintcentral(this->coefficients, this->firstDim, x, this->baseVal); + break; + default: + throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); + } + return polyRes - this->output; +} + +double PolyFracIntCentralResidual::deriv(double x){ + throw CoolProp::NotImplementedError("Derivative of a polynomial frac int is not defined."); +} + double PolyDerResidual::call(double x){ double polyRes = -1; switch (this->dim) { @@ -678,7 +727,7 @@ double PolyDerResidual::deriv(double x){ -/** Implements the same public functions as the +/** Implements the same public functions as the BasePolynomial * but solves the polynomial for the given value * instead of evaluating it. * TODO: This class does not check for bijective @@ -689,7 +738,7 @@ PolynomialSolver::PolynomialSolver(){ this->DEBUG = false; this->macheps = DBL_EPSILON; this->tol = DBL_EPSILON*1e3; - this->maxiter = 100; + this->maxiter = 50; } /** Everything related to the normal polynomials goes in this @@ -781,7 +830,8 @@ double PolynomialSolver::polyfracval(const std::vector< std::vector > &c /// @param coefficients vector containing the ordered coefficients /// @param y double value that represents the current output double PolynomialSolver::polyfracint(const std::vector &coefficients, double y){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function + PolyFracIntResidual residual = PolyFracIntResidual(coefficients, y); + return this->solve(residual); } /// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable @@ -789,7 +839,8 @@ double PolynomialSolver::polyfracint(const std::vector &coefficients, do /// @param x double value that represents the current input in the 1st dimension /// @param z double value that represents the current output double PolynomialSolver::polyfracint(const std::vector< std::vector > &coefficients, double x, double z){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function + PolyFracIntResidual residual = PolyFracIntResidual(coefficients, x, z); + return this->solve(residual); } /// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable @@ -797,7 +848,8 @@ double PolynomialSolver::polyfracint(const std::vector< std::vector > &c /// @param y double value that represents the current output /// @param xbase central x-value for fitted function double PolynomialSolver::polyfracintcentral(const std::vector &coefficients, double y, double xbase){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function + PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, y, xbase); + return this->solve(residual); } /// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable @@ -806,7 +858,8 @@ double PolynomialSolver::polyfracintcentral(const std::vector &coefficie /// @param z double value that represents the current output /// @param ybase central y-value for fitted function double PolynomialSolver::polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function + PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, x, z, ybase); + return this->solve(residual); } @@ -1024,6 +1077,30 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") CAPTURE(val1); CAPTURE(val2); CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyint(cHeat, T); + solver.setGuess(T+100); + val2 = solver.polyint(cHeat, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + +// val1 = poly.polyder(cHeat, T); +// solver.setGuess(T+100); +// val2 = solver.polyder(cHeat, val1); +// CAPTURE(T); +// CAPTURE(val1); +// CAPTURE(val2); +// CHECK(fabs(T-val2) < 1e-1); +// +// val1 = poly.polyfracint(cHeat, T); +// solver.setGuess(T+100); +// val2 = solver.polyfracint(cHeat, val1); +// CAPTURE(T); +// CAPTURE(val1); +// CAPTURE(val2); +// CHECK(fabs(T-val2) < 1e-1); } } SECTION("Solve1DBrent") { @@ -1035,6 +1112,39 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") CAPTURE(val1); CAPTURE(val2); CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyint(cHeat, T); + solver.setLimits(T-300,T+300); + val2 = solver.polyint(cHeat, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyder(cHeat, T); + solver.setLimits(T-300,T+300); + val2 = solver.polyder(cHeat, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyfracint(cHeat, T); + solver.setLimits(T-100,T+100); + val2 = solver.polyfracint(cHeat, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyfracintcentral(cHeat, T, 250.0); + solver.setLimits(T-100,T+100); + val2 = solver.polyfracintcentral(cHeat, val1, 250.0); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + } } @@ -1094,6 +1204,38 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") CAPTURE(val1); CAPTURE(val2); CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyint(cHeat2D, xDim1, T); + solver.setLimits(T-300,T+300); + val2 = solver.polyint(cHeat2D, xDim1, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyder(cHeat2D, xDim1, T); + solver.setLimits(T-300,T+300); + val2 = solver.polyder(cHeat2D, xDim1, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyfracint(cHeat2D, xDim1, T); + solver.setLimits(T-100,T+100); + val2 = solver.polyfracint(cHeat2D, xDim1, val1); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); + + val1 = poly.polyfracintcentral(cHeat2D, xDim1, T, 250); + solver.setLimits(T-100,T+100); + val2 = solver.polyfracintcentral(cHeat2D, xDim1, val1, 250); + CAPTURE(T); + CAPTURE(val1); + CAPTURE(val2); + CHECK(fabs(T-val2) < 1e-1); } }