[cig-commits] [commit] add_thermodynamic_potentials: Attempts to locate tc-burnman discrepancies (9f037a6)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:52:57 PST 2014
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 ''
-"""
More information about the CIG-COMMITS
mailing list