[cig-commits] [commit] inversion: Added Enstatite-Pyrolite Mixing model (e9242b5)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sun Dec 14 15:12:03 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : inversion
Link : https://github.com/geodynamics/burnman/compare/fcc314ffd83374fe1d6a4d79f8bd2f8be734dcd1...e9242b52d6d1ba1cf4d512e6e2c72ba57cc7f26d
>---------------------------------------------------------------
commit e9242b52d6d1ba1cf4d512e6e2c72ba57cc7f26d
Author: Cayman Unterborn <kmanunterborn at gmail.com>
Date: Sun Dec 14 15:10:36 2014 -0800
Added Enstatite-Pyrolite Mixing model
Seems to be missing a set_method() in Mg_fe_PT dependent...
>---------------------------------------------------------------
e9242b52d6d1ba1cf4d512e6e2c72ba57cc7f26d
.../{realdata_adiabat.py => realdata_enstpyro.py} | 64 +++++++++++-----------
1 file changed, 33 insertions(+), 31 deletions(-)
diff --git a/inversion/realdata_adiabat.py b/inversion/realdata_enstpyro.py
similarity index 79%
copy from inversion/realdata_adiabat.py
copy to inversion/realdata_enstpyro.py
index f75e157..e46af5a 100644
--- a/inversion/realdata_adiabat.py
+++ b/inversion/realdata_enstpyro.py
@@ -23,14 +23,31 @@ number_of_points = 5 #set on how many depth slices the computations should be do
# velocity constraints from seismology
+Kd_0 = .5
+
+weight_percents_enst = {'Mg':0.213, 'Fe': 0.0721, 'Si':0.242, 'Ca':0., 'Al':0.}
+weight_percents_pyro = {'Mg':0.228, 'Fe': 0.0626, 'Si':0.21, 'Ca':0., 'Al':0.}
+
+phase_fractions_pyro,relative_molar_percent_pyro = burnman.calculate_phase_percents(weight_percents_pyro)
+phase_fractions_enst,relative_molar_percent_enst = burnman.calculate_phase_percents(weight_percents_enst)
+
+iron_content_enst = lambda p,t: burnman.calculate_partition_coefficient(p,t,relative_molar_percent_enst, Kd_0)
+iron_content_pyro = lambda p,t: burnman.calculate_partition_coefficient(p,t,relative_molar_percent_pyro, Kd_0)
+
+enstatite = burnman.Composite ([phase_fractions_enst['pv'], phase_fractions_enst['fp']],
+ [minerals.SLB_2005.mg_fe_perovskite_pt_dependent(iron_content_enst,0),
+ minerals.SLB_2005.ferropericlase_pt_dependent(iron_content_enst,1)] )
+
+pyrolite = burnman.Composite([phase_fractions_pyro['pv'], phase_fractions_pyro['fp']],
+ [minerals.SLB_2005.mg_fe_perovskite_pt_dependent(iron_content_pyro,0),
+ minerals.SLB_2005.ferropericlase_pt_dependent(iron_content_pyro,1)])
+
seismic_model = burnman.seismic.PREM() # pick from .prem() .slow() .fast() (see code/seismic.py)
depths = np.linspace(1000e3,2500e3, number_of_points)
seis_p, seis_rho, seis_vp, seis_vs, seis_vphi = seismic_model.evaluate_all_at(depths)
+temperature = burnman.geotherm.brown_shankland(seis_p)
-seismic_model2 = burnman.seismic.PREM() # pick from .prem() .slow() .fast() (see code/seismic.py)
-depths = np.linspace(1000e3,2500e3, number_of_points)
-seis_p2, seis_rho2, seis_vp2, seis_vs2, seis_vphi2 = seismic_model2.evaluate_all_at(depths)
'''
print seis_vs-seis_vs2
@@ -74,18 +91,11 @@ print "preparations done"
# Priors on unknown parameters:
-fraction_pv = pymc.Uniform('fraction_pv', 0.0, 1.0)
-fe_pv = pymc.Uniform('fe_pv', 0.0, 0.5)
-fe_pc= pymc.Uniform('fe_pc', 0.0, 0.5)
+fraction_pyro = pymc.Uniform('fraction_pv', 0.0, 1.0)
T0 = pymc.Normal('T0',mu=1500.,tau=1./500.)
-def calc_all_velocities(fraction_pv, fe_pv, fe_pc,T0):
+def calc_all_velocities(fraction_pyro):
method = 'slb3' #slb3|slb2|mgd3|mgd2
- pv=minerals.other.mg_fe_perovskite(fe_pv)
- pc=minerals.other.ferropericlase(fe_pc)
- rock = burnman.Composite( [fraction_pv, 1.0-fraction_pv],[ pv,pc ] )
- temperature = burnman.geotherm.adiabatic(seis_p,T0, rock)
-
-
+ rock = burnman.Composite( [fraction_pyro, 1.0-fraction_pyro],[ pyrolite,enstatite ] )
mat_rho, mat_vp, mat_vs, mat_vphi, mat_K, mat_G = burnman.velocities_from_rock(rock,seis_p, temperature)
return mat_vp, mat_vs, mat_rho, mat_vphi, mat_K, mat_G
@@ -107,10 +117,10 @@ def nrmse(funca,funcb):
nrmse=rmse/(np.max(funcb)-np.min(funcb))
return nrmse
-def error(fraction_pv, fe_pv, fe_pc,T0):
+def error(fraction_pyro):
if True:
- if fraction_pv>0. and fraction_pv<1.0 and fe_pv>0. and fe_pv<1.0 and fe_pc>0. and fe_pc<1.0 and T0>300.:
- mat_vp, mat_vs, mat_rho, mat_vphi, mat_K, mat_G = calc_all_velocities(fraction_pv, fe_pv, fe_pc,T0)
+ if fraction_pyro>0. and fraction_pyro<1.0:
+ mat_vp, mat_vs, mat_rho, mat_vphi, mat_K, mat_G = calc_all_velocities(fraction_pyro)
misfit = nrmse(mat_rho,seis_rho)+nrmse(mat_K,seis_K)+nrmse(mat_G,seis_G)
return misfit
@@ -122,8 +132,8 @@ def error(fraction_pv, fe_pv, fe_pc,T0):
# rewrite as function of f_pv and fe_content
@pymc.deterministic
-def calc_misfit(fraction_pv=fraction_pv, fe_pv=fe_pv, fe_pc=fe_pc,T0=T0):
- return error(fraction_pv,fe_pv,fe_pc,T0)
+def calc_misfit(fraction_pyro=fraction_pyro):
+ return error(fraction_pyro)
#sig = 1e-2
@@ -132,7 +142,7 @@ def calc_misfit(fraction_pv=fraction_pv, fe_pv=fe_pv, fe_pc=fe_pc,T0=T0):
sig=1.e-2 # Some sorts of error
obs = pymc.Normal('d',mu=calc_misfit,tau=1.0/(sig*sig),value=0,observed=True,trace=True)
model = [fraction_pv, fe_pv, fe_pc, T0,obs]
-things = ['fraction_pv', 'fe_pv','fe_pc','T0']
+things = ['fraction_pyro']
whattodo = ""
@@ -206,7 +216,7 @@ if whattodo=="plot":
i=i+1
- plt.savefig("example_inv_big_pv.png")
+ plt.savefig("example_inv_big_pyro.png")
plt.show()
if whattodo=="test":
@@ -220,7 +230,7 @@ if whattodo=="test":
for t in things:
print t,"\t",db.trace(t).stats()['mean']
- print "#samples: %d" % db.trace('fraction_pv').stats()['n']
+ print "#samples: %d" % db.trace('fraction_pyro').stats()['n']
pymc.raftery_lewis(S, q=0.025, r=0.01)
@@ -232,16 +242,8 @@ if whattodo=="test":
pymc.Matplot.trace(db.trace('deviance',chain=None).gettrace(burn=1000,thin=t,chain=None),'deviance',rows=2,columns=4,num=1)
- pymc.Matplot.trace(db.trace('fraction_pv',chain=None).gettrace(thin=t,chain=None),'fraction_pv',rows=2,columns=4,num=2)
- pymc.Matplot.histogram(np.array(db.trace('fraction_pv',chain=None).gettrace(burn=b,chain=None)),'fraction_pv',rows=2,columns=4,num=6)
-
- pymc.Matplot.trace(db.trace('fe_pv',chain=None).gettrace(thin=t,chain=None),'fe_pv',rows=2,columns=4,num=3)
- pymc.Matplot.histogram(np.array(db.trace('fe_pv',chain=None).gettrace(burn=b,chain=None)),'fe_pv',rows=2,columns=4,num=7)
-
- pymc.Matplot.trace(db.trace('fe_pc',chain=None).gettrace(thin=t,chain=None),'fe_pc',rows=2,columns=4,num=4)
- pymc.Matplot.histogram(np.array(db.trace('fe_pc',chain=None).gettrace(burn=b,chain=None)),'fe_pc',rows=2,columns=4,num=8)
-
-
+ pymc.Matplot.trace(db.trace('fraction_pyro',chain=None).gettrace(thin=t,chain=None),'fraction_pyro',rows=2,columns=4,num=2)
+ pymc.Matplot.histogram(np.array(db.trace('fraction_pyro',chain=None).gettrace(burn=b,chain=None)),'fraction_pv',rows=2,columns=4,num=6)
plt.show()
More information about the CIG-COMMITS
mailing list