Added more tests for partial derivatives, and a preliminary class for…

… fitting data to various alpha functions

tests/test_eos.py

@@ -21,8 +21,10 @@
 SOFTWARE.'''
 
 from numpy.testing import assert_allclose
+import numpy as np
 import pytest
 from thermo.eos import *
+from thermo.utils import allclose_variable
 from scipy.misc import derivative
 from scipy.constants import R
 from math import log, exp, sqrt
@@ -678,3 +680,74 @@ def test_TWUSRK_quick():
     assert_allclose(eos.T, 900)
 
 
+@pytest.mark.slow
+def test_fuzz_dV_dT_and_d2V_dT2_derivatives():
+    from thermo import eos
+    eos_list = list(eos.__all__); eos_list.remove('GCEOS')
+    eos_list.remove('ALPHA_FUNCTIONS'); eos_list.remove('VDW')
+
+    phase_extensions = {True: '_l', False: '_g'}
+    derivative_bases_dV_dT = {0:'V', 1:'dV_dT', 2:'d2V_dT2'}
+
+    def dV_dT(T, P, eos, order=0, phase=True, Tc=507.6, Pc=3025000., omega=0.2975):
+        eos = globals()[eos_list[eos]](Tc=Tc, Pc=Pc, omega=omega, T=T, P=P)
+        phase_base = phase_extensions[phase]
+        attr = derivative_bases_dV_dT[order]+phase_base
+        return getattr(eos, attr)
+
+    x, y = [], []
+    for eos in range(len(eos_list)):
+        for T in np.linspace(.1, 1000, 50):
+            for P in np.logspace(np.log10(3E4), np.log10(1E6), 50):
+                T, P = float(T), float(P)
+                for phase in [True, False]:
+                    for order in [1, 2]:
+                        try:
+                            # If dV_dx_phase doesn't exist, will simply abort and continue the loop
+                            numer = derivative(dV_dT, T, dx=1E-4, args=(P, eos, order-1, phase))
+                            ana = dV_dT(T=T, P=P, eos=eos, order=order, phase=phase)
+                        except:
+                            continue
+                        x.append(numer)
+                        y.append(ana)
+    ans = allclose_variable(x, y, limits=[.009, .05, .65, .93],rtols=[1E-5, 1E-6, 1E-9, 1E-10])
+    if not ans:
+        raise Exception('dVdT numerical derivative failing!')
+
+
+@pytest.mark.slow
+def test_fuzz_dV_dP_and_d2V_dP2_derivatives():
+    from thermo import eos
+    eos_list = list(eos.__all__); eos_list.remove('GCEOS')
+    eos_list.remove('ALPHA_FUNCTIONS'); eos_list.remove('VDW')
+
+    phase_extensions = {True: '_l', False: '_g'}
+    derivative_bases_dV_dP = {0:'V', 1:'dV_dP', 2:'d2V_dP2'}
+
+    def dV_dP(P, T, eos, order=0, phase=True, Tc=507.6, Pc=3025000., omega=0.2975):
+        eos = globals()[eos_list[eos]](Tc=Tc, Pc=Pc, omega=omega, T=T, P=P)
+        phase_base = phase_extensions[phase]
+        attr = derivative_bases_dV_dP[order]+phase_base
+        return getattr(eos, attr)
+
+
+    x, y = [], []
+    for eos in range(len(eos_list)):
+        for T in np.linspace(.1, 1000, 50):
+            for P in np.logspace(np.log10(3E4), np.log10(1E6), 50):
+                T, P = float(T), float(P)
+                for phase in [True, False]:
+                    for order in [1, 2]:
+                        try:
+                            # If dV_dx_phase doesn't exist, will simply abort and continue the loop
+                            numer = derivative(dV_dP, P, dx=15., args=(T, eos, order-1, phase))
+                            ana = dV_dP(T=T, P=P, eos=eos, order=order, phase=phase)
+                        except:
+                            continue
+                        x.append(numer)
+                        y.append(ana)
+    ans = allclose_variable(x, y, limits=[.02, .04, .04, .05, .15, .45, .95],
+                            rtols=[1E-2, 1E-3, 1E-4, 1E-5, 1E-6, 1E-7, 1E-9])
+    if not ans:
+        raise Exception('Error')
+

tests/test_utils.py

@@ -437,4 +437,32 @@ def test_mu_JT():
     assert_allclose(mu_JT, 1.621956080529905e-05)
 
     with pytest.raises(Exception):
-        Joule_Thomson(T=390, V=0.00229754, Cp=153.235)
+        Joule_Thomson(T=390, V=0.00229754, Cp=153.235)
+
+
+
+def test_allclose_variable():
+    x = [2.7244322249597719e-08, 3.0105683900110473e-10, 2.7244124924802327e-08, 3.0105259397637556e-10, 2.7243929226310193e-08, 3.0104990272770901e-10, 2.7243666849384451e-08, 3.0104101821236015e-10, 2.7243433745917367e-08, 3.0103707421519949e-10]
+    y = [2.7244328304561904e-08, 3.0105753470546008e-10, 2.724412872417824e-08,  3.0105303055834564e-10, 2.7243914341030203e-08, 3.0104819238021998e-10, 2.7243684057561379e-08, 3.0104299541023674e-10, 2.7243436694839306e-08, 3.010374130526363e-10]
+
+    ans = allclose_variable(x, y, limits=[.0, .5], rtols=[1E-5, 1E-6])
+    if not ans:
+        raise Exception('Failure')
+
+    with pytest.raises(Exception):
+        ans = allclose_variable(x, y, limits=[.0, .1], rtols=[1E-5, 1E-6])
+        if not ans:
+            raise Exception('Failure')
+
+    with pytest.raises(Exception):
+        allclose_variable(x, y[1:], limits=[.0, .5], rtols=[1E-5, 1E-6])
+
+    with pytest.raises(Exception):
+        ans = allclose_variable(x, y, limits=[.0, .1])
+
+
+    x = [1,1,1,1,1,1,1,1,1]
+    y = [.9,.9,.9,.9,.9,.9,.9,.9, .9]
+    ans = allclose_variable(x, y, limits=[.0], atols=[.1])
+    if ans is False:
+        raise Exception('Failure on a unit test')

thermo/__init__.py

@@ -32,6 +32,8 @@
 from . import electrochem
 from . import elements
 from . import environment
+from . import eos
+from . import eos_mix
 from . import heat_capacity
 from . import identifiers
 from . import law
@@ -65,6 +67,8 @@
 from .electrochem import *
 from .elements import *
 from .environment import *
+from .eos import *
+from .eos_mix import *
 from .heat_capacity import *
 from .identifiers import *
 from .law import *
@@ -87,7 +91,7 @@
 
 
 __all__ = ['activity', 'chemical', 'combustion', 'critical',
- 'dipole', 'electrochem', 'elements', 'environment',
+ 'dipole', 'electrochem', 'elements', 'environment', 'eos', 'eos_mix',
  'heat_capacity',  'identifiers', 'law', 'lennard_jones',
  'miscdata',
  'permittivity', 'phase_change', 'reaction',
@@ -109,6 +113,8 @@
 __all__.extend(electrochem.__all__)
 __all__.extend(elements.__all__)
 __all__.extend(environment.__all__)
+__all__.extend(eos.__all__)
+__all__.extend(eos_mix.__all__)
 __all__.extend(heat_capacity.__all__)
 __all__.extend(identifiers.__all__)
 __all__.extend(law.__all__)

thermo/acentric.py

@@ -38,7 +38,7 @@
 omega_methods = ['PSRK', 'PD', 'YAWS', 'LK', 'DEFINITION']
 
 
-def omega(CASRN='', AvailableMethods=False, Method=None, IgnoreMethods=['LK', 'DEFINITION']):
+def omega(CASRN, AvailableMethods=False, Method=None, IgnoreMethods=['LK', 'DEFINITION']):
     r'''This function handles the retrieval of a chemical's acentric factor,
     `omega`, or its calculation from correlations or directly through the
     definition of acentric factor if possible. Requires a known boiling point,

