diff --git a/Web/conf.py b/Web/conf.py index 9df62061..7c7af3e8 100644 --- a/Web/conf.py +++ b/Web/conf.py @@ -26,7 +26,7 @@ except ImportError: print('Unable to import sphinxcontrib.doxylink; try to run "pip install sphinxcontrib-doxylink"') -Release = tags.has('Release') or tags.has('Release') +Release = tags.has('Release') or tags.has('release') if Release: doxylink = { @@ -34,7 +34,7 @@ if Release: } else: doxylink = { - 'cpapi' : ('_static/doxygen/CoolPropDoxyLink.tag', 'http://www.coolprop.dreamhosters.com:8010/sphinx/_static/doxygen/html') + 'cpapi' : ('_static/doxygen/CoolPropDoxyLink.tag', 'http://www.coolprop.dreamhosters.com:8010/binaries/sphinx/_static/doxygen/html') } # -- General configuration ----------------------------------------------------- diff --git a/Web/coolprop/wrappers/SharedLibrary/index.rst b/Web/coolprop/wrappers/SharedLibrary/index.rst index 29dbc45d..610fe32c 100644 --- a/Web/coolprop/wrappers/SharedLibrary/index.rst +++ b/Web/coolprop/wrappers/SharedLibrary/index.rst @@ -106,7 +106,7 @@ For 32-bit compilation:: # Make a build folder mkdir build && cd build # Generate builder - cmake .. -DCOOLPROP_32BIT_STDCALL_SHARED_LIBRARY=ON + cmake .. -DCOOLPROP_32BIT_SHARED_LIBRARY_LINUX_MODULE=ON # Build cmake --build . @@ -122,3 +122,12 @@ For 64-bit compilation:: cmake .. -DCOOLPROP_64BIT_SHARED_LIBRARY=ON # Build cmake --build . + +On Linux, installation could be done by:: + + # Change "32" to match your system bitness + sudo cp libCoolProp.so /usr/local/lib/libCoolProp.so.32.:version: + pushd /usr/local/lib + sudo ln -sf libCoolProp.so.32.:version: libCoolProp.so.5 + sudo ln -sf libCoolProp.so.5 libCoolProp.so + popd diff --git a/Web/develop/code.rst b/Web/develop/code.rst index e0be2bc6..cd4b6540 100644 --- a/Web/develop/code.rst +++ b/Web/develop/code.rst @@ -10,7 +10,8 @@ The source code of CoolProp is stored in a github repository at https://github.c Doxygen formatted documentation of the source files --------------------------------------------------- -Real-time builds of the `doxygen `_ formatted HTML outputs for the development code are at :bbsphinx:`the buildbot development server`. +Builds of the `doxygen `_ formatted HTML outputs for the +current version are integrated this page, :cpapi:`CoolProp`. More information ---------------- diff --git a/Web/develop/index.rst b/Web/develop/index.rst index 82593ee2..c597b9e2 100644 --- a/Web/develop/index.rst +++ b/Web/develop/index.rst @@ -11,6 +11,7 @@ Information for Developers cmake.rst buildbot.rst documentation.rst + release.rst Address Sanitizer ----------------- diff --git a/Web/develop/release.rst b/Web/develop/release.rst new file mode 100644 index 00000000..0c0cddde --- /dev/null +++ b/Web/develop/release.rst @@ -0,0 +1,25 @@ +.. _release: + +****************** +Release Checklist +****************** + +We have made a serious effort to automate the release of new binaries. Even +though things have become easier, there are still many things to remember. +Here is your new best friend, a checklist that helps you to keep track of all +the small things that have to be done when releasing a new version of the CoolProp +library. + +* **Version**: Edit CMakeLists.txt and remove all qualifiers (alpha, dev, ...) from the version number. +* **Changelog**: Update the changelog and generate a list of closed GitHub issues: *HOW?* +* **release branch**: Merge all code from *master* into *release* branch +* **build bots**: Force all buildbots to run on the *release* branch, this will also change the upload folder from *binaries* to *release*. +* **script**: Wait for all bots to finish and run the release script by launching the ``release version`` bot with dry run disabled and the correct version number. This uploads binaries to pypi and sourceforge. +* **clean up**: If everything went well, you can proceed: + - Tag the release branch in GitHub. It is a good idea to include the information on the closed issues here as well. + - Change the default download file on sourceforge to point to the new zipped sources. + - Copy the new Javascript library to the homepage and make a symlink to ``coolprop-latest.js``. *I think I automated this one already* + - Bump the version number in the CMake file and commit. + - Announce the new features if you like... + +That's all folks. \ No newline at end of file diff --git a/Web/fluid_properties/Incompressibles.rst b/Web/fluid_properties/Incompressibles.rst index facc9c46..78bbdcd2 100644 --- a/Web/fluid_properties/Incompressibles.rst +++ b/Web/fluid_properties/Incompressibles.rst @@ -191,67 +191,71 @@ compositions as independent fluids. This should be kept in mind when comparing properties for different compositions. Setting the reference state for one composition will always affect all fluids consisting of the same components. -The approach described in textbooks like Cengel and Boles :cite:`Cengel2007` -is that the internal energy :math:`u` only depends on temperature and does not -change with pressure. - -.. Alternatively, use cancel package with \cancelto{0}{x-d} command - -.. math:: - - du &= \overbrace{ \left( \frac{\partial u}{\partial T} \right)_p}^{=c_p=c_v=c} dT &+ \overbrace{\left( \frac{\partial u}{\partial p} \right)_T}^{\stackrel{!}{=}0} dp \\ - -By using the fourth Maxwell relation, we can extend the simplifications to the -entropy formulation - -.. math:: - - ds &= \left( \frac{\partial s}{\partial T} \right)_p dT &+ \left( \frac{\partial s}{\partial p} \right)_T dp \\ - &= \underbrace{ \left( \frac{\partial h}{\partial T} \right)_p}_{=c_p=c_v=c} T^{-1} dT &-\underbrace{\left( \frac{\partial v}{\partial T} \right)_p}_{\stackrel{!}{=} 0} dp \\ - -As indicated by the braces above, the fluids implemented in CoolProp do also follow -the second common assumption of a constant specific volume :math:`v` that does -change neither with temperature nor with pressure. It should be highlighted, that -this simplification violates the integrity of the implemented equations since there -are changes in density as a function of temperature for all incompressible fluids. - -Employing :math:`h=u+pv`, we can derive the impact on enthalpy as well by -rewriting the equation in terms of our state variables :math:`p` and :math:`T` -as shown by Skovrup :cite:`Skovrup1999`. - -.. dh &= \overbrace{ \left( \frac{\partial h}{\partial T} \right)_p}^{=c_p=c_v=c} dT + \left( \frac{\partial h}{\partial p} \right)_T dp \\ - -.. math:: - dh &= \overbrace{ \left( \frac{\partial h}{\partial T} \right)_p}^{=c_p=c_v=c} dT + \left( \frac{\partial h}{\partial p} \right)_T dp \\ - &= \left( \frac{\partial u}{\partial T} \right)_v dT + \left( v - T \left( \frac{\partial v}{\partial T} \right)_p \right) dp \\ - &= du + \underbrace{p dv}_{\stackrel{!}{=} 0} + v dp \quad \text{ with $v\stackrel{!}{=}v_0=$ const } \\ - -The two assumptions used above :math:`\left( \partial v / \partial T \right)_p \stackrel{!}{=} 0` -and :math:`\left( \partial u / \partial T \right)_p \stackrel{!}{=} \left( \partial u / \partial T \right)_v` -imply that :math:`v` is constant under all circumstances. Hence, we have to use -the specific volume at reference conditions to calculate enthalpy from the -integration in :math:`T` and :math:`p`. Future work could provide a more accurate -formulation of entropy and enthalpy by implementing the term -:math:`\left( \partial v / \partial T \right)_p \neq 0`. - -Using only polynomials for the heat capacity functions, we can derive internal -energy and entropy by integrating the specific heat capacity in temperature. +.. The approach described in textbooks like Cengel and Boles :cite:`Cengel2007` +.. is that the internal energy :math:`u` only depends on temperature and does not +.. change with pressure. +.. +.. .. Alternatively, use cancel package with \cancelto{0}{x-d} command +.. +.. .. math:: +.. +.. du &= \overbrace{ \left( \frac{\partial u}{\partial T} \right)_p}^{=c_p=c_v=c} dT &+ \overbrace{\left( \frac{\partial u}{\partial p} \right)_T}^{\stackrel{!}{=}0} dp \\ +.. +.. By using the fourth Maxwell relation, we can extend the simplifications to the +.. entropy formulation +.. +.. .. math:: +.. +.. ds &= \left( \frac{\partial s}{\partial T} \right)_p dT &+ \left( \frac{\partial s}{\partial p} \right)_T dp \\ +.. &= \underbrace{ \left( \frac{\partial h}{\partial T} \right)_p}_{=c_p=c_v=c} T^{-1} dT &-\underbrace{\left( \frac{\partial v}{\partial T} \right)_p}_{\stackrel{!}{=} 0} dp \\ +.. +.. As indicated by the braces above, the fluids implemented in CoolProp do also follow +.. the second common assumption of a constant specific volume :math:`v` that does +.. change neither with temperature nor with pressure. It should be highlighted, that +.. this simplification violates the integrity of the implemented equations since there +.. are changes in density as a function of temperature for all incompressible fluids. +.. +.. Employing :math:`h=u+pv`, we can derive the impact on enthalpy as well by +.. rewriting the equation in terms of our state variables :math:`p` and :math:`T` +.. as shown by Skovrup :cite:`Skovrup1999`. +.. +.. .. dh &= \overbrace{ \left( \frac{\partial h}{\partial T} \right)_p}^{=c_p=c_v=c} dT + \left( \frac{\partial h}{\partial p} \right)_T dp \\ +.. +.. .. math:: +.. dh &= \overbrace{ \left( \frac{\partial h}{\partial T} \right)_p}^{=c_p=c_v=c} dT + \left( \frac{\partial h}{\partial p} \right)_T dp \\ +.. &= \left( \frac{\partial u}{\partial T} \right)_v dT + \left( v - T \left( \frac{\partial v}{\partial T} \right)_p \right) dp \\ +.. &= du + \underbrace{p dv}_{\stackrel{!}{=} 0} + v dp \quad \text{ with $v\stackrel{!}{=}v_0=$ const } \\ +.. +.. The two assumptions used above :math:`\left( \partial v / \partial T \right)_p \stackrel{!}{=} 0` +.. and :math:`\left( \partial u / \partial T \right)_p \stackrel{!}{=} \left( \partial u / \partial T \right)_v` +.. imply that :math:`v` is constant under all circumstances. Hence, we have to use +.. the specific volume at reference conditions to calculate enthalpy from the +.. integration in :math:`T` and :math:`p`. Future work could provide a more accurate +.. formulation of entropy and enthalpy by implementing the term +.. :math:`\left( \partial v / \partial T \right)_p \neq 0`. +.. +.. Using only polynomials for the heat capacity functions, we can derive internal +.. energy and entropy by integrating the specific heat capacity in temperature. .. _BaseValue: -.. math:: - - c &= \sum_{i=0}^n x^i \cdot \sum_{j=0}^m C_{c}[i,j] \cdot T^j \text{ yielding } \\ - u &= \int_{0}^{1} c\left( x,T \right) dT - = \sum_{i=0}^n x^i \cdot \sum_{j=0}^m \frac{1}{j+1} \cdot C_{c}[i,j] - \cdot \left( T_{1}^{j+1} - T_{0}^{j+1} \right) \text{ and } \\ - s &= \int_{0}^{1} \frac{c\left( x,T \right)}{T} dT - = \sum_{i=0}^n x^i \cdot \left( - C_{c}[i,0] \cdot \ln\left(\frac{T_{1}}{T_{0}}\right) - + \sum_{j=0}^{m-1} \frac{1}{j+1} \cdot C_{c}[i,j+1] \cdot \left( T_{1}^{j+1} - T_{0}^{j+1} \right) - \right) \\ - h &= u + v_{0} \cdot \left( p_{1} - p_{0} \right) +.. note:: + The internal routines for the incompressibles were updated 2015-02-10, the documentation is not fully updated. + We are going to add the new equation as soon as possible, probably mid-March 2015. Please be patient. +.. .. math:: +.. +.. c &= \sum_{i=0}^n x^i \cdot \sum_{j=0}^m C_{c}[i,j] \cdot T^j \text{ yielding } \\ +.. u &= \int_{0}^{1} c\left( x,T \right) dT +.. = \sum_{i=0}^n x^i \cdot \sum_{j=0}^m \frac{1}{j+1} \cdot C_{c}[i,j] +.. \cdot \left( T_{1}^{j+1} - T_{0}^{j+1} \right) \text{ and } \\ +.. s &= \int_{0}^{1} \frac{c\left( x,T \right)}{T} dT +.. = \sum_{i=0}^n x^i \cdot \left( +.. C_{c}[i,0] \cdot \ln\left(\frac{T_{1}}{T_{0}}\right) +.. + \sum_{j=0}^{m-1} \frac{1}{j+1} \cdot C_{c}[i,j+1] \cdot \left( T_{1}^{j+1} - T_{0}^{j+1} \right) +.. \right) \\ +.. h &= u + v_{0} \cdot \left( p_{1} - p_{0} \right) +.. According to Melinder :cite:`Melinder2010` and Skovrup :cite:`Skovrup2013`, using a centred approach for the independent variables enhances the fit quality. @@ -278,10 +282,9 @@ be multiplied with the other coefficients and the concentration. .. math:: - s &= \int_{0}^{1} \frac{c\left( x_\text{in},T_\text{in} \right)}{T_\text{in}} dT_\text{in} = \sum_{i=0}^n x_\text{in}^i \cdot \sum_{j=0}^m C_{c}[i,j] \cdot F(j,T_\text{in,0},T_\text{in,1}) \\ + \int_{0}^{1} \left( \frac{\partial s}{\partial T} \right)_p dT &= \int_{0}^{1} \frac{c\left( x_\text{in},T_\text{in} \right)}{T_\text{in}} dT_\text{in} = \sum_{i=0}^n x_\text{in}^i \cdot \sum_{j=0}^m C_{c}[i,j] \cdot F(j,T_\text{in,0},T_\text{in,1}) \\ F &= (-1)^j \cdot \ln \left( \frac{T_\text{in,1}}{T_\text{in,0}} \right) \cdot T_{base}^j + \sum_{k=0}^{j-1} \binom{j}{k} \cdot \frac{(-1)^k}{j-k} \cdot \left( T_\text{in,1}^{j-k} - T_\text{in,0}^{j-k} \right) \cdot T_{base}^k - .. _Equations: Equations diff --git a/dev/scripts/release.bsh b/dev/scripts/release.bsh index b1a0a88b..a89972da 100755 --- a/dev/scripts/release.bsh +++ b/dev/scripts/release.bsh @@ -142,6 +142,12 @@ rsync $RSYNC_DRY_RUN $RSYNC_OPTS $RSYNC_EXCL "$BINFOLDER/" frs.sf.net-$SFUSER:/h if [ ${CPVERSION:0:7} != "nightly" ]; then printMessage "Publishing the docs on SourceForge" rsync $RSYNC_DRY_RUN $RSYNC_OPTS "$SPHFOLDER/" frs.sf.net-$SFUSER:/home/project-web/coolprop/htdocs + printMessage "Updating the default Javascript library on SourceForge" + pushd "$BINFOLDER/Javascript/" + cp "coolprop.js" "coolprop-${CPVERSION}.js" + ln -sf "coolprop-${CPVERSION}.js" "coolprop-latest.js" + rsync $RSYNC_DRY_RUN $RSYNC_OPTS "coolprop-${CPVERSION}.js" "coolprop-latest.js" frs.sf.net-$SFUSER:/home/project-web/coolprop/htdocs/jscript/ + popd fi # if [[ ("$BINFOLDER" == "release") && ("$DRYRUN" == "false") ]]; then diff --git a/include/IncompressibleFluid.h b/include/IncompressibleFluid.h index 5445b91f..61d29e05 100644 --- a/include/IncompressibleFluid.h +++ b/include/IncompressibleFluid.h @@ -60,15 +60,9 @@ protected: composition_types xid; double TminPsat; - double xref, Tref, pref; - double href, sref; - double uref, rhoref; + double Tref, pref, xref, href, sref; double xbase, Tbase; - double _xref, _Tref, _pref; - double _href, _sref; - double _uref, _rhoref; - /// These are the objects that hold the coefficients /** Note that all polynomials require a 2-dimensional array * of coefficients. This array may have only one row or @@ -201,6 +195,7 @@ protected: double basePolyOffset(IncompressibleData data, double y, double z=0.0); public: + /* All functions need T and p as input. Might not * be necessary, but gives a clearer structure. */ @@ -225,6 +220,31 @@ public: /// Freezing temperature as a function of pressure and composition. double Tfreeze( double p, double x); + + /* Below are direct calculations of the derivatives. Nothing + * special is going on, we simply use the polynomial class to + * derive the different functions with respect to temperature. + */ + /// Partial derivative of density with respect to temperature at constant pressure and composition + double drhodTatPx (double T, double p, double x); + /// Partial derivative of entropy with respect to temperature at constant pressure and composition + double dsdTatPx (double T, double p, double x){return c(T,p,x)/T;}; + /// Partial derivative of internal energy with respect to temperature at constant pressure and composition + double dudTatPx (double T, double p, double x){return c(T,p,x);}; + /// Partial derivative of enthalpy with respect to temperature at constant pressure and composition + double dhdTatPx (double T, double p, double x){return c(T,p,x);}; + + + /* Other useful derivatives + */ + /// Partial derivative of enthalpy with respect to pressure at constant temperature and composition + // \f[ \left( \frac{\partial h}{\partial p} \right)_T = v - T \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-1} \left( 1 + T \rho^{-1} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] + double dhdpatTx (double T, double p, double x); + /// Partial derivative of entropy with respect to pressure at constant temperature and composition + // \f[ \left( \frac{\partial s}{\partial p} \right)_T = - \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-2} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] + double dsdpatTx (double T, double p, double x); + + /// Mass fraction conversion function /** If the fluid type is mass-based, it does not do anything. Otherwise, * it converts the mass fraction to the required input. */ @@ -239,7 +259,6 @@ public: double inputFromMole (double T, double x); - /* Some functions can be inverted directly, those are listed * here. It is also possible to solve for other quantities, but * that involves some more sophisticated processing and is not @@ -250,9 +269,9 @@ public: /// Temperature as a function of heat capacities as a function of temperature, pressure and composition. double T_c (double Cmass, double p, double x); /// Temperature as a function of entropy as a function of temperature, pressure and composition. - double T_s (double Smass, double p, double x); + double T_s (double Smass, double p, double x){throw NotImplementedError(format("%s (%d): T from entropy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} /// Temperature as a function of internal energy as a function of temperature, pressure and composition. - double T_u (double Umass, double p, double x); + double T_u (double Umass, double p, double x){throw NotImplementedError(format("%s (%d): T from internal energy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} /// Temperature as a function of enthalpy, pressure and composition. double T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} /// Viscosity as a function of temperature, pressure and composition. @@ -276,7 +295,7 @@ protected: * 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-_pref)/_rhoref; + return u(T,p,x)+p/rho(T,p,x); }; /// Internal energy from h, p and rho. @@ -284,7 +303,7 @@ protected: * 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-_pref)/_rhoref; + return h(T,p,x)-p/rho(T,p,x); }; diff --git a/include/PolyMath.h b/include/PolyMath.h index 79fc500d..20edc08d 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -288,7 +288,8 @@ public: /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension - double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0); + /// @param ax_val double value that represents the base value for the definite integral on the chosen axis. + double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0, const double &ax_val=0.0); public: /// Returns a vector with ALL the real roots of p(x_in,y_in)-z_in diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp index e353c1b9..0531eaeb 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.cpp +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -90,11 +90,11 @@ void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double valu _T = this->DmassP_flash(value1,value2); break; } - case PUmass_INPUTS: { - _p = value1; - _T = this->PUmass_flash(value1, value2); - break; - } +// case PUmass_INPUTS: { +// _p = value1; +// _T = this->PUmass_flash(value1, value2); +// break; +// } case PSmass_INPUTS: { _p = value1; _T = this->PSmass_flash(value1, value2); @@ -264,7 +264,35 @@ long double IncompressibleBackend::HmassP_flash(long double hmass, long double p @returns T The temperature in K */ long double IncompressibleBackend::PSmass_flash(long double p, long double smass){ - return fluid->T_s(smass, p, _fractions[0]); + + class PSmass_residual : public FuncWrapper1D { + protected: + double p,x,s_in; + IncompressibleFluid* fluid; + protected: + PSmass_residual(){}; + public: + PSmass_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &s_in){ + this->p = p; + this->x = x; + this->s_in = s_in; + this->fluid = fluid; + } + virtual ~PSmass_residual(){}; + double call(double target){ + return fluid->s(target,p,x) - s_in; + } + }; + + PSmass_residual res = PSmass_residual(fluid, p, _fractions[0], smass); + + std::string errstring; + double macheps = DBL_EPSILON; + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Brent(&res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring); + //if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; + return result; } /// Calculate T given pressure and internal energy @@ -274,7 +302,35 @@ long double IncompressibleBackend::PSmass_flash(long double p, long double smass @returns T The temperature in K */ long double IncompressibleBackend::PUmass_flash(long double p, long double umass){ - return fluid->T_u(umass, p, _fractions[0]); + + class PUmass_residual : public FuncWrapper1D { + protected: + double p,x,u_in; + IncompressibleFluid* fluid; + protected: + PUmass_residual(){}; + public: + PUmass_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &u_in){ + this->p = p; + this->x = x; + this->u_in = u_in; + this->fluid = fluid; + } + virtual ~PUmass_residual(){}; + double call(double target){ + return fluid->u(target,p,x) - u_in; + } + }; + + PUmass_residual res = PUmass_residual(fluid, p, _fractions[0], umass); + + std::string errstring; + double macheps = DBL_EPSILON; + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Brent(&res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring); + //if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; + return result; } @@ -436,7 +492,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CAPTURE(res); CHECK( check_abs(val,res,acc) ); } - // Make sure the ruslt does not change -> reference state... + // Make sure the result does not change -> reference state... val = backend.calc_smass(); backend.update( CoolProp::PT_INPUTS, p, T ); res = backend.calc_smass(); diff --git a/src/Backends/Incompressible/IncompressibleFluid.cpp b/src/Backends/Incompressible/IncompressibleFluid.cpp index 5e7a1ea8..9d5b5b8d 100644 --- a/src/Backends/Incompressible/IncompressibleFluid.cpp +++ b/src/Backends/Incompressible/IncompressibleFluid.cpp @@ -19,23 +19,11 @@ and transport properties. //IncompressibleFluid::IncompressibleFluid(); void IncompressibleFluid::set_reference_state(double T0, double p0, double x0, double h0, double s0){ - this->_href = h0; - this->_pref = p0; - this->_Tref = T0; - this->_xref = x0; - this->_rhoref = rho(T0,p0,x0); - this->_uref = h0; // Reset the old reference value - this->_uref = u(T0,p0,x0); // uses _uref - this->_sref = s0; // Reset the old reference value - this->_sref = s(T0,p0,x0); // uses _sref - // Save the values for later calls at composition changes - this->href = h0; - this->pref = p0; - this->Tref = T0; - this->xref = x0; - this->rhoref = this->_rhoref; - this->uref = this->_uref; - this->sref = s0; + this->Tref = T0; + this->pref = p0; + this->xref = x0; + this->href = h0; + this->sref = s0; } void IncompressibleFluid::validate(){ @@ -88,6 +76,7 @@ double IncompressibleFluid::basePolyOffset(IncompressibleData data, double y, do return _HUGE; } + /// Density as a function of temperature, pressure and composition. double IncompressibleFluid::rho (double T, double p, double x){ switch (density.type) { @@ -135,10 +124,15 @@ double IncompressibleFluid::c (double T, double p, double x){ /// Entropy as a function of temperature, pressure and composition. double IncompressibleFluid::s (double T, double p, double x){ + double dsdTatp = _HUGE; + double dsdpatT = _HUGE; switch (specific_heat.type) { case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - //throw NotImplementedError("Here you should implement the polynomial."); - return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - _sref; + //TODO: Optimise this, add a cached value or direct calculation of definite integral + dsdTatp = poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - poly.integral(specific_heat.coeffs, Tref, xref, 0, -1, 0, Tbase, xbase); + dsdpatT = this->rho(T,p,x); + dsdpatT = this->drhodTatPx(T,p,x)/dsdpatT/dsdpatT * (p-pref); + return sref + dsdTatp + dsdpatT; break; case IncompressibleData::INCOMPRESSIBLE_NOT_SET: throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); @@ -151,24 +145,28 @@ double IncompressibleFluid::s (double T, double p, double x){ } /// Internal energy as a function of temperature, pressure and composition. -double IncompressibleFluid::u (double T, double p, double x){ - switch (specific_heat.type) { - case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - //throw NotImplementedError("Here you should implement the polynomial."); - return poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - _uref; - break; - case IncompressibleData::INCOMPRESSIBLE_NOT_SET: - throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); - break; - default: - throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for internal energy.",__FILE__,__LINE__,specific_heat.type)); - break; - } - return _HUGE; -} +double IncompressibleFluid::u (double T, double p, double x){return u_h(T,p,x);}; /// Enthalpy as a function of temperature, pressure and composition. -double IncompressibleFluid::h (double T, double p, double x){return h_u(T,p,x);}; +double IncompressibleFluid::h (double T, double p, double x){ + double dhdTatP = _HUGE; + double dhdpatT = _HUGE; + switch (specific_heat.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + //TODO: Optimise this, add a cached value or direct calculation of definite integral + dhdTatP = poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - poly.integral(specific_heat.coeffs, Tref, xref, 0, 0, 0, Tbase, xbase); + dhdpatT = this->dhdpatTx(T, p, x) * (p-pref); + return href + dhdTatP + dhdpatT; + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for internal energy.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; + } /// Viscosity as a function of temperature, pressure and composition. double IncompressibleFluid::visc(double T, double p, double x){ @@ -280,6 +278,43 @@ double IncompressibleFluid::Tfreeze( double p, double x){ return _HUGE; } + +/* Below are direct calculations of the derivatives. Nothing + * special is going on, we simply use the polynomial class to + * derive the different functions with respect to temperature. + */ +/// Partial derivative of density with respect to temperature at constant pressure and composition +double IncompressibleFluid::drhodTatPx (double T, double p, double x){ + switch (density.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.derivative(density.coeffs, T, x, 0, 0, 0, Tbase, xbase); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,density.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for density.",__FILE__,__LINE__,density.type)); + break; + } + return _HUGE; +} + +/* Other useful derivatives + */ +/// Partial derivative of enthalpy with respect to temperature at constant pressure and composition +double IncompressibleFluid::dhdpatTx (double T, double p, double x){ + double rho = 0; + rho = this->rho(T,p,x); + return 1/rho * ( 1 + T/rho * this->drhodTatPx(T,p,x)); +} +/// Partial derivative of entropy with respect to temperature at constant pressure and composition +double IncompressibleFluid::dsdpatTx (double T, double p, double x){ + double rho = 0; + rho = this->rho(T,p,x); + return 1/rho/rho * this->drhodTatPx(T,p,x); +} + + /// Mass fraction conversion function /** If the fluid type is mass-based, it does not do anything. Otherwise, * it converts the mass fraction to the required input. */ @@ -292,7 +327,7 @@ double IncompressibleFluid::inputFromMass (double T, double x){ throw NotImplementedError("Mass composition conversion has not been implemented."); switch (mass2input.type) { case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - return poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way + return poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way break; case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL: return baseExponential(mass2input, x, 0.0); @@ -301,7 +336,7 @@ double IncompressibleFluid::inputFromMass (double T, double x){ return baseLogexponential(mass2input, x, 0.0); break; case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL: - return exp(poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way + return exp(poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way break; case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET: return basePolyOffset(mass2input, T, x); @@ -329,7 +364,7 @@ double IncompressibleFluid::inputFromVolume (double T, double x){ throw NotImplementedError("Volume composition conversion has not been implemented."); switch (volume2input.type) { case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - return poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way + return poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way break; case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL: return baseExponential(volume2input, x, 0.0); @@ -338,7 +373,7 @@ double IncompressibleFluid::inputFromVolume (double T, double x){ return baseLogexponential(volume2input, x, 0.0); break; case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL: - return exp(poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way + return exp(poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way break; case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET: return basePolyOffset(volume2input, T, x); @@ -366,7 +401,7 @@ double IncompressibleFluid::inputFromMole (double T, double x){ throw NotImplementedError("Mole composition conversion has not been implemented."); switch (mole2input.type) { case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way + return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way break; case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL: return baseExponential(mole2input, x, 0.0); @@ -375,7 +410,7 @@ double IncompressibleFluid::inputFromMole (double T, double x){ return baseLogexponential(mole2input, x, 0.0); break; case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL: - return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way + return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way break; case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET: return basePolyOffset(mole2input, T, x); @@ -428,49 +463,6 @@ double IncompressibleFluid::T_c (double Cmass, double p, double x){ } return _HUGE; } -/// Temperature as a function of entropy as a function of temperature, pressure and composition. -double IncompressibleFluid::T_s (double Smass, double p, double x){ - double s_raw = Smass + _sref; - switch (specific_heat.type) { - case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - return poly.solve_limitsInt(specific_heat.coeffs, x, s_raw, Tmin, Tmax, 0, -1, 0, Tbase, xbase, 0); - break; - case IncompressibleData::INCOMPRESSIBLE_NOT_SET: - throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); - break; - default: - throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type)); - break; - } - return _HUGE; -} -/// Temperature as a function of internal energy as a function of temperature, pressure and composition. -double IncompressibleFluid::T_u (double Umass, double p, double x){ - double u_raw = Umass + _uref; - switch (specific_heat.type) { - case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - return poly.solve_limitsInt(specific_heat.coeffs, x, u_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase, 0); - break; - case IncompressibleData::INCOMPRESSIBLE_NOT_SET: - throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); - break; - default: - throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type)); - break; - } - return _HUGE; -} -///// Temperature as a function of enthalpy, pressure and composition. -//double IncompressibleFluid::T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} -///// Viscosity as a function of temperature, pressure and composition. -//double IncompressibleFluid::T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));} -///// Thermal conductivity as a function of temperature, pressure and composition. -//double IncompressibleFluid::T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));} -///// Saturation pressure as a function of temperature and composition. -//double IncompressibleFluid::T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));} -///// Composition as a function of freezing temperature and pressure. -//double IncompressibleFluid::x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));} - /* * Some more functions to provide a single implementation @@ -737,7 +729,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } // Compare s - val = 145.59157247249246; + val = 144.08; res = XLT.s(T,p,x); { CAPTURE(T); @@ -756,7 +748,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } // Compare u - val = 45212.407309106304; + val = 44724.1; res = XLT.u(T,p,x); { CAPTURE(T); @@ -775,7 +767,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } // Compare h - val = 46388.7; + val = 45937; res = XLT.h(T,p,x); { CAPTURE(T); @@ -863,7 +855,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } // Compare s - expected = -206.62646783739274; + expected = -207.027; actual = CH3OH.s(T,p,x); { CAPTURE(T); @@ -904,7 +896,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } // Compare u - expected = -60043.78429641827; + expected = -60157.1; actual = CH3OH.u(T,p,x); { CAPTURE(T); @@ -927,7 +919,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } // Compare h - expected = -58999.1; + expected = -59119; actual = CH3OH.h(T,p,x); { CAPTURE(T); diff --git a/src/Backends/Incompressible/IncompressibleLibrary.cpp b/src/Backends/Incompressible/IncompressibleLibrary.cpp index b4ae37e9..8fc1a6a0 100644 --- a/src/Backends/Incompressible/IncompressibleLibrary.cpp +++ b/src/Backends/Incompressible/IncompressibleLibrary.cpp @@ -288,8 +288,6 @@ LiBrSolution::LiBrSolution():IncompressibleFluid(){ pref = 0.0; href = 0.0; sref = 0.0; - uref = 0.0; - rhoref = 0.0; xbase = 0.0; Tbase = 0.0; diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 7d667079..6633db9f 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -596,7 +596,8 @@ double Polynomial2DFrac::derivative(const Eigen::MatrixXd &coefficients, const d /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension -double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){ +/// @param ax_val double value that represents the base value for the definite integral on the chosen axis. +double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, const double &ax_val){ Eigen::MatrixXd newCoefficients; int int_exp,other_exp; @@ -628,6 +629,8 @@ double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const dou } 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)); + // TODO: Fix this and allow the direct calculation of integrals + if (std::abs(ax_val)>DBL_EPSILON) throw NotImplementedError(format("%s (%d): This function is only implemented for indefinite integrals, ax_val=%d ",__FILE__,__LINE__,ax_val)); size_t r = newCoefficients.rows(); size_t c = newCoefficients.cols();