ENKI

Source code for equilibrate

'''
The equilibrate module implements a Python interface to the Equilibrate and 
EquilState classes. It also implements a Python class called MELTSmodel that 
wraps the objective-C classes: EquilibrateUsingMELTSv102, 
EquilibrateUsingMELTSv110, EquilibrateUsingMELTSv120, 
EquilibrateUsingpMELTSv561, EquilibrateUsingMELTSwithDEW, 
and EquilibrateUsingStixrude.

The Equilibrate class provides methods to calculate an equilibrium phase 
assemblage given a list of Phase class instances.  

You can calculate equilibrium in closed thermodynamic systems under a 
variety of constraints:

- Temperature and pressure (Gibbs free energy minimization)
 
- Entropy and pressure (enthalpy minimization)

- Temperature and volume (Helmholtz free energy minimization)

- Entropy and volume (internal energy minimization)

You can also calculate equilibrium under constraints of fixed temperature and 
pressure in open thermodynamic systems by specifying one or more fixed elemental 
chemical potentials.

For details of the underlying theory and algorithms implemented in the Equilibrate 
class, see the notebooks in the PublicNotebooks/Equilibrate folder on 
the `ENKI server <https://enki.ofm-research.org/hub/login>`_.
'''

import numpy as np
import scipy as sp

import xml.etree.ElementTree as ET
import locale
import sys

from thermoengine import core
from thermoengine import chem

# Objective-C imports
import ctypes
from ctypes import cdll
from ctypes import util
from collections import OrderedDict
from rubicon.objc import ObjCClass, ObjCProtocol, NSObject, objc_method
cdll.LoadLibrary(util.find_library('phaseobjc'))

# Excel imports
from openpyxl import Workbook
from openpyxl.utils import get_column_letter

__all__ = ['Equilibrate', 'EquilState', 'MELTSmodel']

# array length for elements
NE = 106

