/
utils.py
206 lines (173 loc) · 8.27 KB
/
utils.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
from __future__ import print_function
from operator import pos, neg
import numpy as np
from pycalphad import calculate, variables as v
from pycalphad.core.errors import ConditionError
from pycalphad.core.hyperplane import hyperplane
from pycalphad.core.equilibrium import _adjust_conditions
from pycalphad.core.cartesian import cartesian
from pycalphad.core.constants import MIN_SITE_FRACTION
from .compsets import BinaryCompSet
def v_array(center, distance, step):
"""create a grid like a "V" around a specific value going from inside to outside"""
xs = np.arange(0, distance, step).reshape(1, -1)
v = np.concatenate([xs, -1*xs]).reshape(-1, order='F') + center
return v[1:]
def build_composition_grid(components, conditions):
"""
Create a cartesion grid of compositions, including adding the dependent component.
Parameters
----------
components : list of str
List of component names
conditions : dict
Dictionary of pycalphad conditions
Returns
-------
np.ndarray
2D array of (M compositions, N components)
"""
comp_conds = sorted([x for x in conditions.keys() if isinstance(x, v.X)])
if len(comp_conds) > 0:
comp_values = cartesian([conditions[cond] for cond in comp_conds])
# Insert dependent composition value
# TODO: Handle W(comp) as well as X(comp) here
specified_components = {x.species.name for x in comp_conds}
all_comps = set(components) - {'VA'}
dependent_component = all_comps - specified_components
dependent_component = list(dependent_component)
if len(dependent_component) != 1:
raise ValueError('Number of dependent components is different from one')
else:
dependent_component = dependent_component[0]
insert_idx = sorted(all_comps).index(dependent_component)
comp_values = np.concatenate((comp_values[..., :insert_idx],
1 - np.sum(comp_values, keepdims=True, axis=-1),
comp_values[..., insert_idx:]),
axis=-1)
# Prevent compositions near an edge from going negative
comp_values[np.nonzero(comp_values < MIN_SITE_FRACTION)] = MIN_SITE_FRACTION*10
# TODO: Assumes N=1
comp_values /= comp_values.sum(axis=-1, keepdims=True)
return comp_values
def convex_hull(dbf, comps, phases, conditions, model=None,
calc_opts=None, parameters=None, callables=None):
"""
1D convex hull for fixed potentials.
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
StateVariables and their corresponding value.
model : Model, a dict of phase names to Model, or a seq of both, optional
Model class to use for each phase.
calc_opts : dict, optional
Keyword arguments to pass to `calculate`, the energy/property calculation routine.
parameters : dict, optional
Maps SymPy Symbol to numbers, for overriding the values of parameters in the Database.
callables : dict, optional
Pre-computed callable functions for equilibrium calculation.
Returns
-------
tuple
Tuple of (Gibbs energies, phases, phase fractions, compositions, site fractions, chemical potentials)
Notes
-----
Assumes that potentials are fixed and there is just a 1d composition grid.
Minimizes the use of Dataset objects.
"""
calc_opts = calc_opts or {}
conditions = _adjust_conditions(conditions)
# 'calculate' accepts conditions through its keyword arguments
if 'pdens' not in calc_opts:
calc_opts['pdens'] = 500
grid = calculate(dbf, comps, phases, T=conditions[v.T], P=conditions[v.P],
parameters=parameters, fake_points=True, output='GM',
callables=callables, model=model, **calc_opts)
# assume only one independent component
indep_comp_conds = [c for c in conditions if isinstance(c, v.X)]
num_indep_comp = len(indep_comp_conds)
if num_indep_comp != 1:
raise ConditionError(
"Convex hull independent components different than one.")
max_num_phases = (len(indep_comp_conds) + 1,) # Gibbs phase rule
comp_grid = build_composition_grid(comps, conditions)
calc_grid_shape = comp_grid.shape[:-1]
num_comps = comp_grid.shape[-1:]
grid_energy_values = grid.GM.values.squeeze()
grid_composition_values = grid.X.values.squeeze()
grid_site_frac_values = grid.Y.values.squeeze()
grid_phase_values = grid.Phase.values.squeeze()
# construct the arrays to pass to hyperplane
phase_fractions = np.empty(calc_grid_shape + max_num_phases)
chempots = np.empty(calc_grid_shape + num_comps)
simplex_points = np.empty(calc_grid_shape + max_num_phases, dtype=np.int32)
it = np.nditer(np.empty(calc_grid_shape), flags=['multi_index'])
while not it.finished:
idx = it.multi_index
hyperplane(grid_composition_values, grid_energy_values, comp_grid[idx],
chempots[idx], phase_fractions[idx], simplex_points[idx])
it.iternext()
simplex_phases = grid_phase_values[simplex_points] # shape: (calc_grid_shape + max_num_phases)
GM_values = grid_energy_values[simplex_points] # shape: (calc_grid_shape + max_num_phases)
phase_compositions = grid_composition_values[simplex_points] # shape: (calc_grid_shape + max_num_phases + num_comps)
phase_site_fracs = grid_site_frac_values[simplex_points] # shape: (calc_grid_shape + max_num_phases + num_internal_dof)
return (GM_values, simplex_phases, phase_fractions, phase_compositions, phase_site_fracs, chempots)
def get_num_phases(eq_dataset):
"""Return the number of phases in equilibrium from an equilibrium dataset"""
return int(np.sum(eq_dataset.Phase.values != '', axis=-1, dtype=np.int))
def get_compsets(eq_dataset, indep_comp=None, indep_comp_index=None):
"""Return a list of composition sets in an equilibrium dataset."""
if indep_comp is None:
indep_comp = [c for c in eq_dataset.coords if 'X_' in c][0][2:]
if indep_comp_index is None:
indep_comp_index = eq_dataset.component.values.tolist().index(indep_comp)
return BinaryCompSet.from_dataset_vertices(eq_dataset, indep_comp, indep_comp_index, 3)
def close_zero_or_one(val, tol):
zero = np.isclose(0, val, atol=tol)
one = np.isclose(1, val, atol=tol)
return zero or one
def close_to_same(val_1, val_2, tol):
return np.isclose(val_1, val_2, atol=tol)
def sort_x_by_y(x, y):
"""Sort a list of x in the order of sorting y"""
return [xx for _, xx in sorted(zip(y, x), key=lambda pair: pair[0])]
def opposite_direction(direction):
return neg if direction is pos else pos
def find_two_phase_region_compsets(hull_output, temperature, indep_comp, indep_comp_idx, discrepancy_tol=0.01):
"""
From a dataset at constant T and P, return the composition sets for a two
phase region or that have the smallest index composition coordinate
Parameters
----------
dataset : xr.Dataset
Equilibrium-like from pycalphad that has a `Phase` Data variable.
indep_comp_coord : str
Coordinate name of the independent component
Returns
-------
list
List of two composition sets for different phases
"""
phases, compositions, site_fracs = hull_output[1], hull_output[3], hull_output[4]
grid_shape = phases.shape[:-1]
num_phases = phases.shape[-1]
it = np.nditer(np.empty(grid_shape), flags=['multi_index']) # empty grid for indexing
while not it.finished:
idx = it.multi_index
cs = []
for i in np.arange(num_phases):
compset = BinaryCompSet(str(phases[idx][i]), temperature, indep_comp, compositions[idx][i, indep_comp_idx], site_fracs[idx][i, :])
cs.append(compset)
if len(set([c.phase_name for c in cs])) == 2:
# we found a multiphase region, return them if the discrepancy is
# above the tolerance
if cs[0].xdiscrepancy(cs[1], ignore_phase=True) > discrepancy_tol:
return cs
it.iternext()
return []