[cig-commits] [commit] master: extract debye test, cleanup (0f84018)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/burnman

On branch  : master
Link       : https://github.com/geodynamics/burnman/compare/ab0efb9b1ff17b6d893db6c4d457836c768cdd59...f76fd51b3334a16aae97d4dd7ab1f0204e1d2eb0

>---------------------------------------------------------------

commit 0f840189bd42a04871745d76189555117edf4ed5
Author: Timo Heister <timo.heister at gmail.com>
Date:   Thu Dec 11 00:28:08 2014 -0500

    extract debye test, cleanup


>---------------------------------------------------------------

0f840189bd42a04871745d76189555117edf4ed5
 burnman/debye.py      | 50 ----------------------------------------------
 misc/ref/debye.py.out |  1 +
 test.sh               |  4 ++--
 tests/debye.py        | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 58 insertions(+), 52 deletions(-)

diff --git a/burnman/debye.py b/burnman/debye.py
index 1a337d1..c0f6452 100644
--- a/burnman/debye.py
+++ b/burnman/debye.py
@@ -98,53 +98,3 @@ def heat_capacity_v(T,debye_T,n):
 
 
 
-if __name__ == "__main__":
-    import matplotlib.pyplot as plt
-    import time
-
-    def old_thermal(T, debye_T, n):
-        if T == 0:
-            return 0
-        return 3.*n*constants.R*T * debye_fn(debye_T/T)
-
-    def old_heat(T, debye_T, n):
-        if T ==0:
-            return 0
-        deb = integrate.quad( lambda x : pow(x,4.)*np.exp(x)/pow((np.exp(x)-1.),2.), 0.0, debye_T/T)
-        return 9.*n*constants.gas_constant*deb[0]/pow(debye_T/T,3.)
-
-    temperatures = np.linspace(100,5000, 10000)
-    Debye_T = 1000.
-    old = np.empty_like(temperatures)
-    start = time.clock()
-    for i in range(len(temperatures)):
-        old[i] = old_heat(temperatures[i], Debye_T, 1.0)
-    time_old = time.clock()-start
-
-    new = np.empty_like(temperatures)
-    start = time.clock()
-    for i in range(len(temperatures)):
-        new[i] = heat_capacity_v(temperatures[i], Debye_T, 1.0)
-    time_new = time.clock()-start
-
-    print "error %e"%np.linalg.norm((old-new)/new)
-    print "time old %g, time new %g"%(time_old,time_new)
-
-
-
-    temperatures = np.linspace(0,5000, 200)
-    vibrational_energy = np.empty_like(temperatures)
-    heat_capacity = np.empty_like(temperatures)
-    Debye_T = 1000.
-    for i in range(len(temperatures)):
-      vibrational_energy[i] = thermal_energy(temperatures[i], Debye_T, 1.0)
-      heat_capacity[i] = heat_capacity_v(temperatures[i], Debye_T, 1.0)
-
-    plt.subplot(121)
-    plt.plot(temperatures, vibrational_energy)
-    plt.subplot(122)
-    plt.plot(temperatures, heat_capacity)
-    plt.show()
-
-
-
diff --git a/misc/ref/debye.py.out b/misc/ref/debye.py.out
new file mode 100644
index 0000000..a2a2ce1
--- /dev/null
+++ b/misc/ref/debye.py.out
@@ -0,0 +1 @@
+error 7.013244e-14
diff --git a/test.sh b/test.sh
index 2409ef6..af1c2d4 100755
--- a/test.sh
+++ b/test.sh
@@ -27,6 +27,7 @@ else
   sed -i '/UserWarning: findfont: Font family/d' $t.tmp #remove font warning crap
   sed -i '/tight_layout : falling back to Agg renderer/d' $t.tmp #remove font warning crap
   sed -i '/cannot be converted with the encoding. Glyph may be wrong/d' $t.tmp #remove font warning crap
+  sed -i '/time old .* time new/d' $t.tmp #remove timing from tests/debye.py
 
   (numdiff -r 1e-6 -s ' \t\n[]' -a 1e-6 -q $t.tmp ../misc/ref/$t.out >/dev/null && rm $t.tmp && echo "$t ... ok"
   ) || { 
@@ -45,10 +46,9 @@ echo "*** running test suite..."
 cd tests
 python tests.py || exit 1
 testit "benchmark.py"
+testit "debye.py"
 cd ..
 
-python burnman/partitioning.py || exit 1
-
 cd misc
 echo "gen_doc..."
 python gen_doc.py >/dev/null || exit 1
diff --git a/tests/debye.py b/tests/debye.py
new file mode 100644
index 0000000..8b322a0
--- /dev/null
+++ b/tests/debye.py
@@ -0,0 +1,55 @@
+import os, sys, numpy as np, matplotlib.pyplot as plt
+if not os.path.exists('burnman') and os.path.exists('../burnman'):
+    sys.path.insert(1,os.path.abspath('..')) 
+sys.path.insert(1,os.path.abspath('.')) 
+import burnman
+
+import scipy.integrate
+import time
+
+def old_thermal(T, debye_T, n):
+    if T == 0:
+        return 0
+    return 3.*n*constants.R*T * debye_fn(debye_T/T)
+
+def old_heat(T, debye_T, n):
+    if T ==0:
+        return 0
+    deb = scipy.integrate.quad( lambda x : pow(x,4.)*np.exp(x)/pow((np.exp(x)-1.),2.), 0.0, debye_T/T)
+    return 9.*n*burnman.constants.gas_constant*deb[0]/pow(debye_T/T,3.)
+
+temperatures = np.linspace(100,5000, 10000)
+Debye_T = 1000.
+old = np.empty_like(temperatures)
+start = time.clock()
+for i in range(len(temperatures)):
+    old[i] = old_heat(temperatures[i], Debye_T, 1.0)
+time_old = time.clock()-start
+
+new = np.empty_like(temperatures)
+start = time.clock()
+for i in range(len(temperatures)):
+    new[i] = burnman.debye.heat_capacity_v(temperatures[i], Debye_T, 1.0)
+time_new = time.clock()-start
+
+print "error %e"%np.linalg.norm((old-new)/new)
+print "time old %g, time new %g"%(time_old,time_new)
+
+
+
+temperatures = np.linspace(0,5000, 200)
+vibrational_energy = np.empty_like(temperatures)
+heat_capacity = np.empty_like(temperatures)
+Debye_T = 1000.
+for i in range(len(temperatures)):
+  vibrational_energy[i] = burnman.debye.thermal_energy(temperatures[i], Debye_T, 1.0)
+  heat_capacity[i] = burnman.debye.heat_capacity_v(temperatures[i], Debye_T, 1.0)
+
+plt.subplot(121)
+plt.plot(temperatures, vibrational_energy)
+plt.subplot(122)
+plt.plot(temperatures, heat_capacity)
+plt.show()
+
+
+