Skip to content

Commit

Permalink
added covariance matrices to HP_2011_ds62 and JH_2015
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed Oct 8, 2018
1 parent e12056e commit 759d69c
Show file tree
Hide file tree
Showing 3 changed files with 16,881 additions and 523 deletions.
185 changes: 105 additions & 80 deletions burnman/data/input_raw_endmember_datasets/HPdata_to_burnman.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,15 +44,9 @@ def __init__(self, name, atoms, formula, sites, comp, H, S, V, Cp, a, k, flag, o
self.od = od


def read_dataset(datafile):
f = open(datafile, 'r')
ds = []
for line in f:
ds.append(line.split())
return ds

ds = read_dataset('tc-ds62.txt')

# Read dataset
with open('tc-ds62.txt', 'r') as f:
ds = [line.split() for line in f]

def getmbr(ds, mbr):
mbrarray = []
Expand All @@ -74,74 +68,105 @@ def getmbr(ds, mbr):
ds[i * 4 + 4][1]), float(ds[i * 4 + 4][2]), map(float, ds[i * 4 + 5]), float(ds[i * 4 + 6][0]), list(map(float, ds[i * 4 + 6][1:4])), flag, list(map(float, ds[i * 4 + 6][5:])))
return endmember

print '# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences'
print '# Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU GPL v2 or later.'
print ''
print ''
print '"""'
print 'HP_2011 (ds-62)'
print 'Minerals from Holland and Powell 2011 and references therein'
print 'Update to dataset version 6.2'
print 'The values in this document are all in S.I. units,'
print 'unlike those in the original tc-ds62.txt'
print 'File autogenerated using HPdata_to_burnman.py'
print '"""'
print ''
print 'from ..mineral import Mineral'
print 'from ..solidsolution import SolidSolution'
print 'from ..solutionmodel import *'
print 'from ..processchemistry import dictionarize_formula, formula_mass'
print ''

solutionfile = 'tc-ds62_solutions.txt'
with open(solutionfile, 'r') as fin:
print fin.read()
fin.close()

print '"""'
print 'ENDMEMBERS'
print '"""'
print ''
formula = '0'
for i in range(int(ds[0][0])):
mbr = ds[i * 4 + 3][0]
M = getmbr(ds, mbr)
if mbr == 'and': # change silly abbreviation
mbr = 'andalusite'
if M.flag != -1 and M.flag != -2 and M.k[0] > 0:
print 'class', mbr, '(Mineral):'
print ' def __init__(self):'
print ''.join([' formula=\'', M.formula, '\''])
print ' formula = dictionarize_formula(formula)'
print ' self.params = {'
print ''.join([' \'name\': \'', M.name, '\','])
print ' \'formula\': formula,'
print ' \'equation_of_state\': \'hp_tmt\','
print ' \'H_0\':', M.H * 1e3, ','
print ' \'S_0\':', M.S * 1e3, ','
print ' \'V_0\':', M.V * 1e-5, ','
print ' \'Cp\':', [round(M.Cp[0] * 1e3, 10), round(M.Cp[1] * 1e3, 10), round(M.Cp[2] * 1e3, 10), round(M.Cp[3] * 1e3, 10)], ','
print ' \'a_0\':', M.a, ','
print ' \'K_0\':', M.k[0] * 1e8, ','
print ' \'Kprime_0\':', M.k[1], ','
print ' \'Kdprime_0\':', M.k[2] * 1e-8, ','
print ' \'n\': sum(formula.values()),'
print ' \'molar_mass\': formula_mass(formula)}'
if M.flag != 0:
print ' self.property_modifiers = [['
if M.flag == 1:
print ' \'landau_hp\', {\'P_0\':', 1.e5, ','
print ' \'T_0\':', 298.15, ','
print ' \'Tc_0\':', M.od[0], ','
print ' \'S_D\':', M.od[1] * 1e3, ','
print ' \'V_D\':', M.od[2] * 1e-5, '}]]'
if M.flag == 2:
print ' \'bragg_williams\', {\'deltaH\':', M.od[0] * 1e3, ','
print ' \'deltaV\':', M.od[1] * 1e-5, ','
print ' \'Wh\':', M.od[2] * 1e3, ','
print ' \'Wv\':', M.od[3] * 1e-5, ','
print ' \'n\':', M.od[4], ','
print ' \'factor\':', M.od[5], '}]]'
print ' Mineral.__init__(self)'
print ''
print ''
with open('HP_2011_ds62.py', 'wb') as outfile:
outfile.write('# This file is part of BurnMan - a thermoelastic and '
'thermodynamic toolkit for the Earth and Planetary Sciences\n'
'# Copyright (C) 2012 - 2018 by the BurnMan team, '
'released under the GNU \n# GPL v2 or later.\n\n\n'
'"""\n'
'HP_2011 (ds-62)\n'
'Endmember minerals from Holland and Powell 2011 and references therein\n'
'Update to dataset version 6.2\n'
'The values in this document are all in S.I. units,\n'
'unlike those in the original tc-ds62.txt\n'
'File autogenerated using HPdata_to_burnman.py\n'
'"""\n\n'
'from ..mineral import Mineral\n'
'from ..processchemistry import dictionarize_formula, formula_mass\n\n')

