[cig-commits] [commit] add_thermodynamic_potentials: Attempts to locate tc-burnman discrepancies (9f037a6)

[cig_noreply at geodynamics.org]

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

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51

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

commit 9f037a68aa5a9d0f23ac72fb1de165f9ed0f435b
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Fri Aug 15 23:48:15 2014 +0200

    Attempts to locate tc-burnman discrepancies


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

9f037a68aa5a9d0f23ac72fb1de165f9ed0f435b
 burnman/modified_tait.py | 6 +++---
 burnman/testing.py       | 9 ++++-----
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index a21384a..da58861 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -56,19 +56,19 @@ def entropydisorder(n):
 def equilibrium_Q(Q, deltaS, P, T, params):
     n=params['BW_n']
     W=params['BW_W'] + P*params['BW_Wv']
-    return params['BW_factor']*(params['BW_deltaH'] - T*deltaS + P*params['BW_deltaV']) + R*T*(lnxdisord(n,Q) - lnxord(n,Q)) + (2.*Q - 1.)*W
+    return params['BW_deltaH'] - T*deltaS + P*params['BW_deltaV'] + R*T*(lnxdisord(n,Q) - lnxord(n,Q)) + (2.*Q - 1.)*W
 
 # see Holland and Powell, 1996
 def gibbs_disord_BW(P, T, params):
     n=params['BW_n']
     deltaS=entropydisorder(n)
     Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
-    print P, T, Q, '*** 0 =', equilibrium_Q(Q, deltaS, P, T, params)
+    #print P, T, Q, '*** 0 =', equilibrium_Q(Q, deltaS, P, T, params)
     W=params['BW_W'] + P*params['BW_Wv']
     ideal=(1.-Q)*(params['BW_deltaH'] - T*entropydisorder(n) + P*params['BW_deltaV'] + R*T*lnxdisord(n,Q)) + Q*(R*T*lnxord(n,Q))
     nonideal=(1.-Q)*Q*W
-    print '***', T, ideal/2, (1.-Q)*params['BW_deltaH']/2, (-1.+Q)*T*entropydisorder(n)/2, (1.-Q)*R*T*lnxdisord(n,Q)/2, Q*(R*T*lnxord(n,Q))/2, (1.-Q)*Q*W/2
     Edisord=ideal+nonideal
+    #print '***', T, Q, ideal*(params['BW_factor']), (1.-Q)*params['BW_deltaH']*(params['BW_factor']), (-1.+Q)*T*entropydisorder(n)*(params['BW_factor']), (1.-Q)*R*T*lnxdisord(n,Q)*(params['BW_factor']), Q*(R*T*lnxord(n,Q))*(params['BW_factor']), (1.-Q)*Q*W*(params['BW_factor']), Edisord*(params['BW_factor']), Edisord
     return Edisord
 
 # see Holland and Powell, 2011
diff --git a/burnman/testing.py b/burnman/testing.py
index 6326e64..3dba00d 100644
--- a/burnman/testing.py
+++ b/burnman/testing.py
@@ -127,7 +127,7 @@ print 'Root mean squared error:', rmserror, 'J/mol'
 print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
 print ''
 
-
+"""
 # A mineral with a phase transition described by Bragg-Williams theory: spinel
 print 'Spinel: a phase with order-disorder described by Bragg-Williams'
 sp=minerals.HP_2011.spinel()
@@ -183,7 +183,7 @@ for pi in range(len(Gsill)-1):
         Gburnman=sill.calcgibbs(P,T)
         GHP=Gsill[pi+1][ti+1]*1000. # convert to J/mol
         Gdiff=Gburnman-GHP
-        #print P, T, Gburnman, GHP, Gdiff
+        print P, T, Gburnman, GHP, Gdiff
         serror=serror + pow(Gdiff,2)
         if abs(Gdiff) > abs(maxerror):
             maxerror = Gdiff
@@ -194,7 +194,7 @@ rmserror = np.sqrt(serror/((len(Gsill)-1)*(len(Gsill[0])-1)))
 print 'Root mean squared error:', rmserror, 'J/mol'
 print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
 print ''
-
+"""
 
 # Another mineral with a phase transition described by Bragg-Williams theory: hercynite
 print 'Hercynite: a phase with order-disorder described by Bragg-Williams'
@@ -229,7 +229,7 @@ print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
 print ''
 
 
-
+"""
 # Another mineral with a phase transition described by Bragg-Williams theory: dolomite
 print 'Dolomite: a phase with order-disorder described by Bragg-Williams'
 dol=minerals.HP_2011.dolomite()
@@ -261,4 +261,3 @@ rmserror = np.sqrt(serror/((len(Gdol)-1)*(len(Gdol[0])-1)))
 print 'Root mean squared error:', rmserror, 'J/mol'
 print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
 print ''
-"""