Disagreement with Thermo-Calc for Fe-Ni-Ti on develop

[bocklund]

I noticed some odd behavior where calculating NP vs. T with pycalphad will give 100% liquid until about T=1464 K, where it drops off to NP(liquid)~=0.69 after being NP(liquid)=1 at T=1465 K. I would expect that the phase fraction drops off smoothly and this behavior is shown in Thermo-Calc using the same database.

from pycalphad import Database, equilibrium, variables as v
from pycalphad import __version__ as vers

dbf = Database('FeNiTi.tdb')
comps = ['FE', 'NI', 'TI', 'VA']

out = equilibrium(dbf, comps, ['LIQUID', 'C14'], {v.X('NI'):0.264503, v.X('TI'): 0.26734, v.T: 1675, v.P: 101325, v.N: 1})
# Tested with pdens up to 100_000 without seeing convergence to the C14 and LIQUID phases
#out = equilibrium(dbf, comps, ['LIQUID', 'C14'], {v.X('NI'):0.264503, v.X('TI'): 0.26734, v.T: 1675, v.P: 101325, v.N: 1}, calc_opts={'pdens': 100000})

print(f"pycalphad {vers} Phases={out.Phase.values.squeeze()}, TC=['LIQUID', 'C14']")
print(f"pycalphad {vers} GM={out.GM.values.squeeze()}, TC=-123262.71")
print(f"pycalphad {vers} MU={out.MU.values.squeeze()}, TC=[-106356.02, -133027.43, -143208.06]")

On develop gives:

pycalphad 0.8+8.gabccdf5a Phases=['LIQUID' '' '' ''], TC=['LIQUID', 'C14']
pycalphad 0.8+8.gabccdf5a GM=-123261.07316905519, TC=-123262.71
pycalphad 0.8+8.gabccdf5a MU=[-106447.61455119 -133405.76872427 -142667.20744647], TC=[-106356.02, -133027.43, -143208.06]

FeNiTi.tdb.zip

[bocklund]

Here is the same, using user-starting-point with the starting point from the Thermo-Calc solution in the script. Note that I did not turn global min off in the starting point generation, so it's probably just a grid sampling issue.

from pycalphad import Database, equilibrium, variables as v
# calculating NP vs. T with pycalphad will give 100% liquid until about T=1465, where it drops off to NP(liquid)~=0.69
from pycalphad import __version__ as vers

dbf = Database('FeNiTi.tdb')
comps = ['FE', 'NI', 'TI', 'VA']

# from TC
startpt = [
('LIQUID', [1.0, 101325.0, 1675.0, 0.46239099, 0.29075398, 0.24685503]),
# TC:
('C14', [1.0, 101325.0, 1675.0, 5.0182131E-2, 9.3858474E-3, 0.94043202, 0.78744661, 0.21116113, 1.3922589E-3, 0.42533288, 0.57466605, 1.0614344E-6]),
]

out = equilibrium(dbf, comps, ['LIQUID', 'C14'], {v.X('NI'):0.264503, v.X('TI'): 0.26734, v.T: 1675, v.P: 101325, v.N: 1}, user_starting_point=startpt)

print(f"pycalphad {vers} Phases={out.Phase.values.squeeze()}, TC=['LIQUID', 'C14']")
print(f"pycalphad {vers} GM={out.GM.values.squeeze()}, TC=-123262.71")
print(f"pycalphad {vers} MU={out.MU.values.squeeze()}, TC=[-106356.02, -133027.43, -143208.06]")

pycalphad 0.8+13.g651af073 Phases=['LIQUID' 'C14' '' ''], TC=['LIQUID', 'C14']
pycalphad 0.8+13.g651af073 GM=-123262.69627695215, TC=-123262.71
pycalphad 0.8+13.g651af073 MU=[-106355.99121553 -133027.33770931 -143208.14151373], TC=[-106356.02, -133027.43, -143208.06]

[bocklund]

Going to close this issue for now, since it is a grid sampling/starting point issue. The "correct" solution would be to take a scientifically different approach to sampling high-dimensional energy surfaces, which is more in scope of a larger roadmap task than an issue that could be closed by a simple PR.

This might be a good test case to use when developing a new sampling method, but whatever new method used will probably have different tricky edge cases than the edge cases for the current grid sampler.