Added tests and benchmarks for Vinet and BM4

burnman/eos/__init__.py

@@ -17,6 +17,8 @@
 from .modified_tait import MT
 from .hp import HP_TMT
 from .cork import CORK
+from .vinet import Vinet
+from .birch_murnaghan_4th import BM4
 
 from .property_modifiers import calculate_property_modifications
 

misc/benchmarks/benchmark.py

@@ -10,8 +10,10 @@
 import burnman
 
 import burnman.eos.birch_murnaghan as bm
+import burnman.eos.birch_murnaghan_4th as bm4
 import burnman.eos.mie_grueneisen_debye as mgd
 import burnman.eos.slb as slb
+import burnman.eos.vinet as vinet
 import matplotlib.image as mpimg
 
 
@@ -58,6 +60,82 @@ def check_birch_murnaghan():
 
     plt.show()
 
+def check_birch_murnaghan_4th():
+    """
+    Recreates Stixrude and Lithgow-Bertelloni (2005) Figure 1, bulk and shear modulus without thermal corrections
+    """
+    plt.close()
+
+    #make a test mineral
+    test_mineral = burnman.Mineral()
+    test_mineral.params ={'name':'test',
+                        'V_0': 10.e-6,
+                        'K_0': 72.7e9,
+                        'Kprime_0':  4.14 ,
+                        'Kprime_prime_0':  -0.0484e-9,
+                        }
+
+    test_mineral.set_method('bm4')
+
+    pressure = np.linspace(0.,90.e9, 20)
+    volume = np.empty_like(pressure)
+
+    #calculate its static properties
+    for i in range(len(pressure)):
+        volume[i] = bm4.volume_fourth_order(pressure[i], test_mineral.params)/test_mineral.params.get('V_0')
+
+
+    #compare with figure 1
+    plt.plot(pressure/1.e9,volume)
+    fig1 = mpimg.imread('../../burnman/data/input_figures/Ahmad.png')
+    plt.imshow(fig1, extent=[0.,90.,.65,1.], aspect='auto')
+    plt.plot(pressure/1.e9 , volume,marker='o',color='r',linestyle='',label = 'BM4')
+    plt.legend(loc='lower left')
+    plt.xlim(0.,90.)
+    plt.ylim(.65,1.)
+    plt.xlabel("Volume/V0")
+    plt.ylabel("Pressure (GPa)")
+    plt.title("Comparing with Figure 1 of Ahmad et al., (2012)")
+
+    plt.show()
+
+def check_vinet():
+    """
+    Recreates Stixrude and Lithgow-Bertelloni (2005) Figure 1, bulk and shear modulus without thermal corrections
+    """
+    plt.close()
+
+    #make a test mineral
+    test_mineral = burnman.Mineral()
+    test_mineral.params ={'name':'test',
+                    'V_0': 6.75e-6,
+                    'K_0': 163.4e9,
+                    'Kprime_0': 5.38,
+                        }
+
+    test_mineral.set_method('vinet')
+
+    pressure = np.linspace(17.7e9,300.e9, 20)
+    volume = np.empty_like(pressure)
+
+    #calculate its static properties
+    for i in range(len(pressure)):
+        volume[i] = vinet.volume(pressure[i], test_mineral.params)
+
+    #compare with figure 1
+    plt.plot(pressure/1.e9,volume/6.02e-7)
+    fig1 = mpimg.imread('../../burnman/data/input_figures/Dewaele.png')
+    plt.imshow(fig1, extent=[0.,300.,6.8,11.8], aspect='auto')
+    plt.plot(pressure/1.e9 , volume/6.02e-7,marker='o',color='r',linestyle='',label = 'Vinet Fit')
+    plt.legend(loc='lower left')
+    plt.xlim(0.,300.)
+    plt.ylim(6.8,11.8)
+    plt.ylabel("Volume (Angstroms^3/atom")
+    plt.xlabel("Pressure (GPa)")
+    plt.title("Comparing with Figure 1 of Dewaele et al., (2006)")
+
+    plt.show()
+
 def check_mgd_shim_duffy_kenichi():
     """
     Attemmpts to recreate Shim Duffy Kenichi (2002)
@@ -597,6 +675,8 @@ def check_averaging_3():
     check_averaging_2()
     check_averaging_3()
     check_birch_murnaghan()
+    check_birch_murnaghan_4th()
+    check_vinet()
     check_slb_fig7()
     check_slb_fig3()
     check_mgd_shim_duffy_kenichi()

tests/test_eos.py

@@ -34,6 +34,32 @@ def __init__(self):
             'q_0': 1.5,
             'eta_s_0': 2.8}
 
+class Fe_Dewaele(burnman.Mineral):
+    """
+    Dewaele et al., 2006, Physical Review Letters
+    """
+    def __init__(self):
+        self.params = {
+                'equation_of_state': 'vinet',
+                'V_0': 6.75e-6,
+                'K_0': 163.4e9,
+                'Kprime_0': 5.38,
+                'molar_mass': 0.055845,
+                'n': 1}
+class Liquid_Fe_Anderson(burnman.Mineral):
+    """
+    Anderson & Ahrens, 1994 JGR
+    """
+    def __init__(self):
+        self.params = {
+            'equation_of_state': 'bm4',
+            'V_0': 7.95626e-6,
+            'K_0': 109.7e9,
+            'Kprime_0': 4.66,\
+            'Kprime_prime_0': -0.043e-9,
+            'molar_mass': 0.055845,
+        }
+
 
 class eos(BurnManTest):
     def test_reference_values(self):
@@ -71,6 +97,41 @@ def test_reference_values(self):
             Grun_test = i.grueneisen_parameter(pressure, temperature, rock.params['V_0'], rock.params)
             self.assertFloatEqual(Grun_test, rock.params['grueneisen_0'])
 
+    def test_reference_values_noG(self):
+        #First test Vinet
+        rock = Fe_Dewaele()
+        pressure = 0.
+        temperature = 300.
+        eos = burnman.eos.Vinet()
+
+        Volume_test = eos.volume(pressure, temperature, rock.params)
+        self.assertFloatEqual(Volume_test, rock.params['V_0'])
+        Kt_test = eos.isothermal_bulk_modulus(pressure, 300., rock.params['V_0'], rock.params)
+        self.assertFloatEqual(Kt_test, rock.params['K_0'])
+        # K_S is based on 0 reference temperature:
+        Kt_test = eos.isothermal_bulk_modulus(pressure, 0., rock.params['V_0'], rock.params)
+        K_test = eos.adiabatic_bulk_modulus(pressure, 0., rock.params['V_0'], rock.params)
+        self.assertFloatEqual(K_test, Kt_test)
+        Density_test = eos.density(rock.params['V_0'], rock.params)
+        self.assertFloatEqual(Density_test, rock.params['molar_mass'] / rock.params['V_0'])
+
+        #Now BM4
+        rock = Liquid_Fe_Anderson()
+        pressure = 0.
+        temperature = 300.
+        eos = burnman.eos.BM4()
+
+        Volume_test = eos.volume(pressure, temperature, rock.params)
+        self.assertFloatEqual(Volume_test, rock.params['V_0'])
+        Kt_test = eos.isothermal_bulk_modulus(pressure, 300., rock.params['V_0'], rock.params)
+        self.assertFloatEqual(Kt_test, rock.params['K_0'])
+        # K_S is based on 0 reference temperature:
+        Kt_test = eos.isothermal_bulk_modulus(pressure, 0., rock.params['V_0'], rock.params)
+        K_test = eos.adiabatic_bulk_modulus(pressure, 0., rock.params['V_0'], rock.params)
+        self.assertFloatEqual(K_test, Kt_test)
+        Density_test = eos.density(rock.params['V_0'], rock.params)
+        self.assertFloatEqual(Density_test, rock.params['molar_mass'] / rock.params['V_0'])
+