"""
The equilibrium module defines routines for interacting with
calculated phase equilibria.
"""
import warnings
from collections import OrderedDict
from collections.abc import Mapping
from datetime import datetime
import pycalphad.variables as v
from pycalphad.core.utils import unpack_components, unpack_condition, unpack_phases, filter_phases, instantiate_models, get_state_variables
from pycalphad import calculate
from pycalphad.core.errors import EquilibriumError, ConditionError
from pycalphad.core.starting_point import starting_point
from pycalphad.codegen.callables import build_phase_records
from pycalphad.core.eqsolver import _solve_eq_at_conditions
from pycalphad.core.phase_rec import PhaseRecord
from pycalphad.core.solver import Solver
from pycalphad.core.light_dataset import LightDataset
from pycalphad.model import Model
import numpy as np
def _adjust_conditions(conds):
"Adjust conditions values to be within the numerical limit of the solver."
new_conds = OrderedDict()
minimum_composition = 1e-10
for key, value in sorted(conds.items(), key=str):
if key == str(key):
key = getattr(v, key, key)
if isinstance(key, v.MoleFraction):
vals = unpack_condition(value)
# "Zero" composition is a common pattern. Do not warn for that case.
if np.any(np.logical_and(np.asarray(vals) < minimum_composition, np.asarray(vals) > 0)):
warnings.warn(
f"Some specified compositions are below the minimum allowed composition of {minimum_composition}.")
new_conds[key] = [max(val, minimum_composition) for val in vals]
else:
new_conds[key] = unpack_condition(value)
return new_conds
def _eqcalculate(dbf, comps, phases, conditions, output, data=None, per_phase=False, callables=None, model=None,
parameters=None, **kwargs):
"""
WARNING: API/calling convention not finalized.
Compute the *equilibrium value* of a property.
This function differs from `calculate` in that it computes
thermodynamic equilibrium instead of randomly sampling the
internal degrees of freedom of a phase.
Because of that, it's slower than `calculate`.
This plugs in the equilibrium phase and site fractions
to compute a thermodynamic property defined in a Model.
Parameters
----------
dbf : Database
Thermodynamic database containing the relevant parameters.
comps : list
Names of components to consider in the calculation.
phases : list or dict
Names of phases to consider in the calculation.
conditions : dict or (list of dict)
StateVariables and their corresponding value.
output : str
Equilibrium model property (e.g., CPM, HM, etc.) to compute.
This must be defined as an attribute in the Model class of each phase.
data : Dataset
Previous result of call to `equilibrium`.
Should contain the equilibrium configurations at the conditions of interest.
If the databases are not the same as in the original calculation,
the results may be meaningless.
per_phase : bool, optional
If True, compute and return the property for each phase present.
If False, return the total system value, weighted by the phase fractions.
callables : dict
Callable functions to compute 'output' for each phase.
model : a dict of phase names to Model
Model class to use for each phase.
parameters : dict, optional
Maps SymEngine Symbol to numbers, for overriding the values of parameters in the Database.
kwargs
Passed to `calculate`.
Returns
-------
Dataset of property as a function of equilibrium conditions
"""
if data is None:
raise ValueError('Required kwarg "data" is not specified')
if model is None:
raise ValueError('Required kwarg "model" is not specified')
active_phases = unpack_phases(phases)
conds = _adjust_conditions(conditions)
indep_vars = ['N', 'P', 'T']
# TODO: Rewrite this to use the coord dict from 'data'
str_conds = OrderedDict((str(key), value) for key, value in conds.items())
indep_vals = list([float(x) for x in np.atleast_1d(val)]
for key, val in str_conds.items() if key in indep_vars)
coord_dict = str_conds.copy()
components = [x for x in sorted(comps)]
desired_active_pure_elements = [list(x.constituents.keys()) for x in components]
desired_active_pure_elements = [el.upper() for constituents in desired_active_pure_elements for el in constituents]
pure_elements = sorted(set([x for x in desired_active_pure_elements if x != 'VA']))
coord_dict['vertex'] = np.arange(len(pure_elements) + 1) # +1 is to accommodate the degenerate degree of freedom at the invariant reactions
grid_shape = np.meshgrid(*coord_dict.values(),
indexing='ij', sparse=False)[0].shape
prop_shape = grid_shape
prop_dims = list(str_conds.keys()) + ['vertex']
result = LightDataset({output: (prop_dims, np.full(prop_shape, np.nan))}, coords=coord_dict)
# For each phase select all conditions where that phase exists
# Perform the appropriate calculation and then write the result back
for phase in active_phases:
dof = len(model[phase].site_fractions)
current_phase_indices = (data.Phase == phase)
if ~np.any(current_phase_indices):
continue
points = data.Y[np.nonzero(current_phase_indices)][..., :dof]
statevar_indices = np.nonzero(current_phase_indices)[:len(indep_vals)]
statevars = {key: np.take(np.asarray(vals), idx)
for key, vals, idx in zip(indep_vars, indep_vals, statevar_indices)}
statevars.update(kwargs)
if statevars.get('mode', None) is None:
statevars['mode'] = 'numpy'
calcres = calculate(dbf, comps, [phase], output=output, points=points, broadcast=False,
callables=callables, parameters=parameters, model=model, **statevars)
result[output][np.nonzero(current_phase_indices)] = calcres[output].values
if not per_phase:
out = np.nansum(result[output] * data['NP'], axis=-1)
dv_output = result.data_vars[output]
result.remove(output)
# remove the vertex coordinate because we summed over it
result.add_variable(output, dv_output[0][:-1], out)
else:
dv_phase = data.data_vars['Phase']
dv_np = data.data_vars['NP']
result.add_variable('Phase', dv_phase[0], dv_phase[1])
result.add_variable('NP', dv_np[0], dv_np[1])
return result
[docs]
def equilibrium(dbf, comps, phases, conditions, output=None, model=None,
verbose=False, broadcast=True, calc_opts=None, to_xarray=True,
scheduler='sync', parameters=None, solver=None, callables=None,
phase_records=None, **kwargs):
"""
Calculate the equilibrium state of a system containing the specified
components and phases, under the specified conditions.
Parameters
----------
dbf : Database
Thermodynamic database containing the relevant parameters.
comps : list
Names of components to consider in the calculation.
phases : list or dict
Names of phases to consider in the calculation.
conditions : dict or (list of dict)
StateVariables and their corresponding value.
output : str or list of str, optional
Additional equilibrium model properties (e.g., CPM, HM, etc.) to compute.
These must be defined as attributes in the Model class of each phase.
model : Model, a dict of phase names to Model, or a seq of both, optional
Model class to use for each phase.
verbose : bool, optional
Print details of calculations. Useful for debugging.
broadcast : bool
If True, broadcast conditions against each other. This will compute all combinations.
If False, each condition should be an equal-length list (or single-valued).
Disabling broadcasting is useful for calculating equilibrium at selected conditions,
when those conditions don't comprise a grid.
calc_opts : dict, optional
Keyword arguments to pass to `calculate`, the energy/property calculation routine.
to_xarray : bool
Whether to return an xarray Dataset (True, default) or an EquilibriumResult.
scheduler : Dask scheduler, optional
Job scheduler for performing the computation.
If None, return a Dask graph of the computation instead of actually doing it.
parameters : dict, optional
Maps SymEngine Symbol to numbers, for overriding the values of parameters in the Database.
solver : pycalphad.core.solver.SolverBase
Instance of a solver that is used to calculate local equilibria.
Defaults to a pycalphad.core.solver.Solver.
callables : dict, optional
Pre-computed callable functions for equilibrium calculation.
phase_records : Optional[Mapping[str, PhaseRecord]]
Mapping of phase names to PhaseRecord objects with `'GM'` output. Must include
all active phases. The `model` argument must be a mapping of phase names to
instances of Model objects.
Returns
-------
Structured equilibrium calculation, or Dask graph if scheduler=None.
Examples
--------
None yet.
"""
if not broadcast:
raise NotImplementedError('Broadcasting cannot yet be disabled')
comps = sorted(unpack_components(dbf, comps))
phases = unpack_phases(phases) or sorted(dbf.phases.keys())
list_of_possible_phases = filter_phases(dbf, comps)
if len(list_of_possible_phases) == 0:
raise ConditionError('There are no phases in the Database that can be active with components {0}'.format(comps))
active_phases = filter_phases(dbf, comps, phases)
if len(active_phases) == 0:
raise ConditionError('None of the passed phases ({0}) are active. List of possible phases: {1}.'.format(phases, list_of_possible_phases))
if isinstance(comps, (str, v.Species)):
comps = [comps]
if len(set(comps) - set(dbf.species)) > 0:
raise EquilibriumError('Components not found in database: {}'
.format(','.join([c.name for c in (set(comps) - set(dbf.species))])))
calc_opts = calc_opts if calc_opts is not None else dict()
solver = solver if solver is not None else Solver(verbose=verbose)
parameters = parameters if parameters is not None else dict()
if isinstance(parameters, dict):
parameters = OrderedDict(sorted(parameters.items(), key=str))
# Temporary solution until constraint system improves
if conditions.get(v.N) is None:
conditions[v.N] = 1
if np.any(np.array(conditions[v.N]) != 1):
raise ConditionError('N!=1 is not yet supported, got N={}'.format(conditions[v.N]))
# Modify conditions values to be within numerical limits, e.g., X(AL)=0
# Also wrap single-valued conditions with lists
conds = _adjust_conditions(conditions)
for cond in conds.keys():
if isinstance(cond, (v.MoleFraction, v.ChemicalPotential)) and cond.species not in comps:
raise ConditionError('{} refers to non-existent component'.format(cond))
str_conds = OrderedDict((str(key), value) for key, value in conds.items())
components = [x for x in sorted(comps)]
desired_active_pure_elements = [list(x.constituents.keys()) for x in components]
desired_active_pure_elements = [el.upper() for constituents in desired_active_pure_elements for el in constituents]
pure_elements = sorted(set([x for x in desired_active_pure_elements if x != 'VA']))
if verbose:
print('Components:', ' '.join([str(x) for x in comps]))
print('Phases:', end=' ')
output = output if output is not None else 'GM'
output = output if isinstance(output, (list, tuple, set)) else [output]
output = set(output)
output |= {'GM'}
output = sorted(output)
if phase_records is None:
models = instantiate_models(dbf, comps, active_phases, model=model, parameters=parameters)
phase_records = build_phase_records(dbf, comps, active_phases, conds, models,
output='GM', callables=callables,
parameters=parameters, verbose=verbose,
build_gradients=True, build_hessians=True)
else:
# phase_records were provided, instantiated models must also be provided by the caller
models = model
if not isinstance(models, Mapping):
raise ValueError("A dictionary of instantiated models must be passed to `equilibrium` with the `model` argument if the `phase_records` argument is used.")
active_phases_without_models = [name for name in active_phases if not isinstance(models.get(name), Model)]
active_phases_without_phase_records = [name for name in active_phases if not isinstance(phase_records.get(name), PhaseRecord)]
if len(active_phases_without_phase_records) > 0:
raise ValueError(f"phase_records must contain a PhaseRecord instance for every active phase. Missing PhaseRecord objects for {sorted(active_phases_without_phase_records)}")
if len(active_phases_without_models) > 0:
raise ValueError(f"model must contain a Model instance for every active phase. Missing Model objects for {sorted(active_phases_without_models)}")
if verbose:
print('[done]', end='\n')
state_variables = sorted(get_state_variables(models=models, conds=conds), key=str)
# 'calculate' accepts conditions through its keyword arguments
grid_opts = calc_opts.copy()
statevar_strings = [str(x) for x in state_variables]
grid_opts.update({key: value for key, value in str_conds.items() if key in statevar_strings})
if 'pdens' not in grid_opts:
grid_opts['pdens'] = 60
grid = calculate(dbf, comps, active_phases, model=models, fake_points=True,
phase_records=phase_records, output='GM', parameters=parameters,
to_xarray=False, **grid_opts)
coord_dict = str_conds.copy()
coord_dict['vertex'] = np.arange(len(pure_elements) + 1) # +1 is to accommodate the degenerate degree of freedom at the invariant reactions
coord_dict['component'] = pure_elements
properties = starting_point(conds, state_variables, phase_records, grid)
properties = _solve_eq_at_conditions(properties, phase_records, grid,
list(str_conds.keys()), state_variables,
verbose, solver=solver)
# Compute equilibrium values of any additional user-specified properties
# We already computed these properties so don't recompute them
output = sorted(set(output) - {'GM', 'MU'})
for out in output:
if (out is None) or (len(out) == 0):
continue
# TODO: How do we know if a specified property should be per_phase or not?
# For now, we make a best guess
if (out == 'degree_of_ordering') or (out == 'DOO'):
per_phase = True
else:
per_phase = False
eqcal = _eqcalculate(dbf, comps, active_phases, conditions, out,
data=properties, per_phase=per_phase, model=models,
callables=callables, parameters=parameters, **calc_opts)
properties = properties.merge(eqcal, inplace=True, compat='equals')
if to_xarray:
properties = properties.get_dataset()
properties.attrs['created'] = datetime.utcnow().isoformat()
if len(kwargs) > 0:
warnings.warn('The following equilibrium keyword arguments were passed, but unused:\n{}'.format(kwargs))
return properties