ENKI

Source code for phases

"""
This module provides Python wrappers to individual and collections of phases.

"""

# from thermoengine.core import *
from thermoengine import core
from thermoengine import chem

DATADIR = 'data/phases'

import numpy as np
import pandas as pd
from os import path
from collections import OrderedDict
from abc import ABCMeta, abstractmethod
from importlib import import_module, invalidate_caches
import sys

# specialized numerical imports
from scipy import optimize as optim

# __all__ = ['Rxn','Assemblage','PurePhase','SolutionPhase','get_phaselist']
__all__ = ['Rxn','Assemblage','Phase','PurePhase','SolutionPhase','get_phaselist']

PURE_PHASE_FILE = 'PurePhaseList.csv'
SOLUTION_PHASE_FILE = 'SolutionPhaseList.csv'

# inputs
#
# * P=1, T=1, mol=1
# * P=N_PTX, T=1, mol=1
# * P=1, T=N_PTX, mol=1
# * P=1, T=1, mol = N_PTX
# * P=N_PTX, T=N_PTX, mol=N_PTX
#
# qty
#
# * N_endmem (deriv_type)
# * N_rxn
# * N_phase
#
# outputs
#
# * N_PTX,
# * N_PTX, N_endmem
# * N_PTX, N_endmem, N_endmem
#
# * N_phase,
# * N_rxn, N_phase
# * N_rxn, N_phase, N_PTX
#
# * reduce dimensions if N_PTX=1
#===================================================
def get_phase_info():
    """
    Get acceptable table of basic phase info
        (e.g., name, abbrev, formula, members, etc.).

    Returns
    -------
    phase_table : dict
        Dictionary containing phase information.
            - 'purephases' : pandas dataframe
            - 'purefilenm' : str defining source file for pure phase info

    """

    phase_info = {}
    phase_info['pure'] = _read_phase_info(PURE_PHASE_FILE)
    phase_info['solution'] = _read_phase_info(SOLUTION_PHASE_FILE)

    info_files = {}
    info_files['pure'] = PURE_PHASE_FILE
    info_files['solution'] = SOLUTION_PHASE_FILE
    return phase_info, info_files

def _read_phase_info(filenm):
    """
    Read phase info tables

    Internal method to read phase abbreviations and names from csv files.

    """

    parentpath = path.dirname(__file__)
    pathname = path.join(parentpath,DATADIR,filenm)
    try:
       phases_info_df = pd.read_csv(pathname)
    except:
        assert False,'The '+pathname+' file cannot be found. '\
            'It is needed to define the standard phase abbreviations.'

    return phases_info_df
#===================================================
class FixedRxnSet:
    def __init__(self, phase_symbols, endmember_ids, rxn_coefs,
                 phases_dict, T, P, mols,):

        self._init_rxn_set(phase_symbols, endmember_ids, rxn_coefs,
                           phases_dict)
        self._init_exp_cond(T, P, mols)

    def _init_rxn_set(self, phase_symbols, endmember_ids, rxn_coefs,
                      phases_dict):

        assert np.all(
            [sym in phases_dict.keys() for sym in phase_symbols ]), (
                'phase symbols must appear in phase_obj_dict'
            )

        self._phase_symbols = phase_symbols
        self._endmember_ids = endmember_ids
        self._rxn_coefs = rxn_coefs
        self._phases_dict = phases_dict
        self._rxn_num = rxn_coefs.shape[0]
        self._endmember_num = rxn_coefs.shape[1]

    def _init_exp_cond(self, T, P, mols):
        phase_symbols = self._phase_symbols

        # validate mols
        self._validate_mol_input(mols)

        assert np.all(
            [sym in phase_symbols for sym in mols.keys()]),(
                'molar composition must be defines for every phase'
            )

        self._T = T
        self._P = P
        self._mols = mols

    def _validate_mol_input(self, mols):
        mols = {} if mols is None else mols

        phases = self.phases
        phase_symbols = self.phase_symbols

        for symbol in phase_symbols:
            if symbol not in mols:
                mols[symbol] = None

        assert np.all(np.array(
            [symbol in phase_symbols for symbol in mols.keys()])), (
                'Invalid phase symbol(s) used to define phase '
                'composition in mols. Only use valid phase symbols '
                'found in rxn.phase_symbols.'
            )

        return mols

    def affinity(self):
        rxn_coefs = self._rxn_coefs

        phase_chem_potentials = self.phase_chem_potentials()
        chem_potentials = self.endmem_chem_potentials(
            phase_chem_potentials)

        affinities = np.dot(rxn_coefs, chem_potentials)
        return affinities

    def endmem_chem_potentials(self, phase_chem_potentials):
        endmember_num = self._endmember_num
        phase_symbols = self._phase_symbols
        endmember_ids = self._endmember_ids

        chem_potentials = zeros(endmember_num)
        for ind, (phs_sym, endmem_id) in enumerate(
            zip(phase_symbols, endmember_ids)):
            chem_potentials[ind] = phase_chem_potentials[phs_sym][endmem_id]

        return chem_potentials

    def phase_chem_potentials(self):
        phases_dict = self._phases_dict

        T = self._T
        P = self._P
        mols = self._mols

        phase_chem_potentials = {}

        for phs_sym in phases_dict:
            phs = phases_dict[phs_sym]
            phase_chem_potentials[phs_sym] = phs.chem_pot(T, P, mols=mols)

        return phase_chem_potentials




#===================================================
[docs]class Rxn: """ Class that defines identity/properties of a specific phase reaction. Reactions occur between phases (either pure or solution) and are defined in terms of the participating endmembers, indicating which atoms are exchanged between phases during the reaction. Parameters ---------- phase_objs : array of Phase Objects Defines which phases participate in the reaction. endmember_ids : int array Indicates the endmember of each phase that participates in the reaction. This array must have the same order as the phase array (phase_objs). rxn_coefs : double array Defines the stoichiometric rxn coefficient, where negative values are reactants and positive values are products. The reaction must be balanced (obeying mass conservation). This array must have the same order as the phase array (phase_objs). Attributes ---------- endmember_ids endmember_names phase_num phase_symbols phases product_phases reactant_phases rxn_coefs Notes ----- * The phases themselves may be pure or have realistic intermediate compositions (if they are solution phases). * The reaction is defined in terms of the exchange of endmembers between the participating phases. * Reaction coefficients correspond to a balanced stoichiometric reaction. """ def __init__(self, phase_objs, endmember_ids, rxn_coefs, coefs_per_atom=False): self._validate_inputs(phase_objs, endmember_ids, rxn_coefs) self._init_rxn(phase_objs, endmember_ids, rxn_coefs, coefs_per_atom) # TK # self._validate_rxn_balance pass def _validate_inputs(self, phase_objs, endmember_ids, rxn_coefs): Nphase = len(phase_objs) rxn_coefs = np.array(rxn_coefs) assert len(endmember_ids) == Nphase, ( 'endmember_ids must provide endmember index for ' 'every reaction phase. ' ) assert len(rxn_coefs) == Nphase, ( 'rxn_coefs must provide stoichiometric ' 'coefficients for every reaction phase.' ) assert np.any(rxn_coefs<0), ( 'rxn_coefs must indicate reactants with negative ' 'sign and products with positive sign. Currently, ' 'none of the stoichiometric coefficients are negative.' ) assert np.any(rxn_coefs>0), ( 'rxn_coefs must indicate reactants with negative ' 'sign and products with positive sign. Currently, ' 'none of the stoichiometric coefficients are positive.' ) def _init_rxn(self, phase_objs, endmember_ids, rxn_coefs, coefs_per_atom): phase_objs, endmember_ids, rxn_coefs = ( self._trim_absent_rxn_phases( phase_objs, endmember_ids, rxn_coefs)) phase_objs = np.array(phase_objs) endmember_ids = np.array(endmember_ids) Nphase = len(phase_objs) rxn_coefs = np.array(rxn_coefs) all_atom_nums = [] for phs, endmember_id in zip(phase_objs, endmember_ids): endmember_id = int(endmember_id) iatom_num = phs.props['atom_num'][endmember_id] all_atom_nums.append(iatom_num) all_atom_nums = np.array(all_atom_nums) if coefs_per_atom: rxn_coefs = rxn_coefs/all_atom_nums phase_symbols = [] endmember_names = [] for phase_obj, iendmember_id in zip(phase_objs, endmember_ids): phase_symbols.append( phase_obj.props['abbrev']) endmember_names.append( phase_obj.endmember_names[int(iendmember_id)]) phase_symbols = np.array(phase_symbols) # indsort = np.argsort(phase_symbols) # self._set_phase_assemblage( # phase_objs[indsort], phase_symbols[indsort], # obj_is_classnm=obj_is_classnm) self._phase_num = Nphase self._phase_symbols = phase_symbols self._phases = phase_objs self._endmember_ids = endmember_ids self._endmember_names = endmember_names self._rxn_coefs = rxn_coefs def _trim_absent_rxn_phases(self, phase_objs, endmember_ids, rxn_coefs): remove_set, = np.where(rxn_coefs == 0) # re-create lists for items not in remove set rxn_coefs = [v for i, v in enumerate(rxn_coefs) if i not in remove_set] phase_objs = [v for i, v in enumerate(phase_objs) if i not in remove_set] endmember_ids = [v for i, v in enumerate(endmember_ids) if i not in remove_set] return phase_objs, endmember_ids, rxn_coefs @property def phase_num(self): """ Number of phases Returns ------- Number of phases (int) """ return self._phase_num @property def phases(self): """ Phase objects used in the reaction Returns ------- Array of phase objects used in the reaction """ return self._phases @property def reactant_phases(self): """ Reactant phases Returns ------- Array of reactant phase objects """ return self.phases[self.rxn_coefs<0] @property def product_phases(self): """ Product phases Returns ------- Array of product phase objects """ return self.phases[self.rxn_coefs>0] @property def phase_symbols(self): """ Phase symbols Returns ------- Array of phase symbols used in the reaction (str) """ return self._phase_symbols @property def endmember_ids(self): """ ID number of each endmember in phase Returns ------- Array of ids, [int,...] """ return self._endmember_ids @property def endmember_names(self): """ Name of each endmember Returns ------- List of endmember names for this solution phase, [str,...] """ return self._endmember_names @property def rxn_coefs(self): """ Reaction coefficients Returns ------- Array of reaction coefficients (double) """ return self._rxn_coefs def _validate_mol_input(self, mols): mols = {} if mols is None else mols phases = self.phases phase_symbols = self.phase_symbols for symbol in phase_symbols: if symbol not in mols: mols[symbol] = None assert np.all(np.array( [symbol in phase_symbols for symbol in mols.keys()])), ( 'Invalid phase symbol(s) used to define phase ' 'composition in mols. Only use valid phase symbols ' 'found in rxn.phase_symbols.' ) return mols
[docs] def affinity(self, T, P, mols=None): """ Calculate reaction affinity Parameters ---------- T : array-like Temperature in Kelvin P : array-like Pressure in bars mol : dict of arrays, optional Composition of each phase in terms of mols of endmembers (unneeded for pure phases) Returns ------- value : array-like Reaction affinity in J """ affinity = -self.chem_potential(T, P, mols=mols) return affinity
[docs] def chem_potential(self, T, P, mols=None): """ Calculate net chemical potential change of the reaction Parameters ---------- T : array-like Temperature in Kelvin P : array-like Pressure in bars mol : dict of arrays, optional Composition of each phase in terms of mols of endmembers (unneeded for pure phases) Returns ------- value : array-like Chemical potential in J for the net change of the reaction """ chem_potential = self._calc_net_rxn_values( 'chem_potential', T, P, mols=mols, use_endmember=True) return chem_potential
[docs] def volume(self, T, P, mols=None, peratom=False): """ Calculate net volume change of the reaction Parameters ---------- T : array-like Temperature in Kelvin P : array-like Pressure in bars mol : dict of arrays, optional Composition of each phase in terms of mols of endmembers (unneeded for pure phases) Returns ------- value : array-like Volume in J/bar for the net change of the reaction """ rxn_volume = self._calc_net_rxn_values('volume', T, P, mols=mols, peratom=peratom) return rxn_volume
[docs] def entropy(self, T, P, mols=None, peratom=False): """ Calculate net entropy change of the reaction Parameters ---------- T : array-like Temperature in Kelvin P : array-like Pressure in bars mol : dict of arrays, optional Composition of each phase in terms of mols of endmembers (unneeded for pure phases). Returns ------- value : array-like Entropy in J/K for the net change of the reaction. """ rxn_entropy = self._calc_net_rxn_values('entropy', T, P, mols=mols, peratom=peratom) return rxn_entropy
def boundary(self, T=None, P=None, mols=None, init_guess=None): Tstd = 300.0 Pstd = 1.0 assert (T is not None) or (P is not None), \ 'Both T and P are None; Must define either temperature or pressure.' if T is not None: if init_guess is None: init_guess = Pstd Afun = lambda P, T=T: self.affinity(T, P, mols=mols) dAfun = lambda P, T=T: -self.volume(T, P, mols=mols, peratom=True) else: if init_guess is None: init_guess = Tstd Afun = lambda T, P=P: self.affinity(T, P, mols=mols) dAfun = lambda T, P=P: +self.entropy(T, P, mols=mols, peratom=True) value_bound = optim.newton(Afun, init_guess, fprime=dAfun) return value_bound def clapeyron_slope(self, T, P, mols=None, peratom=False): dV_rxn = self.volume(T, P, mols=mols, peratom=peratom) dS_rxn = self.entropy(T, P, mols=mols, peratom=peratom) dTdP = dV_rxn / dS_rxn return dTdP def simultaneous_rxn_cond(self, other_rxn, Pinit=1.0, TOL=1e-5): P = Pinit while True: Tbnd1 = self.boundary(P=P) Tbnd2 = other_rxn.boundary(P=P) dTdP1 = self.clapeyron_slope(Tbnd1,P) dTdP2 = other_rxn.clapeyron_slope(Tbnd2,Pinit) if np.abs(np.log(Tbnd1/Tbnd2)) < TOL: T = 0.5*(Tbnd1+Tbnd2) break dP = - (Tbnd1 - Tbnd2)/(dTdP1-dTdP2) P += dP T = float(T) P = float(P) return T, P def trace_boundary(self, Tlims=None, Plims=None, init_guess=None, Nsamp=30): assert (Tlims is not None) or (Plims is not None), \ 'Both T and P are None; Must define either temperature or pressure.' if Tlims is not None: xvals = np.linspace(Tlims[0], Tlims[1], Nsamp) bound_fun = lambda T, init_guess: self.boundary( T=T, init_guess=init_guess) dAdx_fun = lambda x, y: +self.entropy(x, y) dAdy_fun = lambda x, y: -self.volume(x, y) else: xvals = np.linspace(Plims[0], Plims[1], Nsamp) bound_fun = lambda P, init_guess: self.boundary( P=P, init_guess=init_guess) dAdx_fun = lambda x, y: -self.volume(y, x) dAdy_fun = lambda x, y: +self.entropy(y, x) dx = xvals[1]-xvals[0] yvals = np.zeros(xvals.shape) for ind, x in enumerate(xvals): y = bound_fun(x, init_guess) dy_guess = dx*dAdx_fun(x=x, y=y)/dAdy_fun(x=x, y=y) yvals[ind] = y init_guess = y + dy_guess if Tlims is not None: T_bnds, P_bnds = xvals, yvals else: T_bnds, P_bnds = yvals, xvals return T_bnds, P_bnds def competing_rxn(self, T_bounds, P_bounds, other_rxn, TOL=1e-6): competing = False A_rxn_other = other_rxn.affinity(T_bounds, P_bounds) is_competing = np.tile(False, len(T_bounds)) if np.all(other_rxn.reactant_phases in self.phases): is_competing[A_rxn_other>TOL] = True if np.all(other_rxn.product_phases in self.phases): is_competing[A_rxn_other<-TOL] = True return is_competing def stability(self, T_bounds, P_bounds, other_rxns, TOL=1e-6): # NOTE: other_rxns must be a list of reactions of the same composition peratom=True peratom=False # self._reac_assemblage # A_rxn_vals = self.affinity(T_bounds, P_bounds, peratom=peratom) A_rxn_vals = self.affinity(T_bounds, P_bounds) assert np.all(np.abs(A_rxn_vals)<TOL), \ 'Rxn affinity must be equal to zero to within TOL.' N_TP = A_rxn_vals.size N_rxn = len(other_rxns) # A_rxn_others = np.zeros((N_rxn, N_TP)) # Calculate gibbs energy of rxn for all possible reactions stable = np.tile(True, N_TP) for ind, irxn in enumerate(other_rxns): is_competing = self.competing_rxn(T_bounds, P_bounds, irxn) stable[is_competing] = False # iA_rxn_other = irxn.affinity(T_bounds, P_bounds, peratom=peratom ) # # if other rxn phase assemblage is a subset of the current # # assemblage, then it is a valid competing reaction # if self.competing_rxn(irxn): # if (rxn._reac_assemblage.issubset(self._reac_assemblage) | \ # rxn._reac_assemblage.issubset(self._prod_assemblage) ): # G_rxn_other_a[ind] = iG_rxn_other_a # elif (rxn._prod_assemblage.issubset(self._reac_assemblage) | \ # rxn._prod_assemblage.issubset(self._prod_assemblage) ): # # Store negative energy as reverse reaction energy change # G_rxn_other_a[ind] = -iG_rxn_other_a # else: # G_rxn_other_a[ind] = 0.0 # print(G_rxn_other_a) # If any rxn is energetically favored, then current reaction is not # stable # stable_a = ~np.any(G_rxn_other_a < -TOL, axis=0) return stable def _validate_state_input(self, T, P, mols): N_T = np.asarray(T).size N_P = np.asarray(P).size N_solution_phases = len(mols) N_X = np.ones(N_solution_phases, dtype=int) for ind, iphs in enumerate(mols): imol = mols[iphs] if imol is None: N_X[ind] = 0 else: imol = np.asarray(imol) if imol.ndim==1: N_X[ind] = 1 else: N_X[ind] = imol.shape[0] assert np.all(N_X<=1), ( 'Multiple compositions not currently supported.' ) N_all = np.hstack((N_T, N_P, N_X)) N_PTX = np.max(N_all) assert np.all(np.isin(N_all, [0, 1, N_PTX])),( 'The number of state points for T, P, mol must be compatible.' ) state = {} state['N_T'] = N_T state['N_P'] = N_P state['N_X'] = N_X state['N_PTX'] = N_PTX # from IPython import embed;embed();import ipdb as pdb;pdb.set_trace() T = np.tile(np.asarray(T), N_PTX) if N_T < N_PTX else T P = np.tile(np.asarray(P), N_PTX) if N_P < N_PTX else P state['T'] = T state['P'] = P state['mols'] = mols return state def _calc_net_rxn_values(self, method_name, T, P, mols=None, peratom=False, use_endmember=False): mols = self._validate_mol_input(mols) phase_num = self.phase_num endmember_ids = self.endmember_ids rxn_coefs = self.rxn_coefs state = self._validate_state_input(T, P, mols) T, P, mols, N_PTX = (state.get(key) for key in ['T', 'P', 'mols', 'N_PTX']) values = np.zeros((phase_num, N_PTX)) for i, (iphs, iendmember) in enumerate(zip( self.phases, endmember_ids)): iphs_abbrev = iphs.abbrev imol = mols[iphs_abbrev] iphs_method = getattr(iphs, method_name) if use_endmember: ival = iphs_method(T, P, mol=imol, endmember=iendmember) else: ival = iphs_method(T, P, mol=imol) # if peratom: # ival /= self.rxn_atomnum # from IPython import embed;embed();import ipdb as pdb;pdb.set_trace() values[i] = ival net_rxn_values = np.dot(rxn_coefs, values) #from IPython import embed;embed();import pdb as pdb;pdb.set_trace() return net_rxn_values ###################### # NOT UPDATED # ###################### def gibbs_energy( self, T_a, P_a, peratom=False ): dG_rxn_a = self._calc_rxn_change('gibbs_energy_all', T_a, P_a, peratom=peratom) return dG_rxn_a def enthalpy( self, T_a, P_a, peratom=False ): dH_rxn_a = self._calc_rxn_change('enthalpy_all', T_a, P_a, peratom=peratom ) return dH_rxn_a def _calc_reac_value( self, method_name, T_a, P_a, peratom=False ): reac_method = getattr(self.reac_assemblage, method_name) val_phs_a = reac_method( T_a, P_a ) val_a = np.dot(self.reac_rxn_coef_a, val_phs_a) # if peratom: # val_a /= self.rxn_atomnum return val_a def _calc_prod_value( self, method_name, T_a, P_a, peratom=False ): prod_method = getattr(self.prod_assemblage, method_name) val_phs_a = prod_method( T_a, P_a ) val_a = np.dot(self.prod_rxn_coef_a, val_phs_a) # # if peratom: # # # val_a /= self.rxn_atomnum return val_a def _calc_rxn_change( self, method_name, T_a, P_a, peratom=False ): val_prod_a = self._calc_prod_value( method_name, T_a, P_a, peratom=peratom ) val_reac_a = self._calc_reac_value( method_name, T_a, P_a, peratom=peratom ) val_rxn_a = val_prod_a-val_reac_a return val_rxn_a
#=================================================== class Assemblage: def __init__(self, phase_objs, obj_is_classnm=False): # Get phase symbol list phase_symbols = [] for phase_obj in phase_objs: phase_symbols.append( phase_obj.props['abbrev']) phase_objs = np.array(phase_objs) phase_symbols = np.array(phase_symbols) indsort = np.argsort(phase_symbols) self._set_phase_assemblage( phase_objs[indsort], phase_symbols[indsort], obj_is_classnm=obj_is_classnm) pass def __eq__(self, other): return self._phases == other._phases def __lt__(self, other): if len(self._phases) < len(other._phases): return True return self._phases[0] < other._phases[0] def __gt__(self, other): if len(self._phases) > len(other._phases): return True return self._phases[0] > other._phases[0] def issubset(self, other): is_member = [phase in other._phases for phase in self._phases] return np.all(is_member) def _set_phase_assemblage(self, phase_objs, phase_symbols, obj_is_classnm=False): if obj_is_classnm: phase_classnms = phase_objs phase_objs = [ PurePhase(phase_classnm, phasesym) for (phase_classnm, phasesym) in zip(phase_classnms, phase_symbols)] self._phase_symbols = phase_symbols self._phases = phase_objs props = {} props['formula_all'] = [phs.props['formula'] for phs in phase_objs] props['name_all'] = [phs.props['phase_name'] for phs in phase_objs] props['abbrev_all'] = [phs.props['abbrev'] for phs in phase_objs ] props['molwt_all'] = np.array([ phs.props['molwt'] for phs in phase_objs]) props['element_comp_all'] = np.vstack([ phs.props['element_comp'] for phs in phase_objs]) # NOTE replace symbols? # props['element_symbols_all'] = phase_objs[0].props['element_symbols'] self._props = props pass def _validate_mol_input(self, mols): phases = self.phases phase_symbols = self.phase_symbols for symbol in phase_symbols: if symbol not in mols: mols[symbol] = None # else: # imol = mols[symbol] # iphase = phases[symbol] # iprops = iphase.props # iendmember_names = iprops['endmember_name'] # iendmember_num = len(iendmember_names) return mols @property def phase_symbols(self): return self._phase_symbols @property def phases(self): return self._phases @property def props(self): return self._props def get_endmember_comp_matrix(self): oxide_num = chem.oxide_props['oxide_num'] all_mol_oxide_comp = np.zeros((0,oxide_num)) all_phase_name = np.zeros((0)) all_endmember_name = np.zeros((0)) all_endmember_ind = np.zeros((0)) all_phase_ind = np.zeros((0)) for iphase in self.phases: iphs_props = iphase.props imol_oxide_comp = iphs_props['mol_oxide_comp'] iendmember_name = iphs_props['endmember_name'] iendmember_num = len(iendmember_name) iendmember_ind = np.arange(iendmember_num) iphase_name_tile = np.tile(np.array([iphs_props['phase_name']]), iendmember_num) all_phase_ind = np.hstack((all_phase_ind, np.tile(ind_phs,(iendmember_num)))) all_mol_oxide_comp= np.vstack((all_mol_oxide_comp, iall_mol_oxide_comp)) all_phase_name = np.hstack((all_phase_name, iphase_name_tile)) all_endmember_name = np.hstack((all_endmember_name, iendmember_name)) all_endmember_ind = np.hstack((all_endmember_ind, iendmember_ind)) return all_mol_oxide_comp def gibbs_energy_all(self, T, P, mols=None): mols = self._validate_mol_input(mols) gibbs_energy_all = [] for iphs in self.phases: iphs_abbrev = iphs.abbrev imol = mols[iphs_abbrev] igibbs_energy = iphs.gibbs_energy( T, P, mol=imol) gibbs_energy_all.append(igibbs_energy) gibbs_energy_all = np.vstack(gibbs_energy_all) return gibbs_energy_all def chem_potential_all(self, T, P, mols=None): mols = self._validate_mol_input(mols) chem_potential_all = [] for iphs in self.phases: iphs_abbrev = iphs.abbrev imol = mols[iphs_abbrev] ichem_potential = iphs.chem_potential( T, P, mol=imol) chem_potential_all.append(ichem_potential) chem_potential_all = np.vstack(chem_potential_all) return chem_potential_all def enthalpy_all(self, T, P): return np.vstack([phs.enthalpy(T, P) for phs in self.phases]) def entropy_all(self, T, P): return np.vstack([phs.entropy(T, P) for phs in self.phases]) def heat_capacity_all(self, T, P): return np.vstack([phs.heat_capacity(T, P) for phs in self.phases]) def dCp_dT_all(self, T, P): return np.vstack([phs.dCp_dT(T, P) for phs in self.phases]) def volume_all(self, T, P): return np.vstack([phs.volume(T, P) for phs in self.phases]) def dV_dT_all(self, T, P): return np.vstack([phs.dV_dT(T, P) for phs in self.phases]) def dV_dP_all(self, T, P): return np.vstack([phs.dV_dP(T, P) for phs in self.phases]) def d2V_dT2_all(self, T, P): return np.vstack([phs.d2V_dT2(T, P) for phs in self.phases]) def d2V_dTdP_all(self, T, P): return np.vstack([phs.d2V_dTdP(T, P) for phs in self.phases]) def d2V_dP2_all(self, T, P): return np.vstack([phs.d2V_dP2(T, P) for phs in self.phases]) #=================================================== class classproperty(property): # https://stackoverflow.com/questions/128573/using-property-on-classmethods def __get__(self, obj, objtype=None): return super(classproperty, self).__get__(objtype) def __set__(self, obj, value): super(classproperty, self).__set__(type(obj), value) def __delete__(self, obj): super(classproperty, self).__delete__(type(obj))
[docs]class Phase(metaclass=ABCMeta): """ Abstract parent class defining generic phase properties. The user must use a subclass, like PurePhase or SolutionPhase, which implements the Phase interface. Parameters ---------- phase_classnm : str Official class name for phase (implemented in src code). String has the form *classname* if source is *objc*, else: - if !calib ['cy', 'phase name', 'module name', ''] - if calib ['cy', 'phase name', 'module name', 'calib', ''] abbrev : str Official abbreviation of phase (regardless of implementation). calib : bool, default True Indicates whether sample phase should be calibration ready. source : str Code source for phase implementation. Default is 'objc' (code in objective-C that is part of the original code base). Alternative is 'coder' (code generated by the coder module). coder_module : str Name of the coder module that contains the phase classes. See documentation for model.Database for additional information and examples. Attributes ---------- abbrev calib class_name endmember_ids endmember_names endmember_num formula identifier module MOLWTS OXIDES param_names param_props phase_name phase_obj phase_type props source Notes ----- * This code is highly dependent on implementation and is likely to change dramatically with changes in the underlying code that calculates phase properties. * The pure water phase, "H2O", is very complex and thus not available for calibration. The water phase will force its calib flag to False, regardless of input value. """ _MINVAL = np.sqrt(np.finfo(float).eps) def __init__(self, phase_classnm, abbrev, calib=True, source='objc', coder_module=None): if abbrev=='H2O': calib=False if source == 'coder': invalidate_caches() self._module = import_module(coder_module) else: self._module = None self._calib = calib self._source = source self.deriv_keys = {'dT','dP','dV','dS','dmol'} self._init_phase_props(phase_classnm, abbrev) self._init_param_props() self._phase_type = None self._abbrev = abbrev self._OXIDES = chem.oxide_props['oxides'] self._MOLWTS = chem.oxide_props['molwt'] #Moved to class property #MINVAL = np.sqrt(np.finfo(float).eps) #self._MINVAL = MINVAL pass def _init_phase_props(self, phase_classnm, abbrev): if self.source == 'objc': phase_obj, phase_cls = core.get_src_object( phase_classnm, return_class=True) self._phase_obj = phase_obj self._phase_cls = phase_cls phase_name = phase_obj.phaseName class_name = phase_cls.name formula = phase_obj.phaseFormula identifier = 'Objective-C-base' self._phase_name = phase_name self._class_name = class_name self._formula = formula self._identifier = identifier elif self.source == 'coder': # if !calib ['cy', 'phase name', 'module name'] # if calib ['cy', 'phase name', 'module name', 'calib'] parts = phase_classnm.split("_") self._phase_obj = phase_classnm self._phase_cls = phase_classnm method = getattr(self.module, phase_classnm+'name') phase_name = method() self._phase_name = phase_name self._class_name = None try: method = getattr(self.module, phase_classnm+'formula') formula = method() self._formula = formula except AttributeError: self._formula = "" print ('Module generated by the coder package ' + 'does not yet provide a formula method.') except TypeError: self._formula = "" method = getattr(self.module, phase_classnm+'identifier') identifier = method() self._identifier = identifier class_name = phase_classnm props = OrderedDict() props['abbrev'] = abbrev props['phase_name'] = phase_name props['class_name'] = class_name props['identifier'] = identifier props['endmember_name'] = [None] props['endmember_ids'] = [None] #TK props['formula'] = [None] props['atom_num'] = [None] props['molwt'] = [None] props['elemental_entropy'] = [None] props['element_comp'] = [None] props['mol_oxide_comp'] = [None] self._props = props pass def _init_param_props(self): if self.source == 'objc': try: phase_obj = self._phase_obj supports_calib = phase_obj.supportsParameterCalibration() param_num = phase_obj.getNumberOfFreeParameters() param_names_NSArray = phase_obj.getArrayOfNamesOfFreeParameters() param_names = [param_names_NSArray.objectAtIndex_(i) for i in range(param_num)] param_units = np.array([ phase_obj.getUnitsForParameterName_(key) for key in param_names]) param0 = np.array([ phase_obj.getValueForParameterName_(key) for key in param_names]) except AttributeError: supports_calib = False param_num = 0 param_names = np.array([]) param_units = np.array([]) param0 = np.array([]) else: if 'calib_' in self._phase_obj: supports_calib = True method_stub = self._phase_obj.replace('calib_', '') method = getattr(self.module, method_stub+'get_param_number') param_num = method() method = getattr(self.module, method_stub+'get_param_names') param_names = method() method = getattr(self.module, method_stub+'get_param_units') param_units = method() method = getattr(self.module, method_stub+'get_param_values') param0 = method() else: supports_calib = False param_num = 0 param_names = np.array([]) param_units = np.array([]) param0 = np.array([]) param_props = OrderedDict() param_props['supports_calib'] = supports_calib param_props['param_num'] = param_num param_props['param_names'] = param_names param_props['param_units'] = param_units param_props['param0'] = param0 self._param_props = param_props pass def __eq__(self, other): return self._props['abbrev'] == other._props['abbrev'] def __lt__(self, other): return self._props['abbrev'] <= other._props['abbrev'] def __gt__(self, other): return self._props['abbrev'] >= other._props['abbrev'] @property def phase_name(self): """ Name of phase Returns ------- Name of phase (str) """ return self._phase_name @property def abbrev(self): """ Official unique abbreviation for phase Returns ------- Abbreviation (str) """ return self._abbrev @property def class_name(self): """ Name of class Returns ------- Name of class (str) """ return self._class_name @property def formula(self): """ Formula of phase Returns ------- Formula of phase (str) """ return self._formula @property def identifier(self): """ Identifier of phase Returns ------- Identifier of phase (str) """ return self._identifier @property def phase_type(self): """ Phase type Returns ------- Phase type (str) Permissible values are 'pure' or 'solution'. """ return self._phase_type @property def calib(self): """ Indicates whether phase calibration is enabled Returns ------- Value of calib (bool) """ return self._calib @property def source(self): """ Indicates origin of source code implementation Returns ------- String indicating origin of source code for implementation Permissible values 'objc' or 'coder'. """ return self._source @property def module(self): """ Python module attribute for coder generated functions Returns ------- module Module attribute returned from importlib.import_module, else None if source is 'objc' """ return self._module @property def props(self): """ Dictionary of phase properties The dictionary defines phase properties with these keys: ``abbrev`` : str Official unique phase abbreviation ``name`` : str Name of phase (implementation dependent) ``class_name`` : str Official class name for phase (implemented in src code) ``formula`` : str Formula of phase ``natom`` : int Number of atoms in formula unit ``molwt`` : double Molecular weight of phase (in g/mol-formula-unit) ``elemental_entropy`` : double Estimated entropy from elemental formula (from Robie et al. 1979) ``element_symbols`` : str array Symbol array string ``element_comp`` : int array Phase formula in terms of number of each element Returns ------- A Python dictionary (dict) Notes ----- Need to update these dictionary values to be vectors for solution phases """ return self._props @property def endmember_names(self): """ Name of each endmember Returns ------- List of endmember names for this solution phase, [str,...] """ return self.props['endmember_name'] @property def endmember_num(self): """ Number of endmembers in phase Returns ------- Number of endmembers in phase (int) """ return len(self.props['endmember_name']) @property def endmember_ids(self): """ ID number of each endmember in phase Returns ------- Array of ids, [int,...] """ return self.props['endmember_ids'] #TK @property def param_props(self): """ Dictionary of phase model parameters This dictionary defines parameter properties for the phase, using these keys: ``supports_calib`` : bool Flag indicating whether phase allows calibration ``param_num`` : int Number of parameters ``param_names`` : str array Name of each parameter ``param_units`` : str array Units for each parameter ``param0`` : double array Initial parameter values Returns ------- Dictionary of phase model parameters (dict) """ return self._param_props @property def phase_obj(self): """ Instance of the phase object Returns ------- Object instance """ return self._phase_obj @property def param_names(self): """ Array of parameter names Returns ------- Array of names for each parameter of the phase model, [str,...] """ return self.param_props['param_names']
[docs] def param_units(self, param_names=[], all_params =False): """ Get units for listed parameters. Parameters ---------- param_names : str array List of parameter names all_params : bool, default False If true, returns units for all parameters Returns ------- units : double array List of units for selected parameters """ if all_params: return self.param_props['param_units'] if type(param_names) == str: param_names = [param_names] units_l = [] for ind,key in enumerate(param_names): unit = self.param_props['param_units'][ind] # unit = self._phase_obj.getUnitsForParameterName_(key) units_l.append(unit) return units_l
@property def OXIDES(self): """ Array of oxide names Returns ------- Array of oxide names, [str,...] """ return self._OXIDES @property def MOLWTS(self): """ Array of molecular weights of oxides Returns ------- Numpy array of molecular weights, (nparray) """ return self._MOLWTS @classproperty def MINVAL(cls): """ Minimum molar concentration of a component in a solution phase Returns ------- Floating point number """ return cls._MINVAL @MINVAL.setter def MINVAL(cls, value): assert value > np.finfo(float).eps, 'MINVAL must have a value ' \ + 'greater than machine precision(' + str(np.finfo(float).eps) + ')' cls._MINVAL = value
[docs] def set_ref_state(self, Tr=298.15, Pr=1.0, Trl=298.15): """ Set reference state P/T conditions. Parameters ---------- Tr : double, default 298.15 Reference temperature in Kelvin Pr : double, default 1.0 Reference pressure in bars Trl : double, default 298.15 Reference temperature for lambda heat capacity correction in Kelvin """ if self.source == 'objc': self._phase_obj.setTr_(Tr) self._phase_obj.setPr_(Pr) self._phase_obj.setTrl_(Trl) pass
[docs] def enable_gibbs_energy_reference_state(self): """ Set Gibbs energy of the reference state. Notes ----- Call method on any phase class, and it automatically applies to all. """ if self.source == 'objc': phase_cls = self._phase_cls phase_cls.enableGibbsFreeEnergyReferenceStateUsed() pass
[docs] def disable_gibbs_energy_reference_state(self): """ Unset Gibbs energy of the reference state. Notes ----- Call method on any phase class, and it automatically applies to all. """ # call method on any phase class (automatically applied to all) if self.source == 'objc': phase_cls = self._phase_cls phase_cls.disableGibbsFreeEnergyReferenceStateUsed() pass
[docs] def get_param_values(self, param_names=[], all_params=False): """ Get current values for listed parameters. Parameters ---------- param_names : str array List of parameter names all_params : bool, default False If true, returns units for all parameters Returns ------- values : double array List of values for selected parameters """ if all_params: param_names = self.param_names elif type(param_names) == str: param_names = [param_names] values_a = np.zeros(len(param_names)) if self.source == 'objc': for ind,key in enumerate(param_names): value = self._phase_obj.getValueForParameterName_(key) values_a[ind] = value elif self.calib: method_stub = self._phase_obj.replace('calib_', '') method = getattr(self.module, method_stub+'get_param_value') for ind,key in enumerate(param_names): value = method(ind) values_a[ind] = value else: values_a = None return values_a
[docs] def set_param_values(self, param_names=[], param_values=[]): """ Set new values for listed parameters. Parameters ---------- param_names : str array List of parameter names param_values : double array List of parameter values """ assert len(param_names)==len(param_values), \ 'param_names and param_values must have the same length' if self.source == 'objc': for name, value in zip(param_names, param_values): self._phase_obj.setParameterName_tovalue_(name, float(value)) elif self.calib: method_stub = self._phase_obj.replace('calib_', '') method = getattr(self.module, method_stub+'set_param_value') for name, value in zip(param_names, param_values): method(name, value) pass
[docs] def get_phase_amount_from_elements(self, elements, kind='mass'): """ Convert list of elements to quantity of phase. Parameters ---------- elements : double array Number of each element kind : {'mass','moles'} Determines how phase amount is determined (mass vs. moles) Returns ------- amount : double Amount of phase (expressed according to kind) """ if self.source == 'objc': elements_vec = core.array_to_double_vector(elements) # NOTE: Need to update to allow for solution phases if kind == 'mass': return self._phase_obj.convertElementsToMassOfPhase_( elements_vec) elif kind == 'moles': return self._phase_obj.convertElementsToMolesOfPhase_( elements_vec) else : raise NotImplemented('The kind "'+kind+'" is not implimented ' + 'for get_phase_amount_from_elements()') else: raise NotImplemented('Function not implemented for coder ' + 'generated modules.')
def _est_site_mult(self, t, p, mu): """ Estimates site multiplicity of a solution phase """ x_ave = 1.0/mu.size s_mid = self.entropy(t, p, mol=np.full(mu.size, x_ave)) s_idl = -8.3143*np.log(x_ave) s = s_mid x_end = np.full(mu.size, 1.0e-8) for i in range(0,mu.size): x_end[i] = 0.99999998 s -= x_ave*self.entropy(t, p, mol=x_end) x_end[i] = 1.0e-8 return np.rint(s/s_idl) def _tailored_affinity_and_comp(self, t, p, mu, debug=False): """ Use a specific affinity and composition solver if available """ if self.identifier == 'Objective-C-base': method = getattr(self.phase_obj, 'affinityAndCompositionFromLiquidChemicalPotentialSum_andT_andP_') if method: cmu = core.array_to_ctype_array(mu) if debug > 0: print ('Calling tailored Affinity and Comp routine for', self.class_name) result = method(cmu, t, p) A = float(result.objectAtIndex_(0)) X = np.zeros(mu.size) for i in range(0,mu.size): X[i] = float(result.objectAtIndex_(i+1)) if debug > 0: print ('... Affinity ', A, 'J/mol') print ('... X', X) print ('... Convergence', result.objectAtIndex_(mu.size+1)) print ('... Iterations', result.objectAtIndex_(mu.size+2)) print ('... Affinity scalar', result.objectAtIndex_(mu.size+3)) print ('... Estimated error on affinity', result.objectAtIndex_(mu.size+4)) return (A, X) else: return None else: return None
[docs] def affinity_and_comp(self, t, p, mu, debug=False, method='generic'): """ Given a temperature, pressure and chemical potential(s), compute and return a chemical affinity and phase composition. Parameters ---------- t,p : float Temperature (K) and pressure (bars) mu : ndarray Chemical potential(s) of endmember components in phase - For a stoichiometric phase, a 1-component 1-D numpy array - For a solution phase, a 1-D numpy array of length equal to the number of endmember components in the phase debug : bool, def False Print iteration details method : str, def 'generic' Algorithm used: 'generic' - generic algorithm on components, 'species' - generic algorithm on species (> = components), 'special' - algorthmic specifically coded in the phase instance, currently, only Objective-C coded methods are utilized, see _tailored_affinity_and_comp() Returns ------- A, X : tuple A is a scalar equal to the chemical affinity of the phase relative to mu. X is a numpy array of mole fractions of endmember components in the phase. X is the same length as mu. Notes ----- The algorithm used is from Ghiorso (2013) and is a globally convergent saturation state algorithm applicable to thermodynamic systems with a stable or metastable omni-component phase. Geochimica et Cosmochimica Acta, 103, 295-300 """ if mu.size == 1: if mu[0] == 0.0: return (999999.0, np.array([1.0])) else: mu0 = self.chem_potential(t, p) if debug: print (mu[0], mu0) return (mu0-mu[0], np.array([1.0])) else: if method == 'special': result = self._tailored_affinity_and_comp(t, p, mu, debug=debug) if result is not None: return result X_mask = np.nonzero(mu)[0] if debug: print ('Mole fraction mask: ', X_mask) X = np.full(mu.size, 1.0e-8) muM = np.zeros(mu.size) phi = np.zeros(mu.size) A_previous = sys.float_info.max loop = True n_loop = 0 site_m = self._est_site_mult(t, p, mu) if debug: print ('Computed site multiplicity:', site_m) while loop: for ind in X_mask: if n_loop == 0: X[ind] = 1.0 # muM = mu0 + RTln X + RT ln gamma (for the endmember) muM[ind] = self.chem_potential(t, p, mol=X, endmember=ind) if n_loop == 0: X[ind] = 1.0e-8 phi[ind] = mu[ind] - muM[ind] if n_loop > 0: # the phi now contains only the mu0 and RT ln gamma phi[ind] += 8.3143*t*site_m*np.log(X[ind]) for i in range(0,X.size): X[i] = 0.0 denom = 1.0 prod = 1.0 for i in range(0, X_mask[:-1].size): ind = X_mask[i] ind_next = X_mask[i+1] r = np.exp((phi[ind_next]-phi[ind])/(8.3143*t*site_m)) prod *= r denom += prod X[ind_next] = prod ind0 = X_mask[0] X[ind0] = 1.0/denom for i in range(1,X_mask.size): ind = X_mask[i] X[ind] *= X[ind0] A = -(phi[ind0] - 8.3143*t*site_m*np.log(X[ind0])) if debug: print('iter: {0:2d}'.format(n_loop),end=' ') print('Anew: {0:10.2f}'.format(A),end=' ') if n_loop > 0: print('Aold: {0:10.2f}'.format(A_previous), end=' ') else: print('Aold: ', end=' ') print ('X: ', X) loop = abs(A - A_previous) > 0.1 n_loop += 1 A_previous = A if n_loop > 50: loop = False return (A, X)
def _validate_deriv(self, deriv): if not deriv: deriv = {} else: assert all([key in self.deriv_keys for key in deriv]), \ 'Specified deriv keys, '+str(deriv)+' are not allowed. ' + \ 'They must belong to the allowable set: '+str(self.deriv_keys) # Check for all zeros (default) if all(deriv[key]==0 for key in deriv): deriv = {} ''' Insert some tests here for values of derivative keys: 'dT' = <3 'dP' = <3 'dmol' = <2 'dw' <= 1 ''' return deriv def _nudge_solution_comp(self, mol): """ Adjust solution composition to ensure that all components are present (at machine precision). """ MINVAL = self.MINVAL ind_row, ind_col = np.where(mol==0) for irow, icol in zip(ind_row, ind_col): mol[ind_row, ind_col] = MINVAL return mol #################################################### # Decorator for validating thermodynamic functions # #################################################### def validate(func): # *args arguments as tuple for a in args # **kwargs dictionary of keyword arguments for a in kwargs def func_wrapper(*args, **kwargs): # execute before function call self = args[0] T = args[1] P = args[2] mol = kwargs['mol'] if ('mol' in kwargs) else None V = kwargs['V'] if ('V' in kwargs) else None deriv = kwargs['deriv'] if ('deriv' in kwargs) else None deriv_param = kwargs['deriv_param'] if ( 'deriv_param' in kwargs) else None mol_deriv_qty = kwargs['mol_deriv_qty'] if ( 'mol_deriv_qty' in kwargs) else False endmember = kwargs['endmember'] if ( 'endmember' in kwargs) else None const = kwargs['const'] if 'const' in kwargs else None species = kwargs['species'] if 'species' in kwargs else None if self.phase_type == 'pure': assert not mol, 'mol must be empty for PurePhase.' deriv = {'dT':0, 'dP':0, 'dV':0, 'dS':0, 'dmol':0} if ( deriv is None) else deriv endmember = -1 if endmember is None else endmember if deriv_param is not None: if not self._param_props['supports_calib']: assert("This module does not support calibration.") if type(deriv_param) is int: deriv_param = [deriv_param] elif type(deriv_param) is str: if deriv_param in self._param_props['param_names']: deriv_param = ( self._param_props['param_names'].index(deriv_param)) else: assert("Parameter name is in param_props list. ") deriv_param = [deriv_param] elif type(deriv_param) is list: for ind,(x) in enumerate(deriv_param): if type(x) is str: if x in self._param_props['param_names']: deriv_param[ind] = ( self._param_props['param_names'].index(x)) else: assert("Parameter name is in param_props list. ") elif type(deriv_param) is np.ndarray: pass else: print ('deriv_param argument must be one of int, ' + 'str, list, or numpy array.') scalar_input = False if V is None: T, P = core.fill_array(T, P) if len(T)==1 & len(P)==1: scalar_input = True else: # overide pressure with input volume P = None T, V = core.fill_array(T, V) if len(T)==1 & len(V)==1: scalar_input = True if mol_deriv_qty and endmember is not None: mol_deriv_qty = False T, endmember = core.fill_array(T, endmember) deriv = self._validate_deriv(deriv) if 'dmol' in deriv and deriv['dmol']>0: scalar_input = False if mol is not None: mol = np.array(mol) num_PT_pts = len(T) if mol.ndim==1: mol = np.tile(mol[np.newaxis, :], (num_PT_pts, 1)) elif mol.ndim==2: assert mol.shape[0]==num_PT_pts, ('molar composition must ' + 'be provided for every PT point, or otherwise be ' + 'fixed to a single value for all PT points.') else: assert False, ('Molar composition must be defined, for ' + 'each PT point. It currently has too many ' + 'dimensions.') mol = self._nudge_solution_comp(mol) assert mol.size > 0, ('mol cannot be initialized to an empty ' + 'numpy array or to an empty list') # call the function that is wrapped if const: result = func(self, T, P, mol=mol, V=V, deriv=deriv, deriv_param=deriv_param, mol_deriv_qty=mol_deriv_qty, endmember=endmember, const=const) elif species: result = func(self, T, P, mol=mol, V=V, deriv=deriv, deriv_param=deriv_param, mol_deriv_qty=mol_deriv_qty, endmember=endmember, species=species) else: result = func(self, T, P, mol=mol, V=V, deriv=deriv, deriv_param=deriv_param, mol_deriv_qty=mol_deriv_qty, endmember=endmember) # execute after function call if result is None: raise NotImplementedError('Cannot find specified derivative: ' +str(deriv)) if scalar_input and not mol_deriv_qty: result = result[0] return result return func_wrapper ################### # State Functions # ################### ########### # order 0 # ########### @validate def gibbs_energy(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate Gibbs energy (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin P : array-like Pressure in bars mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases) V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (Default is zero for all params.) deriv_param : array of strs Parameter names that identify returned derivatives (Default is None.) Returns ------- value : array-like Gibbs energy in J (or derivative in deriv units) """ method = None args_dw = False if (not deriv) and (not deriv_param): method = getattr(self, '_calc_G') else: nT = deriv['dT'] if ('dT' in deriv) else 0 nP = deriv['dP'] if ('dP' in deriv) else 0 nMol = deriv['dmol'] if ('dmol' in deriv) else 0 nW = 0 if not deriv_param else 1 nTot = nT + nP + nMol + nW method_sig = '_calc_' method_sig += 'd' + str(nTot) + 'G_' if (nTot > 1) else ( 'dG_' if (nTot > 0) else 'G') method_sig += 'dT' + str(nT) if (nT > 1) else ( 'dT' if (nT > 0) else '') method_sig += 'dP' + str(nP) if (nP > 1) else ( 'dP' if (nP > 0) else '') method_sig += 'dm' + str(nMol) if (nMol > 1) else ( 'dm' if (nMol > 0) else '') method_sig += 'dw' if (nW > 0) else '' method = getattr(self, method_sig) args_dw = False if nW == 0 else True if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, param=deriv_param ) if args_dw else method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, param=deriv_param ) if args_dw else method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def enthalpy(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate enthalpy (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin P : array-like Pressure in bars mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Enthalpy in J (or derivative in deriv units). """ method = None # Fix me! if (not deriv) and (not deriv_param): method = getattr(self, '_calc_H') else: nT = deriv['dT'] if ('dT' in deriv) else 0 nP = deriv['dP'] if ('dP' in deriv) else 0 nMol = deriv['dmol'] if ('dmol' in deriv) else 0 nW = 0 if not deriv_param else 1 nTot = nT + nP + nMol + nW assert nT == 0 and nP == 0 and nMol == 0 and nW == 1, \ 'Only parameter derivatives of the enthalpy are permitted' method_sig = '_calc_' method_sig += 'd' + str(nTot) + 'H_' if (nTot > 1) else ( 'dH_' if (nTot > 0) else 'H') method_sig += 'dw' if (nW > 0) else '' method = getattr(self, method_sig) if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP, imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def helmholtz_energy(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate helmholtz energy (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Helmholtz energy in J (or derivative in deriv units). """ return None @validate def internal_energy(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate internal energy (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Internal energy in J (or derivative in deriv units). """ return None ########### # order 1 # ########### @validate def volume(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate volume (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Volume in J/bar (or derivative in deriv units). """ method = None #fix me if (not deriv) and (not deriv_param): method = getattr(self, '_calc_V') else: nT = deriv['dT'] if ('dT' in deriv) else 0 nP = deriv['dP'] if ('dP' in deriv) else 0 nMol = deriv['dmol'] if ('dmol' in deriv) else 0 nW = 0 if not deriv_param else 1 nTot = nT + nP + nMol + nW method_sig = '_calc_' method_sig += 'd' + str(nTot) + 'V_' if (nTot > 1) else ( 'dV_' if (nTot > 0) else 'V') method_sig += 'dm' + str(nMol) if (nMol > 1) else ( 'dm' if (nMol > 0) else '') method_sig += 'dT' + str(nT) if (nT > 1) else ( 'dT' if (nT > 0) else '') method_sig += 'dP' + str(nP) if (nP > 1) else ( 'dP' if (nP > 0) else '') method_sig += 'dw' if (nW > 0) else '' method = getattr(self, method_sig) if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def density(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate density (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Density in g*bar/J (or derivative in deriv units). """ method = None if (not deriv) and (not deriv_param): method = getattr(self, '_calc_density') if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def entropy(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=None, endmember=None, const=None, species=None): """ Calculate entropy (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Entropy in J/K (or derivative in deriv units). """ method = None #fix me if (not deriv) and (not deriv_param): method = getattr(self, '_calc_S') else: nT = deriv['dT'] if ('dT' in deriv) else 0 nP = deriv['dP'] if ('dP' in deriv) else 0 nMol = deriv['dmol'] if ('dmol' in deriv) else 0 nW = 0 if not deriv_param else 1 nTot = nT + nP + nMol + nW method_sig = '_calc_' method_sig += 'd' + str(nTot) + 'S_' if (nTot > 1) else ( 'dS_' if (nTot > 0) else 'S') method_sig += 'dm' + str(nMol) if (nMol > 1) else ( 'dm' if (nMol > 0) else '') method_sig += 'dT' + str(nT) if (nT > 1) else ( 'dT' if (nT > 0) else '') method_sig += 'dP' + str(nP) if (nP > 1) else ( 'dP' if (nP > 0) else '') method_sig += 'dw' if (nW > 0) else '' method = getattr(self, method_sig) if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def chem_potential(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, species=False, const=None): """ Calculate chemical potential (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). endmember : None or int scalar or int array If None, retrieve an array of chemical potentials, else chemical potential for endmember index or index set in array species : boolean If False, returned value is for components of the solution. If True, returned value is for species in the solution. Returns ------- value : array-like Chemical potential in J (or derivative in deriv units). """ if not deriv: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(self._calc_mu(iT, iP, V=V, endmember=endmember, species=species)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(self._calc_mu(iT, iP, mol=imol, V=V, endmember=endmember, species=species)) result = np.array(result) else: result = None return result @validate def activity(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, const=None, species=None): """ Calculate activity (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). endmember : None or int scalar or int array If None, retrieve an array of chemical potentials; else chemical potential for endmber index or index set in array Returns ------- value : array-like Activity (or derivatives in deriv units). """ method = None if not deriv: method = getattr(self, '_calc_a') elif deriv=={'dmol':1}: method = getattr(self, '_calc_da_dm') else: pass if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V, endmember=endmember)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V, endmember=endmember)) result = np.array(result) else: result = None return result @validate def fugacity(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, const=None, species=None): """ Calculate activity (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Fugacity in bars (or derivatives in deriv units). """ return None ########### # order 2 # ########### @validate def thermal_exp(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, const=None, species=None): """ Calculate heat capacity (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Thermal expansion in 1/K (or derivatives in deriv units). """ method = None if (not deriv) and (not deriv_param): method = getattr(self, '_calc_alpha') if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def gamma(self, T, P, mol=None, V=None, deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, const=None, species=None): """ Calculate grüneisen parameter (or derivatives) for phase. NOT IMPLEMENTED """ return None @validate def heat_capacity(self, T, P, mol=None, V=None, const='P', deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, species=None): """ Calculate heat capacity (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). const : ['P', 'V'], optional Defines constant path for derivative (yielding C_P vs C_V) V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Thermal expansion in J/K (or derivatives in deriv units). """ method = None if const=='P': if not deriv: method = getattr(self, '_calc_Cp') elif deriv=={'dT':1}: method = getattr(self, '_calc_dCp_dT') elif deriv=={'dmol':1}: method = getattr(self, '_calc_dCp_dm') else: result = None elif const=='V': if not deriv: method = getattr(self, '_calc_Cv') if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None return result @validate def bulk_mod(self, T, P, mol=None, V=None, const='T', deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, species=None): """ Calculate bulk modulus (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). const : ['T', 'S'], optional Defines constant path for derivative (yielding K_T vs K_S) V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Bulk modulus in bars (or derivatives in deriv units). """ if const=='T': method = None if (not deriv) and (not deriv_param): method = getattr(self, '_calc_K') else: nP = deriv['dP'] if ('dP' in deriv) else 0 if nP == 1: method = getattr(self, '_calc_Kp') if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None else: # const=='S' result = None return result @validate def compressibility(self, T, P, mol=None, V=None, const='T', deriv=None, deriv_param=None, mol_deriv_qty=True, endmember=None, species=None): """ Calculate compressibility (or derivatives) for phase. Parameters ---------- T : array-like Temperature in Kelvin. P : array-like Pressure in bars. mol : array-like, optional Composition in terms of mols of endmembers (unneeded for pure phases). const : ['T', 'S'], optional Defines constant path for derivative (yielding Beta_T vs Beta_S) V : array-like, optional (default None) Volume in J/bar. Overrides pressure if not None. deriv : dict of ints Derivative order for each parameter (default is zero for all params). Returns ------- value : array-like Bulk modulus in 1/bars (or derivatives in deriv units). """ if const=='T': method = None if (not deriv) and (not deriv_param): method = getattr(self, '_calc_beta') if method: result = [] if mol is None: for ind,(iT,iP) in enumerate(zip(T,P)): result.append(method(iT, iP, V=V)) else: for ind,(iT,iP,imol) in enumerate(zip(T,P,mol)): result.append(method(iT, iP, mol=imol, V=V)) result = np.array(result) else: result = None else: # const=='S' result = None return result ###################### # Calculator methods # ###################### @abstractmethod def _calc_G(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dG_dT(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dG_dP(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dG_dm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dG_dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d2G_dT2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2G_dTdP(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2G_dTdm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2G_dTdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d2G_dP2(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d2G_dPdm(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d2G_dPdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d2G_dm2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2G_dmdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dT3(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dT2dP(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dT2dm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dT2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dTdP2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dTdPdm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dTdPdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dTdm2(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dTdmdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dP3(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dP2dm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dP2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dPdm2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dPdmdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d3G_dm3(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d3G_dm2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dT2dmdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dTdPdmdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dP2dmdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dm3dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dTdm2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dPdm2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dT3dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dT2dPdw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dTdP2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d4G_dP3dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d5G_dT2dm2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d5G_dTdPdm2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_d5G_dP2dm2dw(self, T, P, mol=None, param=None): pass @abstractmethod def _calc_H(self, T, P, mol=None, V=None): pass # (order=1) @abstractmethod def _calc_S(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_V(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_Cv(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_Cp(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dCp_dT(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dV_dT(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dV_dP(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2V_dT2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2V_dTdP(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2V_dP2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_density(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_alpha(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_beta(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_K(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_Kp(self, T, P, mol=None, V=None): pass # Compositional quantities @abstractmethod def _calc_mu(self, T, P, mol=None, V=None, endmember=None, species=False): pass @abstractmethod def _calc_a(self, T, P, mol=None, V=None, endmember=None): pass @abstractmethod def _calc_dV_dm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_dS_dm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_da_dm(self, T, P, mol=None, V=None, endmember=None): pass # (order=3) @abstractmethod def _calc_dCp_dm(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2V_dmdT(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2V_dmdP(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2S_dm2(self, T, P, mol=None, V=None): pass @abstractmethod def _calc_d2V_dm2(self, T, P, mol=None, V=None): pass
#===================================================
[docs]class PurePhase(Phase): """ Pure stoichiometric phases. Implements the Phase interface. Parameters ---------- phase_classnm : str Official class name for phase (implemented in src code) abbrev : str Official abbreviation of phase (regardless of implementation) calib : bool, default True Indicates whether sample phase should be calibration ready Attributes ---------- Berman_formula Notes ----- * This code is highly dependent on implementation and is likely to change dramatically with changes in the underlying code that calculates phase properties. * The pure water phase, "H2O", is very complex and thus not available for calibration. The water phase will force its calib flag to False, regardless of input value. * In addition to the attributes listed, this class inherits the Phase class attributes. """ def __init__(self, phase_classnm, abbrev, calib=True, source='objc', coder_module=None): super().__init__(phase_classnm, abbrev, calib=calib, source=source, coder_module=coder_module) self._init_pure_phase_props(phase_classnm) self._phase_type = 'pure' self._methods = {} if self.source == 'objc': self._methods['_calc_G'] = getattr(self._phase_obj, 'getGibbsFreeEnergyFromT_andP_') self._methods['_calc_dG_dT'] = getattr(self, 'not_coded') self._methods['_calc_dG_dP'] = getattr(self._phase_obj, 'getVolumeFromT_andP_') self._methods['_calc_dG_dm'] = getattr(self, 'not_coded') self._methods['_calc_dG_dw'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dT2'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dTdP'] = getattr(self._phase_obj, 'getDvDtFromT_andP_') self._methods['_calc_d2G_dTdm'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dTdw'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dP2'] = getattr(self._phase_obj, 'getDvDpFromT_andP_') self._methods['_calc_d2G_dPdm'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dPdw'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dm2'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT3'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT2dP'] = getattr(self._phase_obj, 'getD2vDt2FromT_andP_') self._methods['_calc_d3G_dT2dm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT2dw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdP2'] = getattr(self._phase_obj, 'getD2vDtDpFromT_andP_') self._methods['_calc_d3G_dTdPdm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdPdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdm2'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dP3'] = getattr(self._phase_obj, 'getD2vDp2FromT_andP_') self._methods['_calc_d3G_dP2dm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dP2dw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dPdm2'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dm3'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dP2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dm3dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT3dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT2dPdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdP2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dP3dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dT2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dTdPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dP2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_H'] = getattr(self._phase_obj, 'getEnthalpyFromT_andP_') self._methods['_calc_S'] = getattr(self._phase_obj, 'getEntropyFromT_andP_') self._methods['_calc_dS_dm'] = getattr(self, 'not_coded') self._methods['_calc_d2S_dm2'] = getattr(self, 'not_coded') self._methods['_calc_V'] = getattr(self._phase_obj, 'getVolumeFromT_andP_') self._methods['_calc_dV_dT'] = getattr(self._phase_obj, 'getDvDtFromT_andP_') self._methods['_calc_dV_dP'] = getattr(self._phase_obj, 'getDvDpFromT_andP_') self._methods['_calc_dV_dm'] = getattr(self, 'not_coded') self._methods['_calc_d2V_dT2'] = getattr(self._phase_obj, 'getD2vDt2FromT_andP_') self._methods['_calc_d2V_dTdP'] = getattr(self._phase_obj, 'getD2vDtDpFromT_andP_') self._methods['_calc_d2V_dP2'] = getattr(self._phase_obj, 'getD2vDp2FromT_andP_') self._methods['_calc_d2V_dmdT'] = getattr(self, 'not_coded') self._methods['_calc_d2V_dmdP'] = getattr(self, 'not_coded') self._methods['_calc_d2V_dm2'] = getattr(self, 'not_coded') self._methods['_calc_Cv'] = getattr(self, 'not_coded') self._methods['_calc_Cp'] = getattr(self._phase_obj, 'getHeatCapacityFromT_andP_') self._methods['_calc_dCp_dT'] = getattr(self._phase_obj, 'getDcpDtFromT_andP_') self._methods['_calc_dCp_dm'] = getattr(self, 'not_coded') self._methods['_calc_density'] = getattr(self, 'not_coded') self._methods['_calc_alpha'] = getattr(self, 'not_coded') self._methods['_calc_beta'] = getattr(self, 'not_coded') self._methods['_calc_K'] = getattr(self, 'not_coded') self._methods['_calc_Kp'] = getattr(self, 'not_coded') self._methods['_calc_mu'] = getattr(self._phase_obj, 'getGibbsFreeEnergyFromT_andP_') self._methods['_calc_a'] = getattr(self, 'not_coded') self._methods['_calc_da_dm'] = getattr(self, 'not_coded') else: phase_calibnm = phase_classnm.replace('calib', 'dparam') if ( self.calib) else '' self._methods['_calc_G'] = getattr(self.module, phase_classnm+'g') self._methods['_calc_dG_dT'] = getattr(self.module, phase_classnm+'dgdt') self._methods['_calc_dG_dP'] = getattr(self.module, phase_classnm+'dgdp') self._methods['_calc_dG_dm'] = getattr(self, 'not_coded') self._methods['_calc_dG_dw'] = getattr(self.module, phase_calibnm+'g') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d2G_dT2'] = getattr(self.module, phase_classnm+'d2gdt2') self._methods['_calc_d2G_dTdP'] = getattr(self.module, phase_classnm+'d2gdtdp') self._methods['_calc_d2G_dTdm'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dTdw'] = getattr(self.module, phase_calibnm+'dgdt') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d2G_dP2'] = getattr(self.module, phase_classnm+'d2gdp2') self._methods['_calc_d2G_dPdm'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dPdw'] = getattr(self.module, phase_calibnm+'dgdp') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d2G_dm2'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT3'] = getattr(self.module, phase_classnm+'d3gdt3') self._methods['_calc_d3G_dT2dP'] = getattr(self.module, phase_classnm+'d3gdt2dp') self._methods['_calc_d3G_dT2dm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT2dw'] = getattr(self.module, phase_calibnm+'d2gdt2') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dTdP2'] = getattr(self.module, phase_classnm+'d3gdtdp2') self._methods['_calc_d3G_dTdPdm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdPdw'] = getattr(self.module, phase_calibnm+'d2gdtdp') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dTdm2'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dP3'] = getattr(self.module, phase_classnm+'d3gdp3') self._methods['_calc_d3G_dP2dm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dP2dw'] = getattr(self.module, phase_calibnm+'d2gdp2') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dPdm2'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dm3'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dP2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dm3dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT3dw'] = getattr(self.module, phase_calibnm+'d3gdt3') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dT2dPdw'] = getattr(self.module, phase_calibnm+'d3gdt2dp') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dTdP2dw'] = getattr(self.module, phase_calibnm+'d3gdtdp2') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dP3dw'] = getattr(self.module, phase_calibnm+'d3gdp3') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d5G_dT2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dTdPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dP2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_H'] = getattr(self, 'not_coded') self._methods['_calc_S'] = getattr(self.module, phase_classnm+'s') self._methods['_calc_dS_dm'] = getattr(self, 'not_coded') self._methods['_calc_d2S_dm2'] = getattr(self, 'not_coded') self._methods['_calc_V'] = getattr(self.module, phase_classnm+'v') self._methods['_calc_dV_dT'] = getattr(self.module, phase_classnm+'d2gdtdp') self._methods['_calc_dV_dP'] = getattr(self.module, phase_classnm+'d2gdp2') self._methods['_calc_dV_dm'] = getattr(self, 'not_coded') self._methods['_calc_d2V_dT2'] = getattr(self.module, phase_classnm+'d3gdt2dp') self._methods['_calc_d2V_dTdP'] = getattr(self.module, phase_classnm+'d3gdtdp2') self._methods['_calc_d2V_dP2'] = getattr(self.module, phase_classnm+'d3gdp3') self._methods['_calc_d2V_dmdT'] = getattr(self, 'not_coded') self._methods['_calc_d2V_dmdP'] = getattr(self, 'not_coded') self._methods['_calc_d2V_dm2'] = getattr(self, 'not_coded') self._methods['_calc_Cv'] = getattr(self.module, phase_classnm+'cv') self._methods['_calc_Cp'] = getattr(self.module, phase_classnm+'cp') self._methods['_calc_dCp_dT'] = getattr(self.module, phase_classnm+'dcpdt') self._methods['_calc_dCp_dm'] = getattr(self, 'not_coded') self._methods['_calc_density'] = getattr(self, 'not_coded') self._methods['_calc_alpha'] = getattr(self.module, phase_classnm+'alpha') self._methods['_calc_beta'] = getattr(self.module, phase_classnm+'beta') self._methods['_calc_K'] = getattr(self.module, phase_classnm+'K') self._methods['_calc_Kp'] = getattr(self.module, phase_classnm+'Kp') self._methods['_calc_mu'] = getattr(self.module, phase_classnm+'g') self._methods['_calc_a'] = getattr(self, 'not_coded') self._methods['_calc_da_dm'] = getattr(self, 'not_coded') def _init_pure_phase_props(self, phase_classnm): phase_obj = self.phase_obj if self.source == 'objc': element_comp = core.double_vector_to_array( phase_obj.formulaAsElementArray) else: method = getattr(self.module, phase_classnm+'elements') element_comp = method() atom_num = np.sum(element_comp) elem_num = element_comp.size mol_oxide_comp = chem.calc_mol_oxide_comp(element_comp) props = self.props props['endmember_name'] = np.array([str(props['phase_name'])]) props['endmember_id'] = np.array([0]) #TK if self.source == 'objc': props['formula'] = np.array([phase_obj.phaseFormula]) props['molwt'] = np.array([phase_obj.mw]) props['elemental_entropy'] = np.array( [phase_obj.entropyFromRobieEtAl1979]) else: method = getattr(self.module, phase_classnm+'formula') props['formula'] = np.array([method()]) method = getattr(self.module, phase_classnm+'mw') props['molwt'] = np.array([method()]) props['elemental_entropy'] = np.array([None]) props['atom_num'] = np.array([atom_num]) props['element_comp'] = element_comp[np.newaxis,:] props['mol_oxide_comp'] = mol_oxide_comp[np.newaxis, :] @property def Berman_formula(self): """ Representation of formula using Berman format Returns ------- Chemical formula of phase (str) """ return chem.get_Berman_formula(self.props['element_comp']) def not_coded(self, T, P, mol=None, V=None, param=None): raise AttributeError( 'Function not implemented for this database model.') def _calc_G(self, T, P, mol=None, V=None): return (self._methods['_calc_G'])(T, P) def _calc_dG_dT(self, T, P, mol=None, V=None): return (self._methods['_calc_dG_dT'])(T, P) def _calc_dG_dP(self, T, P, mol=None, V=None): return (self._methods['_calc_dG_dP'])(T, P) def _calc_dG_dm(self, T, P, mol=None, V=None): return (self._methods['_calc_dG_dm'])(T, P) def _calc_dG_dw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_dG_dw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_dG_dw'])(T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d2G_dT2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2G_dT2'])(T, P) def _calc_d2G_dTdP(self, T, P, mol=None, V=None): return (self._methods['_calc_d2G_dTdP'])(T, P) def _calc_d2G_dTdm(self, T, P, mol=None, V=None): return (self._methods['_calc_d2G_dTdm'])(T, P) def _calc_d2G_dTdw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d2G_dTdw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d2G_dTdw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d2G_dP2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2G_dP2'])(T, P) def _calc_d2G_dPdm(self, T, P, mol=None, V=None): return (self._methods['_calc_d2G_dPdm'])(T, P) def _calc_d2G_dPdw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d2G_dPdw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d2G_dPdw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d2G_dm2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2G_dm2'])(T, P) def _calc_d2G_dmdw(self, T, P, mol=None, param=None): return (self._methods['_calc_d2G_dmdw'])(T, P, param) def _calc_d3G_dT3(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dT3'])(T, P) def _calc_d3G_dT2dP(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dT2dP'])(T, P) def _calc_d3G_dT2dm(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dT2dm'])(T, P) def _calc_d3G_dT2dw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d3G_dT2dw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d3G_dT2dw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dTdP2(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dTdP2'])(T, P) def _calc_d3G_dTdPdm(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dTdPdm'])(T, P) def _calc_d3G_dTdPdw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d3G_dTdPdw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d3G_dTdPdw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dTdm2(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dTdm2'])(T, P) def _calc_d3G_dTdmdw(self, T, P, mol=None, param=None): return (self._methods['_calc_d3G_dTdmdw'])(T, P, param) def _calc_d3G_dP3(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dP3'])(T, P) def _calc_d3G_dP2dm(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dP2dm'])(T, P) def _calc_d3G_dP2dw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d3G_dP2dw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d3G_dP2dw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dPdm2(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dPdm2'])(T, P) def _calc_d3G_dPdmdw(self, T, P, mol=None, param=None): return (self._methods['_calc_d3G_dPdmdw'])(T, P, param) def _calc_d3G_dm3(self, T, P, mol=None, V=None): return (self._methods['_calc_d3G_dm3'])(T, P) def _calc_d3G_dm2dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d3G_dm2dw'])(T, P, param) def _calc_d4G_dT2dmdw(self, T, P, mol=None, param=None): return (self._methods['_calc_d4G_dT2dmdw'])(T, P, param) def _calc_d4G_dTdPdmdw(self, T, P, mol=None, param=None): return (self._methods['_calc_d4G_dTdPdmdw'])(T, P, param) def _calc_d4G_dTdm2dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d4G_dTdm2dw'])(T, P, param) def _calc_d4G_dP2dmdw(self, T, P, mol=None, param=None): return (self._methods['_calc_d4G_dP2dmdw'])(T, P, param) def _calc_d4G_dPdm2dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d4G_dPdm2dw'])(T, P, param) def _calc_d4G_dm3dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d4G_dm3dw'])(T, P, param) def _calc_d4G_dT3dw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d4G_dT3dw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dT3dw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dT2dPdw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d4G_dT2dPdw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dT2dPdw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dTdP2dw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d4G_dTdP2dw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dTdP2dw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dP3dw(self, T, P, mol=None, param=None): result = [] if self.source == 'objc': iresult_arr = (self._methods['_calc_d4G_dP3dw'])(T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dP3dw'])( T, P, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d5G_dT2dm2dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d5G_dT2dm2dw'])(T, P, param) def _calc_d5G_dTdPdm2dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d5G_dTdPdm2dw'])(T, P, param) def _calc_d5G_dP2dm2dw(self, T, P, mol=None, param=None): return (self._methods['_calc_d5G_dP2dm2dw'])(T, P, param) def _calc_H(self, T, P, mol=None, V=None): return (self._methods['_calc_H'])(T, P) def _calc_S(self, T, P, mol=None, V=None): return (self._methods['_calc_S'])(T, P) def _calc_dS_dm(self, T, P, mol=None, V=None): return (self._methods['_calc_dS_dm'])(T, P) def _calc_d2S_dm2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2S_dm2'])(T, P) def _calc_Cv(self, T, P, mol=None, V=None): return (self._methods['_calc_Cv'])(T, P) def _calc_Cp(self, T, P, mol=None, V=None): return (self._methods['_calc_Cp'])(T, P) def _calc_dCp_dT(self, T, P, mol=None, V=None): return (self._methods['_calc_dCp_dT'])(T, P) def _calc_dCp_dm(self, T, P, mol=None, V=None): return (self._methods['_calc_dCp_dm'])(T, P) def _calc_V(self, T, P, mol=None, V=None): return (self._methods['_calc_V'])(T, P) def _calc_dV_dT(self, T, P, mol=None, V=None): return (self._methods['_calc_dV_dT'])(T, P) def _calc_d2V_dmdT(self, T, P, mol=None, V=None): return (self._methods['_calc_d2V_dmdT'])(T, P) def _calc_d2V_dmdP(self, T, P, mol=None, V=None): return (self._methods['_calc_d2V_dmdP'])(T, P) def _calc_dV_dP(self, T, P, mol=None, V=None): return (self._methods['_calc_dV_dP'])(T, P) def _calc_dV_dm(self, T, P, mol=None, V=None): return (self._methods['_calc_dV_dm'])(T, P) def _calc_d2V_dT2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2V_dT2'])(T, P) def _calc_d2V_dTdP(self, T, P, mol=None, V=None): return (self._methods['_calc_d2V_dTdP'])(T, P) def _calc_d2V_dP2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2V_dP2'])(T, P) def _calc_d2V_dm2(self, T, P, mol=None, V=None): return (self._methods['_calc_d2V_dm2'])(T, P) def _calc_density(self, T, P, mol=None, V=None): return (self._methods['_calc_density'])(T, P) def _calc_alpha(self, T, P, mol=None, V=None): return (self._methods['_calc_alpha'])(T, P) def _calc_beta(self, T, P, mol=None, V=None): return (self._methods['_calc_beta'])(T, P) def _calc_K(self, T, P, mol=None, V=None): return (self._methods['_calc_K'])(T, P) def _calc_Kp(self, T, P, mol=None, V=None): return (self._methods['_calc_Kp'])(T, P) # Compositional quantities def _calc_mu(self, T, P, mol=None, V=None, endmember=None, species=False): return (self._methods['_calc_G'])(T, P) def _calc_a(self, T, P, mol=None, V=None, endmember=None): return (self._methods['_calc_a'])(T, P) def _calc_da_dm(self, T, P, mol=None, V=None, endmember=None): return (self._methods['_calc_da_dm'])(T, P)
#===================================================
[docs]class SolutionPhase(Phase): """ Solid solution phases. Implements the Phase interface. Parameters ---------- phase_classnm : str Official class name for phase (implemented in src code). abbrev : str Official abbreviation of phase (regardless of implementation). calib : bool, default True Indicates whether sample phase should be calibration ready. Attributes ---------- Attributes for this class are inherited from the Phase class. Notes ----- * This code is highly dependent on implementation and is likely to change dramatically with changes in the underlying code that calculates phase properties. * The pure water phase, "H2O", is very complex and thus not available for calibration. The water phase will force its calib flag to False, regardless of input value. """ def __init__(self, phase_classnm, abbrev, calib=True, source='objc', coder_module=None): super().__init__(phase_classnm, abbrev, calib=calib, source=source, coder_module=coder_module) self._init_endmember_props() self._phase_type = 'solution' self._methods = {} if self.source == 'objc': self._methods['_calc_G'] = getattr(self._phase_obj, 'getGibbsFreeEnergyFromMolesOfComponents_andT_andP_') self._methods['_calc_dG_dT'] = getattr(self, 'not_coded') self._methods['_calc_dG_dP'] = getattr(self._phase_obj, 'getVolumeFromMolesOfComponents_andT_andP_') self._methods['_calc_dG_dm'] = getattr(self._phase_obj, 'getDgDmFromMolesOfComponents_andT_andP_') self._methods['_calc_dG_dw'] = getattr(self._phase_obj, 'getDgDwFromMolesOfComponents_andT_andP_' ) if self.calib else getattr(self, 'not_coded') self._methods['_calc_d2G_dT2'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dTdP'] = getattr(self._phase_obj, 'getDvDtFromMolesOfComponents_andT_andP_') self._methods['_calc_d2G_dTdm'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dTdw'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dP2'] = getattr(self._phase_obj, 'getDvDpFromMolesOfComponents_andT_andP_') self._methods['_calc_d2G_dPdm'] = getattr(self._phase_obj, 'getDvDmFromMolesOfComponents_andT_andP_') self._methods['_calc_d2G_dPdw'] = getattr(self, 'not_coded') self._methods['_calc_d2G_dm2'] = getattr(self._phase_obj, 'getD2gDm2FromMolesOfComponents_andT_andP_') self._methods['_calc_d2G_dmdw'] = getattr(self._phase_obj, 'getChemicalPotentialDerivativesForParameter_usingMolesOfComponents_andT_andP_' ) if self.calib else getattr(self, 'not_coded') self._methods['_calc_d3G_dT3'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT2dP'] = getattr(self._phase_obj, 'getD2vDt2FromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dT2dm'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dT2dw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdP2'] = getattr(self._phase_obj, 'getD2vDtDpFromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dTdPdm'] = getattr(self._phase_obj, 'getD2vDmDtFromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dTdPdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdm2'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dTdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dP3'] = getattr(self._phase_obj, 'getD2vDp2FromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dP2dm'] = getattr(self._phase_obj, 'getD2vDmDpFromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dP2dw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dPdm2'] = getattr(self._phase_obj, 'getD2vDm2FromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dm3'] = getattr(self._phase_obj, 'getD3gDm3FromMolesOfComponents_andT_andP_') self._methods['_calc_d3G_dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dP2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dm3dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT3dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT2dPdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdP2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dP3dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dT2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dTdPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dP2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_H'] = getattr(self._phase_obj, 'getEnthalpyFromMolesOfComponents_andT_andP_') self._methods['_calc_S'] = getattr(self._phase_obj, 'getEntropyFromMolesOfComponents_andT_andP_') self._methods['_calc_dS_dm'] = getattr(self._phase_obj, 'getDsDmFromMolesOfComponents_andT_andP_') self._methods['_calc_d2S_dm2'] = getattr(self._phase_obj, 'getD2sDm2FromMolesOfComponents_andT_andP_') self._methods['_calc_V'] = getattr(self._phase_obj, 'getVolumeFromMolesOfComponents_andT_andP_') self._methods['_calc_dV_dT'] = getattr(self._phase_obj, 'getDvDtFromMolesOfComponents_andT_andP_') self._methods['_calc_dV_dP'] = getattr(self._phase_obj, 'getDvDpFromMolesOfComponents_andT_andP_') self._methods['_calc_dV_dm'] = getattr(self._phase_obj, 'getDvDmFromMolesOfComponents_andT_andP_') self._methods['_calc_d2V_dT2'] = getattr(self._phase_obj, 'getD2vDt2FromMolesOfComponents_andT_andP_') self._methods['_calc_d2V_dTdP'] = getattr(self._phase_obj, 'getD2vDtDpFromMolesOfComponents_andT_andP_') self._methods['_calc_d2V_dP2'] = getattr(self._phase_obj, 'getD2vDp2FromMolesOfComponents_andT_andP_') self._methods['_calc_d2V_dmdT'] = getattr(self._phase_obj, 'getD2vDmDtFromMolesOfComponents_andT_andP_') self._methods['_calc_d2V_dmdP'] = getattr(self._phase_obj, 'getD2vDmDpFromMolesOfComponents_andT_andP_') self._methods['_calc_d2V_dm2'] = getattr(self._phase_obj, 'getD2vDm2FromMolesOfComponents_andT_andP_') self._methods['_calc_Cv'] = getattr(self, 'not_coded') self._methods['_calc_Cp'] = getattr(self._phase_obj, 'getHeatCapacityFromMolesOfComponents_andT_andP_') self._methods['_calc_dCp_dT'] = getattr(self._phase_obj, 'getDcpDtFromMolesOfComponents_andT_andP_') self._methods['_calc_dCp_dm'] = getattr(self._phase_obj, 'getDCpDmFromMolesOfComponents_andT_andP_') self._methods['_calc_density'] = getattr(self, 'not_coded') self._methods['_calc_alpha'] = getattr(self, 'not_coded') self._methods['_calc_beta'] = getattr(self, 'not_coded') self._methods['_calc_K'] = getattr(self, 'not_coded') self._methods['_calc_Kp'] = getattr(self, 'not_coded') self._methods['_calc_mu'] = getattr(self._phase_obj, 'getChemicalPotentialFromMolesOfComponents_andT_andP_') self._methods['_calc_a'] = getattr(self._phase_obj, 'getActivityFromMolesOfComponents_andT_andP_') self._methods['_calc_da_dm'] = getattr(self._phase_obj, 'getDaDmFromMolesOfComponents_andT_andP_') else: phase_calibnm = phase_classnm.replace('calib', 'dparam') if ( self.calib) else '' self._methods['_calc_G'] = getattr(self.module, phase_classnm+'g') self._methods['_calc_dG_dT'] = getattr(self.module, phase_classnm+'dgdt') self._methods['_calc_dG_dP'] = getattr(self.module, phase_classnm+'dgdp') self._methods['_calc_dG_dm'] = getattr(self.module, phase_classnm+'dgdn') self._methods['_calc_dG_dw'] = getattr(self.module, phase_calibnm+'g') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d2G_dT2'] = getattr(self.module, phase_classnm+'d2gdt2') self._methods['_calc_d2G_dTdP'] = getattr(self.module, phase_classnm+'d2gdtdp') self._methods['_calc_d2G_dTdm'] = getattr(self.module, phase_classnm+'d2gdndt') self._methods['_calc_d2G_dTdw'] = getattr(self.module, phase_calibnm+'dgdt') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d2G_dP2'] = getattr(self.module, phase_classnm+'d2gdp2') self._methods['_calc_d2G_dPdm'] = getattr(self.module, phase_classnm+'d2gdndp') self._methods['_calc_d2G_dPdw'] = getattr(self.module, phase_calibnm+'dgdp') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d2G_dm2'] = getattr(self.module, phase_classnm+'d2gdn2') self._methods['_calc_d2G_dmdw'] = getattr(self.module, phase_calibnm+'dgdn') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dT3'] = getattr(self.module, phase_classnm+'d3gdt3') self._methods['_calc_d3G_dT2dP'] = getattr(self.module, phase_classnm+'d3gdt2dp') self._methods['_calc_d3G_dT2dm'] = getattr(self.module, phase_classnm+'d3gdndt2') self._methods['_calc_d3G_dT2dw'] = getattr(self.module, phase_calibnm+'d2gdt2') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dTdP2'] = getattr(self.module, phase_classnm+'d3gdtdp2') self._methods['_calc_d3G_dTdPdm'] = getattr(self.module, phase_classnm+'d3gdndtdp') self._methods['_calc_d3G_dTdPdw'] = getattr(self.module, phase_calibnm+'d2gdtdp') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dTdm2'] = getattr(self.module, phase_classnm+'d3gdn2dt') self._methods['_calc_d3G_dTdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dP3'] = getattr(self.module, phase_classnm+'d3gdp3') self._methods['_calc_d3G_dP2dm'] = getattr(self.module, phase_classnm+'d3gdndp2') self._methods['_calc_d3G_dP2dw'] = getattr(self.module, phase_calibnm+'d2gdp2') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d3G_dPdm2'] = getattr(self.module, phase_classnm+'d3gdn2dp') self._methods['_calc_d3G_dPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d3G_dm3'] = getattr(self.module, phase_classnm+'d3gdn3') self._methods['_calc_d3G_dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dT2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdPdmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dP2dmdw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dTdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d4G_dPdm2dw'] = getattr(self, 'not_coded') # d4gdndt3, d4gdndt2dp, d4gdndtdp2, d4gdndp3 are defined in code # d4gdn2dt2, d4gdn2dtdp, d4gdn2dp2 are defined in code # d4gdn3dt, d4gdn3dp are defined in code self._methods['_calc_d4G_dT3dw'] = getattr(self.module, phase_calibnm+'d3gdt3') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dT2dPdw'] = getattr(self.module, phase_calibnm+'d3gdt2dp') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dTdP2dw'] = getattr(self.module, phase_calibnm+'d3gdtdp2') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dP3dw'] = getattr(self.module, phase_calibnm+'d3gdp3') if self.calib else getattr( self, 'not_coded') self._methods['_calc_d4G_dm3dw'] = getattr(self, 'not_coded') # d5gdn2dt3, d5gdn2dt2dp, d5gdn2dtdp2, d5gdn2dp3 are defined in code # d5gdn3dt2, d5gdn3dtdp, d5gdn3dp2 are defined in code self._methods['_calc_d5G_dT2dm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dTdPdm2dw'] = getattr(self, 'not_coded') self._methods['_calc_d5G_dP2dm2dw'] = getattr(self, 'not_coded') # d6gdn3dt3, d6gdn3dt2dp, d6gdn3dtdp2, d6gdn3dp3 are defined in code self._methods['_calc_H'] = getattr(self, 'not_coded') self._methods['_calc_S'] = getattr(self.module, phase_classnm+'s') self._methods['_calc_dS_dm'] = getattr(self, 'not_coded') self._methods['_calc_d2S_dm2'] = getattr(self, 'not_coded') self._methods['_calc_V'] = getattr(self.module, phase_classnm+'v') self._methods['_calc_dV_dT'] = getattr(self.module, phase_classnm+'d2gdtdp') self._methods['_calc_dV_dP'] = getattr(self.module, phase_classnm+'d2gdp2') self._methods['_calc_dV_dm'] = getattr(self.module, phase_classnm+'d2gdndp') self._methods['_calc_d2V_dT2'] = getattr(self.module, phase_classnm+'d3gdt2dp') self._methods['_calc_d2V_dTdP'] = getattr(self.module, phase_classnm+'d3gdtdp2') self._methods['_calc_d2V_dP2'] = getattr(self.module, phase_classnm+'d3gdp3') self._methods['_calc_d2V_dmdT'] = getattr(self.module, phase_classnm+'d3gdndtdp') self._methods['_calc_d2V_dmdP'] = getattr(self.module, phase_classnm+'d3gdndp2') self._methods['_calc_d2V_dm2'] = getattr(self.module, phase_classnm+'d3gdn2dp') self._methods['_calc_Cv'] = getattr(self.module, phase_classnm+'cv') self._methods['_calc_Cp'] = getattr(self.module, phase_classnm+'cp') self._methods['_calc_dCp_dT'] = getattr(self.module, phase_classnm+'dcpdt') self._methods['_calc_dCp_dm'] = getattr(self, 'not_coded') self._methods['_calc_density'] = getattr(self, 'not_coded') self._methods['_calc_alpha'] = getattr(self.module, phase_classnm+'alpha') self._methods['_calc_beta'] = getattr(self.module, phase_classnm+'beta') self._methods['_calc_K'] = getattr(self.module, phase_classnm+'K') self._methods['_calc_Kp'] = getattr(self.module, phase_classnm+'Kp') self._methods['_calc_mu'] = getattr(self.module, phase_classnm+'dgdn') self._methods['_calc_a'] = getattr(self, 'not_coded') self._methods['_calc_da_dm'] = getattr(self, 'not_coded') def _init_endmember_props(self): phase_obj = self.phase_obj props = self.props if self.source == 'objc': try: result = self.phase_obj.supportsParameterCalibration() self._calib = True except: self._calib = False if self.source == 'objc': try: endmember_num = phase_obj.numberOfSolutionComponents() oxide_num = phase_obj.numberOfSolutionSpecies() species_num = phase_obj.numberOfSolutionSpecies() except: endmember_num = 0 oxide_num = 0 species_num = 0 else: phase_classnm = self._phase_cls method = getattr(self.module, phase_classnm+'endmember_number') endmember_num = method() oxide_num = method() method = getattr(self.module, phase_classnm+'species_number') species_num = method() print ('Solution phase code generated by the coder module ' + 'does not yet provide information on solution species. ' + 'Species are proxied by components.') names = [] molwt = [] mol_oxide_comp = [] element_comp = [] formula = [] atom_num = [] species = [] species_elms = [] if self.source == 'objc': try: for i in range(0, endmember_num): iendmember = phase_obj.componentAtIndex_(i) iendmember_name = str(iendmember.phaseName) iendmember_formula = iendmember.phaseFormula imolwt = iendmember.mw ielems = core.double_vector_to_array( iendmember.formulaAsElementArray) iatom_num = np.sum(ielems) imol_oxide_comp = chem.calc_mol_oxide_comp(ielems) names.append(iendmember_name) molwt.append(imolwt) element_comp.append(ielems) mol_oxide_comp.append(imol_oxide_comp) formula.append(iendmember_formula) atom_num.append(iatom_num) except: print('Error loading endmember properties for ' + self._props['phase_name']) else: for i in range(0, endmember_num): method = getattr(self.module, phase_classnm+'endmember_name') iendmember_name = method(i) method = getattr(self.module, phase_classnm+'endmember_formula') iendmember_formula = method(i) method = getattr(self.module, phase_classnm+'endmember_mw') imolwt = method(i) method = getattr(self.module, phase_classnm+'endmember_elements') ielems = method(i) iatom_num = np.sum(ielems) imol_oxide_comp = chem.calc_mol_oxide_comp(ielems) names.append(iendmember_name) molwt.append(imolwt) element_comp.append(ielems) mol_oxide_comp.append(imol_oxide_comp) formula.append(iendmember_formula) atom_num.append(iatom_num) if self.source == 'objc': try: for i in range(0,species_num): ispecies_name = phase_obj.nameOfSolutionSpeciesAtIndex_(i) ispecies_elms_arr = ( phase_obj.elementalCompositionOfSpeciesAtIndex_(i)) ispecies_elms = core.double_vector_to_array( ispecies_elms_arr) species.append(ispecies_name) species_elms.append(ispecies_elms) except: print('Error loading species properties for ' + self._props['phase_name']) else: for i in range(0,species_num): method = getattr(self.module, phase_classnm+'species_name') ispecies_name = method(i) method = getattr(self.module, phase_classnm+'species_elements') ispecies_elms = method(i) species.append(ispecies_name) species_elms.append(ispecies_elms) print ('Solution phase code generated by the coder module ' + 'does not yet provide information on species properties. ' + 'Species are proxied by components.') names = np.array(names) molwt = np.array(molwt) element_comp = np.array(element_comp) mol_oxide_comp = np.array(mol_oxide_comp) formula = np.array(formula) atom_num = np.array(atom_num) oxide_space = np.any(mol_oxide_comp, axis=0) species = np.array(species) species_elms = np.array(species_elms) # endmember_props = OrderedDict() props['endmember_name'] = names props['endmember_num'] = endmember_num props['endmember_id'] = np.arange(endmember_num) props['species_num'] = species_num props['formula'] = formula props['atom_num'] = atom_num props['molwt'] = molwt props['elemental_entropy'] = [None] props['element_comp'] = element_comp props['mol_oxide_comp'] = mol_oxide_comp props['oxide_space'] = oxide_space props['species_name'] = species props['species_elms'] = species_elms # self._endmember_props = endmember_props # methods:
[docs] def calc_endmember_comp(self, mol_oxide_comp, method='least_squares', output_residual=False, normalize=False, decimals=10): """ Get fraction of each endmember given the composition. Parameters ---------- mol_oxide_comp : double array Amounts of each oxide in standard order (defined in OXIDES) decimals : int, default 10 Number of decimals to round result to method : str, default 'least_squares' Method used to convert oxide composition (in moles) to moles of endmembers 'intrinsic' is alternate method, hardcoded by the solution implementation Returns ------- endmember_comp : double array Best-fit molar composition in terms of endmembers mol_oxide_comp_residual : double array Residual molar oxide composition Notes ----- * Eventually, we may want the ability to calculate endmember comp. using a variety of methods for inputing composition: * kind : ['wt_oxide', 'mol_oxide', 'element'] Identifies how composition is defined. """ assert mol_oxide_comp.size == chem.oxide_props['oxide_num'], ( 'Oxide composition array not standard length') if method == 'least_squares': mol_oxide_comp_endmembers = self.props['mol_oxide_comp'] # NOTE: mol oxide composition is extensive (in absolute mols) # need not sum to 1. output = np.linalg.lstsq( mol_oxide_comp_endmembers.T, mol_oxide_comp, rcond=None) endmember_comp = output[0] endmember_comp = np.round(endmember_comp, decimals=decimals) mol_oxide_comp_model = np.dot(mol_oxide_comp_endmembers.T, endmember_comp) mol_oxide_comp_residual = (mol_oxide_comp - mol_oxide_comp_model) elif method == 'intrinsic': elems = chem.PERIODIC_ORDER monovalent_cations = chem.oxide_props['cations'] monovalent_elem_ind = np.array( [np.where(elems==icat)[0][0] for icat in monovalent_cations]) oxy_elem_ind = np.where(elems=='O')[0][0] moles_elm = np.zeros(chem.PERIODIC_ORDER.size) for i in range(0, mol_oxide_comp.size): ind = monovalent_elem_ind[i] nCat = chem.oxide_props['cat_num'][i] nOx = chem.oxide_props['oxy_num'][i] moles_elm[ind] += mol_oxide_comp[i]*nCat moles_elm[oxy_elem_ind] += mol_oxide_comp[i]*nOx if self.source == 'objc': moles_pot_arr = self._phase_obj.convertElementsToMoles_( core.array_to_ctype_array(moles_elm)) endmember_comp = core.double_vector_to_array(moles_pot_arr) else: method = getattr(self.module, self._phase_cls+'conv_elm_to_moles') endmember_comp = method(moles_elm) mol_oxide_comp_residual = None else: assert False, ( 'Method argument is not valid. Choose either "least_squares" ' + 'or "intrinic".') if normalize: mol_oxide_tot = np.sum(endmember_comp) endmember_comp /= mol_oxide_tot if mol_oxide_comp_residual is not None: mol_oxide_comp_residual /= mol_oxide_tot if output_residual: return endmember_comp, mol_oxide_comp_residual else: return endmember_comp
[docs] def test_endmember_comp(self, mol_comp): """ Tests validity of endmember component moles array Parameters ---------- mol_comp : double array Mole numbers of each component in the solution Returns ------- flag : boolean True is composition is valid, otherwise False. """ if self.source == 'objc': result = self._phase_obj.testPermissibleValuesOfComponents_( core.array_to_ctype_array(mol_comp)) else: method = getattr(self.module, self._phase_cls+'test_moles') result = method(mol_comp) return (True if result else False)
[docs] def compute_formula(self, T, P, mol_comp): """ Converts an input array of moles of endmember components to the chemical formula of the phase Parameters ---------- T : double Temperature in Kelvins P : double Pressure in bars mol_comp : double array Mole numbers of each component in the solution Returns ------- formula : str A string with the formula of the phase """ if self.source == 'objc': return self._phase_obj.getFormulaFromMolesOfComponents_andT_andP_( core.array_to_ctype_array(mol_comp), T, P) else: phase_classnm = self._phase_cls method = getattr(self.module, phase_classnm+'formula') result = method(T, P, mol_comp) return result
[docs] def covert_endmember_comp(self, mol_comp, output='total_moles'): """ Converts an input array of moles of endmember components to the specified quantity Parameters ---------- mol_comp : double array Mole numbers of each component in the solution output : str, default = 'total_moles' Output quantity: - 'total_moles' - double - 'moles_elements' - double array (standard order and length) - 'mole_fraction' - double array (same order and length as input) Returns ------- result : double or double array The computed quantity as double or double array """ if self.source == 'objc': if output == 'total_moles': result = self._phase_obj.totalMolesFromMolesOfComponents_( core.array_to_ctype_array(mol_comp)) return result elif output == 'moles_elements': result = self._phase_obj.convertMolesToElements_( core.array_to_ctype_array(mol_comp)) return core.double_vector_to_array(result) elif output == 'mole_fraction': result = self._phase_obj.convertMolesToMoleFractions_( core.array_to_ctype_array(mol_comp)) return core.double_vector_to_array(result) else: print ("output option must be 'total_moles', " + "'moles_elements', or 'mole_fraction'.") return np.empty(self.props['endmember_num']) else: phase_classnm = self._phase_cls if output == 'total_moles': method = getattr(self.module, self._phase_cls+'conv_moles_to_tot_moles') return method(mol_comp) elif output == 'moles_elements': method = getattr(self.module, self._phase_cls+'conv_moles_to_elm') return method(mol_comp) elif output == 'mole_fraction': method = getattr(self.module, self._phase_cls+'conv_moles_to_mole_frac') return method(mol_comp) else: print ("output option must be 'total_moles', " + "'moles_elements', or 'mole_fraction'.") return np.empty(self.props['endmember_num'])
[docs] def convert_species_to_comp(self, mol_species): """ Converts an input array of moles of species to moles of endmember components Parameters ---------- mol_species : double array Mole numbers of each species in the solution Returns ------- result : double array Moles of endmember components """ if self.source == 'objc': result_arr = self._phase_obj.convertMolesOfSpeciesToMolesOfComponents_( core.array_to_ctype_array(mol_species)) return core.double_vector_to_array(result_arr) else: print ('Solution phase code generated by the coder module ' + 'does not yet provide a conversion method.') return np.empty(self.props['endmember_num'])
[docs] def convert_elements(self, mol_elm, output='moles_end'): """ Converts an array of mole numbers of elements (in the standard order) to the specified output quantity Parameters ---------- mole_elm : double array Mole numbers of elements in the standard order output : str, default = 'moles_end' Output quantity: - 'moles_end' - double array of moles of endmembers - 'total_moles' - double, sum of moles of endmembers - 'total_grams' - double, sum of grams of solution Returns ------- result : double or double array The computed quantity as double or double array """ if self.source == 'objc': if output == 'moles_end': result = self._phase_obj.convertElementsToMoles_( core.array_to_ctype_array(mol_elm)) return core.double_vector_to_array(result) elif output == 'total_moles': result = self._phase_obj.convertElementsToTotalMoles_( core.array_to_ctype_array(mol_elm)) return result elif output == 'total_grams': result = self._phase_obj.convertElementsToTotalMass_( core.array_to_ctype_array(mol_elm)) return result else: print ("output option must be 'moles_end', 'total_moles', " + "or 'total_grams'.") return None else: phase_classnm = self._phase_cls if output == 'moles_end': method = getattr(self.module, phase_classnm+'conv_elm_to_moles') return method(mol_elm) elif output == 'total_moles': method = getattr(self.module, phase_classnm+'conv_elm_to_tot_moles') return method(mol_elm) elif output == 'total_grams': method = getattr(self.module, phase_classnm+'conv_elm_to_tot_grams') return method(mol_elm) else: print ("output option must be 'moles_end', 'total_moles', " + "or 'total_grams'.") return None
[docs] def get_endmember_ind(self, mol_oxide_comp, get_endmember_comp=False, TOL=1e-6): """ Get index of endmember that best matches composition. Parameters ---------- mol_oxide_comp : double array Amounts of each oxide in standard order (defined in OXIDES) TOL : double, default 1e-6 Allowed tolerance for mismatch to defined composition get_endmember_comp : bool, default False If true, also return endmember composition array. Returns ------- endmember_ind : int Index of best-fit endmember endmember_comp : double array, optional Composition array in terms of endmembers. Return if get_endmember_comp==True """ endmember_comp, mol_oxide_comp_residual = ( self.calc_endmember_comp(mol_oxide_comp)) endmember_ind = np.argmax(endmember_comp) endmember_scl_comp = endmember_comp/endmember_comp[endmember_ind] ind_all = np.arange(endmember_comp.size) ind_else = np.delete(ind_all, endmember_ind) assert np.all(np.abs( endmember_scl_comp[ind_else] < TOL)), ( 'Provided composition does not match ' 'an endmember.' ) if get_endmember_comp: return endmember_ind, endmember_comp else: return endmember_ind
def not_coded(self, T, P, mol=None, V=None, param=None): raise AttributeError( 'Function not implemented for this database model.') def _calc_G(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_G'])(cmol, T, P) else: result = (self._methods['_calc_G'])(T, P, mol) return result def _calc_dG_dT(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_dG_dT'])(cmol, T, P) else: result = (self._methods['_calc_dG_dT'])(T, P, mol) return result def _calc_dG_dP(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_dG_dP'])(cmol, T, P) else: result = (self._methods['_calc_dG_dP'])(T, P, mol) return result def _calc_dG_dm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_dG_dm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_dG_dm'])(T, P, mol) return result def _calc_dG_dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_dG_dw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_dG_dw'])(T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d2G_dT2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d2G_dT2'])(cmol, T, P) else: result = (self._methods['_calc_d2G_dT2'])(T, P, mol) return result def _calc_d2G_dTdP(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d2G_dTdP'])(cmol, T, P) else: result = (self._methods['_calc_d2G_dTdP'])(T, P, mol) return result def _calc_d2G_dTdm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2G_dTdm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d2G_dTdm'])(T, P, mol) return result def _calc_d2G_dTdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d2G_dTdw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d2G_dTdw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d2G_dP2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d2G_dP2'])(cmol, T, P) else: result = (self._methods['_calc_d2G_dP2'])(T, P, mol) return result def _calc_d2G_dPdm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2G_dPdm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d2G_dPdm'])(T, P, mol) return result def _calc_d2G_dPdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d2G_dPdw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d2G_dPdw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d2G_dm2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2G_dm2'])(cmol, T, P) result = core.double_matrix_to_array(iresult_mat) else: result = (self._methods['_calc_d2G_dm2'])(T, P, mol) return result def _calc_d2G_dmdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': name = self._param_props['param_names'][iparam] cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d2G_dmdw'])( name, cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) else: iresult = (self._methods['_calc_d2G_dmdw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dT3(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d3G_dT3'])(cmol, T, P) else: result = (self._methods['_calc_d3G_dT3'])(T, P, mol) return result def _calc_d3G_dT2dP(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d3G_dT2dP'])(cmol, T, P) else: result = (self._methods['_calc_d3G_dT2dP'])(T, P, mol) return result def _calc_d3G_dT2dm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dT2dm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d3G_dT2dm'])(T, P, mol) return result def _calc_d3G_dT2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d3G_dT2dw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d3G_dT2dw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dTdP2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d3G_dTdP2'])(cmol, T, P) else: result = (self._methods['_calc_d3G_dTdP2'])(T, P, mol) return result def _calc_d3G_dTdPdm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dTdPdm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d3G_dTdPdm'])(T, P, mol) return result def _calc_d3G_dTdPdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d3G_dTdPdw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d3G_dTdPdw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dTdm2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dTdm2'])(cmol, T, P) result = core.double_matrix_to_array(iresult_mat) else: result = (self._methods['_calc_d3G_dTdm2'])(T, P, mol) return result def _calc_d3G_dTdmdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d3G_dTdmdw'])( iparam, cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) else: iresult = (self._methods['_calc_d3G_dTdmdw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dP3(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d3G_dP3'])(cmol, T, P) else: result = (self._methods['_calc_d3G_dP3'])(T, P, mol) return result def _calc_d3G_dP2dm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dP2dm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d3G_dP2dm'])(T, P, mol) return result def _calc_d3G_dP2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d3G_dP2dw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d3G_dP2dw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dPdm2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dPdm2'])(cmol, T, P) result = core.double_matrix_to_array(iresult_mat) else: result = (self._methods['_calc_d3G_dPdm2'])(T, P, mol) return result def _calc_d3G_dPdmdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d3G_dPdmdw'])( iparam, cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) else: iresult = (self._methods['_calc_d3G_dPdmdw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d3G_dm3(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dm3'])(cmol, T, P) result = core.double_tensor_to_array(iresult_mat) else: result = (self._methods['_calc_d3G_dm3'])(T, P, mol) return result def _calc_d3G_dm2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d3G_dm2dw'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) else: iresult = (self._methods['_calc_d3G_dm2dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dT3dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dT3dw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dT3dw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dT2dPdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dT2dPdw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dT2dPdw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dT2dmdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dT2dmdw'])( iparam, cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) else: iresult = (self._methods['_calc_d4G_dT2dmdw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dTdP2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dTdP2dw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dTdP2dw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dTdPdmdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dTdPdmdw'])( iparam, cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) else: iresult = (self._methods['_calc_d4G_dTdPdmdw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dTdm2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d4G_dTdm2dw'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) else: iresult = (self._methods['_calc_d4G_dTdm2dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dP3dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dP3dw'])(cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) for iparam in param: if self.source == 'objc': iresult_term = iresult[iparam] else: iresult_term = (self._methods['_calc_d4G_dP3dw'])( T, P, mol, iparam) result.append(iresult_term) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dP2dmdw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_arr = (self._methods['_calc_d4G_dP2dmdw'])( iparam, cmol, T, P) iresult = core.double_vector_to_array(iresult_arr) else: iresult = (self._methods['_calc_d4G_dP2dmdw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dPdm2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d4G_dPdm2dw'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) else: iresult = (self._methods['_calc_d4G_dPdm2dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d4G_dm3dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d4G_dm3dw'])(cmol, T, P) iresult = core.double_tensor_to_array(iresult_mat) for iparam in param: if self.source == 'objc': iresult = iresult[iparam] else: iresult = (self._methods['_calc_d4G_dm3dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d5G_dT2dm2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d5G_dT2dm2dw'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) else: iresult = (self._methods['_calc_d5G_dT2dm2dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d5G_dTdPdm2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d5G_dTdPdm2dw'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) else: iresult = (self._methods['_calc_d5G_dTdPdm2dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_d5G_dP2dm2dw(self, T, P, mol=None, param=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) result = [] for iparam in param: if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d5G_dP2dm2dw'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) else: iresult = (self._methods['_calc_d5G_dP2dm2dw'])(T, P, mol, iparam) result.append(iresult) result = np.array(result) if len(result) > 1 else result[0] return result def _calc_H(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_H'])(cmol, T, P) else: result = (self._methods['_calc_G'])(T, P, mol) + T*( self._methods['_calc_S'])(T, P, mol) return result def _calc_S(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_S'])(cmol, T, P) else: result = (self._methods['_calc_S'])(T, P, mol) return result def _calc_V(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_V'])(cmol, T, P) else: result = (self._methods['_calc_V'])(T, P, mol) return result def _calc_mu(self, T, P, mol=None, V=None, endmember=None, species=False): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': chem_pot = [] cmol = core.array_to_ctype_array(mol) for iendmem in endmember: if species: ichem_pot_arr = ( self._phase_obj. chemicalPotentialsOfSpeciesFromMolesOfComponents_andT_andP_( cmol, T, P) ) else: ichem_pot_arr = (self._methods['_calc_mu'])(cmol, T, P) ichem_pot = core.double_vector_to_array(ichem_pot_arr) if iendmem >= 0: ichem_pot = ichem_pot[int(iendmem)] chem_pot.append(ichem_pot) chem_pot = np.array(chem_pot) if endmember[0] == -1 else chem_pot[0] else: chem_pot = (self._methods['_calc_mu'])(T, P, mol) if endmember is not None: result = [] for iendmem in endmember: result.append(chem_pot[iendmem]) chem_pot = np.array(result) return chem_pot def _calc_a(self, T, P, mol=None, V=None, endmember=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': chem_pot = [] cmol = core.array_to_ctype_array(mol) for iendmem in endmember: ichem_pot_arr = ((self._methods['_calc_a'])(cmol, T, P)) ichem_pot = core.double_vector_to_array(ichem_pot_arr) if iendmem >= 0: ichem_pot = ichem_pot[int(iendmem)] chem_pot.append(ichem_pot) chem_pot = np.array(chem_pot) if endmember[0] == -1 else chem_pot[0] else: chem_pot = (self._methods['_calc_a'])(T, P, mol) return chem_pot def _calc_dS_dm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_dS_dm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = -(self._methods['_calc_d2G_dTdm'])(T, P, mol) return result def _calc_dV_dm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_dV_dm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_dV_dm'])(T, P, mol) return result def _calc_da_dm(self, T, P, mol=None, V=None, endmember=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_da_dm'])(cmol, T, P) iresult = core.double_matrix_to_array(iresult_mat) if endmember >= 0: result = [] for iendmem in endmember: iresult = iresult[[int(iendmem)], :][0] result.append(iresult) result = result[0] else: result = iresult else: result = (self._methods['_calc_da_dm'])(T, P, mol) return result def _calc_d2S_dm2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2S_dm2'])(cmol, T, P) result = core.double_matrix_to_array(iresult_mat) else: result = -(self._methods['_calc_d3G_dTdm2'])(T, P, mol) return result def _calc_d2V_dm2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2V_dm2'])(cmol, T, P) result = core.double_matrix_to_array(iresult_mat) else: result = (self._methods['_calc_d2V_dm2'])(T, P, mol) return result def _calc_Cv(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) return (self._methods['_calc_Cv'])(T, P, mol) def _calc_Cp(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for ' 'SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_Cp'])(cmol, T, P) else: result = (self._methods['_calc_Cp'])(T, P, mol) return result def _calc_dV_dT(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_dV_dT'])(cmol, T, P) else: result = (self._methods['_calc_dV_dT'])(T, P, mol) return result def _calc_dV_dP(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_dV_dP'])(cmol, T, P) else: result = (self._methods['_calc_dV_dP'])(T, P, mol) return result def _calc_dCp_dT(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_dCp_dT'])(cmol, T, P) else: result = (self._methods['_calc_dCp_dT'])(T, P, mol) return result def _calc_dCp_dm(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_dCp_dm'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = -(self._methods['_calc_d3G_dT2dm'])(T, P, mol)/T return result def _calc_d2V_dT2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d2V_dT2'])(cmol, T, P) else: result = (self._methods['_calc_d2V_dT2'])(T, P, mol) return result def _calc_d2V_dTdP(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d2V_dTdP'])(cmol, T, P) else: result = (self._methods['_calc_d2V_dTdP'])(T, P, mol) return result def _calc_d2V_dP2(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_d2V_dP2'])(cmol, T, P) else: result = (self._methods['_calc_d2V_dP2'])(T, P, mol) return result def _calc_d2V_dmdT(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2V_dmdT'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d2V_dmdT'])(T, P, mol) return result def _calc_d2V_dmdP(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) iresult_mat = (self._methods['_calc_d2V_dmdP'])(cmol, T, P) result = core.double_vector_to_array(iresult_mat) else: result = (self._methods['_calc_d2V_dmdP'])(T, P, mol) return result def _calc_density(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_density'])(cmol, T, P) else: result = (self._methods['_calc_density'])(T, P, mol) return result def _calc_alpha(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_alpha'])(cmol, T, P) else: result = (self._methods['_calc_alpha'])(T, P, mol) return result def _calc_beta(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_beta'])(cmol, T, P) else: result = (self._methods['_calc_beta'])(T, P, mol) return result def _calc_K(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_K'])(cmol, T, P) else: result = (self._methods['_calc_K'])(T, P, mol) return result def _calc_Kp(self, T, P, mol=None, V=None): assert mol is not None, ( 'Molar composition "mol" must be defined for SolutionPhase.' ) if self.source == 'objc': cmol = core.array_to_ctype_array(mol) result = (self._methods['_calc_Kp'])(cmol, T, P) else: result = (self._methods['_calc_Kp'])(T, P, mol) return result
#=================================================== # class Sample: # """ # Sample Phase with specific properties. # # Parameters # ---------- # phase_classnm : str # Official class name for phase (implemented in src code). # abbrev : str # Official abbreviation of phase (regardless of implementation). # mol : array, optional # Molar composition vector describing amount of each endmember. This # quantity is only meaningful for solution phases, since pure phases have # fixed composition (default is an empty list). # kind : {'Pure','Solution','Aqueous','Gas'} # Defines kind of sample phase. # calib : bool, default True # Indicates whether sample phase should be calibration ready. # # # Attributes # ---------- # mol # kind # phase : Phase object # Stores implementation of phase, including all property parameter values. # # """ # # def __init__(self, phase_classnm, abbrev, mol=None, # kind='Pure', calib=True): # # self.mol = mol # self.kind = kind # # if kind == 'Pure': # phase = PurePhase(phase_classnm, abbrev, calib=calib) # elif kind == 'Solution': # phase = SolutionPhase(phase_classnm, abbrev, calib=calib) # elif kind == 'Aqueous': # phase = AqueousPhase( phase_classnm, abbrev, calib=calib ) # elif kind == 'Gas': # phase = GasPhase( phase_classnm, abbrev, calib=calib ) # else: # raise NotImplemented('Sample kind "'+kind+'" is not yet implimented') # # self.phase = phase # # pass #=================================================== # def get_dG_dw( self, mol_a, T, P ): # assert False, ' # mol_vec = core.array_to_double_vector(mol_a) # print(mol_vec) # dG_dw_vec = self._phase_obj.getDgDwFromMolesOfComponents_andT_andP_( mol_vec, T, P ) # print(dG_dw_vec) # dG_dw_a = core.double_vector_to_array( dG_dw_vec ) # return dG_dw_a # def get_dchempot_dw # getChemicalPotentialDerivativesForParameterArray:(NSArray *)array usingMolesOfComponents:(double *)m andT:(double)t andP: # def get_dchempot_dw( self, T_a, P_a ): # self._phase_obj.getChemicalPotentialDerivativesForParameter_usingMolesOfComponents_andT_andP_('S',m,T,P)