diff --git a/dev/incompressible_liquids/CPIncomp/BaseObjects.py b/dev/incompressible_liquids/CPIncomp/BaseObjects.py index 676d29d9..befdb3ca 100644 --- a/dev/incompressible_liquids/CPIncomp/BaseObjects.py +++ b/dev/incompressible_liquids/CPIncomp/BaseObjects.py @@ -17,6 +17,7 @@ class IncompressibleData(object): INCOMPRESSIBLE_POLYNOMIAL = 'polynomial' INCOMPRESSIBLE_EXPPOLYNOMIAL = 'exppolynomial' INCOMPRESSIBLE_EXPONENTIAL = 'exponential' + INCOMPRESSIBLE_LOGEXPONENTIAL= 'logexponential' INCOMPRESSIBLE_POLYOFFSET = 'polyoffset' INCOMPRESSIBLE_CHEBYSHEV = 'chebyshev' @@ -35,6 +36,7 @@ class IncompressibleData(object): self.data = None # None #np.zeros((10,10)) self.xData = None # In case you need a customised first data set (temperature?) self.yData = None # In case you need a customised second data set (concentration?) + self.rsq = None # Coefficient of determination self.DEBUG = False @staticmethod @@ -53,6 +55,10 @@ class IncompressibleData(object): elif eqnType==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL: #if y!=0.0: raise ValueError("This is 1D only, use x not y.") return IncompressibleData.baseExponential(c, x) + + elif eqnType==IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL: + #if y!=0.0: raise ValueError("This is 1D only, use x not y.") + return IncompressibleData.baseLogexponential(c, x) elif eqnType==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL: return np.exp(np.polynomial.polynomial.polyval2d(x-xbase, y-ybase, c)) @@ -67,6 +73,14 @@ class IncompressibleData(object): raise ValueError("You have to provide a 3,1 matrix of coefficients, not ({0},{1}).".format(r,c)) coeffs_tmp = np.array(coeffs.flat) return np.exp(np.clip( (coeffs_tmp[0]/ ( x+coeffs_tmp[1] ) - coeffs_tmp[2]),IncompressibleData.minLog,IncompressibleData.maxLog)) + + @staticmethod + def baseLogexponential(co, x): + r,c,coeffs = IncompressibleFitter.shapeArray(co) + if not ( (r==4 and c==1) or (r==1 and c==4) ): + raise ValueError("You have to provide a 4,1 matrix of coefficients, not ({0},{1}).".format(r,c)) + coeffs_tmp = np.array(coeffs.flat) + return np.exp(np.clip(np.log(1.0/(x+coeffs_tmp[0]) + 1.0/(x+coeffs_tmp[1])**2)*coeffs_tmp[2]+coeffs_tmp[3],IncompressibleData.minLog,IncompressibleData.maxLog)) @staticmethod def basePolyOffset(co, x): @@ -85,25 +99,51 @@ class IncompressibleData(object): def fitCoeffs(self, xbase, ybase, x=None, y=None): - if (x!=None and self.xData!=None and not np.allclose(x, self.xData)) \ + if (x!=None and self.xData!=None and not IncompressibleFitter.allClose(x, self.xData)) \ or (x==None and self.xData==None): raise ValueError("I do not know which x-value you would like to use. Define either x or self.xData.") - if (y!=None and self.yData!=None and not np.allclose(y, self.yData)) \ + if (y!=None and self.yData!=None and not IncompressibleFitter.allClose(y, self.yData)) \ or (y==None and self.yData==None): raise ValueError("I do not know which y-value you would like to use. Define either y or self.yData.") if x==None and self.xData!=None: x=self.xData if y==None and self.yData!=None: y=self.yData - res = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \ + #res = None + #r2 = None + + res,r2 = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \ xbase=xbase, ybase=ybase, \ eqnType=self.type, \ coeffs=self.coeffs, DEBUG=self.DEBUG) + #count = 0 + #while r2<0.9 and count<1: + ##if self.DEBUG: print("Poor solution found, trying once more with more coefficients.") + #if self.type==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL: + #if self.DEBUG: print("Poor solution found, trying once more with log exponential.") + #self.type=IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL + #self.coeffs = np.array([-250,-250,2,2]) + #res,r2 = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \ + #xbase=xbase, ybase=ybase, \ + #eqnType=self.type, \ + #coeffs=self.coeffs, DEBUG=self.DEBUG) + + #if self.type==IncompressibleData.INCOMPRESSIBLE_POLYNOMIAL or \ + #self.type==IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL: + #if self.DEBUG: print("Poor solution found, trying once more with more coefficients.") + #self.coeffs = np.zeros(np.round(np.array(self.coeffs.shape) * 1.5)) + #res,r2 = IncompressibleFitter.fitter(x=x, y=y, z=self.data, \ + #xbase=xbase, ybase=ybase, \ + #eqnType=self.type, \ + #coeffs=self.coeffs, DEBUG=self.DEBUG) + #count += 1 + if res==None: - raise ValueError("There was a fitting error.") - elif np.allclose(res, self.coeffs): + if self.DEBUG: print("There was a fitting error, no solution found.") + elif IncompressibleFitter.allClose(res, self.coeffs): if self.DEBUG: print("Coefficients did not change.") else: self.coeffs = res + self.rsq = r2 def setxData(self, xData): @@ -138,6 +178,12 @@ class IncompressibleData(object): class IncompressibleFitter(object): + @staticmethod + def allClose(a,b): + if np.array(a).shape==np.array(b).shape: + return np.allclose(a, b) + else: return False + @staticmethod def shapeArray(array, axs=0): """ @@ -239,7 +285,7 @@ class IncompressibleFitter(object): if cr==1 and cc==1: if DEBUG: print("No coefficients left to fit, aborting procedure.") coeffs = np.array([[0.0]]) - return + return coeffs, None if (xr1: # 2D fitting raise ValueError("There are no other 2D fitting functions than polynomials, cannot use {0}.".format(eqnType)) @@ -297,6 +345,8 @@ class IncompressibleFitter(object): if y_order==None: raise ValueError("You did not provide data for y_order.") if DEBUG==None: raise ValueError("You did not provide data for DEBUG.") + r2 = None + x_order += 1 y_order += 1 #To avoid overfitting, we only use the upper left triangle of the coefficient matrix @@ -349,12 +399,18 @@ class IncompressibleFitter(object): if DEBUG: print(rank) if DEBUG: print(singulars) + if resids.size>0: + r2 = 1 - resids / (zz.size * zz.var()) + else: + r2 = 0 + #print("\n r2 2d: ",r2.shape,r2,"\n") + #Rearrange coefficients to a matrix shape C = np.zeros((x_order,y_order)) for i, (xi,yi) in enumerate(xy_exp): # makes columns C[xi][yi] = coeffs[i] - return C + return C, r2 @staticmethod def getCoeffsIterative1D(x_in, z_in, eqnType, coeffs, DEBUG=False): @@ -365,6 +421,8 @@ class IncompressibleFitter(object): if coeffs==None: raise ValueError("You did not provide data for the coefficients.") if DEBUG==None: raise ValueError("You did not provide data for DEBUG.") + r2 = None + #fit = "Powell" # use Powell's algorithm #fit = "BFGS" # use Broyden-Fletcher-Goldfarb-Shanno #fit = "LMA" # use the Levenberg-Marquardt algorithm from curve_fit @@ -377,7 +435,7 @@ class IncompressibleFitter(object): expLog = False # Fitting the logarithm of z_in? - if (eqnType==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL): + if (eqnType==IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL or eqnType==IncompressibleData.INCOMPRESSIBLE_LOGEXPONENTIAL): expLog = True xData = np.array(x_in.flat) @@ -429,16 +487,20 @@ class IncompressibleFitter(object): # print "data: {0}, func: {1}".format(yData[ 2],func(xData[ 2], popt)) # print "data: {0}, func: {1}".format(yData[ 6],func(xData[ 6], popt)) # print "data: {0}, func: {1}".format(yData[-1],func(xData[-1], popt)) - if DEBUG: print("Estimated covariance of parameters: ".format(pcov)) - return popt + #if DEBUG: print("Estimated covariance of parameters: {0}".format(pcov)) + ssErr = np.sqrt(np.diag(pcov)).sum() + ssTot = ((zData-zData.mean())**2).sum() + r2 = 1-(ssErr/ssTot ) + #print("\n r2 FMA: ",r2.shape,r2,"\n") + return popt,r2 else: - print("Fit failed for {0}.".format(algorithm)) - sys.stdout.flush() + if DEBUG: print("Fit failed for {0}.".format(algorithm)) + if DEBUG: sys.stdout.flush() success = False except RuntimeError as e: - print("Exception using "+algorithm+": "+str(e)) - sys.stdout.flush() + if DEBUG: print("Exception using "+algorithm+": "+str(e)) + if DEBUG: sys.stdout.flush() success = False #fit = "MIN" # use a home-made minimisation with Powell and Broyden-Fletcher-Goldfarb-Shanno @@ -452,14 +514,19 @@ class IncompressibleFitter(object): if res.success: success = True if DEBUG: print("Fit succeeded with: {0}".format(algorithm)) - return res.x + if res.has_key('fvec'): + ssErr = (res['fvec']**2).sum() + ssTot = ((zData-zData.mean())**2).sum() + r2 = 1-(ssErr/ssTot ) + #print("\n r2 : ",r2.shape,r2,algorithm,"\n") + return res.x,r2 else: - print("Fit failed for {0}.".format(algorithm)) - sys.stdout.flush() + if DEBUG: print("Fit failed for {0}.".format(algorithm)) + if DEBUG: sys.stdout.flush() success = False except RuntimeError as e: - print("Exception using "+algorithm+": "+str(e)) - sys.stdout.flush() + if DEBUG: print("Exception using "+algorithm+": "+str(e)) + if DEBUG: sys.stdout.flush() success = False # Something went wrong, probably a typo in the algorithm selector @@ -470,16 +537,16 @@ class IncompressibleFitter(object): #print("Fit did not succeed with {0}, reducing tolerance to {1}.".format(algorithm,tol)) success = False counter += 1 - elif tolerance<1e-8: + elif tolerance<1e-10: tolerance *= 1e2 - print("Fit did not succeed, reducing tolerance to {0}.".format(tolerance)) + if DEBUG: print("Fit did not succeed, reducing tolerance to {0}.".format(tolerance)) success = False counter = 0 else: - print("--------------------------------------------------------------") - print("Fit failed for {0}. ".format(fit)) - print("--------------------------------------------------------------") - return None + if DEBUG: print("--------------------------------------------------------------") + if DEBUG: print("Fit failed for {0}. ".format(fit)) + if DEBUG: print("--------------------------------------------------------------") + return None,None diff --git a/dev/incompressible_liquids/CPIncomp/DataObjects.py b/dev/incompressible_liquids/CPIncomp/DataObjects.py index 0ad6a9a8..e3f5565a 100644 --- a/dev/incompressible_liquids/CPIncomp/DataObjects.py +++ b/dev/incompressible_liquids/CPIncomp/DataObjects.py @@ -108,9 +108,45 @@ class SolutionData(object): # objList["viscosity"] = self.viscosity # objList["saturation pressure"] = self.saturation_pressure # return objList + + + def checkT(self, T, p, x): + if self.Tmin <= 0.: raise ValueError("Please specify the minimum temperature.") + if self.Tmax <= 0.: raise ValueError("Please specify the maximum temperature."); + if ((self.Tmin > T) or (T > self.Tmax)): raise ValueError("Your temperature {0} is not between {1} and {2}.".format(T, self.Tmin, self.Tmax)) + TF = 0.0 + if (self.T_freeze.type!=IncompressibleData.INCOMPRESSIBLE_NOT_SET): TF = self.Tfreeze(T, p, x) + if ( T 1.0): raise ValueError("Please specify the minimum concentration between 0 and 1."); + if (self.xmax < 0.0 or self.xmax > 1.0): raise ValueError("Please specify the maximum concentration between 0 and 1."); + if ((self.xmin > x) or (x > self.xmax)): raise ValueError("Your composition {0} is not between {1} and {2}.".format(x, self.xmin, self.xmax)) + else: return True + return False + + def checkTPX(self, T, p, x): + try: + return (self.checkT(T,p,x) and self.checkP(T,p,x) and self.checkX(x)) + except ValueError as ve: + #print("Check failed: {0}".format(ve)) + pass + return False def rho (self, T, p=0.0, x=0.0, c=None): + if not self.checkTPX(T, p, x): return np.NAN if c==None: c=self.density.coeffs if self.density.type==self.density.INCOMPRESSIBLE_POLYNOMIAL: @@ -118,6 +154,7 @@ class SolutionData(object): else: raise ValueError("Unknown function.") def c (self, T, p=0.0, x=0.0, c=None): + if not self.checkTPX(T, p, x): return np.NAN if c==None: c = self.specific_heat.coeffs if self.specific_heat.type==self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL: @@ -131,6 +168,7 @@ class SolutionData(object): return self.c(T,p,x,c) def u (self, T, p=0.0, x=0.0, c=None): + if not self.checkTPX(T, p, x): return np.NAN if c==None: c = self.specific_heat.coeffs if self.specific_heat.type==self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL: @@ -143,15 +181,18 @@ class SolutionData(object): return self.h_u(T,p,x) def visc(self, T, p=0.0, x=0.0, c=None): + if not self.checkTPX(T, p, x): return np.NAN return self.viscosity.baseFunction(T, x, self.Tbase, self.xbase, c=c) def cond(self, T, p=0.0, x=0.0, c=None): + if not self.checkTPX(T, p, x): return np.NAN return self.conductivity.baseFunction(T, x, self.Tbase, self.xbase, c=c) def psat(self, T, p=0.0, x=0.0, c=None): + if (T<=self.TminPsat): return 0.0 return self.saturation_pressure.baseFunction(T, x, self.Tbase, self.xbase, c=c) - def Tfreeze(self, T, p=0.0, x=0.0, c=None): + def Tfreeze(self, T, p=0.0, x=0.0, c=None): if c==None: c = self.T_freeze.coeffs @@ -296,27 +337,31 @@ class DigitalData(SolutionData): if x_in!=None: # Might need update if x!=None: # Both given, check if different mask = np.isfinite(x) - if np.allclose(x[mask], x_in[mask]): + if IncompressibleFitter.allClose(x[mask], x_in[mask]): if DEBUG: print("Both x-arrays are the same, no action required.") updateFile = (updateFile or False) # Do not change a True value to False else: updateFile = True if DEBUG: print("x-arrays do not match. {0} contains \n {1} \n and will be updated with \n {2}".format(self.getFile(dataID),x,x_in)) else: updateFile = True - elif x==None: raise ValueError("Could not load x from file and no x_in provided, aborting.") + elif x==None: + if DEBUG: print("Could not load x from file {0} and no x_in provided, aborting.".format(self.getFile(dataID))) + return None,None,None else: updateFile = (updateFile or False) # Do not change a True value to False if y_in!=None: # Might need update if y!=None: # Both given, check if different mask = np.isfinite(y) - if np.allclose(y[mask], y_in[mask]): + if IncompressibleFitter.allClose(y[mask], y_in[mask]): if DEBUG: print("Both y-arrays are the same, no action required.") updateFile = (updateFile or False) # Do not change a True value to False else: updateFile = True if DEBUG: print("y-arrays do not match. {0} contains \n {1} \n and will be updated with \n {2}".format(self.getFile(dataID),y,y_in)) else: updateFile = True - elif y==None: raise ValueError("Could not load y from file and no y_in provided, aborting.") + elif y==None: + if DEBUG: print("Could not load y from file {0} and no y_in provided, aborting.".format(self.getFile(dataID))) + return None,None,None else: updateFile = (updateFile or False) # Do not change a True value to False if DEBUG: print("Updating data file {0}".format(updateFile)) diff --git a/dev/incompressible_liquids/CPIncomp/DigitalFluids.py b/dev/incompressible_liquids/CPIncomp/DigitalFluids.py index cf576f14..9c6e595e 100644 --- a/dev/incompressible_liquids/CPIncomp/DigitalFluids.py +++ b/dev/incompressible_liquids/CPIncomp/DigitalFluids.py @@ -10,7 +10,8 @@ parameter form. """ from __future__ import division, print_function -from CPIncomp.DataObjects import DigitalData +import numpy as np +from CPIncomp.DataObjects import DigitalData, PureData class LiBrData(DigitalData): """ @@ -63,3 +64,198 @@ class LiBrData(DigitalData): self.saturation_pressure.xData,self.saturation_pressure.yData,self.saturation_pressure.data = self.getArray(dataID=key, func=funcP, x_in=self.temperature.data, y_in=self.concentration.data,DEBUG=self.saturation_pressure.DEBUG) self.saturation_pressure.source = self.saturation_pressure.SOURCE_EQUATION + + +class HyCool20(PureData,DigitalData): + def __init__(self): + DigitalData.__init__(self) + PureData.__init__(self) + + self.name = "HY20" + self.description = "HYCOOL 20, Potassium formate" + self.reference = "Hydro Chemicals" + + self.Tmax = 50 + 273.15 + self.Tmin = -20 + 273.15 + self.TminPsat = self.Tmax + self.Tbase = 0.00 + 273.15 + self.temperature.data = self.getTrange() + + self.density.source = self.density.SOURCE_COEFFS + self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL + self.density.coeffs = np.array([[1202.2],[-0.42918]]) + + self.specific_heat.source = self.specific_heat.SOURCE_COEFFS + self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL + self.specific_heat.coeffs = np.array([[2.955],[0.0023]])*1e3 + + key = 'Cond' + def funcCond(T,x): + T = (T-self.Tbase) + if T <= 20: return 0.001978*T+0.5172 + else: return 0.001005*T+0.5368 + self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(dataID=key,func=funcCond,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.conductivity.DEBUG) + self.conductivity.source = self.conductivity.SOURCE_EQUATION + funcCond = None + + key = 'Mu' + def funcMu(T,x): + T = (T-self.Tbase) + if T <= 20: return 0.07190*np.exp(524.75/(T+142.05)) + else: return T*(0.0005524*T - 0.06281)+2.8536 + self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG) + self.viscosity.source = self.viscosity.SOURCE_EQUATION + funcMu = None + + +class HyCool30(PureData,DigitalData): + def __init__(self): + DigitalData.__init__(self) + PureData.__init__(self) + + self.name = "HY30" + self.description = "HYCOOL 30, Potassium formate" + self.reference = "Hydro Chemicals" + + self.Tmax = 50 + 273.15 + self.Tmin = -30 + 273.15 + self.TminPsat = self.Tmax + self.Tbase = 0.00 + 273.15 + self.temperature.data = self.getTrange() + + self.density.source = self.density.SOURCE_COEFFS + self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL + self.density.coeffs = np.array([[1257.5],[-0.475350]]) + + self.specific_heat.source = self.specific_heat.SOURCE_COEFFS + self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL + self.specific_heat.coeffs = np.array([[2.783],[0.0023]])*1e3 + + key = 'Cond' + def funcCond(T,x): + T = (T-self.Tbase) + if T <= 20: return 0.001840*T+0.4980 + else: return 0.001000*T+0.5140 + self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(dataID=key,func=funcCond,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.conductivity.DEBUG) + self.conductivity.source = self.conductivity.SOURCE_EQUATION + funcCond = None + + key = 'Mu' + def funcMu(T,x): + T = (T-self.Tbase) + if T <= 20: return 0.11100*np.exp(408.17/(T+123.15)) + else: return T*(0.000295*T - 0.0441)+2.6836 + self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG) + self.viscosity.source = self.viscosity.SOURCE_EQUATION + funcMu = None + + +class HyCool40(PureData,DigitalData): + def __init__(self): + DigitalData.__init__(self) + PureData.__init__(self) + + self.name = "HY40" + self.description = "HYCOOL 40, Potassium formate" + self.reference = "Hydro Chemicals" + + self.Tmax = 20 + 273.15 + self.Tmin = -40 + 273.15 + self.TminPsat = self.Tmax + self.Tbase = 0.00 + 273.15 + self.temperature.data = self.getTrange() + + self.density.source = self.density.SOURCE_COEFFS + self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL + self.density.coeffs = np.array([[1304.5],[-0.512290]]) + + self.specific_heat.source = self.specific_heat.SOURCE_COEFFS + self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL + self.specific_heat.coeffs = np.array([[2.646],[0.0023]])*1e3 + + self.conductivity.source = self.conductivity.SOURCE_COEFFS + self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL + self.conductivity.coeffs = np.array([[0.4826],[0.001730]]) + + key = 'Mu' + def funcMu(T,x): + T = (T-self.Tbase) + return 0.07830*np.exp(498.13/(T+130.25)) + self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG) + self.viscosity.source = self.viscosity.SOURCE_EQUATION + funcMu = None + + +class HyCool45(PureData,DigitalData): + def __init__(self): + DigitalData.__init__(self) + PureData.__init__(self) + + self.name = "HY45" + self.description = "HYCOOL 45, Potassium formate" + self.reference = "Hydro Chemicals" + + self.Tmax = 20 + 273.15 + self.Tmin = -45 + 273.15 + self.TminPsat = self.Tmax + self.Tbase = 0.00 + 273.15 + self.temperature.data = self.getTrange() + + self.density.source = self.density.SOURCE_COEFFS + self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL + self.density.coeffs = np.array([[1328.7],[-0.530754]]) + + self.specific_heat.source = self.specific_heat.SOURCE_COEFFS + self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL + self.specific_heat.coeffs = np.array([[2.578],[0.0023]])*1e3 + + self.conductivity.source = self.conductivity.SOURCE_COEFFS + self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL + self.conductivity.coeffs = np.array([[0.4750],[0.001674]]) + + key = 'Mu' + def funcMu(T,x): + T = (T-self.Tbase) + return 0.08990*np.exp(479.09/(T+126.55)) + self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG) + self.viscosity.source = self.viscosity.SOURCE_EQUATION + funcMu = None + + +class HyCool50(PureData,DigitalData): + def __init__(self): + DigitalData.__init__(self) + PureData.__init__(self) + + self.name = "HY50" + self.description = "HYCOOL 50, Potassium formate" + self.reference = "Hydro Chemicals" + + self.Tmax = 20 + 273.15 + self.Tmin = -50 + 273.15 + self.TminPsat = self.Tmax + self.Tbase = 0.00 + 273.15 + self.temperature.data = self.getTrange() + + self.density.source = self.density.SOURCE_COEFFS + self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL + self.density.coeffs = np.array([[1359.0],[-0.552300]]) + + self.specific_heat.source = self.specific_heat.SOURCE_COEFFS + self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL + self.specific_heat.coeffs = np.array([[2.498],[0.0023]])*1e3 + + self.conductivity.source = self.conductivity.SOURCE_COEFFS + self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL + self.conductivity.coeffs = np.array([[0.4660],[0.001610]]) + + key = 'Mu' + def funcMu(T,x): + T = (T-self.Tbase) + res = 0.0491*np.exp(581.12/(T+129.05)) + if T > -10: return res + 0.2 + else: return res + self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID=key,func=funcMu,x_in=self.temperature.data,y_in=self.concentration.data,DEBUG=self.viscosity.DEBUG) + self.viscosity.source = self.viscosity.SOURCE_EQUATION + funcMu = None + diff --git a/dev/incompressible_liquids/CPIncomp/PureFluids.py b/dev/incompressible_liquids/CPIncomp/PureFluids.py index f7ebf239..a0dce377 100644 --- a/dev/incompressible_liquids/CPIncomp/PureFluids.py +++ b/dev/incompressible_liquids/CPIncomp/PureFluids.py @@ -1,6 +1,6 @@ from __future__ import division, print_function import numpy as np -from CPIncomp.DataObjects import PureData +from CPIncomp.DataObjects import PureData , DigitalData class TherminolD12(PureData): """ @@ -336,3 +336,4 @@ class HC10(PureData): self.reference = "Dynalene data sheet" self.reshapeAll() + \ No newline at end of file diff --git a/dev/incompressible_liquids/CPIncomp/SecCoolFluids.py b/dev/incompressible_liquids/CPIncomp/SecCoolFluids.py index ab512598..7f6778c4 100644 --- a/dev/incompressible_liquids/CPIncomp/SecCoolFluids.py +++ b/dev/incompressible_liquids/CPIncomp/SecCoolFluids.py @@ -45,6 +45,41 @@ class SecCoolSolutionData(DigitalData): self.TminPsat = self.Tmax + try: + self.density.xData,self.density.yData,self.density.data = self.getArray(dataID="Rho") + self.density.source = self.density.SOURCE_DATA + except: + if self.density.DEBUG: print("Could not load {}".format(self.getFile("Rho"))) + pass + + try: + self.specific_heat.xData,self.specific_heat.yData,self.specific_heat.data = self.getArray(dataID='Cp') + while np.max(self.specific_heat.data[np.isfinite(self.specific_heat.data)])<500: # Expect values around 1e3 + self.specific_heat.data *= 1e3 + self.specific_heat.source = self.specific_heat.SOURCE_DATA + except: + if self.density.DEBUG: print("Could not load {}".format(self.getFile("Cp"))) + pass + + try: + self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(dataID='Cond') + while np.max(self.conductivity.data[np.isfinite(self.conductivity.data)])>10: # Expect value below 1 + self.conductivity.data *= 1e-3 + self.conductivity.source = self.conductivity.SOURCE_DATA + except: + if self.density.DEBUG: print("Could not load {}".format(self.getFile("Cond"))) + pass + + try: + self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID='Mu') + while np.max(self.viscosity.data[np.isfinite(self.viscosity.data)])>100: # Expect value below 10 + self.viscosity.data *= 1e-3 + self.viscosity.source = self.viscosity.SOURCE_DATA + except: + if self.density.DEBUG: print("Could not load {}".format(self.getFile("Mu"))) + pass + + # def interp(self, tOld, xOld, dataNew): # tCur = self.temperature.data # xCur = self.concentration.data @@ -93,9 +128,7 @@ class SecCoolSolutionData(DigitalData): errList = (ValueError, AttributeError, TypeError, RuntimeError) try: - self.density.xData,self.density.yData,self.density.data = self.getArray(dataID="Rho") self.density.coeffs = np.copy(std_coeffs) - self.density.source = self.density.SOURCE_DATA self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL self.density.fitCoeffs(self.Tbase,self.xbase) except errList as ve: @@ -103,11 +136,7 @@ class SecCoolSolutionData(DigitalData): pass try: - self.specific_heat.xData,self.specific_heat.yData,self.specific_heat.data = self.getArray(dataID='Cp') - while np.max(self.specific_heat.data[np.isfinite(self.specific_heat.data)])<500: # Expect values around 1e3 - self.specific_heat.data *= 1e3 self.specific_heat.coeffs = np.copy(std_coeffs) - self.specific_heat.source = self.specific_heat.SOURCE_DATA self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL self.specific_heat.fitCoeffs(self.Tbase,self.xbase) except errList as ve: @@ -115,84 +144,101 @@ class SecCoolSolutionData(DigitalData): pass try: - self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(data='Cond') - while np.max(self.conductivity.data[np.isfinite(self.conductivity.data)])>10: # Expect value below 1 - self.conductivity.data *= 1e-3 self.conductivity.coeffs = np.copy(std_coeffs) - self.conductivity.source = self.conductivity.SOURCE_DATA self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL self.conductivity.fitCoeffs(self.Tbase,self.xbase) except errList as ve: if self.conductivity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'conductivity',ve)) pass +# try: +# self.viscosity.coeffs = np.copy(std_coeffs) +# self.viscosity.type = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL +# self.viscosity.fitCoeffs(self.Tbase,self.xbase) +# except errList as ve: +# if self.viscosity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'viscosity',ve)) +# pass + try: - self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(data='Mu') - while np.max(self.viscosity.data[np.isfinite(self.viscosity.data)])>100: # Expect value below 10 - self.viscosity.data *= 1e-3 - self.viscosity.coeffs = np.copy(std_coeffs) - self.viscosity.source = self.viscosity.SOURCE_DATA - self.viscosity.type = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL - self.viscosity.fitCoeffs(self.Tbase,self.xbase) + tried = False + if len(self.viscosity.yData)==1:# and np.isfinite(fluidObject.viscosity.data).sum()<10: + self.viscosity.coeffs = np.array([+7e+2, -6e+1, +1e+1]) + self.viscosity.type = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL + self.viscosity.fitCoeffs(self.Tbase,self.xbase) + if self.viscosity.coeffs==None or IncompressibleFitter.allClose(self.viscosity.coeffs, np.array([+7e+2, -6e+1, +1e+1])): # Fit failed + tried = True + if len(self.viscosity.yData)>1 or tried: + self.viscosity.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5)) + self.viscosity.type = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL + self.viscosity.fitCoeffs(self.Tbase,self.xbase) except errList as ve: if self.viscosity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'viscosity',ve)) pass + # reset data for getArray and read special files if self.xid!=self.ifrac_pure and self.xid!=self.ifrac_undefined: allowNegativeData_org = self.allowNegativeData self.allowNegativeData = True - x,_,z = self.getArray(dataID='TFreeze') - self.T_freeze.yData = (x - 273.15) / 100.0 - self.T_freeze.xData = [0.0] - if np.min(z)<150: z += 273.15 - self.T_freeze.data = z.T - try: - self.T_freeze.source = self.T_freeze.SOURCE_DATA - if np.isfinite(self.T_freeze.data).sum()<2: - self.T_freeze.coeffs = np.array([+7e+6, +6e+4, +1e+1]) - self.T_freeze.type = self.T_freeze.INCOMPRESSIBLE_EXPONENTIAL - else: - self.T_freeze.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 2)) - self.T_freeze.type = self.T_freeze.INCOMPRESSIBLE_EXPPOLYNOMIAL - self.T_freeze.fitCoeffs(self.Tbase,self.xbase) + x,_,z = self.getArray(dataID='TFreeze') + self.T_freeze.yData = (x - 273.15) / 100.0 + self.T_freeze.xData = [0.0] + if np.min(z)<150: z += 273.15 + self.T_freeze.data = z.T + try: + self.T_freeze.source = self.T_freeze.SOURCE_DATA + if np.isfinite(self.T_freeze.data).sum()<2: + self.T_freeze.coeffs = np.array([+7e+6, +6e+4, +1e+1]) + self.T_freeze.type = self.T_freeze.INCOMPRESSIBLE_EXPONENTIAL + else: + self.T_freeze.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 2)) + self.T_freeze.type = self.T_freeze.INCOMPRESSIBLE_EXPPOLYNOMIAL + self.T_freeze.fitCoeffs(self.Tbase,self.xbase) + except errList as ve: + if self.T_freeze.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"T_freeze",ve)) + pass except errList as ve: - if self.T_freeze.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"T_freeze",ve)) + if self.T_freeze.DEBUG: print("{0}: Could not load {1} data: {2}".format(self.name,"T_freeze",ve)) pass # reset data for getArray again self.allowNegativeData = allowNegativeData_org - x,_,z = self.getArray(dataID='Vol2Mass') - massData = z[:,0] /100.0 - volData = (x - 273.15) /100.0 - if self.xid==self.ifrac_volume: - _,_,self.mass2input.data = IncompressibleFitter.shapeArray(volData,axs=1) - self.mass2input.xData = [0.0] - self.mass2input.yData = massData - try: - self.mass2input.coeffs = np.copy(std_coeffs) - self.mass2input.source = self.mass2input.SOURCE_DATA - self.mass2input.type = self.mass2input.INCOMPRESSIBLE_POLYNOMIAL - self.mass2input.fitCoeffs(self.Tbase,self.xbase) - except errList as ve: - if self.mass2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"mass2input",ve)) - pass - elif self.xid==self.ifrac_mass: - _,_,self.volume2input.data = IncompressibleFitter.shapeArray(massData,axs=1) - self.volume2input.xData = [0.0] - self.volume2input.yData = volData - try: - self.volume2input.coeffs = np.copy(std_coeffs) - self.volume2input.source = self.volume2input.SOURCE_DATA - self.volume2input.type = self.volume2input.INCOMPRESSIBLE_POLYNOMIAL - self.volume2input.fitCoeffs(self.Tbase,self.xbase) - except errList as ve: - if self.volume2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"volume2input",ve)) - pass - else: - raise ValueError("Unknown xid specified.") + try: + x,_,z = self.getArray(dataID='Vol2Mass') + massData = z[:,0] /100.0 + volData = (x - 273.15) /100.0 + + if self.xid==self.ifrac_volume: + _,_,self.mass2input.data = IncompressibleFitter.shapeArray(volData,axs=1) + self.mass2input.xData = [0.0] + self.mass2input.yData = massData + try: + self.mass2input.coeffs = np.copy(std_coeffs) + self.mass2input.source = self.mass2input.SOURCE_DATA + self.mass2input.type = self.mass2input.INCOMPRESSIBLE_POLYNOMIAL + self.mass2input.fitCoeffs(self.Tbase,self.xbase) + except errList as ve: + if self.mass2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"mass2input",ve)) + pass + elif self.xid==self.ifrac_mass: + _,_,self.volume2input.data = IncompressibleFitter.shapeArray(massData,axs=1) + self.volume2input.xData = [0.0] + self.volume2input.yData = volData + try: + self.volume2input.coeffs = np.copy(std_coeffs) + self.volume2input.source = self.volume2input.SOURCE_DATA + self.volume2input.type = self.volume2input.INCOMPRESSIBLE_POLYNOMIAL + self.volume2input.fitCoeffs(self.Tbase,self.xbase) + except errList as ve: + if self.volume2input.DEBUG: print("{0}: Could not fit {1} coefficients: {2}".format(self.name,"volume2input",ve)) + pass + else: + raise ValueError("Unknown xid specified.") + except errList as ve: + if self.mass2input.DEBUG or self.volume2input.DEBUG: print("{0}: Could not load {1} data: {2}".format(self.name,"Vol2Mass",ve)) + pass # Redefine some functions to avoid data loss @@ -331,6 +377,13 @@ class SecCoolSolutionData(DigitalData): sec += [SecCoolSolutionData(sFile='Dowtherm Q' ,sFolder='xPure',name='DowQ2',desc='Dowtherm Q, Diphenylethane/alkylated aromatics' ,ref='Dow Chemicals, SecCool software')] print(", {0}".format(sec[-1].name), end="") + sec += [SecCoolIceData(sFile='IceEA' ,sFolder='xMass',name='IceEA',desc='Ice slurry with Ethanol' ,ref='Danish Technological Institute, SecCool software')] + print(", {0}".format(sec[-1].name), end="") + sec += [SecCoolIceData(sFile='IceNA' ,sFolder='xMass',name='IceNA',desc='Ice slurry with NaCl' ,ref='Danish Technological Institute, SecCool software')] + print(", {0}".format(sec[-1].name), end="") + sec += [SecCoolIceData(sFile='IcePG' ,sFolder='xMass',name='IcePG',desc='Ice slurry with Propylene Glycol' ,ref='Danish Technological Institute, SecCool software')] + print(", {0}".format(sec[-1].name), end="") + sec += [ThermogenVP1869()] print(", {0}".format(sec[-1].name), end="") sec += [Freezium()] @@ -346,13 +399,107 @@ class SecCoolSolutionData(DigitalData): sec += [AS55()] print(", {0}".format(sec[-1].name), end="") - - print(" ... done") return sec + +class SecCoolIceData(SecCoolSolutionData): + """ + A base class that can be fed with a fluid ID from SecCool + to read data files sitting in data/SecCool/xTables. + """ + def __init__(self,sFile=None,sFolder=None,name=None,desc=None,ref='Danish Technological Institute, SecCool software'): + SecCoolSolutionData.__init__(self,sFile=sFile,sFolder=sFolder,name=name,desc=desc,ref=ref) + + #self.density.xData,self.density.yData,self.density.data = self.getArray(dataID="Rho") + #self.density.source = self.density.SOURCE_DATA + + self.specific_heat.xData,self.specific_heat.yData,self.specific_heat.data = self.getArray(dataID='Hfusion') + self.specific_heat.source = self.specific_heat.SOURCE_DATA + + #self.conductivity.xData,self.conductivity.yData,self.conductivity.data = self.getArray(dataID='Cond') + #self.conductivity.source = self.conductivity.SOURCE_DATA + + #self.viscosity.xData,self.viscosity.yData,self.viscosity.data = self.getArray(dataID='Mu') + #self.viscosity.source = self.viscosity.SOURCE_DATA + + +# def fitFluid(self): +# if self.Tbase==None: +# self.Tbase = (self.Tmin + self.Tmax) / 2.0 +# if self.xbase==None: +# self.xbase = (self.xmin + self.xmax) / 2.0 +# +# std_coeffs = np.zeros((4,6)) +# errList = (ValueError, AttributeError, TypeError, RuntimeError) +# +# try: +# self.density.coeffs = np.copy(std_coeffs) +# self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL +# self.density.fitCoeffs(self.Tbase,self.xbase) +# except errList as ve: +# if self.density.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'density',ve)) +# pass +# +# try: +# self.specific_heat.coeffs = np.copy(std_coeffs) +# self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL +# self.specific_heat.fitCoeffs(self.Tbase,self.xbase) +# except errList as ve: +# if self.specific_heat.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'specific heat',ve)) +# pass +# +# try: +# self.conductivity.coeffs = np.copy(std_coeffs) +# self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL +# self.conductivity.fitCoeffs(self.Tbase,self.xbase) +# except errList as ve: +# if self.conductivity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'conductivity',ve)) +# pass +# +# try: +# self.viscosity.coeffs = np.copy(std_coeffs) +# self.viscosity.type = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL +# self.viscosity.fitCoeffs(self.Tbase,self.xbase) +# except errList as ve: +# if self.viscosity.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(self.name,'viscosity',ve)) +# pass + + + # Redefine some functions to avoid data loss + def getFile(self, data): + return os.path.join(os.path.dirname(__file__), 'data','SecCool','xTables', self.sFolder, self.sFile+"_"+data+".csv") + + + def getFromFile(self, data): + fullPath = self.getFile(data) + content = np.loadtxt(fullPath,dtype=type('string'),delimiter=",",skiprows=3) + r,c,res = IncompressibleFitter.shapeArray(content) +# numbers = res.astype(np.float) + numbers = np.zeros((r,c)) + for i in range(r): + for j in range(c): + nu = np.NAN + try: + nu = np.float(res[i,j]) + if i==0: nu *= 1e-2 # Percent to fraction + if not self.allowNegativeData and nu<0: + nu = np.NAN # invalid entries + except (ValueError, TypeError) as ve: + if False: print("Could not convert entry: {0}".format(ve)) + if i==0: nu = 0.0 # Dummy for tables without concentration (TFreeze and Vol2Mass) + pass + numbers[i,j] = nu + if numbers[1,0]>numbers[-1,0]: + numbers[1:,:] = numbers[1:,:][::-1,:] + if numbers[0,1]>numbers[0,-1]: + numbers[:,1:] = numbers[:,1:][:,::-1] + return numbers + + + class ThermogenVP1869(PureData,DigitalData): """ Source: SecCool Software @@ -731,361 +878,7 @@ class AS55(PureData,DigitalData): - - -if __name__ == '__main__': - - print("The main method converts the AspenTemper coefficients from SecCool") - names = [] - descs = [] - refs = [] - - Tmax = [] - Tmin = [] - TminPsat = [] - Tbase = [] - - cond = [] - cp = [] - mu = [] - rho = [] - - ID = 10 - names.append("AS{0}".format(ID)) - descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID)) - refs.append ("Aspen Petroleum AB, SecCool software") - - Tmax.append("30 + 273.15") - Tmin.append("-{0}".format(ID)) - TminPsat.append("self.Tmax") - Tbase.append("0.00 + 273.15") - - cond.append ("0.001483*(273+T)+0.1094") - cp.append ("1000*(3.54183+T*(0.00201-0.0000132589*T))") - mu.append ("2.80154921+T*(-0.10631376+T*(0.00235381-0.0000211481*T))") - rho.append ("1090+T*(-0.2+T*(1.66533E-16-1.89735E-18*T))") - - ID = 20 - names.append("AS{0}".format(ID)) - descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID)) - refs.append ("Aspen Petroleum AB, SecCool software") - - Tmax.append("30 + 273.15") - Tmin.append("-{0}".format(ID)) - TminPsat.append("self.Tmax") - Tbase.append("0.00 + 273.15") - - cond.append ("0.001342*(273+T)+0.1144") - cp.append ("1000*(3.26252+T*(0.00286-0.0000122673*T))") - mu.append ("(0.97244+(7.14199*exp((-20-T)/18.6014)))") - rho.append ("1147+T*(-0.22142857+T*(-0.00142857+2.98156E-19*T))") - - ID = 30 - names.append("AS{0}".format(ID)) - descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID)) - refs.append ("Aspen Petroleum AB, SecCool software") - - Tmax.append("30 + 273.15") - Tmin.append("-{0}".format(ID)) - TminPsat.append("self.Tmax") - Tbase.append("0.00 + 273.15") - - cond.append ("0.001256*(273+T)+0.1175") - cp.append ("1000*(3.07504+T*(0.00299-0.0000270232*T))") - mu.append ("(1.30143+(16.01368*exp((-30-T)/16.69989)))") - rho.append ("1183.85930736+T*(-0.31085859+T*(-0.00103896+0.0000176768*T))") - - - ID = 40 - names.append("AS{0}".format(ID)) - descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID)) - refs.append ("Aspen Petroleum AB, SecCool software") - - Tmax.append("30 + 273.15") - Tmin.append("-{0}".format(ID)) - - cond.append ("0.001099*(273+T)+0.1433") - cp.append ("1000*(2.97788+T*(0.00228-0.0000387227*T))") - mu.append ("(39.11536*exp((-40-T)/9.99495)+(12.41564*exp((-40-T)/38.45577)))") - rho.append ("1214.83982684+T*(-0.37819865+T*(-0.00109307+0.000016835*T))") - - ID = 55 - names.append("AS{0}".format(ID)) - descs.append("Aspen Temper -{0}, Potassium acetate/formate".format(ID)) - refs.append ("Aspen Petroleum AB, SecCool software") - - Tmax.append("30 + 273.15") - Tmin.append("-{0}".format(ID)) - TminPsat.append("self.Tmax") - Tbase.append("0.00 + 273.15") - - cond.append ("(273+T)*(0.000002287*(273+T)-0.0003108)+0.3402") - cp.append ("1000*(2.83985+T*(0.00229-0.0000248618*T))") - mu.append ("(317.40673*exp((-55-T)/7.24125)+(51.22151*exp((-55-T)/26.28052)))") - rho.append ("1249.7534665+T*(-0.47629615+T*(-0.00117189+0.0000198824*T))") - - spa = " " - - from sympy import symbols, expand - x, T = symbols('x T') - - for i in range(len(names)): - print("class {0}(PureData):\n def __init__(self):".format(names[i])) - print("{0}PureData.__init__(self)".format(spa)) - print("{0}self.name = \"{1}\"".format(spa,names[i])) - print("{0}self.description = \"{1}\"".format(spa,descs[i])) - print("{0}self.reference = \"{1}\"".format(spa,refs[i])) - print("") - print("{0}self.Tmax = {1}".format(spa,Tmax[i])) - print("{0}self.Tmin = {1}".format(spa,Tmin[i])) - print("{0}self.TminPsat = self.Tmax".format(spa)) - print("{0}self.Tbase = 0.00 + 273.15".format(spa)) - print("") - print("{0}self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL".format(spa)) - print("{0}self.density.coeffs = np.array([{1}])[::-1]".format(spa,expand(rho[i]))) - print("") - print("{0}self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL".format(spa)) - print("{0}self.specific_heat.coeffs = np.array([{1}])[::-1]".format(spa,expand(cp[i]))) - print("") - print("{0}self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL".format(spa)) - print("{0}self.conductivity.coeffs = np.array([{1}])[::-1]".format(spa,expand(cond[i]))) - print("") - print("{0}self.viscosity.type = self.viscosity.INCOMPRESSIBLE_EXPPOLYNOMIAL".format(spa)) - print("{0}self.viscosity.coeffs = np.array([{1}])[::-1]".format(spa,expand(mu[i]))) - print("") - print("{0}def fitFluid(self):".format(" ")) - print("{0}key = \"Mu\"".format(spa)) - print("{0}def funcMu(T,x):".format(spa)) - print("") - print("") - print("{0}funcMu=None".format(spa)) - print("") - print("") - print("") - print("") - -# PureData.__init__(self) -# self.name = "AS10" -# self.description = "Aspen Temper -10, Potassium acetate/formate" -# self.reference = "Aspen Petroleum AB, SecCool software" -# -# self.Tmax = 30 + 273.15 -# self.Tmin = -10 + 273.15 -# self.TminPsat = self.Tmax -# -# self.Tbase = 0.00 + 273.15 -# -# self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL -# self.density.coeffs = np.array([[945.5454545],[-1.054545455]]) -# -# self.specific_heat.type = self.specific_heat.INCOMPRESSIBLE_POLYNOMIAL -# self.specific_heat.coeffs = np.array([[2.322218182],[0.003843636]])*1000 -# -# self.conductivity.type = self.conductivity.INCOMPRESSIBLE_POLYNOMIAL -# self.conductivity.coeffs = np.array([[0.001483*273.15+0.1094],[0.001483]]) -# -# self.temperature.data = self.getTrange() -# self.concentration.data = np.array([0]) # mass fraction - - -# -# -#{ THyCool_20 } -# -#constructor THyCool_20.Create; -#begin -# inherited; -# TFreeze := -20; -# FluidType := 'Potassium Formate'; -# TradeName := 'HYCOOL 20'; -# Reference := 'Hydro Chemicals'; -# TMin := -20; -# TMax := 50; -#end; -# -#function THyCool_20.GetCondT(T: Double): Double; -#begin -# if T <= 20 then -# Result := 0.001978*T+0.5172 -# else -# Result := 0.001005*T+0.5368; -#end; -# -#function THyCool_20.GetCpT(T: Double): Double; -#begin -# Result := 1000*(0.0023*T+2.955); -#end; -# -#function THyCool_20.GetMuT(T: Double): Double; -#begin -# if T <= 20 then -# Result := 0.07190*exp(524.75/(T+142.05)) -# else -# Result := T*(0.0005524*T - 0.06281)+2.8536; -#end; -# -#function THyCool_20.GetRhoT(T: Double): Double; -#begin -# Result := -0.429180*T+1202.2; -#end; -# -#{ THyCool_30 } -# -#constructor THyCool_30.Create; -#begin -# inherited; -# TFreeze := -30; -# FluidType := 'Potassium Formate'; -# TradeName := 'HYCOOL 30'; -# Reference := 'Hydro Chemicals'; -# TMin := -30; -# TMax := 50; -#end; -# -#function THyCool_30.GetCondT(T: Double): Double; -#begin -# if T <= 20 then -# Result := 0.001840*T+0.4980 -# else -# Result := 0.001000*T+0.514; -#end; -# -#function THyCool_30.GetCpT(T: Double): Double; -#begin -# Result := 1000*(0.0023*T+2.783); -#end; -# -#function THyCool_30.GetMuT(T: Double): Double; -#begin -# if T <= 20 then -# Result := 0.11100*exp(408.17/(T+123.15)) -# else -# Result := T*(0.000295*T - 0.0441)+2.6836; -#end; -# -#function THyCool_30.GetRhoT(T: Double): Double; -#begin -# Result := -0.475350*T+1257.5; -#end; -# -#{ THyCool_40 } -# -#constructor THyCool_40.Create; -#begin -# inherited; -# TFreeze := -40; -# FluidType := 'Potassium Formate'; -# TradeName := 'HYCOOL 40'; -# Reference := 'Hydro Chemicals'; -# TMin := -40; -# TMax := 20; -#end; -# -#function THyCool_40.GetCondT(T: Double): Double; -#begin -# Result := 0.001730*T+0.4826; -#end; -# -#function THyCool_40.GetCpT(T: Double): Double; -#begin -# Result := 1000*(0.0023*T+2.646); -#end; -# -#function THyCool_40.GetMuT(T: Double): Double; -#begin -# Result := 0.07830*exp(498.13/(T+130.25)); -#end; -# -#function THyCool_40.GetRhoT(T: Double): Double; -#begin -# Result := -0.512290*T+1304.5; -#end; -# -#{ THyCool_45 } -# -#constructor THyCool_45.Create; -#begin -# inherited; -# TFreeze := -45; -# FluidType := 'Potassium Formate'; -# TradeName := 'HYCOOL 45'; -# Reference := 'Hydro Chemicals'; -# TMin := -40; -# TMax := 20; -#end; -# -#function THyCool_45.GetCondT(T: Double): Double; -#begin -# Result := 0.001674*T+0.4750; -#end; -# -#function THyCool_45.GetCpT(T: Double): Double; -#begin -# Result := 1000*(0.0023*T+2.578); -#end; -# -#function THyCool_45.GetMuT(T: Double): Double; -#begin -# Result := 0.08990*exp(479.09/(T+126.55)); -#end; -# -#function THyCool_45.GetRhoT(T: Double): Double; -#begin -# Result := -0.530754*T+1328.7; -#end; -# -#{ THyCool_50 } -# -#constructor THyCool_50.Create; -#begin -# inherited; -# TFreeze := -50; -# FluidType := 'Potassium Formate'; -# TradeName := 'HYCOOL 50'; -# Reference := 'Hydro Chemicals'; -# TMin := -50; -# TMax := 20; -#end; -# -#function THyCool_50.GetCondT(T: Double): Double; -#begin -# Result := 0.001610*T+0.4660; -#end; -# -#function THyCool_50.GetCpT(T: Double): Double; -#begin -# Result := 1000*(0.0023*T+2.498); -#end; -# -#function THyCool_50.GetMuT(T: Double): Double; -#begin -# Result := 0.0491*exp(581.12/(T+129.05)); -# if T > -10 then -# Result := Result+0.2; -#end; -# -#function THyCool_50.GetRhoT(T: Double): Double; -#begin -# Result := -0.552300*T+1359.0; -#end; -# -# -# -# -# -# -# -# -# -# -# -# -# -# - - - diff --git a/dev/incompressible_liquids/CPIncomp/WriterObjects.py b/dev/incompressible_liquids/CPIncomp/WriterObjects.py index 6a2646bf..c01e1571 100644 --- a/dev/incompressible_liquids/CPIncomp/WriterObjects.py +++ b/dev/incompressible_liquids/CPIncomp/WriterObjects.py @@ -5,8 +5,9 @@ import hashlib, os, json, sys import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from CPIncomp.DataObjects import SolutionData -from CPIncomp.BaseObjects import IncompressibleData +from CPIncomp.BaseObjects import IncompressibleData, IncompressibleFitter from matplotlib.patches import Rectangle +from matplotlib.ticker import MaxNLocator class SolutionDataWriter(object): """ @@ -20,8 +21,6 @@ class SolutionDataWriter(object): def fitAll(self, fluidObject=SolutionData()): - - if fluidObject.Tbase==None: fluidObject.Tbase = (fluidObject.Tmin + fluidObject.Tmax) / 2.0 @@ -71,10 +70,10 @@ class SolutionDataWriter(object): fluidObject.viscosity.setxyData(tData,xData) tried = False if len(fluidObject.viscosity.yData)==1:# and np.isfinite(fluidObject.viscosity.data).sum()<10: - fluidObject.viscosity.coeffs = np.array([+7e+2, -6e+1, +1e+1]) + fluidObject.viscosity.coeffs = np.array([+5e+2, -6e+1, +1e+1]) fluidObject.viscosity.type = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL fluidObject.viscosity.fitCoeffs(tBase,xBase) - if fluidObject.viscosity.coeffs==None or np.allclose(fluidObject.viscosity.coeffs, np.array([+7e+2, -6e+1, +1e+1])): # Fit failed + if fluidObject.viscosity.coeffs==None or IncompressibleFitter.allClose(fluidObject.viscosity.coeffs, np.array([+5e+2, -6e+1, +1e+1])): # Fit failed tried = True if len(fluidObject.viscosity.yData)>1 or tried: fluidObject.viscosity.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5)) @@ -86,13 +85,17 @@ class SolutionDataWriter(object): try: fluidObject.saturation_pressure.setxyData(tData,xData) + tried = False if len(fluidObject.saturation_pressure.yData)==1:# and np.isfinite(fluidObject.saturation_pressure.data).sum()<10: fluidObject.saturation_pressure.coeffs = np.array([-5e+3, +6e+1, -1e+1]) fluidObject.saturation_pressure.type = IncompressibleData.INCOMPRESSIBLE_EXPONENTIAL - else: + fluidObject.saturation_pressure.fitCoeffs(tBase,xBase) + if fluidObject.saturation_pressure.coeffs==None or IncompressibleFitter.allClose(fluidObject.saturation_pressure.coeffs, np.array([-5e+3, +6e+1, -1e+1])): # Fit failed + tried = True + if len(fluidObject.saturation_pressure.yData)>1 or tried: fluidObject.saturation_pressure.coeffs = np.zeros(np.round(np.array(std_coeffs.shape) * 1.5)) fluidObject.saturation_pressure.type = IncompressibleData.INCOMPRESSIBLE_EXPPOLYNOMIAL - fluidObject.saturation_pressure.fitCoeffs(tBase,xBase) + fluidObject.saturation_pressure.fitCoeffs(tBase,xBase) except errList as ve: if fluidObject.saturation_pressure.DEBUG: print("{0}: Could not fit polynomial {1} coefficients: {2}".format(fluidObject.name,'saturation pressure',ve)) pass @@ -307,7 +310,7 @@ class SolutionDataWriter(object): pos = np.isfinite(B_a) pos2 = (B_a>eps) - result = np.zeros_like(A_a) + result = np.ones_like(A_a)*np.NAN result[pos & pos2] = (A_a[pos & pos2]-B_a[pos & pos2])/B_a[pos & pos2] @@ -416,9 +419,7 @@ class SolutionDataWriter(object): zError[i,j]= func(tData[i],pData,xData[j]) zError = self.relError(zData, zError) * 1e2 - - if xFunction: axis = 1 - else: axis = 0 + ## Find the column with the largest single error #maxVal = np.amax(zError, axis=0) # largest error per column @@ -426,17 +427,21 @@ class SolutionDataWriter(object): ## Find the column with the largest total error #totVal = np.sum(zError, axis=0) # summed error per column #col2plot = np.argmax(totVal) # largest error in row - # Find the column with the largest average error - avgVal = np.average(zError, axis=axis) # summed error per column - set2plot = np.argmax(avgVal) # largest error in row + # Find the column with the largest average error if xFunction: + #avgVal = np.average(zError, axis=1) # summed error per column + #set2plot = np.argmax(avgVal) # largest error in row + set2plot = int(np.round(r/2.0)) tData = np.array([tData[set2plot]]) zData = zData[set2plot] zError= zError[set2plot] else: + #avgVal = np.average(zError, axis=0) # summed error per column + #set2plot = np.argmax(avgVal) # largest error in row + set2plot = int(np.round(c/2.0)) xData = np.array([xData[set2plot]]) zData = zData.T[set2plot] - zError= zError.T[set2plot] + zError= zError.T[set2plot] else: raise ValueError("You have to provide data and a fitted function.") @@ -472,6 +477,10 @@ class SolutionDataWriter(object): xFunc = xData pFunc = pData zFunc = None + zMiMa = None + xFree = xData + tFree = tData + zFree = None if func!=None: if len(tFunc)1: @@ -488,7 +497,30 @@ class SolutionDataWriter(object): zFunc = np.array(zFunc.flat) else: raise ValueError("Cannot plot non-flat arrays!") - + + if xFunction: + tMiMa = np.array([solObj.Tmin, solObj.Tmax]) + xMiMa = xFunc + else: + tMiMa = tFunc + xMiMa = np.array([solObj.xmin, solObj.xmax]) + + zMiMa = np.zeros((len(tMiMa),len(xMiMa))) + for i in range(len(tMiMa)): + for j in range(len(xMiMa)): + zMiMa[i,j] = func(tMiMa[i],pFunc,xMiMa[j]) + + if not xFunction: # add the freezing front + if solObj.T_freeze.type!=IncompressibleData.INCOMPRESSIBLE_NOT_SET: + cols = len(tMiMa) + conc = np.linspace(solObj.xmin, solObj.xmax, num=cols) + tFree = np.zeros_like(conc) + zFree = np.zeros_like(conc) + for i in range(cols): + tFree[i] = solObj.Tfreeze(10.0, p=pFunc, x=conc[i]) + zFree[i] = func(tFree[i],pFunc,conc[i]) + #zMiMa = np.hstack((zMiMa,temp.reshape((len(conc),1)))) + fitFormatter = {} fitFormatter['color'] = 'red' fitFormatter['ls'] = 'solid' @@ -500,12 +532,20 @@ class SolutionDataWriter(object): errorFormatter['alpha'] = 0.25 pData = None + pFree = None if xFunction: pData = xData pFunc = xFunc + pMiMa = xMiMa + zMiMa = zMiMa.T + #pFree = xFree + #zFree = zFree.T else: pData = tData - 273.15 - pFunc = tFunc - 273.15 + pFunc = tFunc - 273.15 + pMiMa = tMiMa - 273.15 + #zMiMa = zMiMa + pFree = tFree - 273.15 if zData!=None and axVal!=None: axVal.plot(pData, zData, label='data', **dataFormatter) @@ -513,6 +553,12 @@ class SolutionDataWriter(object): if zFunc!=None and axVal!=None: axVal.plot(pFunc, zFunc, label='function' , **fitFormatter) + if zMiMa!=None and axVal!=None: + axVal.plot(pMiMa, zMiMa, alpha=0.25, ls=':', color=fitFormatter["color"]) + + if zFree!=None and axVal!=None: + axVal.plot(pFree, zFree, alpha=0.25, ls=':', color=fitFormatter["color"]) + if zError!=None and axErr!=None: axErr.plot(pData, zError, label='error' , **errorFormatter) if solObj.xid!=solObj.ifrac_pure and not xFunction: @@ -671,7 +717,7 @@ class SolutionDataWriter(object): density_axis.set_ylabel(r'Density [$\mathdefault{kg/m^3\!}$]') density_axis.set_xlabel(tempLabel) density_error.set_ylabel(errLabel) - if solObj.density.source!=solObj.density.SOURCE_NOT_SET: + if solObj.density.source!=solObj.density.SOURCE_NOT_SET: self.plotValues(density_axis,density_error,solObj=solObj,dataObj=solObj.density,func=solObj.rho) else: raise ValueError("Density data has to be provided!") @@ -751,13 +797,15 @@ class SolutionDataWriter(object): #mass2input_axis = plt.subplot2grid((3,2), (2,0)) #volume2input_axis = plt.subplot2grid((3,2), (2,0)) - # Set a minimum error level + # Set a minimum error level and do some more formatting minAbsErrorScale = 0.05 # in per cent for a in fig.axes: if a.get_ylabel()==errLabel: mi,ma = a.get_ylim() if mi>-minAbsErrorScale: a.set_ylim(bottom=-minAbsErrorScale) if ma< minAbsErrorScale: a.set_ylim( top= minAbsErrorScale) + a.xaxis.set_major_locator(MaxNLocator(5)) + #a.yaxis.set_major_locator(MaxNLocator(7)) @@ -779,7 +827,7 @@ class SolutionDataWriter(object): table_axis.legend( legVal, legKey, - bbox_to_anchor=(0.0, 0.015, 1., 0.015), + bbox_to_anchor=(0.0, -0.025, 1., -0.025), ncol=len(legKey), mode="expand", borderaxespad=0., numpoints=1) #table_axis.legend(handles, labels, bbox_to_anchor=(0.0, -0.1), loc=2, ncol=3) diff --git a/dev/incompressible_liquids/all_incompressibles.py b/dev/incompressible_liquids/all_incompressibles.py index 6177e02e..f9f0d46a 100644 --- a/dev/incompressible_liquids/all_incompressibles.py +++ b/dev/incompressible_liquids/all_incompressibles.py @@ -10,22 +10,28 @@ if __name__ == '__main__': writer = SolutionDataWriter() doneObjs = [] -# # To debug single fluids -# from CPIncomp.SecCoolFluids import Freezium -# #solObjs = [SecCoolSolutionData(sFile='Melinder, Ammonia' ,sFolder='xMass',name='MAM2',desc='Melinder, Ammonia' ,ref='Melinder-BOOK-2010, SecCool software')] -# solObjs = [Freezium()] -# #solObjs[0].density.DEBUG = True -# #solObjs[0].specific_heat.DEBUG = True -# #solObjs[0].conductivity.DEBUG = True -# #solObjs[0].viscosity.DEBUG = True -# solObjs[0].T_freeze.DEBUG = True -# writer.fitSecCoolList(solObjs) -# # -# ##from CPIncomp.ExampleObjects import SecCoolExample -# ##solObjs = [SecCoolExample()] -# writer.writeFluidList(solObjs) -# writer.writeReportList(solObjs) -# sys.exit(0) + # To debug single fluids + from CPIncomp.SecCoolFluids import SecCoolSolutionData,SecCoolIceData + from CPIncomp.PureFluids import Therminol72 + #solObjs = [SecCoolSolutionData(sFile='Melinder, Ammonia' ,sFolder='xMass',name='MAM2',desc='Melinder, Ammonia' ,ref='Melinder-BOOK-2010, SecCool software')] + #solObjs += [SecCoolIceData(sFile='IceNA' ,sFolder='xMass',name='IceNA',desc='Ice slurry with NaCl' ,ref='Danish Technological Institute, SecCool software')] + #solObjs = [Freezium()] + #solObjs[0].density.DEBUG = True + #solObjs[0].specific_heat.DEBUG = True + #solObjs[0].conductivity.DEBUG = True + #solObjs[0].viscosity.DEBUG = True + #solObjs[0].T_freeze.DEBUG = True + #writer.fitSecCoolList(solObjs) + solObjs = [Therminol72()] + solObjs[0].viscosity.DEBUG=True + #solObjs[0].saturation_pressure.DEBUG=True + writer.fitFluidList([solObjs[-1]]) + # + ##from CPIncomp.ExampleObjects import SecCoolExample + ##solObjs = [SecCoolExample()] + writer.writeFluidList(solObjs) + writer.writeReportList(solObjs) + sys.exit(0) fluidObjs = getExampleNames(obj=True) examplesToFit = ["ExamplePure","ExampleSolution","ExampleDigital","ExampleDigitalPure"] diff --git a/src/Backends/Incompressible/IncompressibleFluid.cpp b/src/Backends/Incompressible/IncompressibleFluid.cpp index 1f111ddf..5e9d2d88 100644 --- a/src/Backends/Incompressible/IncompressibleFluid.cpp +++ b/src/Backends/Incompressible/IncompressibleFluid.cpp @@ -439,26 +439,21 @@ double IncompressibleFluid::T_u (double Umass, double p, double x){ * freezing point calculation. This is not necessarily * defined for all fluids, default values do not cause errors. */ bool IncompressibleFluid::checkT(double T, double p, double x) { - if (Tmin <= 0.) { - throw ValueError("Please specify the minimum temperature."); - } else if (Tmax <= 0.) { - throw ValueError("Please specify the maximum temperature."); - } else if ((Tmin > T) || (T > Tmax)) { - throw ValueError( + if (Tmin <= 0.) throw ValueError("Please specify the minimum temperature."); + if (Tmax <= 0.) throw ValueError("Please specify the maximum temperature."); + if ((Tmin > T) || (T > Tmax)) throw ValueError( format("Your temperature %f is not between %f and %f.", T, Tmin, Tmax)); + double TF = 0.0; + if (T_freeze.type!=IncompressibleData::INCOMPRESSIBLE_NOT_SET) { + TF = Tfreeze(p, x); + } + if ( T