Solver: convergence failure in Cr-Ti-V

[bocklund]

This test is a convergence failure for Cr-Ti-V in a three phase region. It seems to be representative of a class of convergence failures where a starting point of 3 or more phases sees one phase removed inside the solver, which gets added back by the global grid minimization, and this loop continues until the maximum number of global iterations are reached.

The test is below. I attached the TDB and the diff to patch this test and TDB directly in.

@pytest.mark.solver
@select_database("crtiv_ghosh.tdb")
def test_ternary_three_phase_dilute(load_database):
    components = ["CR", "TI", "V"]
    dbf = load_database()
    phases = sorted(dbf.phases.keys())
    conditions = {v.N: 1.0, v.P: 101325.0, v.T: 500.0}

    # Got the MG to be stable with pdens=10
    eq_res = equilibrium(dbf, components+["VA"], set(phases), {v.X("CR"): 0.20, v.X("TI"): 0.20, **conditions}, verbose=True, calc_opts=dict(pdens=60), to_xarray=False)
    print(eq_res.Phase)
    assert np.all(np.isclose(eq_res.MU.squeeze(), [-26103.482, -16869.033, -19606.779]))
    # Sorting is pretty hacky, but gets the job done
    assert np.all(np.isclose(sorted(eq_res.NP.squeeze()[:3]), [8.9018E-02, 1.6608E-01, 7.4490E-01]))
    assert sorted(eq_res.Phase.squeeze()[:3]) == ["BCC_A2", "HCP_A3", "LAVES_C15"]
    raise ValueError()

The output I'm seeing on current develop:

Components: CR TI V VA
Phases: BCC_A2 
HCP_A3 
LAVES_C14 
LAVES_C15 
LAVES_C36 
LIQUID 
[done]
Adding metastable CompositionSet(HCP_A3, [1.e-14 1.e+00 1.e-14], NP=0.0, GM=-16842.677493938863) Driving force: 27605.287077657165
Adding metastable CompositionSet(LAVES_C14, [1.e-14 1.e+00 0.e+00], NP=0.0, GM=-11842.677493937581) Driving force: 22605.287077656212
Adding metastable CompositionSet(LAVES_C36, [1.e-14 1.e+00 1.e-14], NP=0.0, GM=-11842.67749393893) Driving force: 22605.28707765723
Adding metastable CompositionSet(LIQUID, [1.e-14 1.e+00 1.e-14], NP=0.0, GM=-8138.7314939390635) Driving force: 18901.341077657366
Chemical Potentials [-27994.07356931 -13295.12487247 -19132.76557249]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.56905861e-01
 2.54514624e-02 1.76426761e-02 5.31736011e-12 9.99984073e-01
 1.59269213e-05 9.35921618e-02 1.63389500e-01 7.43018338e-01
 1.00000000e+00 1.95531261e-01 8.04468739e-01]
Adding CompositionSet(HCP_A3, [1.e-14 1.e+00 1.e-14], NP=1e-06, GM=-16842.677493938863) Driving force: 3547.552621463874
Chemical Potentials [-26122.22086588 -16876.38809401 -17672.27700642]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.81582156e-01
 5.19687470e-03 1.32209695e-02 8.66336786e-12 9.99980790e-01
 1.92103497e-05 4.03931007e-02 3.18662062e-01 6.40944837e-01
 1.00000000e+00 2.25009643e-05 9.95498905e-01 4.47859373e-03
 1.00000000e+00 1.94330105e-01 8.05669895e-01 1.00000000e-16]
Adding CompositionSet(BCC_A2, [1.69491525e-01 1.00000000e-14 8.30508475e-01], NP=1e-06, GM=-20655.233593235593) Driving force: 1550.7627123296224
Chemical Potentials [-28010.1421399  -13260.35279569 -19159.16314892]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.48236755e-01
 3.25239289e-02 1.92393157e-02 4.47296745e-12 9.99984632e-01
 1.53683916e-05 1.11657593e-01 1.14404810e-01 7.73937597e-01
 1.00000000e+00 1.94330105e-01 8.05669895e-01]