outfile.write('"""\n'
'ENDMEMBERS\n'
'"""\n\n')

formula = '0'
for i in range(int(ds[0][0])):
mbr = ds[i * 4 + 3][0]
M = getmbr(ds, mbr)
if mbr == 'and': # change silly abbreviation
mbr = 'andalusite'
if M.flag != -1 and M.flag != -2 and M.k[0] > 0:
outfile.write('class {0} (Mineral):\n'.format(mbr)+
' def __init__(self):\n'
' formula = \'{0}\''.format(M.formula)+'\n'
' formula = dictionarize_formula(formula)\n'
' self.params = {{\n'
' \'name\': \'{0}\','.format(M.name)+'\n'
' \'formula\': formula,\n'
' \'equation_of_state\': \'hp_tmt\',\n'
' \'H_0\': {0},'.format(M.H * 1e3)+'\n'
' \'S_0\': {0},'.format(M.S * 1e3)+'\n'
' \'V_0\': {0},'.format(M.V * 1e-5)+'\n'
' \'Cp\': {0},'.format([round(M.Cp[0] * 1e3, 10), round(M.Cp[1] * 1e3, 10), round(M.Cp[2] * 1e3, 10), round(M.Cp[3] * 1e3, 10)])+'\n'
' \'a_0\': {0},'.format(M.a)+'\n'
' \'K_0\': {0},'.format(M.k[0] * 1e8)+'\n'
' \'Kprime_0\': {0},'.format(M.k[1])+'\n'
' \'Kdprime_0\': {0},'.format(M.k[2] * 1e-8)+'\n'
' \'n\': sum(formula.values()),\n'
' \'molar_mass\': formula_mass(formula)}\n')
if M.flag != 0:
outfile.write(' self.property_modifiers = [[\n')
if M.flag == 1:
outfile.write(' \'landau_hp\', {{\'P_0\': 100000.0,\n'
' \'T_0\': 298.15,\n'
' \'Tc_0\': {0},'.format(M.od[0])+'\n'
' \'S_D\': {0},'.format(M.od[1] * 1e3)+'\n'
' \'V_D\': {0}}}]]'.format(M.od[2] * 1e-5)+'\n\n')
if M.flag == 2:
outfile.write(' \'bragg_williams\', {{\'deltaH\': {0},'.format(M.od[0] * 1e3)+'\n'
' \'deltaV\': {0},'.format(M.od[1] * 1e-5)+'\n'
' \'Wh\': {0},'.format(M.od[2] * 1e3)+'\n'
' \'Wv\': {0},'.format(M.od[3] * 1e-5)+'\n'
' \'n\': {0},'.format(M.od[4])+'\n'
' \'factor\': {0}}}]]'.format(M.od[5])+'\n')
outfile.write(' Mineral.__init__(self)\n\n')


# Process uncertainties
import numpy as np
n_mbrs = int(ds[0][0])

names = []
for i in range(n_mbrs):
names.append(ds[i*4+3][0])

cov = []
for i in range(n_mbrs*4+4, len(ds)-2):
cov.extend(map(float, ds[i]))

i_utr = np.triu_indices(n_mbrs)
i_ltr = np.tril_indices(n_mbrs)
M = np.zeros((n_mbrs, n_mbrs))

M[i_utr] = cov[1:]
M[i_ltr] = M.T[i_ltr]

M = M*1.e6 # (kJ/mol)^2 -> (J/mol)^2


d = {'endmember_names':names,
'covariance_matrix':M}

np.set_printoptions(threshold='nan')

import pprint
pp = pprint.PrettyPrinter(indent=0, width=160, depth=3, stream=outfile)

outfile.write('\n# Variance-covariance matrix\n'
'# cov is a dictionary containing:\n'
'# - names: a list of endmember names\n'
'# - covariance_matrix: a 2D variance-covariance array for the endmember enthalpies of formation\n\n'
'from numpy import array\n'
'cov = ')
pp.pprint(d)

outfile.write('\n')

0 comments on commit 759d69c

Please sign in to comment.