thermo/combustion.py

@@ -25,7 +25,7 @@
 __all__ = ['Hcombustion']
 
 
-def Hcombustion(atoms=None, Hf=None, HfH2O=-285825, HfCO2=-393474,
+def Hcombustion(atoms, Hf=None, HfH2O=-285825, HfCO2=-393474,
                 HfSO2=-296800, HfBr2=30880, HfI2=62417, HfHCl=-92173,
                 HfHF=-272711, HfP4O10=-3009940, HfO2=0, HfN2=0):
     '''Calculates the heat of combustion, in J/mol.

thermo/critical.py

@@ -90,7 +90,7 @@
 Tc_methods = [IUPAC, MATTHEWS, CRC, PSRK, PD, YAWS, SURF]
 
 
-def Tc(CASRN='', AvailableMethods=False, Method=None, IgnoreMethods=[SURF]):
+def Tc(CASRN, AvailableMethods=False, Method=None, IgnoreMethods=[SURF]):
     r'''This function handles the retrieval of a chemical's critical
     temperature. Lookup is based on CASRNs. Will automatically select a data
     source to use if no Method is provided; returns None if the data is not
@@ -275,7 +275,7 @@ def list_methods():
 Pc_methods = [IUPAC, MATTHEWS, CRC, PSRK, PD, YAWS, SURF]
 
 
-def Pc(CASRN='', AvailableMethods=False, Method=None, IgnoreMethods=[SURF]):
+def Pc(CASRN, AvailableMethods=False, Method=None, IgnoreMethods=[SURF]):
     r'''This function handles the retrieval of a chemical's critical
     pressure. Lookup is based on CASRNs. Will automatically select a data
     source to use if no Method is provided; returns None if the data is not
@@ -459,7 +459,7 @@ def list_methods():
 Vc_methods = [IUPAC, MATTHEWS, CRC, PSRK, YAWS, SURF]
 
 
-def Vc(CASRN='', AvailableMethods=False, Method=None, IgnoreMethods=[SURF]):
+def Vc(CASRN, AvailableMethods=False, Method=None, IgnoreMethods=[SURF]):
     r'''This function handles the retrieval of a chemical's critical
     volume. Lookup is based on CASRNs. Will automatically select a data
     source to use if no Method is provided; returns None if the data is not
@@ -634,7 +634,7 @@ def list_methods():
 Zc_methods = [IUPAC, MATTHEWS, CRC, PSRK, YAWS, COMBINED]
 
 
-def Zc(CASRN='', AvailableMethods=False, Method=None, IgnoreMethods=[COMBINED]):
+def Zc(CASRN, AvailableMethods=False, Method=None, IgnoreMethods=[COMBINED]):
     r'''This function handles the retrieval of a chemical's critical
     compressibility. Lookup is based on CASRNs. Will automatically select a
     data source to use if no Method is provided; returns None if the data is

thermo/dipole.py

@@ -48,7 +48,7 @@
 dipole_methods = [CCCBDB, MULLER, POLING]
 
 
-def dipole_moment(CASRN='', AvailableMethods=False, Method=None):
+def dipole_moment(CASRN, AvailableMethods=False, Method=None):
     r'''This function handles the retrieval of a chemical's dipole moment.
     Lookup is based on CASRNs. Will automatically select a data source to use
     if no Method is provided; returns None if the data is not available.

thermo/environment.py

@@ -47,7 +47,7 @@
 GWP_methods = [IPCC100, IPCC100SAR, IPCC20, IPCC500]
 
 
-def GWP(CASRN='', AvailableMethods=False, Method=None):
+def GWP(CASRN, AvailableMethods=False, Method=None):
     r'''This function handles the retrieval of a chemical's Global Warming
     Potential, relative to CO2. Lookup is based on CASRNs. Will automatically
     select a data source to use if no Method is provided; returns None if the
@@ -148,7 +148,7 @@ def list_methods():
 ODP_methods = [ODP2MAX, ODP1MAX, ODP2LOG, ODP1LOG, ODP2MIN, ODP1MIN, ODP2STR, ODP1STR]
 
 
-def ODP(CASRN='', AvailableMethods=False, Method=None):
+def ODP(CASRN, AvailableMethods=False, Method=None):
     r'''This function handles the retrieval of a chemical's Ozone Depletion
     Potential, relative to CFC-11 (trichlorofluoromethane). Lookup is based on
     CASRNs. Will automatically select a data source to use if no Method is
@@ -285,7 +285,7 @@ def list_methods():
 logP_methods = [SYRRES, CRC]
 
 
-def logP(CASRN='', AvailableMethods=False, Method=None):
+def logP(CASRN, AvailableMethods=False, Method=None):
     r'''This function handles the retrieval of a chemical's octanol-water
     partition coefficient. Lookup is based on CASRNs. Will automatically
     select a data source to use if no Method is provided; returns None if the