[docs]class EquilState: """Class for holding the phase state for an equilibrium model determined by Equilibrate Parameters ---------- elements_l : [] A list of strings identifying the chemical elements present in the system phases_l : [] A list of class instances for phases that are in the system. These instances must derive from the *PurePhase* or *SolutionPhase* classes, which are defined in the *phases* module. Attributes ---------- c_matrix c_omni c_omni_qr c_qr_decomp element_l phase_d pressure temperature """ def __init__(self, elements_l, phases_l): elements_present = [False for i in range (0,NE)] for elm in elements_l: index = core.chem.PERIODIC_ORDER.tolist().index(elm) elements_present[index] = True self._element_l = [i for i,e in enumerate(elements_present) if e] self._phase_d = OrderedDict() for phase in phases_l: entry = dict() entry['obj'] = phase nc = int(phase.props['endmember_num']) \ if 'endmember_num' in phase.props else 1 entry['name'] = phase.props['phase_name'] entry['abbrev'] = phase.props['abbrev'] entry['nc'] = nc entry['end_name'] = phase.props['endmember_name'] entry['end_formula'] = phase.props['formula'] entry['end_molwts'] = phase.props['molwt'] entry['moles'] = np.zeros(nc) entry['mole_frac'] = np.zeros(nc) entry['grams'] = 0.0 entry['affinity'] = 0.0 # create a list of length nc where each entry is a numpy array # that contains stoichiometric coefficients converting each # endmember to system elements (in the order _element_l) conv_to_elm = [] conv_to_elm_sum = np.zeros(len(self._element_l)) conv_to_oxd = [] for i in range(0,nc): mol = np.zeros(nc) mol[i] = 1.0 if nc == 1: mol_elm = phase.props['element_comp'][0] else: mol_elm = entry['obj'].covert_endmember_comp( mol,output='moles_elements') conv_entry = np.zeros(len(self._element_l)) for i in range(0,NE): if mol_elm[i] != 0.0: conv_entry[self._element_l.index(i)] = mol_elm[i] conv_to_elm_sum = np.add(conv_to_elm_sum, conv_entry) conv_to_elm.append(conv_entry) conv_to_oxd.append(core.chem.calc_mol_oxide_comp(mol_elm)) entry['conv_to_elm'] = conv_to_elm entry['conv_to_oxd'] = conv_to_oxd entry['omnicomponent'] = True if np.count_nonzero( conv_to_elm_sum) == len(conv_to_elm_sum) else False name = phase.props['phase_name'] self._phase_d[name] = entry first = True for entry in self._phase_d.values(): nc = entry['nc'] for row in entry['conv_to_elm']: if first: c_matrix = np.reshape(row,(1,row.size)) first = False else: c_matrix = np.vstack((c_matrix,np.reshape( row,(1,row.size)))) self._c_matrix = c_matrix self._c_qr_decomp = np.linalg.qr(self.c_matrix) self._c_omni = None self._c_omni_qr = None self._temperature = 1000.0 self._pressure = 1000.0 ############################ # Properties of the system # ############################ @property def phase_d(self): """ A dictionary of dictionaries that holds properties of system phases Returns ------- Dictionary of phases in the system (dict) """ return self._phase_d @property def element_l(self): """ List of atomic numbers of elements in the system Returns ------- A Python list """ return self._element_l @property def c_matrix(self): """ Numpy matrix that maps elements to moles of endmembers Moles of endmembers are indexed by row, and moles of elements are indexed by column. Returns ------- Numpy 2-D array """ return self._c_matrix @property def c_qr_decomp(self): """ Tuple of Numpy matrices of the Q,R decomposition of the c_matrix Returns ------- tuple """ return self._c_qr_decomp @property def c_omni(self): """ Numpy conversion matrix for the omnicomponent phase (if one exists) Moles of endmembers are indexed by row, and moles of elements are indexed by column. Returns ------- np 2-D array """ return self._c_omni @property def c_omni_qr(self): """ Tuple of Numpy matrices of the Q,R decomposition of c_omni Returns ------- tuple """ return self._c_omni_qr @property def temperature(self): """ Temperature of the assemblage, in Kelvins Returns ------- float """ return self._temperature @temperature.setter def temperature(self, value): assert isinstance(value, (int, float)) assert value > 0 self._temperature = value @property def pressure(self): """ Pressure of the assemblage, in bars Returns ------- float """ return self._pressure @pressure.setter def pressure(self, value): assert isinstance(value, (int, float)) assert value > 0 self._pressure = value ################################ # Methods for phases in system # ################################
[docs] def omni_phase(self): """ Returns first omnicomponent phase in the phase dictionary """ for entry in self.phase_d.values(): if entry['omnicomponent']: return entry['name'] return None
[docs] def c_matrix_for_phase(self, phase_name): """ Returns a slice of c_matrix relevant to the specified phase_name Parameters ---------- phase_name : str Name of a system phase """ if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return None row_min = 0 for key,entry in self.phase_d.items(): if key == phase_name: row_max = row_min + entry['nc'] break else: row_min += entry['nc'] return self.c_matrix[row_min:row_max,:]
[docs] def set_phase_comp(self, phase_name, cmp, input_as_elements=False): """ Sets the endmember moles and total moles a phase in the system Parameters ---------- phase_name : str Name of a system phase cmp : ndarray 1-D Numpy array with compositional data input_as_elements : bool, def False If True, convert input array from moles of elements to moles of endmember components. Returns ------- valid : bool True if input composition is valid for phase """ if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return False assert(type(cmp) is np.ndarray, 'bulk_comp must be set to an numpy array') if not input_as_elements: assert(cmp.size == self.phase_d[phase_name]['nc'], 'cmp array must have length ' + str(self.phase_d[phase_name]['nc'])) else: assert(cmp.size == len(self.element_l), 'cmp array must have length ' + str(len(self.element_l))) cmp = np.matmul( np.linalg.inv(self.c_matrix_for_phase(phase_name).T), cmp) self.phase_d[phase_name]['moles'] = cmp return self.phase_d[phase_name]['obj'].test_endmember_comp(cmp)
[docs] def tot_moles_phase(self, phase_name): """ Returns total moles of phase Parameters ---------- phase_name : str Name of a system phase """ if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return 0.0 return np.sum(self.phase_d[phase_name]['moles'])
[docs] def tot_grams_phase(self, phase_name): """ Returns total grams of phase Parameters ---------- phase_name : str Name of a system phase """ if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return 0.0 grams = np.matmul(self.phase_d[phase_name]['moles'], self.phase_d[phase_name]['end_molwts']) return np.sum(grams)
[docs] def moles_elements(self, phase_name): """ Returns array of moles of elements in phase Parameters ---------- phase_name : str Name of a system phase Returns ------- result : np_array Numpy array of mole numbers of the reduced set of elements in the phase """ if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return None elm = np.zeros(len(self.element_l)) for moles,coeff in zip(self.phase_d[phase_name]['moles'], self.phase_d[phase_name]['conv_to_elm']): if moles != 0: elm = np.add(elm,moles*coeff) return elm
[docs] def tot_moles_elements(self): """ Returns array of moles of elements in system Returns ------- result : np_array Numpy array of mole numbers of the reduced set of elements in the system """ elm = np.zeros(len(self.element_l)) for phase in self.phase_d.keys(): if self.tot_moles_phase(phase) > 0: elm += self.moles_elements(phase) return elm
[docs] def oxide_comp(self, phase_name, output='wt%', prune_list=True): """ Returns a dictionary of concentrations of oxides in phase Parameters ---------- phase_name : str Name of a system phase output : str, default "wt%" Units of output: 'wt%' (default), 'moles', 'grams' prune_list : bool, default True Remove zeroed valued entries Returns ------- result : collections.OrderedDict Dictionary of oxide concentration values, keyed by oxide """ if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return None oxd = np.zeros(core.chem.oxide_props['oxide_num']) for moles,coeff in zip(self.phase_d[phase_name]['moles'], self.phase_d[phase_name]['conv_to_oxd']): if moles != 0: oxd = np.add(oxd,moles*coeff) result = OrderedDict() for key,value in zip(core.chem.oxide_props['oxides'],oxd): result[key] = value if output == 'wt%' or output == 'grams': mod_result = OrderedDict() sum_oxd = 0.0 for index,(key,value) in enumerate(result.items()): mod_result[key] = value*chem.oxide_props['molwt'][index] sum_oxd += mod_result[key] result = mod_result if output == 'wt%' and sum_oxd > 0: mod_result = OrderedDict() for (key,value) in result.items(): mod_result[key] = 100.0*value/sum_oxd result = mod_result if prune_list: mod_result = OrderedDict() for (key,value) in result.items(): if value != 0: mod_result[key] = value result = mod_result return result
[docs] def properties(self, phase_name=None, props=None, units=False): """ Returns properties of phases in the system Parameters ---------- phase_name : str, default None Name of phase in system. If None, returns a dictionary of system phases, with names as keys and 'stable' or 'unstable' as values. props : str, default None Name of property to be retrieved. If None, a list of valid properties is returned. units : bool, default False Property units are provided as the second entry of a returned tuple. """ prop_d = OrderedDict([ ('Mass','g'), ('GibbsFreeEnergy','J'), ('Enthalpy','J'), ('Entropy','J/K'), ('HeatCapacity','J/K'), ('DcpDt','J/K^2'), ('Volume','J/bar'), ('DvDt','J/bar-K'), ('DvDp','J/bar^2'), ('D2vDt2','J/bar-K^2'), ('D2vDtDp','J/bar^2-K'), ('D2vDp2','J/bar^3'), ('Density','g/cm^3'), ('Alpha','1/K'), ('Beta','1/bar'), ('K','GPa'), ("K'",'none'), ('Gamma','none') ]) if phase_name is None: phase_d = {} for key in self.phase_d.keys(): phase_d[key] = 'stable' if self.tot_moles_phase( key) > 0 else 'unstable' phase_d['System'] = 'stable' return phase_d if phase_name not in self.phase_d and phase_name != 'System': print('Error: ' + phase_name + ' is not in phase dictionary.') return None if props is None: if units: return set(zip(list(prop_d.keys()), list(prop_d.values()))) else: return list(prop_d.keys()) if props not in prop_d: print('Error: ' + props + ' is not a valid property.') return None t = self.temperature p = self.pressure if phase_name != 'System': entry = self.phase_d[phase_name] moles = entry['moles'] if np.sum(moles) == 0: result = 0.0 if units: result = (0.0, None) return result nc = entry['nc'] sys = False else: sys = True result = None if units: result = (None, None) if props == 'Mass': if sys: result = 0.0 for phs in self.phase_d.keys(): grams = np.matmul(self.phase_d[phs]['moles'], self.phase_d[phs]['end_molwts']) result += np.sum(grams) else: grams = np.matmul(self.phase_d[phase_name]['moles'], self.phase_d[phase_name]['end_molwts']) result = np.sum(grams) if units: result = (result, prop_d['Mass']) elif props == 'GibbsFreeEnergy': if sys: result = self.G(t, p) else: if nc == 1: result = moles[0]*entry['obj'].gibbs_energy(t,p) else: result = entry['obj'].gibbs_energy(t,p,mol=moles) if units: result = (result, prop_d['GibbsFreeEnergy']) elif props == 'Enthalpy': if sys: result = self.G(t, p) - t*self.dGdT(t, p) else: if nc == 1: result = moles[0]*(entry['obj'].gibbs_energy(t,p) +t*entry['obj'].entropy(t,p)) else: result = (entry['obj'].gibbs_energy(t,p,mol=moles) +t*entry['obj'].entropy(t,p,mol=moles)) if units: result = (result, prop_d['Enthalpy']) elif props == 'Entropy': if sys: result = -self.dGdT(t, p) else: if nc == 1: result = moles[0]*entry['obj'].entropy(t,p) else: result = entry['obj'].entropy(t,p,mol=moles) if units: result = (result, prop_d['Entropy']) elif props == 'HeatCapacity': if sys: result = -t*self.d2GdT2(t, p) else: if nc == 1: result = moles[0]*entry['obj'].heat_capacity(t,p) else: result = entry['obj'].heat_capacity(t,p,mol=moles) if units: result = (result, prop_d['HeatCapacity']) elif props == 'DcpDt': if sys: result = -t*self.d3GdT3(t, p) - self.d2GdT2(t, p) else: if nc == 1: result = moles[0]*entry['obj'].heat_capacity( t,p,deriv={'dT':1}) else: result = entry['obj'].heat_capacity( t,p,mol=moles,deriv={'dT':1}) if units: result = (result, prop_d['DcpDt']) elif props == 'Volume': if sys: result = self.dGdP(t, p) else: if nc == 1: result = moles[0]*entry['obj'].volume(t,p) else: result = entry['obj'].volume(t,p,mol=moles) if units: result = (result, prop_d['Volume']) elif props == 'DvDt': if sys: result = self.d2GdTdP(t, p) else: if nc == 1: result = moles[0]*entry['obj'].volume(t,p,deriv={'dT':1}) else: result = entry['obj'].volume(t,p,mol=moles,deriv={'dT':1}) if units: result = (result, prop_d['DvDt']) elif props == 'DvDp': if sys: result = self.d2GdP2(t, p) else: if nc == 1: result = moles[0]*entry['obj'].volume(t,p,deriv={'dP':1}) else: result = entry['obj'].volume(t,p,mol=moles,deriv={'dP':1}) if units: result = (result, prop_d['DvDp']) elif props == 'D2vDt2': if sys: result = self.d3GdT2dP(t, p) else: if nc == 1: result = moles[0]*entry['obj'].volume(t,p,deriv={'dT':2}) else: result = entry['obj'].volume(t,p,mol=moles,deriv={'dT':2}) if units: result = (result, prop_d['D2vDt2']) elif props == 'D2vDtDp': if sys: result = self.d3GdTdP2(t, p) else: if nc == 1: result = moles[0]*entry['obj'].volume( t,p,deriv={'dT':1,'dP':1}) else: result = entry['obj'].volume(t,p,mol=moles, deriv={'dT':1,'dP':1}) if units: result = (result, prop_d['D2vDtDp']) elif props == 'D2vDp2': if sys: result = self.d3GdP3(t, p) else: if nc == 1: result = moles[0]*entry['obj'].volume(t,p,deriv={'dP':2}) else: result = entry['obj'].volume(t,p,mol=moles,deriv={'dP':2}) if units: result = (result, prop_d['D2vDp2']) elif props == 'Density': v = self.properties(phase_name=phase_name, props='Volume', units=False) m = self.properties(phase_name=phase_name, props='Mass', units=False) result = m/v/10.0 if units: result = (result, prop_d['Density']) elif props == 'Alpha': v = self.properties(phase_name=phase_name, props='Volume', units=False) dvdt = self.properties(phase_name=phase_name, props='DvDt', units=False) result = 10.0*dvdt/v if units: result = (result, prop_d['Alpha']) elif props == 'Beta': v = self.properties(phase_name=phase_name, props='Volume', units=False) dvdp = self.properties(phase_name=phase_name, props='DvDp', units=False) result = -10.0*dvdp/v if units: result = (result, prop_d['Beta']) elif props == 'K': beta = self.properties(phase_name=phase_name, props='Beta', units=False) result = 1.0/beta/10000.0 if units: result = (result, prop_d['K']) elif props == "K'": v = self.properties(phase_name=phase_name, props='Volume', units=False) dvdp = self.properties(phase_name=phase_name, props='DvDp', units=False) d2vdp2 = self.properties(phase_name=phase_name, props='D2vDp2', units=False) result = v*d2vdp2/dvdp/dvdp - 1 if units: result = (result, prop_d["K'"]) elif props == 'Gamma': v = self.properties(phase_name=phase_name, props='Volume', units=False) dvdp = self.properties(phase_name=phase_name, props='DvDp', units=False) dvdt = self.properties(phase_name=phase_name, props='DvDt', units=False) cp = self.properties(phase_name=phase_name, props='HeatCapacity', units=False) result = (-v*dvdt/dvdp)/(cp+t*dvdt*dvdt/dvdp) if units: result = (result, prop_d['Gamma']) return result
[docs] def compositions(self, phase_name=None, ctype='components', units='moles'): """ Returns compositions of phases in the system Parameters ---------- phase_name : str, default None Name of phase in system ctype : str, default 'components' Compositional descriptor. Permitted: 'components', 'oxides', 'elements' units : str, default 'moles' Units of composition. Permitted: 'moles', 'mole_frac', 'wt%' Returns ------- result : numpy array One dimensional vector Notes ----- If units is set to 'oxides' or 'elements', results are returned in an array of standard length and order. """ if phase_name is None: print ('Error: Method requires a phase_name.') return None if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return None if ctype not in ['components', 'oxides', 'elements']: print ('Error: ' + ctype + ' is not permitted. See doc string.') return None if units not in ['moles', 'mole_frac', 'wt%']: print ('Error: ' + units + ' is not permitted. See doc string.') return None entry = self.phase_d[phase_name] if ctype == 'components': if units == 'moles' or units == 'mole_frac': result = entry['moles'] if units == 'mole_frac': result = result/np.sum(result) elif units == 'wt%': result = entry['moles']*entry['end_molwts'] result = 100.0*result/np.sum(result) elif ctype == 'oxides': if units == 'moles' or units == 'mole_frac': result = self.oxide_comp(phase_name, output='moles', prune_list=False) result = np.array(list(result.values())) if units == 'mole_frac': result = result/np.sum(result) elif units == 'wt%': result = self.oxide_comp(phase_name, output='wt%', prune_list=False) result = np.array(list(result.values())) elif ctype == 'elements': moles = entry['moles'] if entry['nc'] == 1: result = moles[0]*entry['obj'].props['element_comp'][0] else: result = np.zeros(NE) for i in range(0,entry['nc']): mol = np.zeros(entry['nc']) mol[i] = 1.0 mol_elm = entry['obj'].covert_endmember_comp( mol,output='moles_elements') for j in range(1,NE): result[j] += moles[i]*mol_elm[j] if units == 'mole_frac': result = result/np.sum(result) elif units == 'wt%': for i in range(1,NE): result[i] *=chem.PERIODIC_WEIGHTS[i] result = 100.0*result/np.sum(result) return result
[docs] def affinities(self, phase_name=None): """ Returns affinity of phases in the system Parameters ---------- phase_name : str, default None Name of phase in system. If None, returns a dictionary of system phases, with names as keys and 'stable' or 'unstable' as values. """ if phase_name is None: print ('Error: Please specify a phase name.') return None if phase_name not in self.phase_d: print('Error: ' + phase_name + ' is not in phase dictionary.') return None tot_mol = self.tot_moles_phase(phase_name) if tot_mol > 0: return 0.0 else: return self.phase_d[phase_name]['affinity']
[docs] def print_state(self, level='summary', wt_as_oxides=True): """ Prints results about the system state Parameters ---------- level : str Level of detail to be printed wt_as_oxides : bool, default True Print wt% values on an oxide basis; otherwise print wt% of endmember components. """ print (' ') print ('T = {0:10.2f} °C, P = {1:10.1f} MPa'.format( self.temperature-273.15, self.pressure/10.0)) for key,entry in self.phase_d.items(): tot_mol = self.tot_moles_phase(key) if tot_mol > 0: tot_grm = self.tot_grams_phase(key) print ('{0:<15s}'.format(entry['name']), end=' ') print ('moles: {0:10.6f}'.format(tot_mol), end=' ') print ('grams: {0:7.3f}'.format(tot_grm)) if wt_as_oxides: oxd = self.oxide_comp(key) oxd_s = list(oxd.keys()) oxd_v = list(oxd.values()) n_oxd = len(oxd) nc = entry['nc'] if nc > 1: for i in range(0,nc): gram = entry['moles'][i]*entry['end_molwts'][i] print ('{0:>15s}'.format( entry['end_name'][i]), end=' ') print ('form: {0:<10s}'.format( entry['end_formula'][i]), end=' ') print ('X: {0:7.4f}'.format( entry['moles'][i]/tot_mol), end=' ') if wt_as_oxides: if i < n_oxd: print ('wt% {0:>7s} {1:7.2f}'.format( oxd_s[i], oxd_v[i])) else: print ('') else: print ('wt% {0:7.2f}'.format(100.0*gram/tot_grm)) if wt_as_oxides and nc < n_oxd: for i in range(nc,n_oxd): print ('{0:49s} wt% {1:>7s} {2:7.2f}'.format( ' ', oxd_s[i], oxd_v[i])) else: print ('{0:<15s}'.format(entry['name']), end=' ') print ('affn: {0:10.2f}'.format(entry['affinity'])) nc = entry['nc'] if nc > 1: for i in range(0,nc): print ('{0:>15s}'.format(entry['end_name'][i]), end=' ') print ('form: {0:<10s}'.format( entry['end_formula'][i]), end=' ') print ('X: {0:7.4f}'.format(entry['mole_frac'][i]))
################################################### # System potentials and compositional derivatives # ###################################################
[docs] def moles_v(self, reaction_v): """ Computes of a reaction product moles (scalar) = (mole vector of all endmembers of all active phases in the system) x (input mole vector of reaction coefficients) Parameters ---------- reaction_v : ndarray Array of reaction coefficients Returns ------- moles : float Moles of reaction product """ moles = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: for x in entry['moles']: moles.append(x) moles = np.array(moles) moles = np.matmul(reaction_v, moles) return moles
[docs] def G(self, t=1000.0, p=1000.0): """ Returns Gibbs free energy of the system; a scalar """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy(t,p) else: result += entry['obj'].gibbs_energy(t,p,mol=moles) return result
[docs] def dGdT(self, t=1000.0, p=1000.0): """ Returns the temperature derivative of the Gibbs free energy of the system; a scalar, and the negative of the system entropy """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += -moles[0]*entry['obj'].entropy(t,p) else: result += -entry['obj'].entropy(t,p,mol=moles) return result
[docs] def dGdP(self, t=1000.0, p=1000.0): """ Returns the pressure derivative of the Gibbs free energy of the system; a scalar, and the system volume """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy( t,p,deriv={'dP':1}) else: result += entry['obj'].gibbs_energy( t,p,mol=moles,deriv={'dP':1}) return result
[docs] def dGdn(self, t=1000.0, p=1000.0, element_basis=True, full_output=False, use_omni_phase=False): """ Returns the first molar derivative of G of the system; a 1-D numpy array Parameters ---------- element_basis : bool, def True If True, returns a projected array of element chemical potentials full_output : bool, def False If True, returns a tuple with full output from numpy lstlq method use_omni_phase : bool, def False If True, returns a projected array of element chemical potentials using the omnicomponent phase as a basis. The system must have an omnicomponent phase, and element_basis must be set to True for this keyword to have effect. """ if use_omni_phase: omni_phase_name = self.omni_phase() if omni_phase_name is None: assert(True, 'ERROR in dGdn method: use omni_phase keyword specified' + ' when the system has no omnicomponent phase') element_basis = True if use_omni_phase and element_basis: # computes chemical potentials for all component endmembers entry = self.phase_d[omni_phase_name] result = [] nc = entry['nc'] if nc == 1: row = np.array([entry['obj'].gibbs_energy(t, p)]) else: moles = entry['moles'] row = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dmol':1})[0] result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) else: # computes chemical potentials for all component endmembers in a # system phase if the total moles of that phase is > 0 result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: nc = entry['nc'] if nc == 1: row = np.array([entry['obj'].gibbs_energy(t, p)]) else: moles = entry['moles'] row = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dmol':1})[0] result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) if not element_basis: return result else: if use_omni_phase: if self.c_omni is None: self._c_omni = self.c_matrix_for_phase(omni_phase_name) self._c_omni_qr = np.linalg.qr(self.c_omni) Q, R = self.c_omni_qr else: first = True for key,entry in self._phase_d.items(): if self.tot_moles_phase(key) > 0.0: nc = entry['nc'] for row in entry['conv_to_elm']: if first: c_matrix = np.reshape(row,(1,row.size)) first = False else: c_matrix = np.vstack((c_matrix,np.reshape( row,(1,row.size)))) Q, R = np.linalg.qr(c_matrix) Qb = np.matmul(Q.T, result) result, residuals, rank, S = np.linalg.lstsq(R, Qb, rcond=None) return result if not full_output else (result, residuals, rank, S)
[docs] def d2GdT2(self, t=1000.0, p=1000.0): """ Returns the second temperature derivative of the Gibbs free energy of the system; a scalar, and the negative of the temperature times the system heat capacity """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += -moles[0]*entry['obj'].heat_capacity(t,p)/t else: result += -entry['obj'].heat_capacity(t,p,mol=moles)/t return result
[docs] def d2GdTdP(self, t=1000.0, p=1000.0): """ Returns the cross temperature-pressure derivative of the Gibbs free energy of the system; a scalar, and the volume times the system's coefficient of isothermal expansion """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy( t,p,deriv={'dT':1, 'dP':1}) else: result += entry['obj'].gibbs_energy( t,p,mol=moles,deriv={'dT':1, 'dP':1}) return result
[docs] def d2GdP2(self, t=1000.0, p=1000.0): """ Returns the second pressure derivative of the Gibbs free energy of the system; a scalar, and the negative of the volume times the system's coefficient of isothermal compressibility """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy( t,p,deriv={'dP':2}) else: result += entry['obj'].gibbs_energy( t,p,mol=moles,deriv={'dP':2}) return result
[docs] def d2GdTdn(self, t=1000.0, p=1000.0): """ Returns the second derivative of G of the system with respect to temperature and mole numbers; a 1-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: nc = entry['nc'] if nc == 1: row = np.array([-entry['obj'].entropy(t, p)]) else: moles = entry['moles'] row = -entry['obj'].entropy(t, p, mol=moles, deriv={'dmol':1})[0] result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) return result
[docs] def d2GdPdn(self, t=1000.0, p=1000.0): """ Returns the second derivative of G of the system with respect to pressure and mole numbers; a 1-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: nc = entry['nc'] if nc == 1: row = np.array([entry['obj'].gibbs_energy( t, p, deriv={'dP':1})]) else: moles = entry['moles'] row = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dP':1, 'dmol':1})[0] result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) return result
[docs] def d2Gdn2(self, t=1000.0, p=1000.0, element_basis=True): """ Returns the second molar derivative of G of system; a 2-D numpy array Parameters ---------- element_basis : bool, def True If True, returns a projected Hessian on an element basis """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: if entry['nc'] == 1: mat = np.array([[0.0]]) else: moles = entry['moles'] mat = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dmol':2})[0] result.append(mat) result = sp.linalg.block_diag(*result) if not element_basis: return result else: first = True for key,entry in self._phase_d.items(): if self.tot_moles_phase(key) > 0.0: nc = entry['nc'] for row in entry['conv_to_elm']: if first: c_matrix = np.reshape(row,(1,row.size)) first = False else: c_matrix = np.vstack((c_matrix,np.reshape( row,(1,row.size)))) Q, R = np.linalg.qr(c_matrix) #Qb = np.matmul(Q.T, np.matmul(result, Q)) Qb = np.matmul(Q.T, result) result, residuals, rank, S = np.linalg.lstsq(R, Qb, rcond=None) return result
[docs] def d3GdT2dn(self, t=1000.0, p=1000.0): """ Returns the third derivative of G of the system twice with respect to temperature and once with respect to mole numbers; a 1-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: nc = entry['nc'] if nc == 1: row = np.array([-entry['obj'].heat_capacity(t, p)/t]) else: moles = entry['moles'] row = -entry['obj'].heat_capacity(t, p, mol=moles, deriv={'dmol':1})[0]/t result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) return result
[docs] def d3GdP2dn(self, t=1000.0, p=1000.0): """ Returns the third derivative of G of the system with respect to pressure twice and mole numbers once; a 1-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: nc = entry['nc'] if nc == 1: row = np.array([entry['obj'].gibbs_energy( t, p, deriv={'dP':2})]) else: moles = entry['moles'] row = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dP':2, 'dmol':1})[0] result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) return result
[docs] def d3GdTdPdn(self, t=1000.0, p=1000.0): """ Returns the third derivative of G of the system with respect to temperature, pressure and mole numbers; a 1-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: nc = entry['nc'] if nc == 1: row = np.array([entry['obj'].gibbs_energy( t, p, deriv={'dT':1, 'dP':1})]) else: moles = entry['moles'] row = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dT':1, 'dP':1, 'dmol':1})[0] result.append(np.reshape(row,(nc,1))) result = np.vstack(tuple(result)) return result
[docs] def d3GdTdn2(self, t=1000.0, p=1000.0): """ Returns the first temperature and second molar derivative of G of system; a 2-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: if entry['nc'] == 1: mat = np.array([[0.0]]) else: moles = entry['moles'] mat = -entry['obj'].entropy(t, p, mol=moles, deriv={'dmol':2})[0] result.append(mat) result = sp.linalg.block_diag(*result) return result
[docs] def d3GdPdn2(self, t=1000.0, p=1000.0): """ Returns the first pressure and second molar derivative of G of system; a 2-D numpy array """ result = [] for entry in self.phase_d.values(): if self.tot_moles_phase(entry['name']) > 0.0: if entry['nc'] == 1: mat = np.array([[0.0]]) else: moles = entry['moles'] mat = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dP':1, 'dmol':2})[0] result.append(mat) result = sp.linalg.block_diag(*result) return result
[docs] def d3GdT3(self, t=1000.0, p=1000.0): """ Returns the third temperature derivative of the Gibbs free energy of the system; a scalar """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: cp = entry['obj'].heat_capacity(t,p) dcpdt = entry['obj'].heat_capacity(t,p,deriv={'dT':1}) result += moles[0]*(cp/t/t - dcpdt/t) else: cp = entry['obj'].heat_capacity(t,p,mol=moles) dcpdt = entry['obj'].heat_capacity(t,p,mol=moles, deriv={'dT':1}) result += cp/t/t - dcpdt/t return result
[docs] def d3GdT2dP(self, t=1000.0, p=1000.0): """ Returns the third pressure derivative of the Gibbs free energy of the system; a scalar """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy( t,p,deriv={'dT':2, 'dP':1}) else: result += entry['obj'].gibbs_energy( t,p,mol=moles,deriv={'dT':2, 'dP':1}) return result
[docs] def d3GdTdP2(self, t=1000.0, p=1000.0): """ Returns the third pressure derivative of the Gibbs free energy of the system; a scalar """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy( t,p,deriv={'dT':1, 'dP':2}) else: result += entry['obj'].gibbs_energy( t,p,mol=moles,deriv={'dT':1, 'dP':2}) return result
[docs] def d3GdP3(self, t=1000.0, p=1000.0): """ Returns the third pressure derivative of the Gibbs free energy of the system; a scalar """ result = 0.0 for entry in self.phase_d.values(): moles = entry['moles'] if np.sum(moles) > 0.0: if entry['nc'] == 1: result += moles[0]*entry['obj'].gibbs_energy( t,p,deriv={'dP':3}) else: result += entry['obj'].gibbs_energy( t,p,mol=moles,deriv={'dP':3}) return result
[docs] def d3Gdn3(self, t=1000.0, p=1000.0, phase=None, cmp=None): """ Returns the third molar derivative of G system for the cmp component of the phase named phase as a 2-D numpy array Parameters ---------- phase : str, def None cmp : int, def None for the phase named phase d3Gdn3[cmp][*][*] is returned """ if phase is None or cmp is None: return None if self.tot_moles_phase(phase) <= 0.0: return None entry = self.phase_d[phase] if entry['nc'] == 1: mat = np.array([[0.0]]) else: moles = entry['moles'] mat = entry['obj'].gibbs_energy(t, p, mol=moles, deriv={'dmol':3})[0] nc = mat.shape[0] result = np.zeros((nc, nc)) for i in range(0,nc): for j in range(i,nc): for k in range(j,nc): if i == cmp: result[j][k] = mat[i][j][k] result[k][j] = mat[i][j][k] elif j == cmp: result[i][k] = mat[i][j][k] result[k][i] = mat[i][j][k] elif k == cmp: result[i][j] = mat[i][j][k] result[j][i] = mat[i][j][k] return result
[docs]class Equilibrate: """Class for minimizing a generic thermodynamic potential in order to calculate an equilibrium phase assemblage. The default potential is the Gibbs free energy. Parameters ---------- element_l : [], default None See documentation for element_list attribute. phase_l : [], default None See documentation for phase_list attribute. lagrange_l : [], default None See documentation for lagrange_list attribute. Attributes ---------- A_omni_inv bulk_comp CTf element_list entropy lagrange_list lagrange_moles lagrange_no_mol_deriv lagrange_use_omni max_linear max_quad moles_in moles_out phase_list phase_separ_threshold reactions rotate_orthog_proj volume VT_null Notes ----- Alternate potentials are specified by stipulating appropriate Lagrange transformations of the Gibbs potential. """ def __init__(self, element_l=None, phase_l=None, lagrange_l=None): self._element_list = element_l self._phase_list = phase_l self._lagrange_list = lagrange_l self._bulk_comp = None self._moles_in = 1.0e-5 self._moles_out = 1.0e-6 self._max_quad = 100 self._max_linear = 200 self._phase_separ_threshold = -0.1 self._entropy = False self._volume = False self._CTf = None self._VT_null = None self._reactions = None self._lagrange_moles = None self._lagrange_use_omni = True self._rotate_orthog_proj = False self._A_omni_inv = None self._lagrange_no_mol_deriv = False @property def element_list(self): """ A list of strings that identify the element basis for the system For example,['H','C','O','Na','Mg','Al','Si','P','K','Ca','Ti','Cr', 'Mn','Fe','Co','Ni'] readonly Returns ------- A Python list """ return self._element_list @property def phase_list(self): """ A list of class instances for phases that are permitted to form These instances must derive from the *PurePhase* or *SolutionPhase* classes, which are defined in the *phases* module. readonly Returns ------- A Python list """ return self._phase_list @property def lagrange_list(self): """ A list of tuples characterizing the Lagrange transformation of the Gibbs free energy to form the operative thermodynamic potential Each tuple has two terms: 1) Either of the following: - A string with content 'T' or 'P' - A dictionary keyed on element symbols with values corresponding to stoichiometric coefficients of element chemical potentials. These linear combinations of potentials are fixed by the constraint. 2) A function to compute values for the constraints specified in (1) The function has a signature func(T, P, state), where *T* is temperature in K, *P* is pressure in bars, and *state* is an instance of the EquilState class. readonly Returns ------- A Python list of tuples """ return self._lagrange_list @property def bulk_comp(self): """ Bulk composition of the system as moles of elements readwrite Returns ------- numpy array """ return self._bulk_comp @bulk_comp.setter def bulk_comp(self, value): assert(type(value) is np.ndarray, 'bulk_comp must be set to an numpy array') assert(value.size == len(self.element_list), 'bulk_comp array must have same length as element_list') self._bulk_comp = value @property def moles_in(self): """ Moles of phase added to the system on detection of saturation readwrite Returns ------- Number of moles added (number) """ return self._moles_in @moles_in.setter def moles_in(self, value): assert(isinstance(value, (int, long, float)), 'ERROR: moles_in must a a number.') self._moles_in = value @property def moles_out(self): """ Minimum total moles of phase allowed in system If there is less than this quantity, the phase is discarded. readwrite Returns ------- Minimum moles allowed (number) """ return self._moles_out @moles_out.setter def moles_out(self, value): assert(isinstance(value, (int, long, float)), 'ERROR: moles_in must a a number.') self._moles_out = value @property def max_quad(self): """ Maximum number of Newton quadratic minimization steps allowed in computation readwrite Returns ------- Number of steps allowed (int) """ return self._max_quad @max_quad.setter def max_quad(self, value): assert(isinstance(value, int), 'ERROR: max_quad must be an integer.') self._max_quad = value @property def max_linear(self): """ Maximum number of linear search steps associated with each quadratic iteration readwrite Returns ------- Number of linear search steps (int) """ return self._max_linear @max_linear.setter def max_linear(self, value): assert(isinstance(value, int), 'ERROR: max_quad must be an integer.') self._max_linear = value @property def phase_separ_threshold(self): """ Minimum energy in Joules for phase separation Gibbs free energy threshold that must be attained in order to add another instance of a phase to the equilibrium assemblage. In the phase unmixing (immiscibility) algorithm, candidate compositions are tested and if one is found with a Gibbs free energy below the projected tangent plane to the input composition by this amount, then the new phase instance is added to the system. The value is in Joules and should always be negative. Values larger than the default will likely encourage the detection of false unmixing events; values more negative will likely prevent the detection of unmixing and result in metastable assemblages. Use caution when altering the default value. Default is -0.1 readwrite Returns ------- Threshold value (number) """ return self._phase_separ_threshold @phase_separ_threshold.setter def phase_separ_threshold(self, value): self._phase_separ_threshold = value @property def entropy(self): """ Indicates if entropy is an independent variable of the minimal potential and if temperature is a dependent variable True if yes; False if reverse. readonly Returns ------- Flag (bool) """ return self._entropy @property def volume(self): """ Indicates if volume is an independent variable of the minimal potential and if pressure is a dependent variable True if yes; False if reverse. readonly Returns ------- Flag (bool) """ return self._volume @property def CTf(self): """ Stoichiometric vectors of element concentrations that embody the Lagrange constraints Default is None readonly Returns ------- np_array or None """ return self._CTf @property def VT_null(self): """ Orthogonal projection operator obtained from constraints on chemical potentials derived from the lagrange_list entries Default is None readonly Returns ------- np array or None """ return self._VT_null @property def reactions(self): """ Balanced reactions pertaining to the Lagrange constraints Default is None readonly Returns ------- np.array or None """ return self._reactions @property def lagrange_moles(self): """ Moles of chemical potential entities stipulated in dictionaries, which are contained in the lagrange_list constraints Default is None readonly Returns ------- np.array or None """ return self._lagrange_moles @property def lagrange_use_omni(self): """ If set to True, the algorithm uses the omnicomponent phase exclusively to balance non-linear constraint reactions. Default is True readwrite Returns ------- Flag (bool) """ return self._lagrange_use_omni @lagrange_use_omni.setter def lagrange_use_omni(self, value): self._lagrange_use_omni = value @property def rotate_orthog_proj(self): """ Prevents creation of a search direction for potential minimization that has unwanted coupling of oxygen-bearing components The basis vectors of the orthogonal projection matrix, VT_null, generated from the lagrange_list attribute, are rotated to zero to minimize the contribution of oxygen to the null space. This flag is enabled to avoid creating a search direction for potential minimization that has unwanted coupling of oxygen-bearing components (for more information, see method _compute_null_space(...)). This option is seldom needed and can be applied only if no omnicomponent phase is present in the system. Default is False readwrite Returns ------- Flag (bool) """ return self._rotate_orthog_proj @rotate_orthog_proj.setter def rotate_orthog_proj(self, value): self._rotate_orthog_proj = value @property def A_omni_inv(self): """ A matrix that transforms an array/matrix of mole numbers of elements to an array/matrix of mole numbers of components of the omnicomponent phase This matrix is the inverse of a matrix that maps elemental abundances to component mole numbers for the omnicomponent phase. This property is available only if chemical potential constraints are specified in lagrange_list and if there is an omnicomponent phase in the system. Default is None readonly Returns ------- 2-d nd_array """ return self._A_omni_inv @property def lagrange_no_mol_deriv(self): """ Produces a simplified gradient and hessian of the generalized Khorzhinskii potential Assume that the generalized Khorzhinkii potential is constructed so that the imposed chemical potential is not a function of mole numbers of any phase in the system. Alternatively, the imposed potential is stoichiometrically assembled using reaction coefficients involving equilibrium phases in the assemblage. This option produces a simplified gradient and hessian of the generalized Khorzhinskii potential, but does not affect construction of the equality constraint matrix nor the Lagrange multiplier terms in the hessian of the Lagrangian function. It is seldom necessary to set this flag to True. Default is False readwrite Returns ------- Flag (bool) """ return self._lagrange_no_mol_deriv @lagrange_no_mol_deriv.setter def lagrange_no_mol_deriv(self, value): self._lagrange_no_mol_deriv = value ################################ # Class private helper methods # ################################ def _print_matrix_repr(self, M): if M.shape[0] == 0 or M.shape[1] == 0: return for i in range(0,M.shape[0]): print (' ', end=' ') for j in range(0,M.shape[1]): print ('-', end=' ') if M[i][j] == 0 else print ('x', end=' ') print (' ') def _calc_sat_state_affinities(self, state, t, p, mu_elm, omni_phase=None, debug=0): """ Determine saturation state affinities for all system phases with zero moles. """ A_min = sys.float_info.max A_min_name = None for entry in state.phase_d.values(): name = entry['name'] if (name != omni_phase) and (state.tot_moles_phase(name) == 0.0): if debug > 0: print ('******************************** ') print ('Calculating saturation state for ' + name) c_mat = state.c_matrix_for_phase(name) nc = entry['nc'] mu_end = np.zeros(nc) for i in range(0,nc): c_row = c_mat[i,:] c_mask = np.not_equal(c_row,0) valid = np.count_nonzero(c_row) == np.count_nonzero( np.extract(c_mask, self.bulk_comp)) mu_end[i] = np.matmul(c_row, mu_elm) if valid else 0.0 if debug > 1: print ('Target chemical potentials: ', mu_end) # Now call a phase method that converts mu_elm to A and X (A,X) = entry['obj'].affinity_and_comp(t, p, mu_end, debug=(True if debug > 2 else False)) if debug > 0: print ('Affinity, mole fraction', A, X) print (' ') entry['affinity'] = A entry['mole_frac'] = X if A < A_min: A_min = A A_min_name = name return (A_min, A_min_name) def _determine_phase_stability(self, state, t, p, phase_name, debug=0): """ Calculate if the specified phase is stable with respect to unmixing. Returns: -------- result: tuple bool, ndarray - unmixing detected; None or composition of second phase """ if phase_name == state.omni_phase(): return False, None phase = state.phase_d[phase_name] moles = phase['moles'] if moles.size == 1: return False, None tot_moles = np.sum(moles) if tot_moles <= 0: return False, None if debug > 0: print ('') print ('Unmixing calculation for', phase_name, moles/tot_moles) mu0 = [] nc = phase['nc'] for i in range(0,nc): x = np.full(nc,0.000001) x[i] = 1.0-np.sum(x) mu0.append(phase['obj'].gibbs_energy(t, p, mol=x)) gHat = phase['obj'].gibbs_energy(t, p, mol=moles) dgdn = phase['obj'].gibbs_energy(t, p, mol=moles, deriv={'dmol':1})[0] for i in range(0,nc): gHat -= mu0[i]*moles[i] dgdn[i] -= mu0[i] gHat /= tot_moles def deltaG(x, *args): state, t, p, phase, moles, gHat, dgdn, mu0 = args f = phase['obj'].gibbs_energy(t, p, mol=x) df = phase['obj'].gibbs_energy(t, p, mol=x, deriv={'dmol':1})[0] for i in range(0,phase['nc']): f -= mu0[i]*x[i] df[i] -= mu0[i] f /= np.sum(x) f -= gHat + np.matmul(x/np.sum(x)-moles/np.sum(moles), dgdn) df -= dgdn return f, df bound_set = [] for i in range(0,nc): bound_set.append((0.01,0.99)) lowest_fun = 0.0 lowest_comp = None debug_flag = True if debug > 1 else False for i in range(0,nc): x = np.full(nc,0.05) x[i] = 1.0-np.sum(x) result = sp.optimize.minimize(deltaG, x, args=(state, t, p, phase, moles, gHat, dgdn, mu0), method='L-BFGS-B', jac=True, bounds=bound_set, options={'maxiter':250, 'disp':debug_flag}) if debug > 0: print (result) print ('') if result.fun < lowest_fun: lowest_fun = result.fun lowest_comp = result.x if lowest_fun < self.phase_separ_threshold: return True, lowest_comp else: return False, None def _add_phase_to_system(self, state, name): """ Add a small amount of phase 'name' to the system and maintain mass balance by substracting amount added from the omnicomponent phase. """ omni_phase_name = state.omni_phase() omni_phase = state.phase_d[omni_phase_name] new_phase = state.phase_d[name] for i in range(0,new_phase['nc']): new_phase['moles'][i] = new_phase['mole_frac'][i]*self.moles_in c_new_phase = state.c_matrix_for_phase(name) # determine elements associated with new_phase elm = np.matmul(c_new_phase.T, new_phase['moles']) # cast these elements into moles of the omnicomponent phase comp = np.matmul( np.linalg.inv(state.c_matrix_for_phase(omni_phase_name).T), elm) # now subtract this quanity from the omnicomponent phase omni_phase['moles'] = np.subtract(omni_phase['moles'], comp) return omni_phase['obj'].test_endmember_comp(omni_phase['moles']) def _remove_phase_from_system(self, state, name): """ Remove a phase 'name' from the system and maintain mass balance by adding amount removed to the omnicomponent phase. """ omni_phase_name = state.omni_phase() omni_phase = state.phase_d[omni_phase_name] old_phase = state.phase_d[name] c_old_phase = state.c_matrix_for_phase(name) # determine elements associated with old_phase elm = np.matmul(c_old_phase.T, old_phase['moles']) # cast these elements into moles of the omnicomponent phase comp = np.matmul( np.linalg.inv(state.c_matrix_for_phase(omni_phase_name).T), elm) # now add this quanity to the omnicomponent phase omni_phase['moles'] = np.add(omni_phase['moles'], comp) # finally, zero mole numbers of the phase removed from the system for i in range(0,old_phase['nc']): old_phase['moles'][i] = 0.0 return omni_phase['obj'].test_endmember_comp(omni_phase['moles']) def _apply_non_linear_constraints(self, tp_l, state, debug): """ Enforce the non-linear constraints (entropy, volume, chemical potential(s)) on the system. """ def linear_search(step, state, t, p, n_ref, row, f): index = 0 for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: for i in range(0,entry['nc']): corr = step*self.reactions[row][index] entry['moles'][i] = n_ref[index] + corr index += 1 dgdn = state.dGdn(t, p, element_basis=False) result = np.matmul(self.reactions[row], dgdn)[0] result -= f(t, p, state) return np.sqrt(result*result) def test_bounds_for_search(step, state, t, p, row): index = 0 for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: nc = entry['nc'] moles = np.empty(nc) for i in range(0,nc): corr = step*self.reactions[row][index] moles[i] = entry['moles'][i] + corr index += 1 if nc == 1: if moles[0] < 0.0: return False else: if not entry['obj'].test_endmember_comp(moles): return False return True def root_target(x, state, t, p, c, f): if c == 'T': result = state.dGdT(x, p) - f(t, p, state) Dresult = state.d2GdT2(x, p) return result, Dresult elif c == 'P': result = state.dGdP(t, x) - f(t, p, state) Dresult = state.d2GdP2(t, x) return result, Dresult else: return None, None t = tp_l[0] p = tp_l[1] row = 0 func_for_row = [] for c,f in self.lagrange_list: if type(c) is dict: if self.reactions.shape[0] > 1: func_for_row.append(f) continue step_max = 0.5 while not test_bounds_for_search(step_max, state, t, p, row): step_max *= 0.9 if step_max < sys.float_info.epsilon: break if debug > 1: print ('Maximum linear step search length; ', step_max) step_min = -0.5 while not test_bounds_for_search(step_min, state, t, p, row): step_min *= 0.9 if step_min > -sys.float_info.epsilon: break if debug > 1: print ('Minimum linear step search length; ', step_min) n_ref = [] for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: for i in range(0,entry['nc']): n_ref.append(entry['moles'][i]) n_ref = np.array(n_ref) result = sp.optimize.minimize_scalar(linear_search, bracket=(0,1), bounds=(step_min,step_max), args=(state, t, p, n_ref, row, f), method='Bounded', options={'maxiter':250, 'disp':debug}) if debug > 1: print ('Optimal linear search length:', result) if np.absolute(result.fun) > 1.0: print ('Imposed chemical potential constraints', end=' ') print ('cannot be reconciled within 1 Joule.') print ('Consider setting the "lagrange_use_omni"', end=' ') print ('property to True.') self._lagrange_moles[row] += result.x row += 1 elif type(c) is str: if self.entropy and self.volume: if c == 'T': tFun = f elif c == 'P': pFun = f continue result = sp.optimize.root_scalar(root_target, args=(state, t, p, c, f), method='newton', x0=t if c == 'T' else p, fprime=True, maxiter=50) if debug > 0: print (result) if c == 'T': tp_l[0] = result.root elif c == 'P': tp_l[1] = result.root if self.entropy and self.volume: # tFun and pFun from above # TBD def roots_fun(x, state, t, p, tFun, pFun): resT = state.dGdT(x[0], x[1]) - tFun(t, p, state) DresTT = state.d2GdT2(x[0], x[1]) resP = state.dGdP(x[0], x[1]) - pFun(t, p, state) DresPP = state.d2GdP2(x[0], x[1]) DresTP = state.d2GdTdP(x[0], x[1]) return [resT,resP], np.array([[DresTT,DresTP],[DresTP,DresPP]]) result = sp.optimize.root(roots_fun, [t,p], args=(state, t, p, tFun, pFun), method='hybr', jac=True) if debug > 0: print (result) tp_l[0] = result.x[0] tp_l[1] = result.x[1] if len(func_for_row) > 0: if debug > 0: print ('Solving with multiple chemical potential constraints.') def multi_search(step, state, t, p, n_ref, func_for_row): index = 0 for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: for i in range(0,entry['nc']): corr = 0.0 for j in range(0,len(func_for_row)): corr += step[j]*self.reactions[j][index] entry['moles'][i] = n_ref[index] + corr index += 1 dgdn = state.dGdn(t, p, element_basis=False) result = [] for row,f in enumerate(func_for_row): result.append(np.matmul(self.reactions[row], dgdn)[0] - f(t, p, state)) return result def test_bounds_for_multi_search(step, state, t, p, func_for_row): index = 0 for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: nc = entry['nc'] moles = np.empty(nc) for i in range(0,nc): corr = 0.0 for j in range(0,len(func_for_row)): corr += step[j]*self.reactions[j][index] moles[i] = entry['moles'][i] + corr index += 1 if nc == 1: if moles[0] < 0.0: return False else: if not entry['obj'].test_endmember_comp(moles): return False return True n_ref = [] for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: for i in range(0,entry['nc']): n_ref.append(entry['moles'][i]) n_ref = np.array(n_ref) if debug > 0: comp = state.oxide_comp('Liquid') for key,value in comp.items(): print ("{0:<6s} {1:6.2f}".format(key, value), end=' ') print ('') result = sp.optimize.root(multi_search, np.zeros(len(func_for_row)), args=(state, t, p, n_ref, func_for_row), method='hybr', jac=False) if debug > 0: comp = state.oxide_comp('Liquid') for key,value in comp.items(): print ("{0:<6s} {1:6.2f}".format(key, value), end=' ') print ('') print (result) state.print_state() for i in range(0,len(func_for_row)): self._lagrange_moles[i] += result.x[i] def _compute_a_and_qr(self, t, p, state, debug): """ From the system state, calculate the mass balance constraint matrix and its QR decomposition. """ filtr = lambda x : x if abs(x) > float(1000*sys.float_info.epsilon) else 0 vfiltr = np.vectorize(filtr, otypes=[float]) first = True for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0.0: if first: c_mat = state.c_matrix_for_phase(name) first = False else: c_mat = np.vstack((c_mat, state.c_matrix_for_phase(name))) A = c_mat.T if self.lagrange_list is not None: if debug > 1: print ('A before projection and augmentation', A.shape) if debug > 2: print (A) else: self._print_matrix_repr(A) A = np.matmul(self.A_omni_inv, A) if debug > 1: print ('A after projection and before augmentation', A.shape) if debug > 2: print (A) else: self._print_matrix_repr(A) if self.CTf.shape[0] > 0: ne = len(self.element_list) omni_phase_name = state.omni_phase() if self.lagrange_use_omni and (omni_phase_name is not None): first = True for phase in state.phase_d.values(): if state.tot_moles_phase(phase['name']) > 0: if phase['name'] == omni_phase_name: entry = np.zeros((ne,phase['nc'])) else: entry = np.ones((ne,phase['nc'])) A_mask = entry if first else np.hstack((A_mask, entry)) first = False A_react = np.ma.array(A, mask=A_mask, fill_value=0) A_react = A_react.filled() else: A_react = np.copy(A) if debug > 1: print ('A reaction matrix', A_react.shape) if debug > 2: print (A_react) else: self._print_matrix_repr(A_react) # Create a mass balance constraint on constrained entity react, res, rank, s = np.linalg.lstsq(A_react, np.matmul(self.A_omni_inv,self.CTf.T), rcond=None) self._reactions = vfiltr(react.T) if not self.entropy and not self.volume: if debug > 0: print (".....") print ('Enforcing chemical potential constraints.') self._apply_non_linear_constraints([t,p], state, debug) if debug > 0: print (".....") if debug > 1: print ('Balanced Lagrange reactions', self._reactions.shape) print (self._reactions) if debug > 2: print ('... residuals', res.shape) print (res) print ('... rank', rank) print ('... singular values', s.shape) print (s) if debug > 1: dgdn = state.dGdn(t, p, element_basis=False) deltaG = np.matmul(self.reactions, dgdn) row = 0 for c,f in self.lagrange_list: if type(c) is dict: deltaG[row] -= f(t, p, state) row += 1 print ('Lagrange chemical potential constraints', deltaG.shape) print (deltaG) # Project the A matrix into the NULL space Aproj = np.matmul(self.VT_null, A) if debug > 1: print ('Aproj', Aproj.shape) if debug > 2: print (Aproj) else: self._print_matrix_repr(Aproj) # Add rows for the constraint ###### d2gdn2 = state.d2Gdn2(t, p, element_basis=True) ###### con_rows = np.matmul(self.CTf, d2gdn2) d2gdn2 = state.d2Gdn2(t, p, element_basis=False) con_rows = np.matmul(self._reactions, d2gdn2) if debug > 2: print ('Constraint rows', con_rows.shape) print (con_rows) Aproj = np.vstack((Aproj, con_rows)) if debug > 1: print ('A proj augmented with constraints', Aproj.shape) if debug > 2: print (Aproj) else: self._print_matrix_repr(Aproj) A = Aproj if self.entropy and self.volume: d2GdTdn = state.d2GdTdn(t, p) d2GdT2 = state.d2GdT2(t, p) d2GdPdn = state.d2GdPdn(t, p) d2GdP2 = state.d2GdP2(t, p) d2GdTdP = state.d2GdTdP(t, p) A = np.vstack((A, d2GdTdn.T, d2GdPdn.T)) rows,cols = A.shape col = np.zeros((rows,2)) col[-2][0] = d2GdT2 col[-1][1] = d2GdP2 col[-2][1] = d2GdTdP col[-1][0] = d2GdTdP A =np.hstack((A,col)) elif self.entropy: d2GdTdn = state.d2GdTdn(t, p) d2GdT2 = state.d2GdT2(t, p) A = np.vstack((A, d2GdTdn.T)) rows,cols = A.shape col = np.zeros((rows,1)) col[-1][0] = d2GdT2 A =np.hstack((A,col)) elif self.volume: d2GdPdn = state.d2GdPdn(t, p) d2GdP2 = state.d2GdP2(t, p) A = np.vstack((A, d2GdPdn.T)) rows,cols = A.shape col = np.zeros((rows,1)) col[-1][0] = d2GdP2 A = np.hstack((A,col)) if debug > 1: print ('A augmented with constraints', A.shape) if debug > 2: print (A) else: self._print_matrix_repr(A) row,col = A.shape df = col - row R, Q = sp.linalg.rq(A, mode='full') if debug > 1: print ('RQ = A, R matrix', R.shape) if debug > 2: print (vfiltr(R)) else: self._print_matrix_repr(vfiltr(R)) if debug > 1: print ('RQ = A, Q matrix', Q.shape) if debug > 2: print (vfiltr(Q)) else: self._print_matrix_repr(vfiltr(Q)) R11 = vfiltr(R[:,df:]) Q1 = vfiltr(Q[df:,:]) Q2 = vfiltr(Q[0:df,:]) return (A, df, Q1, Q2, R11) def _compute_null_space(self, state, debug): """ Construct a projection operator to account for open system constraints """ ne = len(self.element_list) CTf = np.empty((0,ne)) for c,f in self.lagrange_list: if type(c) is str: if c == 'T': self._entropy = True elif c == 'P': self._volume = True else: assert(True,'Legendre transform must specify "T", "P", ' + 'or a dictionary of chemical potential constraints.') elif type(c) is dict: row = np.zeros(ne) for key,value in c.items(): index = self.element_list.index(key) row[index] = value CTf = np.vstack((CTf, row)) else: assert(True,'Legendre transform must specify "T", "P", ' + 'or a dictionary of chemical potential constraints.') self._CTf = CTf self._lagrange_moles = np.zeros(len(self.lagrange_list)) if state.omni_phase() is not None: self._A_omni_inv = np.linalg.inv( state.c_matrix_for_phase(state.omni_phase()).T) else: self._A_omni_inv = np.eye(ne) if debug > 1: print ('CTf', CTf.shape) print (CTf) print ('A_omni_inv', self.A_omni_inv.shape) if debug > 2: print (self.A_omni_inv) else: self._print_matrix_repr(self.A_omni_inv) CTf_pad = np.pad(CTf,((0,CTf.shape[1]-CTf.shape[0]),(0,0)),mode='constant') U,S,VT = np.linalg.svd(np.matmul(CTf_pad,self.A_omni_inv.T)) if debug > 1: print ('U', U.shape) if debug > 2: print (U) else: self._print_matrix_repr(U) print ('S', S.shape) if debug > 2: print (S) print ('VT', VT.shape) if debug > 2: print (VT) else: self._print_matrix_repr(VT) rank = np.count_nonzero(S) self._VT_null = VT[rank:,:] if debug > 1: print ('VT null space', self._VT_null.shape) if debug > 2: print (self._VT_null) else: self._print_matrix_repr(self._VT_null) if self.rotate_orthog_proj and state.omni_phase() is None: try: row_ind = self.element_list.index("O") except ValueError: if debug > 1: print ("ERROR! Oxygen is not present in the system and ") print (" a null-space rotation on O is specified.") return if debug > 1: print ("----------") print ("Adjusting the null space basis for oxygen.") ns = self._VT_null.T filtr = lambda x : x if abs(x) > float( 1000*sys.float_info.epsilon) else 0 vfiltr = np.vectorize(filtr, otypes=[float]) # zero the null space bassis coefficient at col_ind and row_ind def zeros(deg_a, v_a, u, ns, col_ind_a, row_ind, ret_matrix): ns_p = ns for (deg,v,col_ind) in zip(deg_a,v_a,col_ind_a): theta = (deg/180.0)*np.pi P = np.eye(v.shape[0]) + np.sin(theta)*(np.outer(u,v) -np.outer(v,u)) + (np.cos(theta)-1)*(np.outer(u,u) -np.outer(v,v)) ns_p = np.matmul(P,ns_p) if ret_matrix: return ns_p else: result = [] for col_ind in col_ind_a: result.append(ns_p[row_ind][col_ind]) return np.array(result) u = np.zeros(ns.shape[0]) u[row_ind] = 1 v_a = [] col_ind_a = [] guess = [] for col in range(0,ns.shape[1]): if ns[row_ind][col] != 0.0: v_a.append(ns[:,col]) col_ind_a.append(col) guess.append(0) result = sp.optimize.root(zeros, args=(v_a, u, ns, col_ind_a, row_ind, False), x0=np.array(guess)) if debug > 1: print ("... rotation algoritm output:") print (result) ns_p = vfiltr(zeros(result.x, v_a, u, ns, col_ind_a, row_ind, True)) self._VT_null = ns_p.T if debug > 1: print ('rotated null space basis', self._VT_null.shape) if debug > 2: print (self._VT_null) else: self._print_matrix_repr(self._VT_null) print ("----------") def _augment_gradient_hessian(self, t, p, g, H, state, debug): """ Ammend the gradient and Hessian of the Gibbs free energy to account for the Khorzhinskii corections. Returns ------- result: tuple corrected (g, H) """ dgdn = np.copy(g) d2gdn2 = np.copy(H) rank_H = H.shape[0] idx = 0 for c,f in self.lagrange_list: if type(c) is dict: # determine molar equivalent of open system reactant moles = state.moles_v(self.reactions[idx]) # determine the first order correction to the gradient gAdd_1 = np.matmul(self.reactions[idx], dgdn)[0] gAdd_1 *= self.reactions[idx] gAdd_1 = np.reshape(gAdd_1, (gAdd_1.shape[0],1)) if debug > 2: print ('moles', moles) print ('1st gradient addition', gAdd_1.shape) print (-gAdd_1) # determine the second order correction to the gradient if not self.lagrange_no_mol_deriv: gAdd_2 = np.matmul(self.reactions[idx], d2gdn2) gAdd_2 = np.reshape(gAdd_2, (gAdd_2.shape[0],1)) if debug > 2: print ('2nd gradient addition', gAdd_2.shape) print (-moles*gAdd_2) else: gAdd_2 = np.zeros(g.shape) # add contributions to gradient g = g - gAdd_1 - moles*gAdd_2 # add contributions to Hessian # ... first, index the structure of the hessian phase_info = [] offset = 0 for entry in state.phase_d.values(): if state.tot_moles_phase(entry['name']) > 0.0: for idx2 in range(0,entry['nc']): phase_info.append((entry['name'],offset,idx2)) offset += entry['nc'] if debug > 1: print ('phase_info', len(phase_info), 'reactions', self.reactions.shape) print (phase_info) # ... second, build the Hessian additions if not self.lagrange_no_mol_deriv: hAdd_1 = np.outer(self.reactions[idx],gAdd_2) hAdd_2 = np.outer(gAdd_2, self.reactions[idx]) hAdd_3 = np.zeros(H.shape) for idx2,coeff in enumerate(self.reactions[idx]): if coeff != 0: name,offset,col = phase_info[idx2] HAdd = coeff*state.d3Gdn3(t, p, name, cmp=col) idx3 = rank_H - HAdd.shape[0] - offset hAdd_3 += np.pad(HAdd,((offset,idx3),(offset,idx3)), 'constant',constant_values=0) if debug > 2: print ('1st hessian addition', hAdd_1.shape) print (-hAdd_1) print ('2nd hessian addition', hAdd_2.shape) print (-hAdd_2) print ('3rd hessian addition', hAdd_3.shape) print (-moles*hAdd_3) H += -hAdd_1 - hAdd_2 - moles*hAdd_3 idx += 1 if self.entropy and self.volume: g = g - t*state.d2GdTdn(t, p) g = g - p*state.d2GdPdn(t, p) d2gdtdp = state.d2GdTdP(t, p) g = np.vstack((g, -t*state.d2GdT2(t, p) - p*d2gdtdp)) g = np.vstack((g, -p*state.d2GdP2(t, p) - t*d2gdtdp)) H = H - t*state.d3GdTdn2(t, p) - p*state.d3GdPdn2(t, p) d3gdt2dn = state.d3GdT2dn(t, p) d3gdp2dn = state.d3GdP2dn(t, p) d3gdtdpdn = state.d3GdTdPdn(t, p) colT = -t*d3gdt2dn - p*d3gdtdpdn colP = -p*d3gdp2dn - t*d3gdtdpdn H = np.hstack((H, colT, colP)) d3gdt2dp = state.d3GdT2dP(t,p) d3gdtdp2 = state.d3GdTdP2(t,p) termTT = -state.d2GdT2(t,p) - t*state.d3GdT3(t,p) - p*d3gdt2dp termPP = -state.d2GdP2(t,p) - p*state.d3GdP3(t,p) - t*d3gdtdp2 termTP = -state.d2GdTdP(t,p) - t*d3gdt2dp - p*d3gdtdp2 rowT = np.hstack((colT.T, np.array([termTT],ndmin=2), np.array([termTP],ndmin=2))) rowP = np.hstack((colP.T, np.array([termTP],ndmin=2), np.array([termPP],ndmin=2))) H = np.vstack((H, rowT, rowP)) elif self.entropy: g = g - t*state.d2GdTdn(t, p) g = np.vstack((g, -t*state.d2GdT2(t, p))) H = H - t*state.d3GdTdn2(t, p) col = -t*state.d3GdT2dn(t, p) H = np.hstack((H, col)) term = -state.d2GdT2(t,p) - t*state.d3GdT3(t,p) row = np.hstack((col.T, np.array([term],ndmin=2))) H = np.vstack((H,row)) elif self.volume: g = g - p*state.d2GdPdn(t, p) g = np.vstack((g, -p*state.d2GdP2(t, p))) H = H - p*state.d3GdPdn2(t, p) col = -p*state.d3GdP2dn(t, p) H = np.hstack((H, col)) term = -state.d2GdP2(t,p) - p*state.d3GdP3(t,p) row = np.hstack((col.T, np.array([term],ndmin=2))) H = np.vstack((H,row)) if debug > 1: print ('Augmented g', g.shape) if debug > 2: print (g) print ('Augmented H', H.shape) if debug > 2: print (H) else: self._print_matrix_repr(H) return (g,H) def _compute_lagrange_multipliers(self, g, A, debug): """ Estimate Lagrange multipliers from the gradient and equality constraint derivative matrix. """ soln, res, rank, s = np.linalg.lstsq(A.T, g, rcond=None) if debug > 1: print ('Lagrange multiplier solution', rank) print ('... solution', soln.shape) print (soln) if debug > 2: print ('... residuals', res.shape) print (res) print ('... singular values', s.shape) print (s) return soln def _augment_hessian_return_wronskian(self, t, p, H, l_mult, state, debug): """ Augment the Hessian matrix with Lagrange multiplier terms for the non-linear constraints. """ if self.reactions is not None: rows = self.reactions.shape[0] row = 0 for i in range(0,rows): lagrange = l_mult[l_mult.size-rows+row] react = self.reactions[i] col = 0 result = [] if debug > 1: print ('Lagrange multiplier:', lagrange) for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0.0: nc = entry['nc'] if entry['nc'] == 1: coeff = react[col] mat = coeff*np.array([[0.0]]) col += 1 else: mat = np.zeros((nc,nc)) for j in range(0,nc): if react[col] != 0.0: coeff = react[col] mat_r = state.d3Gdn3(t, p, phase=name, cmp=j) mat = np.add(mat, coeff*mat_r) col += 1 result.append(mat) result = -lagrange*sp.linalg.block_diag(*result) if debug > 1: print ('Hessian addition:', result.shape) if debug > 2: print (result) H = np.add(H, result) row += 1 if self.entropy and self.volume: lam_T = l_mult[-2] lam_P = l_mult[-1] result = -lam_T*state.d3GdTdn2(t, p) - lam_P*state.d3GdPdn2(t, p) d3gdtdpdn = state.d3GdTdPdn(t,p) col_T = -lam_T*state.d3GdT2dn(t,p) - lam_P*d3gdtdpdn col_P = -lam_T*d3gdtdpdn - lam_P*state.d3GdP2dn(t,p) result = np.hstack((result, col_T, col_P)) d3gdt2dp = state.d3GdT2dP(t,p) d3gdtdp2 = state.d3GdTdP2(t,p) term_TT = -lam_T*state.d3GdT3(t,p) - lam_P*d3gdt2dp term_TP = -lam_T*d3gdt2dp - lam_P*d3gdtdp2 term_PP = -lam_T*d3gdtdp2 - lam_P*state.d3GdP3(t,p) row_T = np.hstack((col_T.T, np.array([term_TT],ndmin=2), np.array([term_TP],ndmin=2))) row_P = np.hstack((col_P.T, np.array([term_TP],ndmin=2), np.array([term_PP],ndmin=2))) result = np.vstack((result, row_T, row_P)) if debug > 1: print ('Hessian addition (S,V):', result.shape) if debug > 2: print (result) H = np.add(H, result) elif self.entropy: lam_T = l_mult[-1] result = -lam_T*state.d3GdTdn2(t, p) col_T = -lam_T*state.d3GdT2dn(t,p) result = np.hstack((result, col_T)) term_TT = -lam_T*state.d3GdT3(t,p) row_T = np.hstack((col_T.T, np.array([term_TT],ndmin=2))) result = np.vstack((result, row_T)) if debug > 1: print ('Hessian addition (S):', result.shape) if debug > 2: print (result) H = np.add(H, result) elif self.volume: lam_P = l_mult[-1] result = -lam_P*state.d3GdPdn2(t, p) col_P = -lam_P*state.d3GdP2dn(t,p) result = np.hstack((result, col_P)) term_PP = -lam_P*state.d3GdP3(t,p) row_P = np.hstack((col_P.T, np.array([term_PP],ndmin=2))) result = np.vstack((result, row_P)) if debug > 1: print ('Hessian addition (V):', result.shape) if debug > 2: print (result) H = np.add(H, result) return H ######################## # Class public methods # ########################
[docs] def execute(self, t=1000.0, p=1000.0, bulk_comp=None, state=None, debug=0): """ Calculates an equilibrium assemblage, returning the results as an instance of the EquilState class Parameters ---------- t : float, default 1000.0 Temperature in Kelvins p : float, default 1000.0 Pressure in bars bulk_comp : numpy array, default None Bulk composition of the system. Array of element concentrations, of the length and order of self.element_list state : EquilState, default None An instance of the EquilState class generated by this method. This parameter is only specified for a sequence of equilibrium calculations where the state of the system is initialized from a previous computation. debug : int, default 0 Level of detail printed by the method: - 0, no information - 1, minimal progress information - 2, normal debugint output level - 3, verbose debuging output level Returns ------- state: EquilState An instance of the EquilState class that contains the equilibrium state of the system Notes ----- Both the bulk_comp and state parameters cannot be set simultaneously. """ if state is not None and bulk_comp is not None: assert(True, 'Both "bulk_comp" and "state" parameters cannot be specified.') if bulk_comp is not None: self.bulk_comp = bulk_comp if self.bulk_comp is None and state is None: print ('Bulk composition of the system must be set.') return None # Initialize the EquilState class instance if one is not specified require_quad_pass_before_phase_addition = False if state is None: state = EquilState(self.element_list, self.phase_list) # Currently, one phase must be omnicomponent omni_phase = state.omni_phase() if omni_phase is None: print ("There is no omnicomponent phase in the system.") return None # Assume the system inititally has just the omnicomponent phase state.set_phase_comp(omni_phase, self.bulk_comp, input_as_elements=True) force_once_quad_pass = False state.temperature = t state.pressure = p if self.lagrange_list is not None: self._compute_null_space(state, debug) force_once_quad_pass = True if not self.entropy and not self.volume: #imposed chemical potential constraints require_quad_pass_before_phase_addition = True else: omni_phase = state.omni_phase() force_once_quad_pass = True state.temperature = t state.pressure = p if debug > 0: print (omni_phase + ' is the omnicomponent phase.') ne = len(self.element_list) if self.entropy or self.volume: tp_l = [t, p] self._apply_non_linear_constraints(tp_l, state, debug) (t,p) = tuple(tp_l) state.temperature = t state.pressure = p if debug > 0 and self.entropy: print ('Final corrected temperature', t) if debug > 0 and self.volume: print ('Final corrected pressure', p) ########################### # Calculation starts here # ########################### while True: # Obtain chemical potentials of elements in the system mu_elm = state.dGdn(t, p, element_basis=True, use_omni_phase=True) # If a phase has zero moles, calculate its chemical affinity A_min,A_min_name = self._calc_sat_state_affinities(state, t, p, mu_elm, omni_phase=omni_phase, debug=debug) # if a phase is supersaturated, add it to the system if A_min >= 0.0 and not force_once_quad_pass: if debug > 0: print ('Stable phase assemblage computed. Exiting.') break force_once_quad_pass = False if A_min < 0.0 and not require_quad_pass_before_phase_addition: if debug > 0: print ('Adding phase ', A_min_name, ' to system. with affinity = ', A_min) okay = self._add_phase_to_system(state, A_min_name) if okay and debug > 1: print ('... phase successfully added.') elif not okay: print ('... phase addition unsuccessful. Exiting.') break require_quad_pass_before_phase_addition = False quad_loop = True quad_iter = 0 quad_energies = [] while quad_loop: # Deal with equality/inequality constraints if (quad_iter == 0) or self.lagrange_list is not None: A, df, Q1, Q2, R11 = self._compute_a_and_qr(t, p, state, debug) if debug > 0: print ('Exit constraint matrix and RQ decomposition.') if debug > 1: print ('... df', df, 'A', A.shape) if debug > 2: print (A) else: self._print_matrix_repr(A) if debug > 1: print ('... Q1', Q1.shape) if debug > 2: print (Q1) else: self._print_matrix_repr(Q1) if debug > 1: print ('... Q2', Q2.shape) if debug > 2: print (Q2) else: self._print_matrix_repr(Q2) if debug > 1: print ('... R11', R11.shape) if debug > 2: print (R11) else: self._print_matrix_repr(R11) # Compute gradient and Hessian and their projection into the # null space g = state.dGdn(t, p, element_basis=False) H = state.d2Gdn2(t, p, element_basis=False) if debug > 1: print ('Gradient', g.shape) if debug > 2: print (g) if debug > 1: print ('Hessian', H.shape) if debug > 2: print (H) else: self._print_matrix_repr(H) if self.lagrange_list is not None: g,H = self._augment_gradient_hessian(t, p, g, H, state, debug) l_mult = self._compute_lagrange_multipliers(g, A, debug) H = self._augment_hessian_return_wronskian(t, p, H, l_mult, state, debug) if debug > 1: print ('Final Hessian', H.shape) if debug > 2: print (H) else: self._print_matrix_repr(H) g_p = np.matmul(Q2, g) H_p = np.matmul(np.matmul(Q2, H), Q2.T) if debug > 1: print ('Projected gradient:', g_p.shape) if debug > 2: print (g_p) if debug > 1: print ('Projected Hessian:', H_p.shape) if debug > 2: print (H_p) else: self._print_matrix_repr(H_p) # Solve the null space problem result, residuals, rank, S = np.linalg.lstsq(H_p, -g_p, rcond=None) if debug > 0: print ('') print ('Quadratic solution, rank', rank) if debug > 1: print (result) if debug > 2: print ('... residuals:') print (residuals) print ('... S:') print (S) n2 = np.matmul(Q2.T, result) if debug > 0: print ('Re-inflated delta n solution vector:', n2.shape) for ind in range(0,n2.shape[0]): if n2[ind][0] == 0.0: print ('0.0',end=' ') else: print ("{0:10.3e}".format(n2[ind][0]),end=' ') print ('') norm = np.linalg.norm(n2) if debug > 0: print ('Norm: ', norm) if norm < 1.0e-6: quad_normal_exit = True break # Perform a linear search along the quadratic search direction n_ref = [] for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: for i in range(0,entry['nc']): n_ref.append(entry['moles'][i]) if self.lagrange_list is not None: row = 0 for c,f in self.lagrange_list: if type(c) is dict: n_ref.append(self.lagrange_moles[row]) elif type(c) is str: n_ref.append(t if c == 'T' else p) row += 1 n_ref = np.array(n_ref) if debug > 2: print ('Reference phase/component moles', n_ref.shape) print (n_ref) # Potential function to be minimized def linear_search(step, state, t, p, n_ref, n2): index = 0 for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: for i in range(0,entry['nc']): entry['moles'][i] = n_ref[index] + step*n2[index] index += 1 if self.entropy: t = n_ref[index] + step*n2[index] index += 1 if self.volume: p = n_ref[index] + step*n2[index] index += 1 if self.entropy or self.volume: tp_l = [t, p] self._apply_non_linear_constraints(tp_l, state, debug) (t,p) = tuple(tp_l) result = state.G(t, p) if self.lagrange_list is not None: dgdn = state.dGdn(t, p, element_basis=False) row = 0 for c,f in self.lagrange_list: if type(c) is dict: moles = state.moles_v(self.reactions[row]) contrib = moles*( np.matmul(self.reactions[row], dgdn)[0]) if debug > 1: print ('... linear search. ', end=' ') print ('ref = ', n_ref[index], end=' ') print ('index = ', index, 'i = ',row, end=' ') print ('moles = ', moles, end=' ') print ('contrib ', contrib) result -= contrib self.lagrange_moles[row] = moles index += 1 row += 1 if self.entropy: result -= t*state.dGdT(t,p) if self.volume: result -= p*state.dGdP(t,p) return result # Determine domain bounds on the potential linear search def test_bounds_for_search(step, state, t, p, n_ref, n2): index = 0 for name,entry in state.phase_d.items(): if state.tot_moles_phase(name) > 0: nc = entry['nc'] moles = np.empty(nc) for i in range(0,nc): moles[i] = n_ref[index] + step*n2[index] index += 1 if nc == 1: if moles[0] < 0.0: return False else: if not entry['obj'].test_endmember_comp(moles): return False if self.entropy: t_est = n_ref[index] + step*n2[index] if (t_est < 773.15) or (t_est > 2773.15): return False index += 1 if self.volume: p_est = n_ref[index] + step*n2[index] if (p_est < 1.0) or (p_est > 100000.0): return False index += 1 return True step_max = 2.0 while not test_bounds_for_search(step_max, state, t, p, n_ref, n2): step_max *= 0.9 if step_max < sys.float_info.epsilon: break if debug > 1: print ('Maximum linear step search length; ', step_max) step_min = -2.0 while not test_bounds_for_search(step_min, state, t, p, n_ref, n2): step_min *= 0.9 if step_min > -sys.float_info.epsilon: break if debug > 1: print ('Minimum linear step search length; ', step_min) result = sp.optimize.minimize_scalar(linear_search, bounds=(step_min,step_max), args=(state, t, p, n_ref, n2), method='Bounded', options={'maxiter':250, 'disp':debug, 'xatol':1e-5}) step_val = result.x[0] if isinstance(result.x, np.ndarray ) else result.x func_val = result.fun[0] if isinstance(result.fun, np.ndarray ) else result.fun quad_energies.append(func_val) if self.entropy and self.volume: t_est = n_ref[-2] + step_val*n2[-2] p_est = n_ref[-1] + step_val*n2[-1] t = t_est[0] if isinstance(t_est, np.ndarray) else t_est p = p_est[0] if isinstance(p_est, np.ndarray) else p_est elif self.entropy: t_est = n_ref[-1] + step_val*n2[-1] t = t_est[0] if isinstance(t_est, np.ndarray) else t_est elif self.volume: p_est = n_ref[-1] + step_val*n2[-1] p = p_est[0] if isinstance(p_est, np.ndarray) else p_est if debug > 0: print ('Optimal linear search length:', step_val) print ('Optimal function value', func_val) state.print_state() if debug > 1: print (result) if self.entropy: print ('Approximate corrected temperature', t) if self.volume: print ('Approximate corrected pressure', p) if self.reactions is not None: blk_cmp = state.tot_moles_elements() diff = blk_cmp - self.bulk_comp self._lagrange_moles += np.matmul(self.CTf, diff) if debug > 1: print ('Change in bulk composition:', diff.shape) print (diff) print ('lagrange_moles', self.lagrange_moles.shape) print (self.lagrange_moles) if self.entropy or self.volume: tp_l = [t, p] self._apply_non_linear_constraints(tp_l, state, debug) (t,p) = tuple(tp_l) state.temperature = t state.pressure = p if debug > 0 and self.entropy: print ('Final corrected temperature', t) if debug > 0 and self.volume: print ('Final corrected pressure', p) quad_iter += 1 if quad_iter > self.max_quad: quad_normal_exit = False break if len(quad_energies) > 4: last_three = quad_energies[-3:] average = (last_three[0]+last_three[1]+last_three[2])/3.0 sumsqr = (last_three[0]-average)**2 + (last_three[1]- average)**2 + (last_three[2]-average)**2 if np.sqrt(sumsqr) < 0.1: quad_normal_exit = False print ('Minimal energy termination of quad loop.') break for name,entry in state.phase_d.items(): mol_p = state.tot_moles_phase(name) if mol_p > 0.0 and mol_p < self.moles_out: self._remove_phase_from_system(state, name) if debug > 0: print ('----------') print ('--> Phase', name, 'removed from the system.') print ('----------') quad_iter = 0 # termination of quad_loop if quad_normal_exit: for phase_name in state.phase_d.keys(): result, result_comp = self._determine_phase_stability(state, t, p, phase_name, debug=debug) if result: break if result: if debug > 0: print ('Phase separation found for', phase_name) print (' composition:', result_comp) new_phase = phase_name + '_1' ident = 1 while new_phase in state.phase_d: ident += 1 new_phase = phase_name + '_' + str(ident) state.phase_d[new_phase] = state.phase_d[phase_name].copy() entry = state.phase_d[new_phase] entry['moles'] = np.zeros(entry['nc']) entry['mole_frac'] = np.zeros(entry['nc']) for i in range(0,entry['nc']): entry['mole_frac'][i] = result_comp[i] state._c_matrix = np.vstack( (state._c_matrix, state.c_matrix_for_phase(phase_name))) if debug > 2: print ('Old phase dictionary ...') print (state.phase_d[phase_name]) print ('New phase dictionary ...') print (state.phase_d[new_phase]) okay = self._add_phase_to_system(state, new_phase) if okay and debug > 1: print ('... phase successfully added.') elif not okay: print ('... phase addition unsuccessful. Exiting.') force_once_quad_pass = True else: print ('No check made for phase separation.') # termination of phase saturation loop return state
[docs]class MELTSmodel: """Class for creating an instance of the Equilibrate PhaseObjC class that is tailored to calculate equilibrium phase assemblages using one of the MELTS calibrations. Valid initializers are version='1.0.2', '1.1.0', '1.2.0', '5.6.1', 'DEW', 'OnlyDEW' """ def __init__(self, version='1.0.2'): if version == '1.0.2': MELTSclass = ObjCClass('EquilibrateUsingMELTSv102') elif version == '1.1.0': MELTSclass = ObjCClass('EquilibrateUsingMELTSv110') elif version == '1.2.0': MELTSclass = ObjCClass('EquilibrateUsingMELTSv120') elif version == '5.6.1': MELTSclass = ObjCClass('EquilibrateUsingpMELTSv561') elif version == 'DEW': MELTSclass = ObjCClass('EquilibrateUsingMELTSwithDEW') elif version == 'OnlyDEW': MELTSclass = ObjCClass('EquilibrateUsingMELTSandOnlyDEW') else: assert False, 'Unknown version of MELTS stipulated' self.melts = MELTSclass.alloc().init() self.melts.setDebugS_(0) self.melts.setDebugV_(0) oxide_names_NSArray = self.melts.oxideNames() self.noxides = len(oxide_names_NSArray) # was .count converted to len() in going from 0.2.7 -> 0.2.10 self.oxide_names_a = [oxide_names_NSArray.objectAtIndex_(i) for i in range(self.noxides)] phase_names_NSArray = self.melts.phaseNames() self.nphases = len(phase_names_NSArray) # was .count converted to len() in going from 0.2.7 -> 0.2.10 self.phase_names_a = [phase_names_NSArray.objectAtIndex_(i) for i in range(self.nphases)] locale.setlocale(locale.LC_NUMERIC, '')
[docs] def print_default_settings(self): """Prints list of options and tolerance settings for algorithm. """ option_settings = { 'PPOPTIONS_DISTRIBUTION': ('Make failure in element redistribution non-fatal', self.melts.PPOPTIONS_DISTRIBUTION.boolValue), 'PPOPTIONS_CHEM_POTENTIAL': ('Make chemical potential recast errors non-fatal', self.melts.PPOPTIONS_CHEM_POTENTIAL.boolValue), 'PPOPTIONS_RANK': ('Make rank deficient quadratic solutions fatal', self.melts.PPOPTIONS_RANK.boolValue), 'PPOPTIONS_ZERO_LINEAR': ('Make "zero" steplength linear searches non-fatal', self.melts.PPOPTIONS_ZERO_LINEAR.boolValue), 'PPOPTIONS_QUAD_ITERS': ('Terminate if quadratic iterations ever exceed maximum', self.melts.PPOPTIONS_QUAD_ITERS.boolValue), 'PPOPTIONS_FORCE_COMPONENTS': ('Never zero a component concentration in a system phase', self.melts.PPOPTIONS_FORCE_COMPONENTS.boolValue), 'PPOPTIONS_DETECT_MINIMAL': ('Prevent minimal energy test for convergence', self.melts.PPOPTIONS_DETECT_MINIMAL.boolValue) } for key in option_settings: value = option_settings[key] print ("{0:<60.60s} {1:<30.30s} {2:1d}".format(value[0], key, value[1])) tolerance_settings = { 'PPPARAMETERS_CHEM_POTENTIAL': \ ('Residual tolerance for acceptable recasting of phase chemical potentials to those of the elements', \ float(self.melts.PPPARAMETERS_CHEM_POTENTIAL)), 'PPPARAMETERS_RANK': \ ('Tolerance (relative to Hessian norm) for computation of pseudorank during quadratic search', \ self.melts.PPPARAMETERS_RANK), 'PPPARAMETERS_QUAD_OPTIMAL': \ ('Residual norm relative tolerance for optimal quadratic convergence', \ self.melts.PPPARAMETERS_QUAD_OPTIMAL), 'PPPARAMETERS_QUAD_SUBOPTIMAL': \ ('Residual norm relative tolerance for sub-optimal quadratic convergence', \ self.melts.PPPARAMETERS_QUAD_SUBOPTIMAL), 'PPPARAMETERS_QUAD_MAX_ITERS': \ ('Maximal number of quadratic iterations', \ self.melts.PPPARAMETERS_QUAD_MAX_ITERS), 'PPPARAMETERS_LINEAR_MAX_ITERS': \ ('Maximal number of linear search iterations', \ self.melts.PPPARAMETERS_LINEAR_MAX_ITERS), 'PPPARAMETERS_LINEAR_RELTEST': \ ('Relative convergence criteria for estimation of the minimum during a linear search step', \ self.melts.PPPARAMETERS_LINEAR_RELTEST), 'PPPARAMETERS_LINEAR_MIN_STEPLENGTH': \ ('Minimum allowed steplength in a linear search step', \ self.melts.PPPARAMETERS_LINEAR_MIN_STEPLENGTH), 'PPPARAMETERS_COMPONENT_MINIMUM': \ ('Minimum molar concentration of a component in a phase (absolute value)', \ self.melts.PPPARAMETERS_COMPONENT_MINIMUM), 'PPPARAMETERS_MASSOUT': \ ('Phase mass that triggers removal of that phase from the system', \ self.melts.PPPARAMETERS_MASSOUT), 'PPPARAMETERS_MOLESIN': \ ('Moles of phase added to system on detection of phase saturation', \ self.melts.PPPARAMETERS_MOLESIN), 'PPPARAMETERS_LINEAR_MINIMAL': \ ('Number of quadratic iterations over which to average potential for minimal energy convergence test', \ self.melts.PPPARAMETERS_LINEAR_MINIMAL), 'PPPARAMETERS_CYCLES_MAX': \ ('Maximal number of phase inclusion/phase removal cycles permitted', \ self.melts.PPPARAMETERS_CYCLES_MAX) } for key in tolerance_settings: value = tolerance_settings[key] print ("{0:<70.70s} {1:<30.30s} {2:13.6e}".format(value[0], key, value[1]))
[docs] def set_debug_state(self, debugS=False, debugV=False): """Sets debug output level for Equilibrate class and subclasses Parameters ========== debugS : boolean, optional Sets on or off low level debug output. Default is off (False). debugV : boolean, optional Sets on or off high level debug output. Default is off (False). """ if debugS == False: self.melts.setDebugS_(False) else: self.melts.setDebugS_(True) if debugV == False: self.melts.setDebugV_(False) else: self.melts.setDebugV_(True)
[docs] def get_oxide_names(self): """Retrieves a list of system oxides Composition of the system can be expressed in terms of these oxides. Returns ------- array : strings An array of strings naming system components in terms of oxides """ return self.oxide_names_a
[docs] def get_phase_names(self): """Retrieves a list of system phases Names of phases known to the system Returns ------- array : strings An array of strings naming system phases """ return self.phase_names_a
[docs] def get_phase_inclusion_status(self): """Retrieves a dictionary of the inclusion status of each phase Returns ------- dict : dictionary A dictionary of boolean values indicating the inclusion status of each phase (key) known to the system """ dict = {} state_NSdict = self.melts.getPermissablePhasesState() for phase in self.phase_names_a: value = state_NSdict.valueForKey_(phase) #.boolValue removed in going from 0.2.7 -> 0.2.10 if value == 1: value = True else: value = False dict[phase] = value return dict
[docs] def set_phase_inclusion_status(self, status_d): """Sets the inclusion status of specified phases Parameters ---------- status_d : dictionary A dictionary of phase name keys and boolean values. True sets inclusion, and False prevents inclusion of a phase in the equilibrium assemblage. Note that the chemical affinity of the phase will still be calculated even if the inclusion level is set to False. """ state_NSdict = self.melts.getPermissablePhasesState() nsarray_cls = ObjCClass('NSMutableArray') nsarray = nsarray_cls.arrayWithCapacity_(self.nphases) for phase in self.phase_names_a: value = state_NSdict.valueForKey_(phase) #.boolValue removed in going from 0.2.7 -> 0.2.10 if phase in status_d: if status_d[phase] == True: value = 1 else: value = 0 nsarray.addObject_(ObjCClass('NSNumber').numberWithBool_(value)) self.melts.setPermissablePhasesState_(nsarray)
[docs] def set_bulk_composition(self, oxide_d={}): """Sets the bulk composition of the system This function first tests if the composition is feasible before setting the bulk composition of the system. You should check to make sure that the composition is feasible before proceeding. Parameters ---------- oxide_d : a python dictionary A dictionary of oxide names and values, e.g. {'SiO2':77.8, 'Al2O3':12.0, ..., 'H2O':3.74} Returns ------- boolean : True or False True if the composition is feasible, in which case the composition of the system is defined. False if the composition is infeasible, in which case the composition of the system is undefined. Notes ----- Feasibility call has yet to be implemented: (Objective-C method call:) -(BOOL)compositionIsFeasible:(NSArray \*)compositionInWtPercentOxides forSolution:(id <SolutionPhaseProtocol>)omniComponentPhase; """ wt = (ctypes.c_double*self.noxides)() ctypes.cast(wt, ctypes.POINTER(ctypes.c_double)) for i in range(0, self.noxides): wt[i] = 0.0 for key,value in oxide_d.items(): index = self.oxide_names_a.index(key) wt[index] = value self.melts.setComposition_(wt)
[docs] def set_temperature(self, t_c=800.0): """Sets the temperature of the system and reinitializes a calculation sequence Parameters ---------- t_c : float optional Float value to set the system temperature in degrees centigrade. Default is 800 °C. """ self.t = t_c + 273.15 self.melts.setTemperature_(self.t)
[docs] def set_entropy(self, s): """Sets the entropy of the system and reinitializes a calculation sequence Parameters ---------- s : float Float value to set the system entropy in J/K. """ self.entropy = s self.melts.setEntropy_(self.entropy)
[docs] def set_pressure(self, p_mpa=200.0): """Sets the pressure of the system and reinitializes a calculation sequence Parameters ---------- p_mpa : float optional Float value to set the system pressure in mega-Pascals. Default is 200 MPa. """ self.p = p_mpa*10.0 self.melts.setPressure_(self.p)
[docs] def set_volume(self, v): """Sets the volume of the system and reinitializes a calculation sequence Parameters ---------- v : float Float value to set the system volume in J/bar """ self.volume = v self.melts.setVolume_(self.volume)
[docs] def get_object_for_phase(self, phase_name): """Retrieve an object instance for the named phase. Parameters ---------- phase_name: string Name of phase Returns ------- object: instance of a phase class null if the phase is not in the stable assemblage """ phase_list = self.melts.dictionaryOfPhasesInSystem # dictionary of EquilibrateStatePhase instances if phase_name in phase_list: return phase_list[phase_name].phaseClassInstance else: return None
[docs] def get_properties_of_DEWFluid(self, property_name='species', T=1000, P=1000): """Retrieves a dictionary of properties of the specified type DEWFluid must exist in the equilibrium assemblage; otherwise an empty dictionary is returned. Parameters ---------- property_name: string, optional 'species' (default) returns a dictionary of species mole fractions in the equilibrated solution. 'mu' returns a dictionary of species chemical potentials in the equilibrated solution. 'activity' returns a dictionary of species activities in the equilibrated solution. T : float optional Temperature in degrees centigrade (default is 1000°C) P : float optional Pressure in bars (default is 1000 bars) Returns ------- dictionary: a python dictionary Keys are species names (strings). Values are species concentrations in mole fraction, or chemical potential, or thermodynamic activity, depending on the value of ``property_name``. """ phase_list = self.melts.dictionaryOfPhasesInSystem # dictionary of EquilibrateStatePhase instances if 'DEWFluid' in phase_list: DEW_object = phase_list['DEWFluid'].phaseClassInstance # instantiated class elements = phase_list['DEWFluid'].bulkCompositionInElements moles = DEW_object.convertElementsToMoles_(elements.pointerToDouble()).pointerToDouble() if property_name == 'species': return DEW_object.getSpeciesMoleFractionsForBulkComposition_aT_andP_(moles, T+273.15, P) elif property_name == 'mu': ns = DEW_object.numberOfSolutionSpecies() mu = DEW_object.chemicalPotentialsOfSpeciesFromMolesOfComponents_andT_andP_(moles, T+273.15, P).pointerToDouble() result = dict() for i in range(0,ns): name = DEW_object.nameOfSolutionSpeciesAtIndex_(i) result[name] = mu[i] return result elif property_name == 'activity': ns = DEW_object.numberOfSolutionSpecies() activity = DEW_object.activitiesOfSpeciesFromMolesOfComponents_andT_andP_(moles, T+273.15, P).pointerToDouble() result = dict() for i in range(0,ns): name = DEW_object.nameOfSolutionSpeciesAtIndex_(i) result[name] = activity[i] return result else: return dict() else: return dict()
[docs] def equilibrate_tp(self, T_a, P_a, initialize=False): """Determines the equilibrium phase assemblage at a temperature-pressure point or along a series of T-P points. The bulk composition of the system must first be set by calling the function: ``set_bulk_composition`` Parameters ---------- t_a : float or numpy array of floats Temperature in degrees centigrade. Either a scaler values or a numpy array of float values must be provided. p_a : float or numpy array of floats Pressure in mega-Pascals. Either a scaler values or a numpy array of float values must be provided. NOTE: If both ``t_a`` and ``p_a`` are arrays, then they must both be the same length. initialize : bool, optional True if this is a T-, P-point that starts a sequence of calculations. False if this is a continuation T-,P-pair or series of pairs. Returns ------- output_a : an array of tuples tuple = (status, T, P, xmlout): * **status** is a string indicating the status of the calculation: success/failiure, Reason for success/failure. * **T** is a float value corresponding to the temperature in degrees centigrade. * **P** is a float value corresponding to the pressure in mega-Pascals. * **xmlout** is an xml document tree of the type xml.etree.ElementTree. The xml tree contains information on the masses and abundances of all phases in the system. ``xmlout`` is utilized as input for a number of functions in this package that retrieve properties of the equilibrium assemblage. Notes ----- The ``xmlout`` document tree will be expanded to include thermodynamic properties of the phases and chemical affinities of phases not present in the equilibrium assemblage. """ T_a, P_a = core.fill_array(T_a, P_a) output_a = [] for ind,(T,P) in enumerate(zip(T_a,P_a)): if initialize: self.melts.setTemperature_(T+273.15) self.melts.setPressure_(P*10.0) nsarray_cls = ObjCClass('NSMutableArray') nsarraykeys = nsarray_cls.arrayWithCapacity_(5) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('ordinate')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('abscissa')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('imposeBuffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('buffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('bufferOffset')) nsarrayvalues = nsarray_cls.arrayWithCapacity_(5) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Temperature (°C)')) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Pressure (MPa)')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool_(0)) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('QFM')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool(0)) # Double initialization does not work (?) nsdict_cls = ObjCClass('NSDictionary') nsdict = nsdict_cls.dictionaryWithObjects_forKeys_(nsarrayvalues, nsarraykeys) self.melts.setCalculationOptions_(nsdict) else: self.melts.incrementTemperature_(T+273.15) self.melts.incrementPressure_(P*10.0) output_NSDictionary = self.melts.execute() xmlout = ET.fromstring(self.melts.equilibrateResultsAsXML()) output_a.append((output_NSDictionary.objectForKey_('status'), T, P, xmlout)) return output_a
[docs] def equilibrate_sp(self, S_a, P_a, initialize=False): """Determines the equilibrium phase assemblage at an entropy-pressure point or along a series of S-P points. The bulk composition of the system must first be set by calling the function: ``set_bulk_composition`` Parameters ---------- S_a : float or numpy array of floats Entropy in Joules per Kelvins. Either a scaler values or a numpy array of float values must be provided. P_a : float or numpy array of floats Pressure in mega-Pascals. Either a scaler values or a numpy array of float values must be provided. NOTE: If both ``t_a`` and ``p_a`` are arrays, then they must both be the same length. initialize : bool, optional True if this is a S-, P-point that starts a sequence of calculations. False if this is a continuation S-,P-pair or series of pairs. Returns ------- output_a : an array of tuples tuple = (status, T, P, xmlout): * **status** is a string indicating the status of the calculation: success/failiure, Reason for success/failure. * **T** is a float value corresponding to the temperature in degrees centigrade. * **P** is a float value corresponding to the pressure in mega-Pascals. * **xmlout** is an xml document tree of the type xml.etree.ElementTree. The xml tree contains information on the masses and abundances of all phases in the system. ``xmlout`` is utilized as input for a number of functions in this package that retrieve properties of the equilibrium assemblage. Notes ----- The ``xmlout`` document tree will be expanded to include thermodynamic properties of the phases and chemical affinities of phases not present in the equilibrium assemblage. """ S_a, P_a =core.fill_array(S_a, P_a) output_a = [] for ind,(S,P) in enumerate(zip(S_a,P_a)): if initialize: self.melts.setEntropy_(S) self.melts.setPressure_(P*10.0) nsarray_cls = ObjCClass('NSMutableArray') nsarraykeys = nsarray_cls.arrayWithCapacity_(5) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('ordinate')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('abscissa')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('imposeBuffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('buffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('bufferOffset')) nsarrayvalues = nsarray_cls.arrayWithCapacity_(5) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Entropy (J/K-kg)')) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Pressure (MPa)')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool_(0)) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('QFM')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool(0)) # Double initialization does not work (?) nsdict_cls = ObjCClass('NSDictionary') nsdict = nsdict_cls.dictionaryWithObjects_forKeys_(nsarrayvalues, nsarraykeys) self.melts.setCalculationOptions_(nsdict) else: self.melts.incrementEntropy_(S) self.melts.incrementPressure_(P*10.0) output_NSDictionary = self.melts.execute() xmlout = ET.fromstring(self.melts.equilibrateResultsAsXML()) T = locale.atof(xmlout.find(".//Temperature").text) output_a.append((output_NSDictionary.objectForKey_('status'), T, P, xmlout)) return output_a
[docs] def equilibrate_tv(self, T_a, V_a, initialize=False): """Determines the equilibrium phase assemblage at a temperature-volume point or along a series of T-V points. The bulk composition of the system must first be set by calling the function: ``set_bulk_composition`` Parameters ---------- T_a : float or numpy array of floats Temperature in degrees centigrade. Either a scaler values or a numpy array of float values must be provided. V_a : float or numpy array of floats Volume in Joules per bar. Either a scaler values or a numpy array of float values must be provided. NOTE: If both ``t_a`` and ``p_a`` are arrays, then they must both be the same length. initialize : bool, optional True if this is a T-, V-point that starts a sequence of calculations. False if this is a continuation T-,V-pair or series of pairs. Returns ------- output_a : an array of tuples tuple = (status, T, P, xmlout): * **status** is a string indicating the status of the calculation: success/failiure, Reason for success/failure. * **T** is a float value corresponding to the temperature in degrees centigrade. * **P** is a float value correcponding to the pressure in mega-Pascals. * **xmlout** is an xml document tree of the type xml.etree.ElementTree. The xml tree contains information on the masses and abundances of all phases in the system. ``xmlout`` is utilized as input for a number of functions in this package that retrieve properties of the equilibrium assemblage. Notes ----- The ``xmlout`` document tree will be expanded to include thermodynamic properties of the phases and chemical affinities of phases not present in the equilibrium assemblage. """ T_a, V_a = core.fill_array(T_a, V_a) output_a = [] for ind,(T,V) in enumerate(zip(T_a,V_a)): if initialize: self.melts.setTemperature_(T+273.15) self.melts.setVolume_(V) nsarray_cls = ObjCClass('NSMutableArray') nsarraykeys = nsarray_cls.arrayWithCapacity_(5) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('ordinate')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('abscissa')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('imposeBuffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('buffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('bufferOffset')) nsarrayvalues = nsarray_cls.arrayWithCapacity_(5) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Temperature (°C)')) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Volume (cc/kg)')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool_(0)) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('QFM')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool(0)) # Double initialization does not work (?) nsdict_cls = ObjCClass('NSDictionary') nsdict = nsdict_cls.dictionaryWithObjects_forKeys_(nsarrayvalues, nsarraykeys) self.melts.setCalculationOptions_(nsdict) else: self.melts.incrementTemperature_(T+273.15) self.melts.incrementVolume_(V) output_NSDictionary = self.melts.execute() xmlout = ET.fromstring(self.melts.equilibrateResultsAsXML()) output_a.append((output_NSDictionary.objectForKey_('status'), T, P, xmlout)) return output_a
[docs] def equilibrate_sv(self, S_a, V_a, initialize=False): """Determines the equilibrium phase assemblage at an entropy-volume point or along a series of S-V points. The bulk composition of the system must first be set by calling the function: ``set_bulk_composition`` Parameters ---------- S_a : float or numpy array of floats Entropy in Joules per Kelvins. Either a scaler values or a numpy array of float values must be provided. V_a : float or numpy array of floats Volume in Joules per bar. Either a scaler values or a numpy array of float values must be provided. NOTE: If both ``t_a`` and ``p_a`` are arrays, then they must both be the same length. initialize : bool, optional True if this is a S-, V-point that starts a sequence of calculations. False if this is a continuation S-,V-pair or series of pairs. Returns ------- output_a : an array of tuples tuple = (status, T, P, xmlout): * **status** is a string indicating the status of the calculation: success/failiure, Reason for success/failure. * **T** is a float value corresponding to the temperature in degrees centigrade. * **P** is a float value corresponding to the pressure in mega-Pascals. * **xmlout** is an xml document tree of the type xml.etree.ElementTree. The xml tree contains information on the masses and abundances of all phases in the system. ``xmlout`` is utilized as input for a number of functions in this package that retrieve properties of the equilibrium assemblage. Notes ----- The ``xmlout`` document tree will be expanded to include thermodynamic properties of the phases and chemical affinities of phases not present in the equilibrium assemblage. """ S_a, V_a = core.fill_array(S_a, V_a) output_a = [] for ind,(S,V) in enumerate(zip(S_a,V_a)): if initialize: self.melts.setEntropy_(S) self.melts.setVolume_(V) nsarray_cls = ObjCClass('NSMutableArray') nsarraykeys = nsarray_cls.arrayWithCapacity_(5) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('ordinate')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('abscissa')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('imposeBuffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('buffer')) nsarraykeys.addObject_(ObjCClass('NSString').stringWithString_('bufferOffset')) nsarrayvalues = nsarray_cls.arrayWithCapacity_(5) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Entropy (J/K-kg)')) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('Volume (cc/kg)')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool_(0)) nsarrayvalues.addObject_(ObjCClass('NSString').stringWithString_('QFM')) nsarrayvalues.addObject_(ObjCClass('NSNumber').numberWithBool(0)) # Double initialization does not work (?) nsdict_cls = ObjCClass('NSDictionary') nsdict = nsdict_cls.dictionaryWithObjects_forKeys_(nsarrayvalues, nsarraykeys) self.melts.setCalculationOptions_(nsdict) else: self.melts.incrementEntropy_(S) self.melts.incrementVolume_(V) output_NSDictionary = self.melts.execute() xmlout = ET.fromstring(self.melts.equilibrateResultsAsXML()) output_a.append((output_NSDictionary.objectForKey_('status'), T, P, xmlout)) return output_a
[docs] def get_list_of_phases_in_assemblage(self, root): """Returns a list of phases in the specified equilibrium assemblage. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx`` Returns ------- list : list A Python list of all phases in the equilibrium assemblage """ list = [] for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: list.append(phase) return list
[docs] def get_dictionary_of_default_fractionation_coefficients(self, fracLiq=False, fracSolid=True, fracFluid=True, fracCoeff=1.0): """Returns a dictionary of default coefficients for phase fractionation. These coefficients may be modified by the user. They are used when setting the extent to which a phase will fractionate from the equilibrium assemblage. Parameters ---------- fracLiq : bool, optional Flag to indicate if liquid phases should be fractionated from the system. Default is False. fracSolids : bool, optional Flag to indicate if solid phases should be fractionated from the system. Default is True. fracFluids : bool, optional Flag to indicate if fluid phases should be fractionated from the system. Default is True. fracCoeff : float, optional Fractionation coefficient, which gives the fractional extend to which the mass of the phase is extracted during phase fractionation. Default is 1.0 (100%). Returns ------- dict : dictionary A Python dictionary keyed on phases with values giving the extent (in fractional units) that a phase will fractionation mass. Notes ----- This dictionary is provided as input to the function ``fractionate_phases()``. The default configuration fractionates all solid/fluid phases and retains liquid. """ dict ={} for phase in self.phase_names_a: if phase == 'Liquid': if fracLiq: dict[phase] = fracCoeff else: dict[phase] = 0.0 elif phase == 'Water' or phase == 'Fluid': if fracFluid: dict[phase] = fracCoeff else: dict[phase] = 0.00 else: if fracSolid: dict[phase] = fracCoeff else: dict[phase] = 0.0 return dict
[docs] def get_mass_of_phase(self, root, phase_name='System'): """Returns the mass of a phase in the specified equilibrium assemblage. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx`` phase_name : string, optional The name of the phase whose abundance is requested, or the string 'System', which returns the combined mass of all phases in the system. Default value is 'System'. Returns ------- value : float The mass of the phase in the equilibrium assemblage specified by ``root``, in grams. If the specified phase is not in the equilibrium assemblage, a value of zero is retuned. """ if phase_name == 'System': value = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: value += float(root.find(".//System/Phase[@Type='" + phase + "']/Mass").text) else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: value = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Mass").text) else: value = 0.0 return value
[docs] def get_composition_of_phase(self, root, phase_name='System', mode='oxide_wt'): """Returns the composition of a phase in the specified equilibrium assemblage as a dictionary, with composition tabluated in the specified mode. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx`` phase_name : string, optional The name of the phase whose abundance is requested, or the string 'System', which returns the combined mass of all phases in the system. Default value is 'System'. mode : string, optional Controls the contents of the returned dictionary. * 'formula' returns a dictionary with the string 'formula' as key and value set to a string representation of the phase formula. For pure component phases, this is the standard phase formula. For solutions, this is the actual formula constructed by weighting the endmember components by their mole fractions. * 'oxide_wt' returns a dictionary of oxide string keys with values in wt%. This is a valid ``mode`` for all ``phase_name`` entries. * 'component' returns a dictionary of endmember component keys with values in mole fraction. The length of this dictionary will vary dependening on the number of components that describe the solution. Pure phases return an empty dictionary, as does ``phase_name`` set to 'System'. The default value of ``mode`` is 'oxide_wt'. Returns ------- dict : dictionary A dictionary describing the composition of ``phase_name`` according to the ``mode`` specified. The dictionary will be empty if ``phase_name`` is not present in the equilibrium assemblage. It will also be empty for certain cases described above under ``mode``. """ dict = {} if phase_name == 'System': if mode == 'oxide_wt': oxides = list(root.findall(".//Composition/Oxide")) for oxide in oxides: key = oxide.attrib['Type'] value = float(oxide.text) dict[key] = value else: phase = root.find(".//System/Phase[@Type='" + phase_name + "']") if phase != None: if mode == 'formula': dict['formula'] = phase.find("Formula").text elif mode == 'oxide_wt': oxides = list(phase.findall("Oxide")) for oxide in oxides: key = oxide.attrib['Type'] value = float(oxide.text) dict[key] = value elif mode == 'component': components = list(phase.findall("Component")) if not components: dict['formula'] = phase.find("Formula").text else: for component in components: key = component.attrib['Name'] value = float(component.text) dict[key] = value return dict
[docs] def fractionate_phases(self, root, frac_coeff): """Fractionates phases from the system. Partitions and maintains an internal dictionary of fractionates and automatically modifies system bulk composition to reflect fractionation. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx`` frac_coeff : dictionary A dictionary keyed on phase names with values that indicate the fraction of each phase that should fractionate. See get_dictionary_of_default_fractionation_coefficients(). Returns ------- dict : dictionary A dictionary keyed on phase names with values corresponding to a dictionary of phase properties. Keys are property names. """ dict = {} bc = {} bcFactor = locale.atof(root.find(".//Mass").text)/100.0 oxides = list(root.findall(".//Composition/Oxide")) for oxide in oxides: key = oxide.attrib['Type'] value = locale.atof(oxide.text) bc[key] = value phases = list(root.findall(".//System/Phase")) for phase in phases: phase_name = phase.attrib['Type'] dict[phase_name] = {} fraction = frac_coeff[phase_name] if fraction > 0.0: mass = locale.atof(phase.find("Mass").text) oxides = list(phase.findall("Oxide")) for oxide in oxides: key = oxide.attrib['Type'] value = locale.atof(oxide.text) dict[phase_name][key] = value*fraction*mass/100.0 bc[key] -= value*fraction*mass/100.0 if bc[key] < 0.0: bc[key] = 0.0 print ("Zeroed: " + key) self.set_bulk_composition(bc) return dict
[docs] def get_thermo_properties_of_phase_components(self, root, phase_name, mode='mu'): """Returns a dictionary of the specified component thermodynamic properties of the designated phase. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx``. phase_name : string The name of the phase whose abundance is requested. mode : string, optional Controls the contents of the returned dictionary. * 'mu' returns a dictionary of endmember component keys with values of chemical potential. The length of this dictionary will vary depending on the number of components that describe the solution. Pure phases return a dictionary of unit length, with the phase name as key and their specific Gibbs free energy as value (J/g). * 'excess' returns a dictionary of endmember component keys with values of excess chemical potential. The length of this dictionary will vary depending on the number of components that describe the solution. Pure phases return a dictionary of unit length, with the phase name as key and zero as value. * 'activity' returns a dictionary of endmember component keys with values of component activity. The length of this dictionary will vary depending on the number of components that describe the solution. Pure phases return a dictionary of unit length, with the phase name as key and unity as value. Returns ------- dict : dictionary A dictionary describing the thermodynamic properties of components in ``phase_name`` according to the ``mode`` specified. The dictionary will be empty if ``phase_name`` is not present in the equilibrium assemblage. """ dict = {} phase = root.find(".//System/Phase[@Type='" + phase_name + "']") if phase != None: if mode == 'mu': mus = list(phase.findall("ChemicalPotential")) if not mus: value = locale.atof(phase.find("GibbsFreeEnergy").text) mass = float(phase.find("Mass").text) dict[phase_name] = value/mass else: for mu in mus: key = mu.attrib['Name'] value = locale.atof(mu.text) dict[key] = value elif mode == 'excess': mus = list(phase.findall("ExcessChemicalPotential")) if not mus: dict[phase_name] = 0.0 else: for mu in mus: key = mu.attrib['Name'] value = locale.atof(mu.text) dict[key] = value elif mode == 'activity': mus = list(phase.findall("ExcessChemicalPotential")) if not mus: dict[phase_name] = 1.0 else: t = locale.atof(root.find(".//Temperature").text) + 273.15 for mu in mus: key = mu.attrib['Name'] value = locale.atof(mu.text) dict[key] = np.exp(value/8.3143/t) return dict
[docs] def get_list_of_properties(self): """Returns a list of properties reported for each phase in an equilibrium assemblage. Returns ------- list : list A Python list of all properties of phases in an equilibrium assemblage """ list = ['Mass','GibbsFreeEnergy','Enthalpy','Entropy','HeatCapacity','DcpDt', 'Volume','DvDt','DvDp','D2vDt2','D2vDtDp','D2vDp2', \ 'Density','Alpha','Beta','K',"K'",'Gamma'] return list
[docs] def get_units_of_property(self, prop='Mass'): """Returns the units of a specified property. Returns ------- string : string The units of the specified property. Returns 'none' if property is invalid. """ dict = {'Mass':'g','GibbsFreeEnergy':'J','Enthalpy':'J','Entropy':'J/K','HeatCapacity':'J/K','DcpDt':'J/K^2', \ 'Volume':'J/bar','DvDt':'J/bar-K','DvDp':'J/bar^2','D2vDt2':'J/bar-K^2','D2vDtDp':'J/bar^2-K','D2vDp2':'J/bar^3', \ 'Density':'g/cm^3','Alpha':'1/K','Beta':'1/bar','K':'GPa',"K'":'none','Gamma':'none'} return dict.get(prop)
[docs] def get_property_of_phase(self, root, phase_name='System', property_name='Mass'): """Returns the specified property of a phase in the specified equilibrium assemblage. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx`` phase_name : string, optional The name of the phase whose property is requested, or the string 'System', which returns the combined property of all phases in the system. Default value is 'System'. property_name : string, optional The name of the property to be returned. Default value is 'Mass'. Returns ------- value : float The property of the phase in the equilibrium assemblage specified by ``root``, in standard units. If the specified phase is not in the equilibrium assemblage, a value of zero is returned. If the property is not in the standard list, a value of zero is returned. """ standard = ['Mass','GibbsFreeEnergy','Enthalpy','Entropy','HeatCapacity','DcpDt', 'Volume','DvDt','DvDp','D2vDt2','D2vDtDp','D2vDp2'] derived = ['Density','Alpha','Beta','K',"K'",'Gamma'] if property_name in standard: if phase_name == 'System': value = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/" + property_name) != None: value += locale.atof(root.find(".//System/Phase[@Type='" + phase + "']/" + property_name).text) else: if root.find(".//System/Phase[@Type='" + phase_name + "']/" + property_name) != None: value = locale.atof(root.find(".//System/Phase[@Type='" + phase_name + "']/" + property_name).text) else: value = 0.0 elif property_name in derived: if property_name == 'Density': if phase_name == 'System': volume = 0.0 mass = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: volume += float(root.find(".//System/Phase[@Type='" + phase + "']/Volume").text) mass += float(root.find(".//System/Phase[@Type='" + phase + "']/Mass").text) value = mass/volume/10.0 # g/cc else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: volume = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Volume").text) mass = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Mass").text) value = mass/volume/10.0 # g/cc else: value = 0.0 elif property_name == 'Alpha': if phase_name == 'System': volume = 0.0 dvdt = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: volume += float(root.find(".//System/Phase[@Type='" + phase + "']/Volume").text) dvdt += float(root.find(".//System/Phase[@Type='" + phase + "']/DvDt").text) value = dvdt/volume else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: volume = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Volume").text) dvdt = float(root.find(".//System/Phase[@Type='" + phase_name + "']/DvDt").text) value = dvdt/volume else: value = 0.0 elif property_name == 'Beta': if phase_name == 'System': volume = 0.0 dvdp = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: volume += float(root.find(".//System/Phase[@Type='" + phase + "']/Volume").text) dvdp += float(root.find(".//System/Phase[@Type='" + phase + "']/DvDp").text) value = -dvdp/volume else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: volume = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Volume").text) dvdp = float(root.find(".//System/Phase[@Type='" + phase_name + "']/DvDp").text) value = -dvdp/volume else: value = 0.0 elif property_name == 'K': if phase_name == 'System': volume = 0.0 dvdp = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: volume += float(root.find(".//System/Phase[@Type='" + phase + "']/Volume").text) dvdp += float(root.find(".//System/Phase[@Type='" + phase + "']/DvDp").text) value = -volume/dvdp/10000.0 else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: volume = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Volume").text) dvdp = float(root.find(".//System/Phase[@Type='" + phase_name + "']/DvDp").text) value = -volume/dvdp/10000.0 else: value = 0.0 elif property_name == "K'": if phase_name == 'System': volume = 0.0 dvdp = 0.0 d2vdp2 = 0.0 for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: volume += float(root.find(".//System/Phase[@Type='" + phase + "']/Volume").text) dvdp += float(root.find(".//System/Phase[@Type='" + phase + "']/DvDp").text) d2vdp2 += float(root.find(".//System/Phase[@Type='" + phase + "']/D2vDp2").text) value = (volume*d2vdp2/dvdp/dvdp-1.0) else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: volume = float(root.find(".//System/Phase[@Type='" + phase_name + "']/Volume").text) dvdp = float(root.find(".//System/Phase[@Type='" + phase_name + "']/DvDp").text) d2vdp2 = float(root.find(".//System/Phase[@Type='" + phase_name + "']/D2vDp2").text) value = (volume*d2vdp2/dvdp/dvdp-1.0) else: value = 0.0 elif property_name == 'Gamma': dvdp = 0.0 dvdt = 0.0 cp = 0.0 t = locale.atof(root.find(".//Temperature").text)+273.15 if phase_name == 'System': for phase in self.phase_names_a: if root.find(".//System/Phase[@Type='" + phase + "']/Mass") != None: dvdp += float(root.find(".//System/Phase[@Type='" + phase + "']/DvDp").text) dvdt += float(root.find(".//System/Phase[@Type='" + phase + "']/DvDt").text) cp += float(root.find(".//System/Phase[@Type='" + phase + "']/HeatCapacity").text) value = -dvdt/(cp*dvdp + t*dvdt*dvdt) else: if root.find(".//System/Phase[@Type='" + phase_name + "']/Mass") != None: dvdp = float(root.find(".//System/Phase[@Type='" + phase_name + "']/DvDp").text) dvdt = float(root.find(".//System/Phase[@Type='" + phase_name + "']/DvDt").text) cp = float(root.find(".//System/Phase[@Type='" + phase_name + "']/HeatCapacity").text) value = -dvdt/(cp*dvdp + t*dvdt*dvdt) else: value = 0.0 else: value = 0.0 return value
[docs] def get_dictionary_of_affinities(self, root, sort=True): """Returns an ordered dictionary of tuples of chemical affinity and phase formulae for undersaturated phases in the system. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_xx`` sort : boolean A flag when set to sort the dictionary in order of ascending affinities Returns ------- dict : OrderedDict A Python ordered dictionary. Dictionary keys are strings naming the phases not in the equilibrium assemblage but known to the system. These phases are by definition undersaturated. Dictionary values are tuples consisting of a float value and a string: (affinity, formula). For a solution phase, the formula is the composition closest to equilibrium with the reported phase assemblage. Dictionary ordering corresponds to array order in get_phase_names(), unless ``sorted`` is set to True; then entries are ordered by ascending chemical affinity. """ dict = OrderedDict() for phase in self.phase_names_a: if root.find(".//Potential/Phase[@Type='" + phase + "']") != None: affinity = locale.atof(root.find(".//Potential/Phase[@Type='" + phase + "']/Affinity").text) formulae = root.find(".//Potential/Phase[@Type='" + phase + "']/Formula").text if affinity == 0.0: affinity = 999999.0 dict[phase] = (affinity, formulae) if sort == True: return OrderedDict(sorted(dict.items(), key=lambda t:t[1])) else: return dict
[docs] def output_summary(self, root, printT=True, printP=True, printMass=False, printSysWt=False, printSysM=False, printPhs=True, printPhsWt=False, printPhsM=False): """Prints information about the specified equilibrium phase assemblage. Parameters ---------- root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_tp`` printT : bool, optional, default=True Print the system temperature in degrees centigrade printP : bool, optional, default=True Print the system pressure in mega-Pascals printMass ; bool, optional, default=False Print the mass of the system in grams printSysWt : bool, optional, default=False Print the composition of the system in wt% oxides printSysM : bool, optional, default=False Print the composition of the system in moles of elements printPhs : bool, optional, default=True Print the phases present in the system, their masses (in grams) and their chemical formulas printPhsWt : bool, optional, default=False Print the composition of each phase in wt% oxides (most often used in conjunction with printPhs=True) printPhsM : bool, optional, default=False Print the composition of each phase in moles of endmember components (most often used in conjunction with printPhs=True) """ if printT: print ("{0:<10s} {1:>8.2f}".format("T (°C)", locale.atof(root.find(".//Temperature").text))) if printP: print ("{0:<10s} {1:>8.2f}".format("P (MPa)", locale.atof(root.find(".//Pressure").text)*1000.0)) if printMass: print ("{0:<10s} {1:>8.3f}".format("Mass (g)", float(root.find(".//Mass").text))) if printSysM: print ("Bulk composition in elemental abundances (moles):") bcElements = list(root.findall(".//Composition/Element")) for element in bcElements: print (" {0:>2s} {1:>8.5f}".format(element.attrib['Type'], float(element.text))) if printSysWt: print ("Bulk composition in oxide abundances (wt %):") bcOxides = list(root.findall(".//Composition/Oxide")) for oxide in bcOxides: print (" {0:>5s} {1:>8.4f}".format(oxide.attrib['Type'], float(oxide.text))) phases = list(root.findall(".//System/Phase")) for phase in phases: if printPhs: formula = phase.find("Formula").text.split(" ") if len(formula) == 1: print ("{0:<15.15s} {1:>8.4f} (g) {2:<60s}".format(phase.attrib['Type'], float(phase.find("Mass").text), formula[0])) else: # this formula consists of oxide, value pairs (liquid phase) n_oxides = int(len(formula)/2) n_oxides_1 = int(n_oxides/2) n_oxides_2 = int(n_oxides - n_oxides_1) if n_oxides <= 8: n_oxides_1 = n_oxides n_oxides_2 = 0 line_1 = "" for i in range(n_oxides_1): line_1 = line_1 + ' ' + formula[2*i] + ' ' + formula[2*i+1] print ("{0:<15.15s} {1:>8.4f} (g) {2:<s}".format(phase.attrib['Type'], float(phase.find("Mass").text), line_1)) line_2 = "" for i in range(n_oxides_2): line_2 = line_2 + ' ' + formula[n_oxides_1*2+2*i] + ' ' + formula[n_oxides_1*2+2*i+1] print ("{0:<33.33s}{1:<s}".format(" ", line_2)) if printPhsWt: oxides = list(phase.findall("Oxide")) for oxide in oxides: value = float(oxide.text) if (value != 0.0): print (" {0:<5s} {1:>8.4f}".format(oxide.attrib['Type'], float(oxide.text))) if printPhsM: components = list(phase.findall("Component")) for component in components: value = float(component.text) if (value != 0.0): print (" {0:<15s} {1:>9.6f}".format(component.attrib['Name'], float(component.text)))
# ================================================================= # Block of Excel workbook functions for output # =================================================================
[docs] def start_excel_workbook_with_sheet_name(self, sheetName="Summary"): """Create an Excel workbook with one named sheet. Parameters ---------- sheetName : string Sheet name in the new empty Excel workbook Returns ------- wb : type Workbook A pointer to an Excel workbook Notes ----- Part of the Excel workbook functions subpackage """ wb = Workbook() ws = wb.active ws.title = sheetName self.row = 0 return wb
[docs] def add_sheet_to_workbook_named(self, wb, sheetName): """Creates a new sheet in an existing Excel workbook. Parameters ---------- wb : type Workbook A pointer to an existing Excel workbook sheetName : string New sheet name in the specified Excel workbook Returns ------- ws : type Worksheet A pointer to an Excel worksheet in the specified workbook, ``wb`` Notes ----- Part of the Excel workbook functions subpackage """ ws = wb.create_sheet(title=sheetName) return ws
[docs] def write_to_cell_in_sheet(self, ws, col, row, value, format='general'): """Writes information into the specified row, col on the specified worksheet. Parameters ---------- ws : type Worksheet A pointer to a previously created Excel worksheet (see add_sheet_to_workbook_named) col : int Column number to write entry to. Numbering starts at column one. row : int Row number to write entry to. Numbering starts at row one. value : float Value to place at entry format : string, optional Format to use for value * 'general' is the default; no formatting applied * 'number' formats as '0.00' * 'scientific' formats as '0.00E+00' Notes ----- Part of the Excel workbook functions subpackage """ if format == 'number': ws.cell(column=col, row=row, value=float(value)).number_format = '0.00' elif format == 'scientific': ws.cell(column=col, row=row, value=float(value)).number_format = '0.00E+00' else: ws.cell(column=col, row=row, value=value)
[docs] def write_excel_workbook(self, wb, fileName="junk.xlsx"): """Writes the specified Excel workbook to disk. Parameters ---------- wb : type Workbook A pointer to an existing Excel workbook fileName : string, optional Name of the file that will contain the specified Excel notebook. Default file name is 'junk.xlsx'. Notes ----- Part of the Excel workbook functions subpackage """ wb.save(filename = fileName)
[docs] def update_excel_workbook(self, wb, root): """Writes the specified equilibrium system state to the specified Excel workbook. Parameters ---------- wb : type Workbook A pointer to an existing Excel workbook root : type xml.etree.ElementTree An xml document tree returned by the function ``equilibrate_tp`` Notes ----- The workbook is structured with a Summary worksheet and one worksheet for each equilibrium phase. The evolution of the system is recorded in successive rows. Row integrity is maintained across all sheets. This function may be called repeatedly using different ``root`` objects. Part of the Excel workbook functions subpackage. """ t = locale.atof(root.find(".//Temperature").text) p = locale.atof(root.find(".//Pressure").text)*1000.0 bcElements = list(root.findall(".//Composition/Element")) bcOxides = list(root.findall(".//Composition/Oxide")) wsSummary = wb["Summary"] # .get_sheet_by_name("Summary") upgraded if (self.row == 0): col = 1 self.row = 1 self.write_to_cell_in_sheet(wsSummary, col, self.row, "T °C") col += 1 self.write_to_cell_in_sheet(wsSummary, col, self.row, "P MPa") col += 1 self.write_to_cell_in_sheet(wsSummary, col, self.row, "Mass g") col += 1 for element in bcElements: self.write_to_cell_in_sheet(wsSummary, col, self.row, element.attrib['Type']) col += 1 for oxide in bcOxides: self.write_to_cell_in_sheet(wsSummary, col, self.row, oxide.attrib['Type']) col += 1 self.row += 1 col = 1 self.write_to_cell_in_sheet(wsSummary, col, self.row, t, format='number') col += 1 self.write_to_cell_in_sheet(wsSummary, col, self.row, p, format='number') col += 1 self.write_to_cell_in_sheet(wsSummary, col, self.row, root.find(".//Mass").text, format='number') col += 1 for element in bcElements: self.write_to_cell_in_sheet(wsSummary, col, self.row, element.text, format='scientific') col += 1 for oxide in bcOxides: self.write_to_cell_in_sheet(wsSummary, col, self.row, oxide.text, format='number') col += 1 phases = list(root.findall(".//System/Phase")) for phase in phases: phaseType = phase.attrib['Type'] oxides = list(phase.findall("Oxide")) components = list(phase.findall("Component")) try: wsPhase = wb[phaseType] # .get_sheet_by_name(phaseType) upgraded except KeyError: wsPhase = wb.create_sheet(phaseType) col = 1 self.write_to_cell_in_sheet(wsPhase, col, 1, "T °C") col += 1 self.write_to_cell_in_sheet(wsPhase, col, 1, "P MPa") col += 1 self.write_to_cell_in_sheet(wsPhase, col, 1, "Mass g") col += 1 self.write_to_cell_in_sheet(wsPhase, col, 1, "Formula") col += 1 for oxide in oxides: self.write_to_cell_in_sheet(wsPhase, col, 1, oxide.attrib['Type']) col += 1 for component in components: self.write_to_cell_in_sheet(wsPhase, col, 1, component.attrib['Name']) col += 1 col = 1 self.write_to_cell_in_sheet(wsPhase, col, self.row, t, format='number') col += 1 self.write_to_cell_in_sheet(wsPhase, col, self.row, p, format='number') col += 1 self.write_to_cell_in_sheet(wsPhase, col, self.row, phase.find("Mass").text, format='number') col += 1 self.write_to_cell_in_sheet(wsPhase, col, self.row, phase.find("Formula").text) col += 1 for oxide in oxides: self.write_to_cell_in_sheet(wsPhase, col, self.row, oxide.text, format='number') col += 1 for component in components: self.write_to_cell_in_sheet(wsPhase, col, self.row, component.text, format='scientific') col += 1