Skip to content

Instantly share code, notes, and snippets.

@bocklund
Last active July 6, 2024 19:54
Show Gist options
  • Select an option

  • Save bocklund/45adb94036a2ee567adc0cffc53a3ded to your computer and use it in GitHub Desktop.

Select an option

Save bocklund/45adb94036a2ee567adc0cffc53a3ded to your computer and use it in GitHub Desktop.
Compute properties per phase in PyCalphad from the result of equilibrium
# coding: utf-8
import matplotlib.pyplot as plt
import numpy as np
from pycalphad import Database, calculate, equilibrium, variables as v
import pycalphad
dbf = Database('Al-Mg_Zhong.tdb')
comps = ['AL', 'MG', 'VA']
phases = pycalphad.core.utils.filter_phases(dbf, pycalphad.core.utils.unpack_components(dbf, comps), candidate_phases=None)
mg_composition = 0.08 # mole fraction
# Compute equilibrium as normal
eq = equilibrium(dbf, ['AL', 'MG', 'VA'], phases, {v.N: 1, v.P: 1e5, v.T: 880, v.X('MG'): mg_composition})
# Loop over N, P, T and compute the target property
output = "HM"
# Whatever per-phase output should have same dimension as NP
output_phase = np.full_like(eq.NP.values, np.nan)
it = np.nditer(output_phase, flags=['multi_index'])
while not it.finished:
# Assumes coordinates for the NP data variable are (N, P, T, ...)
N = float(eq["N"][it.multi_index[0]])
P = float(eq["P"][it.multi_index[1]])
T = float(eq["T"][it.multi_index[2]])
phase_name = str(eq.Phase.values[it.multi_index])
# assumes internal_dof is the last dimension
eq_sitefracs = np.atleast_2d(eq.Y.values[it.multi_index +np.index_exp[:]])
dof = len(np.nonzero(~np.isnan(eq_sitefracs))[1])
if phase_name != '':
cr = calculate(dbf, ["AL", "MG", "VA"], phase_name, N=N, P=P, T=T, points=eq_sitefracs[:, :dof], output=output)
output_phase[it.multi_index] = cr[output].values.squeeze()
it.iternext()
eq[f"{output}_phase"] = (eq.NP.dims, output_phase) # again, using same dimension as NP
print(eq)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment