A commit for the REFPROP backend to include the functions of satu… by dongkeun-oh · Pull Request #2016 · CoolProp/CoolProp

dongkeun-oh · 2021-04-13T05:07:06Z

Purpose

As a follow up of the issue #2013, I tried to include the missing functions which is particularly related to the saturated states.
They are critical to associated the REFPROP with the recently updated ExternalMedia library (for the Modelica platform).

calc_hmolar
calc_smolar
calc_first_saturation_deriv
calc_first_two_phase_deriv
calc_first_two_phase_deriv_splined (not tested!)

Description of the Change

protected part:

Instead of pointers of saturated state, as it is different from the HEOS backend, some additional cached variables
(CachedElement) are declared.

// Cached values of saturated state
    CachedElement _hLmolar, _hVmolar, _sLmolar, _sVmolar,_uLmolar, _uVmolar;
    CachedElement _wL, _wV, _cvLmolar, _cvVmolar, _cpLmolar, _cpVmolar;
    CachedElement _viscosityL, _viscosityV, _conductivityL, _conductivityV; //transport properties of saturated state

Whenever the saturated states (new cached variables) are referred by _calc_saturated_liquid_output( CachedElement&) or _calc_saturated_vapor_output( CachedElement&), the variable is checked whether it is already cached or not.
If not cached yet, the saturated (liquid or vapor) state will be updated by the current temperature _T and the saturated density _rhoLmolar (liquid) or _rhoVmolar (vapor) using the THERMdll subroutine in the REFPROP DLL.
To make sure the validity of saturated state, i.e. whether it is "two-phase" or not, an inline function _RPcheckTwophase() is introduced.
inline bool _RPcheckTwophase( void) { if (ValidNumber(_Q) && (_Q >= 0.00) && (_Q <= 1.00)) return true; }
The virtual function clear() in the base class (AbstractState) is redefined to chain the original function, and call the local (protected) one _RPclearSat() to erase the new cached variables

public part:

The functions of saturated properties, calc_saturated_liquid_keyed_output(parameters key) and calc_saturated_vapor_keyed_output(parameters key) are extended for many key parameters to cope with the reqiurment of coolpropsolver in the ExternalMedia library.
calc_hmolar and calc_smolar is implemented to specialize the calculation, particularly, in case of two-phase fluid.
Some derivatives of CoolPropDbl output - see the Purpose section - in the saturated and two-phase state are implemented.
calc_first_saturation_deriv(parameters Of1, parameters Wrt1);
calc_first_saturation_deriv(parameters Of1, parameters Wrt1, CoolPropDbl rhoL_mol, CoolPropDbl rhoV_mol);
calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant);
calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end);
update_DmolarT_direct() is implemented, as it is just required by the splined derivative, using THERMdll subroutine.

Verification Process

This is a preliminary implementation just to know whether it works, and will be developed based on the discussions.
Thus, a simple CO2 thermodynamic state model (TestCoolProp.mo) is tested successfully in OpenModelica
[modelica-3rdparty/ExternalMedia/pull/27].
I will be bring some test runs after clarifying the following issues, mainly, about some decisions to make.

Discussion for Open Issues

In mixture, the REFPROP backend just seems to prohibit the well-defined phase indicator of pure fluid. So, I cannot help making a local function of _RPcheckTwophase(), which looks redundant however. Is there any possibility to actualize the GetRPphase() function also in case of mixture? Moreover, the ExternalMedia library also requires such a strict indicator of phase as a basis. So, how can we manage that?
In case of two-phase, the HEOS backend evaluate some thermodynamic quantities (like h or s) by means of weighted sum with respect to the vapor quality Q. Am I right to apply such a manner also to the REFPROP backend's calc_hmolar() and so on?
In the functions of two-phase derivatives, there are some copy of backends (instantiations of REFPROP backend) just to duplicate the existing routines for the HEOS backend. However, I'm not sure such an instantiation is safe. So, do we need to flatten the logic to avoid a risk?
As it is designed, the cached variables are always supposed to be referred by _calc_saturated_liquid_output( CachedElement&) or _calc_saturated_vapor_output( CachedElement&) for an update in case of empty state. However, it's error-prone a little bit and lengthy, because the cached element is assigned as an argument to the functions with respect to the state (liquid or vapor). So, it looks better to design a derived class of CachedElement to have a shared pointer to the host class (REFPROP backend class), and to implement an automatic checking and updating. Is it plausible?

Plan of Completion

I'll revise the commit according to the conclusions of following discussions. Moreover, there will be updates based on the bug fixing with test runs.

Additional Information

In addition, I fixed the saturated quantity calculation in PQ_INPUTS and QT_INPUTS in update(), on the issue [/issues/1502], and a mistake of the factor 1000 with kPa unt.

…rated states

ibell · 2021-04-16T02:39:06Z

Why are you doing this in C++? What are you catching:

try {
some call..
    }
    catch(std::logic_error) {
	clear();
	throw ValueError(format("Error in THERMdll: during the calculation of vapor state.").c_str());
    }

ibell · 2021-04-16T02:39:59Z

Please use consistent indentation. Recommend you use 4 spaces per tab, no tabs.

ibell · 2021-04-16T02:37:11Z

+       if (_RPcheckTwophase() == false)
+	 throw ValueError("Not two phase fluid: saturation derivative is not possible");
+
+       shared_ptr<REFPROPMixtureBackend> SatL(new REFPROPMixtureBackend(this->calc_fluid_names())), SatV(new REFPROPMixtureBackend(this->calc_fluid_names()));


Definitely this won't work. You need to just call methods with given temperature and density

ibell · 2021-04-16T02:37:30Z

+       if (_RPcheckTwophase() == false)
+	 throw ValueError("Not two phase fluid: saturation derivative is not possible");
+
+       shared_ptr<REFPROPMixtureBackend> SatL(new REFPROPMixtureBackend(this->calc_fluid_names())), SatV(new REFPROPMixtureBackend(this->calc_fluid_names()));


Definitely this won't work. You need to just call methods with given temperature and density

ibell · 2021-04-16T02:37:34Z

+       SatL->update(QT_INPUTS, 0.0, _T); // <--- This looks never effecient,
+       SatV->update(QT_INPUTS, 1.0, _T); //      but just mimics the saturate states of HEOS..
+
+       shared_ptr<REFPROPMixtureBackend> Liq(new REFPROPMixtureBackend(this->calc_fluid_names())), End(new REFPROPMixtureBackend(this->calc_fluid_names()));


Definitely this won't work. You need to just call methods with given temperature and density

dongkeun-oh · 2021-04-16T08:26:28Z

Why are you doing this in C++? What are you catching:

try {
some call..
    }
    catch(std::logic_error) {
	clear();
	throw ValueError(format("Error in THERMdll: during the calculation of vapor state.").c_str());
    }

In my intention, it means "converting standard exceptions of logical (mathematical) error to the ValueError".
Uncontrolled exit by exception, in the ExternalMedia library in Modelica in my case, has been a painful matter, which makes the application hang, whenever any bold attempt (extremely large or small inputs, for instance) usually happens with respect to the direction of iterative step to seek the solution depending on the properties of fluid. So, I would like to make "all computations out of sense" just throw the ValueError, to manage them, for instance, simply showing warning messages.
However, it matters whether such a purpose is implemented as intended; am I right?

dongkeun-oh · 2021-04-16T08:28:38Z

Please use consistent indentation. Recommend you use 4 spaces per tab, no tabs.

Oops, sorry!!
Thank you for reminding me.

dongkeun-oh · 2021-04-16T08:58:39Z

#2016 (comment)

Definitely this won't work. You need to just call methods with given temperature and density

Strange to say, it seems work so far, but, of course, I will dig more..
So, could you please to check the attached test which I tried for verification?
At least in such a pure fluid, I couldn't find anything wrong yet.
I didn't do the mixture of "Methane&Ethane", because it is slightly inconsistent with the outputs between the REFPROP backend and the HEOS one; is there any good model of mixture to be an appropriate benchmarking?

from pylab import *
import CoolProp.CoolProp as CP

Medium = "Methane"
# two abstract states based on the REFPROP backend and HEOS backend for each..
REFP = CP.AbstractState("REFPROP", Medium)
HEOS = CP.AbstractState("HEOS", Medium)

T = linspace(91, 189, 50)
dTdPsat_RP = []  #dTdP along saturation boundary by REFPROP
dTdPsat_HS = []  #dTdP along saturation boundary by HEOS
dDdTsat_RP = []  #dDdT along saturation boundary by REFPROP
dDdTsat_HS = []  #dDdT along saturation boundary by HEOS
dHdTsat_RP = []  #dHdT along saturation boundary by REFPROP
dHdTsat_HS = []  #dHdT along saturation boundary by HEOS
dDdPsat_RP = []  #dDdP along saturation boundary by REFPROP
dDdPsat_HS = []  #dDdP along saturation boundary by HEOS
dHdPsat_RP = []  #dHdP along saturation boundary by REFPROP
dHdPsat_HS = []  #dHdP along saturation boundary by HEOS

for Tv in T:
    HEOS.update(CP.QT_INPUTS, 0.1, Tv)
    REFP.update(CP.QT_INPUTS, 0.1, Tv)
    dTdPsat_RP.append( REFP.first_saturation_deriv(CP.iT, CP.iP))
    dTdPsat_HS.append( HEOS.first_saturation_deriv(CP.iT, CP.iP))
    dDdTsat_RP.append( REFP.first_saturation_deriv(CP.iDmolar, CP.iT))
    dDdTsat_HS.append( HEOS.first_saturation_deriv(CP.iDmolar, CP.iT))
    dHdTsat_RP.append( REFP.first_saturation_deriv(CP.iHmolar, CP.iT))
    dHdTsat_HS.append( HEOS.first_saturation_deriv(CP.iHmolar, CP.iT))
    dDdPsat_RP.append( REFP.first_saturation_deriv(CP.iDmolar, CP.iP))
    dDdPsat_HS.append( HEOS.first_saturation_deriv(CP.iDmolar, CP.iP))
    dHdPsat_RP.append( REFP.first_saturation_deriv(CP.iHmolar, CP.iP))
    dHdPsat_HS.append( HEOS.first_saturation_deriv(CP.iHmolar, CP.iP))

_, (ax1,ax2) = subplots(2,1)
ax1.plot( T, dTdPsat_RP)
ax1.plot( T, dTdPsat_HS)
ax1.set_ylabel('$\\left(\\frac{dT}{dp}\\right)_{\\rm sat.}$')
ax2.plot(T, array(dTdPsat_RP)-array(dTdPsat_HS))
ax2.set_ylabel('error of $\\left(\\frac{dT}{dp}\\right)_{\\rm sat.}$')
ax2.set_xlabel('$T~(K)$')

_, (ax1,ax2) = subplots(2,1)
ax1.plot( T, dDdTsat_RP)
ax1.plot( T, dDdTsat_HS)
ax1.set_ylabel('$\\left(\\frac{d\\rho}{dT}\\right)_{\\rm sat.}$')
ax2.plot(T, array(dDdTsat_RP)-array(dDdTsat_HS))
ax2.set_ylabel('error of $\\left(\\frac{d\\rho}{dT}\\right)_{\\rm sat.}$')
ax2.set_xlabel('$T~(K)$')

_, (ax1,ax2) = subplots(2,1)
ax1.plot( T, dHdTsat_RP)
ax1.plot( T, dHdTsat_HS)
ax1.set_ylabel('$\\left(\\frac{dH}{dT}\\right)_{\\rm sat.}$')
ax2.plot(T, array(dHdTsat_RP)-array(dHdTsat_HS))
ax2.set_ylabel('error of $\\left(\\frac{dH}{dT}\\right)_{\\rm sat.}$')
ax2.set_xlabel('$T~(K)$')

