A commit for the REFPROP backend to include the functions of satu…

[dongkeun-oh]

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.

[ibell]

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]

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

[ibell]

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

[ibell]

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

[ibell]

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

[dongkeun-oh]

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]

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

Oops, sorry!!
Thank you for reminding me.

[dongkeun-oh]

#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]

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]

#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') ..