Adding CompositionSet(HCP_A3, [1.e-14 1.e+00 1.e-14], NP=1e-06, GM=-16842.677493938863) Driving force: 3582.3246982484816
Chemical Potentials [-26122.22086588 -16876.38809401 -17672.27700642]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.81582156e-01
 5.19687470e-03 1.32209695e-02 8.66336786e-12 9.99980790e-01
 1.92103497e-05 4.03931007e-02 3.18662062e-01 6.40944837e-01
 1.00000000e+00 2.25009643e-05 9.95498905e-01 4.47859373e-03
 1.00000000e+00 1.94330105e-01 8.05669895e-01 1.00000000e-16]
Adding CompositionSet(BCC_A2, [1.69491525e-01 1.00000000e-14 8.30508475e-01], NP=1e-06, GM=-20655.233593235593) Driving force: 1550.762712329626
Chemical Potentials [-28010.1421399  -13260.35279569 -19159.16314892]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.48236755e-01
 3.25239289e-02 1.92393157e-02 4.47296745e-12 9.99984632e-01
 1.53683916e-05 1.11657593e-01 1.14404810e-01 7.73937597e-01
 1.00000000e+00 1.94330105e-01 8.05669895e-01]
Adding CompositionSet(HCP_A3, [1.e-14 1.e+00 1.e-14], NP=1e-06, GM=-16842.677493938863) Driving force: 3582.3246982467517
Chemical Potentials [-26122.22086588 -16876.38809401 -17672.27700642]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.81582156e-01
 5.19687470e-03 1.32209695e-02 8.66336786e-12 9.99980790e-01
 1.92103497e-05 4.03931007e-02 3.18662062e-01 6.40944837e-01
 1.00000000e+00 2.25009643e-05 9.95498905e-01 4.47859373e-03
 1.00000000e+00 1.94330105e-01 8.05669895e-01 1.00000000e-16]
Adding CompositionSet(BCC_A2, [1.69491525e-01 1.00000000e-14 8.30508475e-01], NP=1e-06, GM=-20655.233593235593) Driving force: 1550.7627123296807
Chemical Potentials [-28010.1421399  -13260.35279569 -19159.16314892]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.48236755e-01
 3.25239289e-02 1.92393157e-02 4.47296745e-12 9.99984632e-01
 1.53683916e-05 1.11657593e-01 1.14404810e-01 7.73937597e-01
 1.00000000e+00 1.94330105e-01 8.05669895e-01]
Adding CompositionSet(HCP_A3, [1.e-14 1.e+00 1.e-14], NP=1e-06, GM=-16842.677493938863) Driving force: 3582.324698248958
Chemical Potentials [-26122.22086588 -16876.38809401 -17672.27700642]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.81582156e-01
 5.19687470e-03 1.32209695e-02 8.66336786e-12 9.99980790e-01
 1.92103497e-05 4.03931007e-02 3.18662062e-01 6.40944837e-01
 1.00000000e+00 2.25009643e-05 9.95498905e-01 4.47859373e-03
 1.00000000e+00 1.94330105e-01 8.05669895e-01 1.00000000e-16]
Adding CompositionSet(BCC_A2, [1.69491525e-01 1.00000000e-14 8.30508475e-01], NP=1e-06, GM=-20655.233593235593) Driving force: 1550.7627123296588
Chemical Potentials [-28010.1421399  -13260.35279569 -19159.16314892]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.48236755e-01
 3.25239289e-02 1.92393157e-02 4.47296745e-12 9.99984632e-01
 1.53683916e-05 1.11657593e-01 1.14404810e-01 7.73937597e-01
 1.00000000e+00 1.94330105e-01 8.05669895e-01]
