10
$\begingroup$

What's the most efficient way to get the equilibrium phases and phase fractions for a set of alloy compositions? For context, I'm coming from Thermo-Calc/TC-Python where I would just loop through the alloy compositions and do a single equilbrium calculation. However, when I attempt to do this with pycalphad it seems quite slow. For example, here is what my code for the "one at a time" looping:

from pycalphad import Database, calculate, equilibrium, variables as v
dbf = Database("Cr-Fe-Ni.tdb")
components = ['FE','CR', 'NI', 'VA'] 
phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'SIGMA']

#Example set of alloy compositions, normally I would want to sample more finely but I'll keep it simple for demonstration purposes
compositions = array([[0. , 0. ],
                      [0. , 0.2],
                      [0. , 0.4],
                      [0. , 0.6],
                      [0. , 0.8],
                      [0.2, 0. ],
                      [0.2, 0.2],
                      [0.2, 0.4],
                      [0.2, 0.6],
                      [0.2, 0.8],
                      [0.4, 0. ],
                      [0.4, 0.2],
                      [0.4, 0.4],
                      [0.4, 0.6],
                      [0.6, 0. ],
                      [0.6, 0.2],
                      [0.6, 0.4],
                      [0.8, 0. ],
                      [0.8, 0.2]])

phase_data = []
for i in tqdm.tqdm(range(len(compositions)), total=len(compositions), position=0, leave=True, smoothing=0.3):
    eq_result = equilibrium(dbf, components, phases, {v.X('CR'):compositions[i,0], v.X('NI'):compositions[i,1], v.T:1000, v.P:101325})
    phase_data.append([*eq_result.NP.values.squeeze().tolist(), *eq_result.Phase.values.squeeze().tolist()])

This takes about 15 seconds, or about 1.2 alloys evaluated per second. This seems quite slow for a ternary, and I suspect it's because I'm calling the equilibrium function each time instead of just changing the composition.

I tried feeding the whole column at once (see code below), but this seemed to take the Cartesian product of the compositions, which is not what I want.

#database, phase, and composition setup is same as in above code

eq_result = equilibrium(dbf, components, phases, {v.X('CR'):compositions[:,0], v.X('NI'):compositions[:,1], v.T:1000, v.P:101325})
 

So, please let me know what the better approach is. Thanks!

$\endgroup$
2
  • 2
    $\begingroup$ Can I ask for some clarification on what your specific goal is (particularly why the Cartesian product is non-viable)? There is some overhead in calling equilibrium. There are some possible workarounds today and the project is in the process of stabilizing a lower-level API for tight loops of calculations like this (including a new, faster minimizer). For calculations today, a Cartesian product may be faster even if there are more overall compositions being computed because of the amortized overhead. $\endgroup$ – Brandon Bocklund May 31 at 5:42
  • $\begingroup$ @BrandonBocklund Usually I like to have specific alloy compositions, which are often subject to constraints (e.g. at 40at% of element A, etc.) and I usually want to test many more alloys (100s of thousands). I guess for the time being I can just test the whole range of compositions. For that should I feed v.X('CR'):[0.0, 0.2, 0.4, 0.6, 0.8] and same for NI? If so, won't the cartesian product result in alloys with atomic fractions whose total is greater than one (e.g. CR=0.6, NI=0.6)? And if that's the case how can I parse the data from eq_result to remove/ignore these entries? $\endgroup$ – sgp45 May 31 at 15:22
6
$\begingroup$

In your comment, you mentioned that you may be testing many alloys (100,000s). I will assume that interested in taking more of a screening approach, rather than specifying precisely 100,000 alloys compositions.

In that case, it is easy and faster to let pycalphad broadcast the conditions for you. It will skip any conditions that have independent mole fractions that sum to greater than unity, e.g. X(CR)=0.8 and X(NI)=0.8. The result that you get from equilibrium is an xarray Dataset object and uses padded NumPy arrays. Floating point values (like NP) will be padded with nan and string values (like Phase) will be padded with '' when no phases are present.

xarray and most of the tools in the ecosystem (like matplotlib) play nicely with these padded values. Below, I give an example where equilibrium calculations are performed on a dense grid and I plot the phase fractions of FCC over the whole grid and under a constraint where the composition of Cr is under 40 at%.

from pycalphad import Database, calculate, equilibrium, variables as v
import matplotlib.pyplot as plt

dbf = Database("Cr-Fe-Ni.tdb")
components = ['FE','CR', 'NI', 'VA'] 
phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'SIGMA']

eq_result = equilibrium(dbf, components, phases, {v.X('CR'): (0, 1, 0.1), v.X('NI'): (0, 1, 0.1), v.T:1000, v.P:101325})


# verify what the Phase and NP look like by printing
# note these are multi-dimensional arrays and are padded with '' and nan, respectively
print(eq_result.Phase.values)
print(eq_result.NP.values)


# Plot FCC phase fractions over the whole grid
plt.figure()
plt.subplot(projection="triangular")  # projection provided by pycalphad

mask = eq_result.Phase == "FCC_A1" # choose only FCC, others will be NaN
plt.scatter(
    eq_result.X_CR.where(mask).broadcast_like(mask),
    eq_result.X_NI.where(mask).broadcast_like(mask), 
    c=eq_result.NP.where(mask)
)
plt.title("Cr-Fe-Ni FCC Phase Fractions")
plt.xlabel("X(CR)")
plt.ylabel("X(NI)")
plt.colorbar(label="NP(FCC_A1)")
plt.show()


# Plot FCC phase fractions only at compositions under 40 at% Cr
plt.figure()
plt.subplot(projection='triangular')  # projection provided by pycalphad

# choose only FCC AND where the composition of Cr is less than 40%, others will be NaN
mask = (eq_result.Phase == 'FCC_A1') & (eq_result.X_CR < 0.4)
plt.scatter(
    eq_result.X_CR.where(mask).broadcast_like(mask),
    eq_result.X_NI.where(mask).broadcast_like(mask), 
    c=eq_result.NP.where(mask)
)
plt.title("Cr-Fe-Ni FCC Phase Fractions")
plt.xlabel("X(CR)")
plt.ylabel("X(NI)")
plt.colorbar(label="NP(FCC_A1)")
plt.show()

This produces the following two plots:

Cr-Fe-Ni full grid

Cr-Fe-Ni constrained grid

$\endgroup$
1
  • 1
    $\begingroup$ Fantastic, thanks! $\endgroup$ – sgp45 Jun 1 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.