calculating activities of ternary system

[dpttw]

Hi,

I followed the example to calculate the activities of Ni-Al-Cr
setup: ref_eq = equilibrium(TDB, elements, AllPhases, {v.X('AL'): 1, v.X('CR'): 0, v.T: 674, v.P: 101325})
result:
(N: 1, P: 1, T: 1, X_AL: 1, X_CR: 1, vertex: 4, component: 3)
refeq.MU.values.ravel(): [nan nan nan]
engine: pycalphad 0.8.3

Is there anything I am missing?

[bocklund]

Hi @dpttw, if all of your output chemical potentials are nan, that is usually indicating a convergence failure. Can you create a minimal, complete and verifiable example that we can use to reproduce these values and track down the source of the issue?

[bocklund]

This might be caused by v.X('CR'): 0, setting the composition condition to 1e-14 or another small, near zero number should fix it.

[richardotis]

Very small or zero conditions will be silently changed to be within numerical limits:

pycalphad/pycalphad/core/equilibrium.py

Lines 27 to 28 in d217b4b

	if isinstance(key, v.MoleFraction):
	new_conds[key] = [max(val, MIN_SITE_FRACTION*1000) for val in unpack_condition(value)]

However, even that change will give you 1e-12, which may still be too small for convergence. At one time we had the minimum composition condition set at 1e-8, but this was changed as the solver was updated. In my experience achieving convergence for the very dilute case is still challenging. Another thing to try is increasing the pdens setting to get more points near the edge of composition space.

[bocklund]

@richardotis you're right. Maybe this is another example case of #271, but instead of converging to a wrong answer, it doesn't converge at all.

Anecdotally, I've seen this come up several times.

[richardotis]

Checking this with the development version (373d1212c), you do still fail to get convergence. I suspect this has to do with the solver literally setting the constraint that a mole fraction should be 1, which is incompatible with having non-zero amounts of the other components. If you modify the X(AL) condition to be a bit less than one, you do get convergence as expected:

dbf = Database('mc_ni_v2.034.pycalphad.tdb')
comps = ['NI', 'AL', 'CR', 'VA']
phases = sorted(set(dbf.phases.keys()) - {'BCC_B2'}) # workaround for https://github.com/pycalphad/pycalphad/issues/345
ref_eq = equilibrium(dbf, comps, phases, {v.X('AL'): 1-2e-10, v.X('CR'): 0, v.T: 674, v.P: 101325})

One possible solution to this issue is analogous to how we modify X(CR): 0 to be X(CR): 1e-10. We could force X(*): 1 to be 1-n*1e-10, where n is the number of active components. That would make it easier for users to set "from zero to one" conditions without having to be pedantic about actually meaning "a little more than" zero and "a little less than" one.

[bocklund]

It looks like it's not solver related, but it's because the solver doesn't even run because of this check:

pycalphad/pycalphad/core/eqsolver.pyx

Line 197 in 373d121

if indep_sum > 1:

[bocklund]

Not that this doesn't affect the solver as well - I commented out that check and it still fails to converge. It looks like the allowed_mass_residual gets set to -1.000000082740371e-11 and the best we can do with X('AL')=1 is obviously 1-2e-10 without also breaking the N=1 condition or other composition conditions