Adding CompositionSet(HCP_A3, [1.e-14 1.e+00 1.e-14], NP=1e-06, GM=-16842.677493938863) Driving force: 3582.3246982507662
Chemical Potentials [-26122.22086588 -16876.38809401 -17672.27700642]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.81582156e-01
 5.19687470e-03 1.32209695e-02 8.66336786e-12 9.99980790e-01
 1.92103497e-05 4.03931007e-02 3.18662062e-01 6.40944837e-01
 1.00000000e+00 2.25009643e-05 9.95498905e-01 4.47859373e-03
 1.00000000e+00 1.94330105e-01 8.05669895e-01 1.00000000e-16]
Adding CompositionSet(BCC_A2, [1.69491525e-01 1.00000000e-14 8.30508475e-01], NP=1e-06, GM=-20655.233593235593) Driving force: 1550.7627123296043
Chemical Potentials [-28010.1421399  -13260.35279569 -19159.16314892]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 9.48236755e-01
 3.25239289e-02 1.92393157e-02 4.47296745e-12 9.99984632e-01
 1.53683916e-05 1.11657593e-01 1.14404810e-01 7.73937597e-01
 1.00000000e+00 1.94330105e-01 8.05669895e-01]
[[[[[['' '' '' '']]]]]]

crtiv_fail.patch.txt
crtiv_ghosh.tdb.txt

[bocklund]

On closer inspection, it's not clear whether this is a solver failure or starting point failure. I noticed that we are adding composition sets with positive driving force before even running the solver. I printed out the chemical potentials just before add nearly stable and got: [-39857.862145 10762.60958372 -22565.20996254]

The starting point composition sets are:

[
CompositionSet(LAVES_C15, [6.66666667e-01 3.33333333e-01 1.00000000e-14], NP=0.2795026017216821, GM=-22984.371568762246),
CompositionSet(BCC_A2, [1.00000000e-14 1.52542373e-01 8.47457627e-01], NP=0.2740378620490976, GM=-17481.305285989976), 
CompositionSet(BCC_A2, [0.03060733 0.14565728 0.8237354 ], NP=0.44645953622922036, GM=-18240.052345507993)
]

These chemical potentials do seem to be correct w.r.t. the composition sets (i.e. sum(MU(i) * X(i)) = GM and mass balance is satisfied w.r.t. the conditions), but the composition sets do not seem to be the global minimum amongst the grid.

[richardotis]

Reproducing on Windows.

[richardotis]

Looks like hyperplane is hitting max_iterations for this case and silently returning a metastable result.

best_guess_simplex [0 1 2]
best_guess_simplex [0 1 2]
best_guess_simplex [1091    1    2]
best_guess_simplex [93  1  2]
best_guess_simplex [93  1 68]
best_guess_simplex [ 93 265  68]
best_guess_simplex [ 93 230  68]
best_guess_simplex [ 93 174  68]
best_guess_simplex [122 174  68]
best_guess_simplex [358 174  68]
best_guess_simplex [358 506  68]
best_guess_simplex [643 506  68]
best_guess_simplex [578 506  68]
best_guess_simplex [583 506  68]
best_guess_simplex [602 506  68]
best_guess_simplex [602 506 601]
best_guess_simplex [519 506 601]
best_guess_simplex [524 506 601]
best_guess_simplex [524 539 601]
best_guess_simplex [524 541 601]
best_guess_simplex [524 541 723]
best_guess_simplex [524 541 643]
best_guess_simplex [524 541 650]
best_guess_simplex [692 541 650]
best_guess_simplex [297 541 650]
best_guess_simplex [ 34 541 650]
best_guess_simplex [ 34 541  57]
best_guess_simplex [ 40 541  57]
best_guess_simplex [ 34 541  57]
best_guess_simplex [ 40 541  57]
best_guess_simplex [ 34 541  57]
best_guess_simplex [ 40 541  57]
best_guess_simplex [ 34 541  57]
[snip]

There's some kind of oscillation happening. Details:

best_guess_simplex [ 34 541  57]
compositions [[4.74576271e-01 5.25423729e-01 1.00000000e-14]
 [3.06073251e-02 1.45657277e-01 8.23735398e-01]
 [8.64406780e-01 1.35593220e-01 1.00000000e-14]]
result_fractions [0.14641294 0.72838924 0.12519783]
candidate_potentials [-14195.14560945 -17018.56111746 -11985.99718834]
candidate_energy (NP*GM) -13434.339658386109