_, (ax1,ax2) = subplots(2,1)
ax1.plot( T, dDdPsat_RP)
ax1.plot( T, dDdPsat_HS)
ax1.set_ylabel('$\\left(\\frac{d\\rho}{dp}\\right)_{\\rm sat.}$')
ax2.plot(T, array(dDdPsat_RP)-array(dDdPsat_HS))
ax2.set_ylabel('error of $\\left(\\frac{d\\rho}{dp}\\right)_{\\rm sat.}$')
ax2.set_xlabel('$T~(K)$')

_, (ax1,ax2) = subplots(2,1)
ax1.plot( T, dHdPsat_RP)
ax1.plot( T, dHdPsat_HS)
ax1.set_ylabel('$\\left(\\frac{dH}{dp}\\right)_{\\rm sat.}$')
ax2.plot(T, array(dHdPsat_RP)-array(dHdPsat_HS))
ax2.set_ylabel('error of $\\left(\\frac{dH}{dp}\\right)_{\\rm sat.}$')
ax2.set_xlabel('$T~(K)$')

show()

dongkeun-oh · 2021-04-16T12:56:28Z

By the way, Ian(@ibell). I'm also a bit leery; please see the 3rd issue item in my original posting, and clarify whether your mention implies such a risk. So, I guess that it is possible to flatten the routines, as I mentioned in the same lines, without such a copy of state. If you still recommend a change, I'm able to try such a direction.

dongkeun-oh · 2021-04-17T04:33:03Z

#2016 (comment)

Definitely this won't work. You need to just call methods with given temperature and density

I just found many trivial (but serious) mistakes not only in my commit, but also with the previous example which was nothing to do with your comment, but a test for first_saturation_deriv;
I fixed the codes with the following tests of first_two_phase_derive based on shared pointers to the backend

from pylab import *
import CoolProp.CoolProp as CP
rcParams['mathtext.fontset'] = 'stix'
rcParams['font.family'] = 'StixGeneral'

Medium = "Methane"
# two abstract states based on the REFPROP backend and HEOS backend for each..
REFP = CP.AbstractState("REFPROP", Medium)
HEOS = CP.AbstractState("HEOS", Medium)

T = linspace(91, 189, 50)
Q = linspace(0.1, 1.0, 5)

_, ax = subplots(4,1, figsize=(8, 10))
subplots_adjust(hspace=0.2)
ax[0].set_ylabel(r'$\left(\frac{\partial \rho}{\partial H}\right)_{_p}$', fontsize=18)
ax[1].set_ylabel(r'error of $\left(\frac{\partial \rho}{\partial H}\right)_{_p}$', fontsize=18)
ax[2].set_ylabel(r'$\left(\frac{\partial \rho}{\partial p}\right)_{_H}$', fontsize=18)
ax[3].set_ylabel(r'error of $\left(\frac{\partial \rho}{\partial p}\right)_{_H}$', fontsize=18)            

for _Q in Q:
    dDdH_p_RP = []  #(dDdH)_p in two phase by REFPROP
    dDdH_p_HS = []  #(dDdH)_p in two phase by  by HEOS
    dDdp_H_RP =[]  #(dDdp)_H in two phase by REFPROP
    dDdp_H_HS = []  #(dDdp)_H in two phase by  by HEOS
    for _T in T:
        HEOS.update(CP.QT_INPUTS, _Q, _T) 
        REFP.update(CP.QT_INPUTS, _Q, _T) 
        dDdH_p_RP.append( REFP.first_two_phase_deriv(CP.iDmolar, CP.iHmolar, CP.iP))
        dDdH_p_HS.append( HEOS.first_two_phase_deriv(CP.iDmolar, CP.iHmolar, CP.iP))
        dDdp_H_RP.append( REFP.first_two_phase_deriv(CP.iDmolar, CP.iP, CP.iHmolar))
        dDdp_H_HS.append( HEOS.first_two_phase_deriv(CP.iDmolar, CP.iP, CP.iHmolar))
    
    ax[0].plot( T, dDdH_p_RP, '-')
    ax[0].plot( T, dDdH_p_HS, '.', color='gray')
    ax[1].plot(T, array(dDdH_p_RP)-array(dDdH_p_HS))
    
    ax[2].plot( T, dDdp_H_RP)
    ax[2].plot( T, dDdp_H_HS, '.', color='gray')
    ax[3].plot(T, array(dDdp_H_RP)-array(dDdp_H_HS))
