ENH: Minimizer: refactor and convergence improvements

[bocklund]

This PR contains 3 changes to the minimizer:

1. The take_step function is replaced with two functions to solve the equilibrium matrix and extract the chemical potentials (solve_state) and to take the step with the corrections in the phase amounts, state variables, and site fractions (advance_state). The extract_equilibrium_solution function is also now a part of advance_state. This refactor keeps the logic for advancing the state separate and contained in one function.
2. The convergence criterion has been changed to converge based on changes in phase amounts, state variables, and site fractions, instead of depending on mass residual and/or chemical potential changes. The rationale being that if the phase amounts, state variables, and site fractions aren’t changing, the other quantities that could be convergence criteria are guaranteed to not be changing.
3. Steps in phase amounts are now scaled such that no phase amount can ever become negative. This change helps convergence by preventing phases having negative amount, then being removed, which messes up the mass balance. If a phase amount can never become negative, the moles in the system will not change if phases with zero amounts are not removed. Scaling in phase amounts is valid, because the phase energies are always in agreement with the chemical potentials.

These changes (particularly the third change) helps convergence as demonstrated here:

from pycalphad import Database, equilibrium, variables as v
from pycalphad.plot.utils import phase_legend
import matplotlib.pyplot as plt
import numpy as np

db = Database('scratch/mc_ni_v2.034.pycalphad.tdb')
phases = ['FCC_A1','GAMMA_PRIME', 'LIQUID']
components = ['AL', 'CO', 'CR', 'MO', 'NB', 'NI', 'TI', 'W','VA']
conditions = {v.X('AL'): 0.073, v.X('CO'): 0.077, v.X('CR'): 0.151,
               v.X('MO'): 0.021, v.X('NB'): 0.022, v.X('TI'): 0.030, v.X('W'): 0.011}
conditions[v.T] = (1550, 1650, 1) # Kelvin
conditions[v.P] = 1e5 # Pa
conditions[v.N] = 1
eq = equilibrium(db, components, phases, conditions)

phases = sorted(set(eq.Phase.values.flatten()) - {''})
phase_handles, phasemap = phase_legend(phases)

plt.gca().set_title('Phase fractions vs T')
plt.gca().set_xlabel('Temperature (K)')
plt.gca().set_ylabel('Phase Fraction')
plt.gca().set_ylim((0,1.1))
plt.gca().set_xlim((1550, 1650))

for name in phases:
    plt.scatter(eq.T.broadcast_like(eq.NP), eq.NP.where(eq.Phase == name), color=phasemap[name])
plt.gca().legend(phase_handles, phases, loc='upper left')
plt.savefig('NP.png')


plt.figure()
plt.gca().scatter(eq.T.broadcast_like(eq.GM), eq.GM)
plt.gca().set_title('Gibbs energy')
plt.gca().set_xlabel('Temperature (K)')
plt.gca().set_ylabel('Gibbs energy')
plt.gca().set_ylim(-120000, -107000)
plt.savefig('GM.png')

Thanks to @james-saal for sharing the motivating example on our Gitter channel!

Develop

this PR (e57eb04)

Also note that if I run this same script using phases = list(db.phases.keys()), then develop has a bunch of output lines containing ** On entry to DORMBRP parameter number 13 had an illegal value, finally ending with terminated by signal SIGSEGV (Address boundary error) (catastrophic failure). It works well on this PR.

[codecov]

Codecov Report

Merging #373 (e57eb04) into develop (e36e369) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff            @@
##           develop     #373   +/-   ##
========================================
  Coverage    89.92%   89.92%           
========================================
  Files           44       44           
  Lines         4368     4368           
========================================
  Hits          3928     3928           
  Misses         440      440           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e36e369...e57eb04. Read the comment docs.