best_guess_simplex [ 40 541  57]
compositions [[5.76271186e-01 4.23728814e-01 1.00000000e-14]
 [3.06073251e-02 1.45657277e-01 8.23735398e-01]
 [8.64406780e-01 1.35593220e-01 1.00000000e-14]]
result_fractions [0.19808809 0.72838924 0.07352267]
candidate_potentials [-14099.24186462 -17629.94749077 -11881.45206683]
candidate_energy (NP*GM) -13474.709111175876

[richardotis]

Here's a visualization of the candidate simplices (triangles here) with the target composition point marked. I don't see any problem numerically with estimating the barycentric coordinates (phase fractions) here. Both triangles are feasible, but one is lower energy than the other at the target point, so we shouldn't be seeing this oscillation.

[richardotis]

There was no oscillation problem in the develop branch. I had made a mistake in some debugging code that accidentally changed the behavior.

The issue may be a problem in some code that is supposed to accelerate the hyperplane calculation by throwing out metastable points:

pycalphad/pycalphad/core/hyperplane.pyx

Lines 241 to 247 in 04b870e

	for i in range(num_points):
	if driving_forces[i] < 1.0:
	remaining_point_indices[ici] = remaining_point_indices[i]
	if driving_forces[i] < lowest_df:
	lowest_df = driving_forces[i]
	min_df = ici
	ici += 1

If I force that first branch to always be taken (if True:), the starting point/potentials are significantly improved and the test passes. The minimizer still adds some composition sets, but they are all metastable with respect to the starting point. The number of iterations in hyperplane also does not increase significantly for this case (from 7 to 10).

Components: CR TI V VA
Phases: BCC_A2
HCP_A3
LAVES_C14
LAVES_C15
LAVES_C36
LIQUID
[done]
candidate_potentials [-26080.7644499  -16842.67749394 -19577.50191689]
candidate_energy (MU*X) -20331.189538901348
smallest_fractions [  0.09676002 -34.50172414 -16.09090909]
best_guess_simplex [ 486   77 1937]
lowest_df -9.379164112033322e-13
trial_simplices [[1937   77 1937]
 [ 486 1937 1937]
 [ 486   77 1937]]
num_iterations 10
Adding metastable CompositionSet(LAVES_C14, [0.66666667 0.33333333 0.        ], NP=0.0, GM=-22373.026568761084) Driving force: -628.3564780091801
Adding metastable CompositionSet(LAVES_C36, [0.65536723 0.33333333 0.01129944], NP=0.0, GM=-22739.43040771062) Driving force: -188.48853107354807
Chemical Potentials [-26103.48124821 -16869.03299259 -19606.77561676]
[1.00000000e+00 1.01325000e+05 5.00000000e+02 5.11950690e-05
 9.93548306e-01 6.40049853e-03 1.00000000e+00 1.90022112e-01
 6.65953816e-03 8.03318350e-01 1.00000000e+00 9.84813523e-01
 6.02520547e-03 9.16127166e-03 1.73905853e-11 9.99975704e-01
 2.42964143e-05 1.66081420e-01 7.44900946e-01 8.90176345e-02]
Composition Sets [CompositionSet(HCP_A3, [5.11950690e-05 9.93548306e-01 6.40049853e-03], NP=0.1660814198904083, GM=-16887.02866843402), CompositionSet(BCC_A2, [0.19002211 0.00665954 0.80331835], NP=0.7449009455715806, GM=-20823.061237758367), CompositionSet(LAVES_C15, [0.65654235 0.33734204 0.00611561], NP=0.08901763453668607, GM=-22948.582313288694)]
[[[[[['HCP_A3' 'BCC_A2' 'LAVES_C15' '']]]]]]

The assumption made in the accelerator code (that any point metastable at any particular iteration must be metastable at every subsequent iteration, and at the solution) may be incorrect, or maybe the assumption is being implemented incorrectly. Either way, it may be worth it to turn it off since it's counterproductive here.