ax[1].legend([*map(lambda x : "Q = " +str(x), Q)])
ax[3].legend([*map(lambda x : "Q = " +str(x), Q)])

Eventually, I made it work.
Here is the result, and I'm terribly sorry for such a mass!!

And another example of splined derivative is also successful, which is the same to the famous example http://www.coolprop.org/coolprop/LowLevelAPI-1.py except for the backend; I just applied AS = CoolProp.AbstractState('REFPROP','Water') ..

thorade · 2025-02-26T08:48:04Z

How much work would it be to finish this PR?
Came here from this ticket: modelica-3rdparty/ExternalMedia#48

ibell · 2025-02-28T02:33:14Z

Rather a lot I think. From a cursory glance quite a few functions should have been implemented in the base class, rather than in the derived class. In fact I think many of them should already work (the splined ones for instance).

The first commit for REFPROP backend to include the functions of satu…

02b0b6a

…rated states

jfulem reviewed Apr 15, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

jfulem reviewed Apr 15, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

jfulem reviewed Apr 15, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

jfulem reviewed Apr 15, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

dongkeun-oh added 4 commits April 15, 2021 17:39

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

9395f21

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

6e5954c

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

90a51b9

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

418f8c9

jfulem mentioned this pull request Apr 15, 2021

Testing CoolProp 6 (6.4.2) modelica-3rdparty/ExternalMedia#27

Closed

ibell requested changes Apr 16, 2021

View reviewed changes

dongkeun-oh commented Apr 16, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp

dongkeun-oh added 2 commits April 16, 2021 18:06

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

6ec1b14

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

59f60bd

dongkeun-oh requested a review from ibell April 16, 2021 12:25

dongkeun-oh commented Apr 17, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

dongkeun-oh added 2 commits April 17, 2021 12:34

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

8548d15

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

dd1ff51

dongkeun-oh commented Apr 17, 2021

View reviewed changes

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

Comment thread src/Backends/REFPROP/REFPROPMixtureBackend.cpp Outdated

dongkeun-oh added 2 commits April 17, 2021 13:16

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

7888be2

Update src/Backends/REFPROP/REFPROPMixtureBackend.cpp

25d1344

ianhbell mentioned this pull request Apr 28, 2021

Excel: dH/dT at constant T, or at constant P, wrong values at odd input values usnistgov/REFPROP-wrappers#369

Closed

Conversation

dongkeun-oh commented Apr 13, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Purpose

Description of the Change

protected part:

public part:

Verification Process

Discussion for Open Issues

Plan of Completion

Additional Information

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ibell commented Apr 16, 2021

Uh oh!

ibell commented Apr 16, 2021

Uh oh!

ibell Apr 16, 2021

Choose a reason for hiding this comment

Uh oh!

ibell Apr 16, 2021

Choose a reason for hiding this comment

Uh oh!

ibell Apr 16, 2021

Choose a reason for hiding this comment

Uh oh!

dongkeun-oh commented Apr 16, 2021

Uh oh!

dongkeun-oh commented Apr 16, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dongkeun-oh commented Apr 16, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

dongkeun-oh commented Apr 16, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

dongkeun-oh commented Apr 17, 2021

Uh oh!

thorade commented Feb 26, 2025

Uh oh!

ibell commented Feb 28, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

dongkeun-oh commented Apr 13, 2021 •

edited

Loading

dongkeun-oh commented Apr 16, 2021 •

edited

Loading

dongkeun-oh commented Apr 16, 2021 •

edited

Loading

dongkeun-oh commented Apr 16, 2021 •

edited

Loading