CoolProp/dev/TTSE/validate_TTSE.py

import matplotlib
matplotlib.use('WXAgg')
import CoolProp
from CoolProp.Plots import Ph
from CoolProp.Plots.Plots import Trho,Ps,PT,Prho
import CoolProp.CoolProp as CP
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import random
import numpy as np
from math import log,exp
random.seed()

def check_Trho(N=5000,param='P',fluid='R245fa'):
    values = []
    CP.enable_TTSE_LUT(fluid)
    try:
        CP.Props('D','P',CP.Props(fluid,'ptriple')+1,'Q',1,fluid)
    except:
        return []
    #CP.set_TTSESinglePhase_LUT_size(fluid,500,500)
    hmin,hmax,pmin,pmax = CP.get_TTSESinglePhase_LUT_range(fluid)
    for i in range(N):
        x1 = random.random()
        h = x1*hmin+(1-x1)*hmax
        x2 = random.random()
        logp = x2*log(pmin)+(1-x2)*log(pmax)
        p = exp(logp)
        try:
            try:
                #Get the T,rho from the EOS directly without the LUT
                CP.disable_TTSE_LUT(fluid)
                s = CP.Props('S','P',p,'H',h,fluid)   
                T = CP.Props('T','P',p,'H',h,fluid)
                rho = CP.Props('D','P',p,'H',h,fluid)
            except:
                print 'EOS failed: ', p,h
                raise
            #Now get p,h from the T,rho
            CP.enable_TTSE_LUT(fluid)
            val = CP.Props(param,'T',T,'D',rho,fluid)
            CP.disable_TTSE_LUT(fluid)
            valREFPROP = CP.Props(param,'T',T,'D',rho,fluid)
            #print T,rho,val,valREFPROP,(val/valREFPROP-1)*100
            if abs(val-valREFPROP)>0.00001:
                raise ValueError
        except ValueError:
            print 'TTSE failed: ', T,rho
            values.append((T,rho,0,0))
            pass
    return values


fluid = 'R245fa'
values = check_Trho(param='P',fluid=fluid)
if len(values) == 0:
    print 'good'

T,rho,values_withTTSE,values_noTTSE = zip(*values)
    
CP.disable_TTSE_LUT(fluid)
Trho(fluid)

plt.plot(rho,T,'.',mfc='none')
plt.savefig('Trho_TTSE_Validation.png',dpi=300)
plt.gca().set_xscale('log')
plt.show()
    
quit()

def check_Pother(N=5000,param='T',other='S',fluid='R245fa'):
    values = []
    CP.enable_TTSE_LUT(fluid)
    try:
        CP.Props('D','P',CP.Props(fluid,'ptriple')+1,'Q',1,fluid)
    except:
        return []
    #CP.set_TTSESinglePhase_LUT_size(fluid,500,500)
    hmin,hmax,pmin,pmax = CP.get_TTSESinglePhase_LUT_range(fluid)
    for i in range(N):
        x1 = random.random()
        h = x1*hmin+(1-x1)*hmax
        x2 = random.random()
        logp = x2*log(pmin)+(1-x2)*log(pmax)
        p = exp(logp)
        try:
            try:
                #Get the T,rho from the EOS directly without the LUT
                CP.disable_TTSE_LUT(fluid)
                s = CP.Props('S','P',p,'H',h,fluid)   
                T = CP.Props('T','P',p,'H',h,fluid)
                rho = CP.Props('D','P',p,'H',h,fluid)
            except:
                print 'EOS failed: ', p,h
                raise
            #Now get p,h from the T,rho
            CP.enable_TTSE_LUT(fluid)
            if other =='S':
                other_val = s
            elif other =='T':
                other_val = T
            elif other == 'D':
                other_val = rho
            else:
                raise ValueError
            val = CP.Props(param,'P',p,other,other_val,fluid)
        except ValueError:
            print 'TTSE failed: ', p,other_val
            values.append((p,other_val,0,0))
            pass
    return values

fluid = 'R245fa'
for other in ['D','T','S']:
    if other == 'D':
        param = 'T'
    elif other == 'T':
        param = 'D'
    elif other == 'S':
        param = 'T'
    values = check_Pother(other=other)
    if len(values) == 0:
        print 'good: ',other
        continue
    
    p,othervals,values_withTTSE,values_noTTSE = zip(*values)
    
    CP.disable_TTSE_LUT(fluid)
    if other =='D':
        Prho(fluid)
    elif other =='S':
        Ps(fluid)
    elif other =='T':
        PT(fluid)
    
    plt.plot(othervals,p,'.',mfc='none')
    plt.gca().set_yscale('log')
    plt.savefig('P'+other+'_TTSE_Validation.png',dpi=300)
    plt.show()

quit()

def check(N=5000,param='D',fluid = 'R245fa'):
    values = []
    CP.enable_TTSE_LUT(fluid)
    
    try:
        CP.Props('D','P',CP.Props(fluid,'ptriple')+1,'Q',1,fluid)
    except:
        return []
    #CP.set_TTSESinglePhase_LUT_size(fluid,500,500)
    hmin,hmax,pmin,pmax = CP.get_TTSESinglePhase_LUT_range(fluid)
    for i in range(N):
        x1 = random.random()
        h = x1*hmin+(1-x1)*hmax
        x2 = random.random()
        logp = x2*log(pmin)+(1-x2)*log(pmax)
        p = exp(logp)
        try:
            CP.enable_TTSE_LUT(fluid)
            value_withTTSE = CP.Props(param,'P',p,'H',h,fluid)
            CP.disable_TTSE_LUT(fluid)
            value_noTTSE = CP.Props(param,'P',p,'H',h,fluid)
            values.append((h,p,value_withTTSE,value_noTTSE))
        except ValueError:
            pass
        
    return values

Ncols = 10
Nrows = 10
for parameter in ['D','T','S','C']:
    fig = plt.figure(figsize=(40,40))  
    for Index,fluid in enumerate(sorted(CoolProp.__fluids__)):
        print fluid
        ax = fig.add_subplot(Ncols,Nrows,Index+1)
        ax.set_title(fluid)
        values = check(fluid = fluid, param = parameter)
        if len(values) == 0:
            continue
        h,p,values_withTTSE,values_noTTSE = zip(*values)
        ax = fig.add_subplot(Ncols,Nrows,Index+1)
        CP.disable_TTSE_LUT(fluid)
        Ph(fluid)
        CP.enable_TTSE_LUT(fluid)
        error = (np.array(values_withTTSE)/np.array(values_noTTSE)-1)*100
        plt.scatter(h,p,s=8,c=np.abs(error),norm = LogNorm(),edgecolor='none',vmin = 1e-16, vmax = 10)
        plt.gca().set_yscale('log')
        plt.colorbar()
    
    plt.savefig(parameter+'_TTSE.png',dpi=200)
    plt.tight_layout()
    plt.close()