Source code for pycalphad.core.lower_convex_hull
"""
The lower_convex_hull module handles geometric calculations associated with
equilibrium calculation.
"""
from pycalphad.property_framework.computed_property import LinearCombination
from .hyperplane import hyperplane
from pycalphad.variables import ChemicalPotential, MassFraction, MoleFraction, IndependentPotential, SiteFraction, SystemMolesType
import numpy as np
[docs]
def lower_convex_hull(global_grid, state_variables, conds_keys, phase_record_factory, result_array):
"""
Find the simplices on the lower convex hull satisfying the specified
conditions in the result array.
Parameters
----------
global_grid : Dataset
A sample of the energy surface of the system.
state_variables : List[v.StateVariable]
A list of the state variables (e.g., P, T) used in this calculation.
conds_keys : List
A list of the keys of the conditions used in this calculation.
phase_record_factory : PhaseRecordFactory
PhaseRecordFactory object corresponding to this calculation.
result_array : Dataset
This object will be modified!
Coordinates correspond to conditions axes.
Returns
-------
None. Results are written to result_array.
Notes
-----
This routine will not check if any simplex is degenerate.
Degenerate simplices will manifest with duplicate or NaN indices.
Examples
--------
None yet.
"""
state_variables = sorted(state_variables, key=str)
local_conds_keys = [c for c in conds_keys if getattr(c, 'phase_name', None) is not None]
str_conds_keys = [str(c) for c in conds_keys]
# factored out via profiling
result_array_GM_values = result_array.GM
result_array_points_values = result_array.points
result_array_MU_values = result_array.MU
result_array_NP_values = result_array.NP
result_array_X_values = result_array.X
result_array_Y_values = result_array.Y
result_array_Phase_values = result_array.Phase
global_grid_GM_values = global_grid.GM
global_grid_X_values = global_grid.X
global_grid_Y_values = global_grid.Y
global_grid_Phase_values = global_grid.Phase
num_comps = len(result_array.coords['component'])
it = np.nditer(result_array_GM_values, flags=['multi_index'])
while not it.finished:
primary_index = it.multi_index
grid_index = []
# Relies on being ordered
for lc in local_conds_keys:
coord_idx = conds_keys.index(lc)
grid_index.append(primary_index[coord_idx])
for sv in state_variables:
if sv in conds_keys:
coord_idx = conds_keys.index(sv)
grid_index.append(primary_index[coord_idx])
else:
# free state variable
grid_index.append(0)
grid_index = tuple(grid_index)
idx_global_grid_X_values = global_grid_X_values[grid_index]
idx_global_grid_GM_values = global_grid_GM_values[grid_index]
idx_result_array_MU_values = result_array_MU_values[it.multi_index]
idx_result_array_MU_values[:] = 0
idx_fixed_lincomb_molefrac_coefs = []
idx_fixed_lincomb_molefrac_rhs = []
idx_fixed_chempot_indices = []
for coord_idx, str_cond_key in enumerate(sorted(result_array.coords.keys())):
try:
cond_key = conds_keys[str_conds_keys.index(str_cond_key)]
except ValueError:
continue
rhs = result_array.coords[str_cond_key][primary_index[coord_idx]]
if isinstance(cond_key, IndependentPotential):
# Already handled above in construction of grid_index
continue
if isinstance(cond_key, SiteFraction):
# Already handled above in construction of grid_index
continue
elif isinstance(cond_key, ChemicalPotential):
component_idx = result_array.coords['component'].index(str(cond_key.species))
idx_fixed_chempot_indices.append(component_idx)
idx_result_array_MU_values[component_idx] = rhs
elif isinstance(cond_key, MassFraction):
# wA = k -> (1-k)*MWA*xA - k*MWB*xB - k*MWC*xC = 0
component_idx = result_array.coords['component'].index(str(cond_key.species))
coef_vector = np.zeros(num_comps)
coef_vector -= rhs
coef_vector[component_idx] += 1
# multiply coef_vector times a vector of molecular weights
coef_vector = np.multiply(coef_vector, phase_record_factory.molar_masses)
idx_fixed_lincomb_molefrac_coefs.append(coef_vector)
idx_fixed_lincomb_molefrac_rhs.append(0.)
elif isinstance(cond_key, MoleFraction):
if cond_key.phase_name is not None:
# Phase-local condition already handled in construction of grid
continue
component_idx = result_array.coords['component'].index(str(cond_key.species))
coef_vector = np.zeros(num_comps)
coef_vector[component_idx] = 1
idx_fixed_lincomb_molefrac_coefs.append(coef_vector)
idx_fixed_lincomb_molefrac_rhs.append(rhs)
elif isinstance(cond_key, SystemMolesType):
coef_vector = np.ones(num_comps)
idx_fixed_lincomb_molefrac_coefs.append(coef_vector)
idx_fixed_lincomb_molefrac_rhs.append(rhs)
elif isinstance(cond_key, LinearCombination):
coef_vector = np.zeros(num_comps)
if cond_key.denominator == 1:
for symbol_idx, symbol in enumerate(cond_key.symbols):
if symbol != 1:
coef_idx = result_array.coords['component'].index(str(symbol.species))
coef_vector[coef_idx] = cond_key.coefs[symbol_idx]
else:
idx_fixed_lincomb_molefrac_rhs.append(rhs-cond_key.coefs[symbol_idx])
else:
# This is a molar ratio
denominator_idx = cond_key.symbols.index(cond_key.denominator)
for symbol_idx, symbol in enumerate(cond_key.symbols):
if symbol_idx == denominator_idx:
coef_idx = result_array.coords['component'].index(str(symbol.species))
coef_vector[coef_idx] = cond_key.coefs[symbol_idx] - rhs
elif symbol != 1:
coef_idx = result_array.coords['component'].index(str(symbol.species))
coef_vector[coef_idx] = cond_key.coefs[symbol_idx]
else:
if cond_key.coefs[symbol_idx] != 0:
# Constant term for molar ratio should be zero
raise ValueError(f'Unsupported condition {cond_key}')
idx_fixed_lincomb_molefrac_rhs.append(-cond_key.coefs[symbol_idx])
idx_fixed_lincomb_molefrac_coefs.append(coef_vector)
else:
raise ValueError(f'Unsupported condition {cond_key}')
idx_fixed_lincomb_molefrac_coefs = np.atleast_2d(idx_fixed_lincomb_molefrac_coefs)
idx_fixed_lincomb_molefrac_rhs = np.atleast_1d(idx_fixed_lincomb_molefrac_rhs)
idx_fixed_chempot_indices = np.array(idx_fixed_chempot_indices, dtype=np.uintp)
idx_result_array_NP_values = result_array_NP_values[it.multi_index]
idx_result_array_points_values = result_array_points_values[it.multi_index]
result_array_GM_values[it.multi_index] = \
hyperplane(idx_global_grid_X_values, idx_global_grid_GM_values,
idx_result_array_MU_values, idx_fixed_chempot_indices, idx_fixed_lincomb_molefrac_coefs, idx_fixed_lincomb_molefrac_rhs,
idx_result_array_NP_values, idx_result_array_points_values)
# Copy phase values out
points = result_array_points_values[it.multi_index]
result_array_Phase_values[it.multi_index][:num_comps] = global_grid_Phase_values[grid_index].take(points, axis=0)[:num_comps]
result_array_X_values[it.multi_index][:num_comps] = global_grid_X_values[grid_index].take(points, axis=0)[:num_comps]
result_array_Y_values[it.multi_index][:num_comps,:global_grid_Y_values.shape[-1]] = \
global_grid_Y_values[grid_index].take(points, axis=0)[:num_comps]
# Special case: Sometimes fictitious points slip into the result
if '_FAKE_' in result_array_Phase_values[it.multi_index]:
new_energy = 0.
molesum = 0.
for idx in range(len(result_array_Phase_values[it.multi_index])):
midx = it.multi_index + (idx,)
if result_array_Phase_values[midx] == '_FAKE_':
result_array_Phase_values[midx] = ''
result_array_X_values[midx] = np.nan
result_array_Y_values[midx] = np.nan
idx_result_array_NP_values[idx] = np.nan
else:
new_energy += idx_result_array_NP_values[idx] * global_grid.GM[np.index_exp[grid_index + (points[idx],)]]
molesum += idx_result_array_NP_values[idx]
if molesum != 0:
result_array_GM_values[it.multi_index] = new_energy / molesum
it.iternext()
result_array.remove('points')
return result_array