CoolProp/dev/scripts/fit_avoid_ECS.py

from CoolProp import CoolProp as CP
from PDSim.misc.datatypes import Collector
import numpy as np
import matplotlib.pyplot as plt
from scipy.odr import *
import textwrap

#
#fluid = 'Propane'
#Rfluid = 'REFPROP-propane'
#e_k = 263.88
#sigma = 0.49748
#
fluid = 'DimethylEther'
Rfluid = 'REFPROP-DME'
e_k = 329.72
sigma = 0.5529
molemass = CP.Props(fluid,'molemass')

Ttriple = CP.Props(fluid,'Ttriple')
Tcrit = CP.Props(fluid,'Tcrit')
rhocrit = CP.Props(fluid,'rhocrit')
n = 6
m = 3
NP = 1
Nb = 0
N = (n-1)*(m+1)+3+Nb

mu,mu_dilute,RHO,TTT = Collector(),Collector(),Collector(),Collector()

rhomax = CP.Props('D','T',Ttriple,'Q',0,fluid)
#Build a database of "experimental" data
for T in np.linspace(Ttriple,Tcrit+30,400):
    for rho in np.linspace(1e-10,rhomax,400):
        muval = CP.Props('V','T',T,'D',rho,Rfluid)
        mudilute = CP.viscosity_dilute(fluid,T,rho,e_k,sigma)
        
        #Want positive value, and single-phase
        if (muval > 0 and T > Tcrit or rho > CP.rhosatL_anc(fluid,T) or rho < CP.rhosatV_anc(fluid,T)):
            mu << muval
            mu_dilute << mudilute
            TTT << T
            RHO << rho

from CoolProp.Plots.Plots import Trho
Trho(fluid)
plt.plot(RHO.vec,TTT.vec,'.')
plt.show()

#tau = np.array(TTT.vec)/Tcrit
tau = np.array(TTT.vec)/Tcrit
delta = np.array(RHO.vec)/rhocrit
Tstar = np.array(TTT.vec)/e_k

#Define the objective function
def OBJECTIVE_fit(c,x):
    tau = x[0,:]
    delta = x[1,:]
    #Unpack the inputs into e matrix and f vector
    e = np.zeros((n+1,m+1))
    
    sum = 0
    k = 0
    for i in range(2,n+1):
        for j in range(0,m+1):
            e[i][j] = c[k]
            sum += e[i][j]*delta**i/tau**j
            k += 1
            
    for o in range(0,NP):
        f1 = c[k+o*3]
        g1 = c[k+1+o*3]
        g2 = c[k+2+o*3]
    
    delta_0 = g1*(1+g2*tau**0.5)
    sum += f1*(delta/(delta_0-delta)-delta/delta_0)
    return sum + np.array(mu_dilute.vec)
    
print 'starting fit'
XXX = np.r_[np.array(tau,ndmin = 2), np.array(delta,ndmin=2)]
mod = Model(OBJECTIVE_fit)
mydata = Data(XXX.copy(), np.array(mu.vec))
beta0  = [1 for _ in range(N)]
myodr = ODR(mydata, mod, beta0=beta0)
myoutput = myodr.run()
myoutput.pprint()
print myoutput.sum_square
YFIT = OBJECTIVE_fit(myoutput.beta,XXX)
plt.plot(np.array(mu.vec),YFIT/np.array(mu.vec),'o')
plt.show()

rel_error = (YFIT)/np.array(mu.vec)-1
MAE = np.mean(np.abs(rel_error))*100
SSE = np.sum(np.power(YFIT-np.array(mu.vec),2))
print SSE

def write_output(c):
    e = np.zeros((n+1,m+1))
    k = 0
    edata = ''
    for i in range(2,n+1):
        erow = ''
        for j in range(0,m+1):
            e[i][j] = c[k]
            erow += 'e[{i:d}][{j:d}] = {val:0.16g}; '.format(val = e[i][j],i=i,j=j)
            k += 1
        edata += erow + '\n'
    f1 = c[k]
    g1 = c[k+1]
    g2 = c[k+2]
        
    template = textwrap.dedent(
    """
    double {name:s}Class::viscosity_Trho(double T, double rho)
    {{
        // This function was generated by fitting REFPROP ECS data
        // to the functional form of Vogel, 1998 (propane viscosity)
        // The script entitled dev/fit_avoid_ECS.py was used to make this 
        // function.  The mean absolute error of the fit is equal to
        // {MAE:g} % 
        
        double delta_0, sum, DELTA_H_eta, e_k, sigma, tau, delta;
        double e[{n:d}+1][{m:d}+1];
        
        tau = T/reduce.T; //[Opposite to normal definition]
        delta = rho/reduce.rho;
        
        // Load the coefficients
        double f1 = {f1:0.16g}, g1 = {g1:0.16g}, g2 = {g2:0.16g};
        for (int i=0;i<={n:d};i++){{ for(int j=0;j<={m:d};j++){{ e[i][j]=0.0; }} }}
        {edata:s}
        delta_0=g1*(1+g2*sqrt(tau)); //[no units]
        sum=0;
        for (int i=2;i<={n:d};i++){{
            for (int j=0;j<={m:d};j++){{
                sum += e[i][j]*pow(delta,i)/pow(tau,j);
            }}
        }}
        DELTA_H_eta = sum + f1*(delta/(delta_0-delta)-delta/delta_0); //[Pa-s]
        
        try{{
            // Get the ECS params for the fluid if it has them
            ECSParams(&e_k,&sigma);
        }}
        catch(NotImplementedError)
        {{
            throw ValueError(format("Your fluid does not implement ECSParams"));
        }}
        
        return viscosity_dilute(T,e_k,sigma) + DELTA_H_eta;
    }} 
    """
    )
    
    values = dict(f1 = f1,
                  g1 = g1,
                  g2 = g2,
                  n = n,
                  m = m,
                  name = fluid,
                  edata = edata,
                  MAE = MAE)
    
    print template.format(**values)

write_output(myoutput.beta)