pycalphad / scheil

A Scheil-Gulliver simulation tool using pycalphad.
https://scheil.readthedocs.io/
MIT License
18 stars 6 forks source link

Multicomponent systems stop prematurely #14

Closed AndRicci closed 3 years ago

AndRicci commented 3 years ago

Hi Brandon, I am trying to use pycalphad to simulate solidification behaviour of alloys 625 and Monel. I noticed that simulation ends at some temperature dependent on composition and temperature step, even if solidification is far from complete. I attach an extreme example for Monel, comparing temperature steps of 5K and 2K. The simulation with 2K stops prematurely. I noted a similar behavior with alloys 625, even if not so dramatic. Also changes in composition affects the end point.

image image

I attach the code and the database I am using, derived from mc_ni_v2.034.tdb Thank you very much for your help and all you did with pycalphad!

MonelScheil.txt ARNi.txt

bocklund commented 3 years ago

Hi @AndRicci, thanks for reporting this and especially for providing a script. I am able to reproduce something similar to these two plots with pycalphad 0.8.5. Can you share which version of pycalphad you are using?

The expected behavior is that:

If I use the verbose flag to simulate_scheil_solidification, I get the output:

building callables... done
T=1650.000, X_CU=0.28, X_FE=0.03, X_C=0.01, X_SI=0.02, X_MN=0.01, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1645.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1640.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1635.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1630.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1625.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1620.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
No liquid phase found at T=1615.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
T=1620.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=4.167, NL: 1.000, NP(LIQUID)=1.000
No liquid phase found at T=1615.833, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
T=1620.000, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 1.000, NP(LIQUID)=1.000
New phases seen: {'FCC_A1'}. T=1616.528, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 1.000, NP(LIQUID)=1.000 NP(FCC_A1)=0.000
T=1613.056, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 1.000, NP(LIQUID)=1.000
T=1609.583, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.937, NP(LIQUID)=0.937 NP(FCC_A1)=0.063
T=1606.111, X_C=0.01, X_CU=0.28, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.848, NP(LIQUID)=0.905 NP(FCC_A1)=0.095
T=1602.639, X_C=0.01, X_CU=0.29, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.768, NP(LIQUID)=0.906 NP(FCC_A1)=0.094
T=1599.167, X_C=0.01, X_CU=0.30, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.697, NP(LIQUID)=0.908 NP(FCC_A1)=0.092
T=1595.694, X_C=0.01, X_CU=0.31, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.634, NP(LIQUID)=0.909 NP(FCC_A1)=0.091
T=1592.222, X_C=0.01, X_CU=0.32, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.577, NP(LIQUID)=0.910 NP(FCC_A1)=0.090
T=1588.750, X_C=0.01, X_CU=0.33, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.525, NP(LIQUID)=0.911 NP(FCC_A1)=0.089
T=1585.278, X_C=0.01, X_CU=0.33, X_FE=0.03, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.479, NP(LIQUID)=0.912 NP(FCC_A1)=0.088
T=1581.806, X_C=0.01, X_CU=0.34, X_FE=0.02, X_MN=0.01, X_SI=0.02, ΔT=3.472, NL: 0.438, NP(LIQUID)=0.914 NP(FCC_A1)=0.086
T=1578.333, X_C=0.01, X_CU=0.35, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.401, NP(LIQUID)=0.915 NP(FCC_A1)=0.085
T=1574.861, X_C=0.01, X_CU=0.36, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.368, NP(LIQUID)=0.917 NP(FCC_A1)=0.083
T=1571.389, X_C=0.01, X_CU=0.37, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.338, NP(LIQUID)=0.919 NP(FCC_A1)=0.081
T=1567.917, X_C=0.01, X_CU=0.38, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.311, NP(LIQUID)=0.921 NP(FCC_A1)=0.079
T=1564.444, X_C=0.01, X_CU=0.40, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.288, NP(LIQUID)=0.924 NP(FCC_A1)=0.076
T=1560.972, X_C=0.01, X_CU=0.41, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.267, NP(LIQUID)=0.928 NP(FCC_A1)=0.072
T=1557.500, X_C=0.00, X_CU=0.42, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.249, NP(LIQUID)=0.932 NP(FCC_A1)=0.068
T=1554.028, X_C=0.00, X_CU=0.43, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.233, NP(LIQUID)=0.936 NP(FCC_A1)=0.064
T=1550.556, X_C=0.00, X_CU=0.44, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.219, NP(LIQUID)=0.940 NP(FCC_A1)=0.060
T=1547.083, X_C=0.00, X_CU=0.45, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.206, NP(LIQUID)=0.944 NP(FCC_A1)=0.056
T=1543.611, X_C=0.00, X_CU=0.46, X_FE=0.02, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.195, NP(LIQUID)=0.947 NP(FCC_A1)=0.053
T=1540.139, X_C=0.00, X_CU=0.47, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.186, NP(LIQUID)=0.950 NP(FCC_A1)=0.050
T=1536.667, X_C=0.00, X_CU=0.48, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.177, NP(LIQUID)=0.953 NP(FCC_A1)=0.047
T=1533.194, X_C=0.00, X_CU=0.49, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.169, NP(LIQUID)=0.956 NP(FCC_A1)=0.044
T=1529.722, X_C=0.00, X_CU=0.50, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.162, NP(LIQUID)=0.957 NP(FCC_A1)=0.043
T=1526.250, X_C=0.00, X_CU=0.50, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.155, NP(LIQUID)=0.959 NP(FCC_A1)=0.041
T=1522.778, X_C=0.00, X_CU=0.51, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.149, NP(LIQUID)=0.960 NP(FCC_A1)=0.040
T=1519.306, X_C=0.00, X_CU=0.52, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.143, NP(LIQUID)=0.961 NP(FCC_A1)=0.039
T=1515.833, X_C=0.00, X_CU=0.53, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.138, NP(LIQUID)=0.962 NP(FCC_A1)=0.038
T=1512.361, X_C=0.00, X_CU=0.53, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.133, NP(LIQUID)=0.962 NP(FCC_A1)=0.038
T=1508.889, X_C=0.00, X_CU=0.54, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.128, NP(LIQUID)=0.962 NP(FCC_A1)=0.038
T=1505.417, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=3.472, NL: 0.123, NP(LIQUID)=0.962 NP(FCC_A1)=0.038
No liquid phase found at T=1501.944, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
No liquid phase found at T=1505.417, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
No liquid phase found at T=1508.310, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
T=1510.721, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=2.009, NL: 0.123, NP(LIQUID)=1.000
T=1508.712, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=2.009, NL: 0.123, NP(LIQUID)=1.000
No liquid phase found at T=1506.703, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
T=1508.712, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=1.674, NL: 0.123, NP(LIQUID)=1.000
No liquid phase found at T=1507.038, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
T=1508.712, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=1.395, NL: 0.123, NP(LIQUID)=1.000
No liquid phase found at T=1507.317, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
T=1508.712, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02, ΔT=1.163, NL: 0.123, NP(LIQUID)=1.000
No liquid phase found at T=1507.549, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Stepping back and reducing step size.
No liquid phase found at T=1508.712, X_C=0.00, X_CU=0.55, X_FE=0.01, X_MN=0.02, X_SI=0.02. (Found set()) (Convergence failure) Maximum step size reduction exceeded. Stopping.

So it looks like convergence failures are probably causing the early termination because every failure it tries to step backward and reduce the step size up to a minimum allowed step size.

Maybe the correct behavior here should be that the step size gets reset if an equilibrium is successfully found?

AndRicci commented 3 years ago

Hi Brandon, I am using pycalphad 0.8.5 too. I agree that maybe resetting the step size could avoid some early terminations caused by convergence failures. In my case there should be no need to continue the simulation with such small step size. I share another example with alloy 625, where a change in composition makes the simulation stop earlier. The first two graphs are obtained with the code below and the same database already shared. The third one setting Si and Mn at 1w%. This system should reach an eutectic (NbC carbide), as opposed to the Monel one. Is this the reason for the stop?

image imageimage

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

# Load database
dbf = Database('mc_ni_v2.034.pycalphad.tdb')
# Define the components and phases
comps = ['NI', 'CR', 'MO', 'NB', 'FE', 'SI', 'MN', 'C', 'VA']
phases = ['LIQUID','BCC_A2','FCC_A1', 'LAVES']
# Create the dictionary of conditions
mass_fracs = {v.W('CR'): 0.21, 
              v.W('MO'): 0.08, 
              v.W('NB'): 0.035, 
              v.W('FE'): 0.04, 
              v.W('SI'): 0.5/100, 
              v.W('MN'):0.5/100, 
              v.W('C'): 0.03/100}
conds = v.get_mole_fractions(mass_fracs, 'NI', dbf)
conds[v.P] = 1e5
conds[v.N] = 1

phase_handles, phasemap = phase_legend(phases)

from scheil import simulate_scheil_solidification

# setup the simulation parameters
initial_composition = v.get_mole_fractions(mass_fracs, 'NI', dbf)
liquid_phase_name = 'LIQUID'
start_temperature = 1650

# perform the simulation
sol_res = simulate_scheil_solidification(dbf, comps, phases, initial_composition, start_temperature, step_temperature=5)

# plot the result
plt.figure(figsize=(9,6))
for phase_name, amounts in sol_res.cum_phase_amounts.items():
    plt.plot(sol_res.temperatures, amounts, label=phase_name)
plt.plot(sol_res.temperatures, sol_res.fraction_liquid, label='LIQUID')
plt.ylabel('Phase Fraction')
plt.xlabel('Temperature (K)')
plt.title('Alloy 625 Scheil simulation, phase fractions')
plt.legend(loc='best')
plt.show()

How do I use the verbose flag?

Thank you

Andrea

bocklund commented 3 years ago

This system should reach an eutectic (NbC carbide), as opposed to the Monel one. Is this the reason for the stop?

Using the verbose flag should tell you the reason why it stops, so you'll be able to check that way. To enable it, replace the line:

sol_res = simulate_scheil_solidification(dbf, comps, phases, initial_composition, start_temperature, step_temperature=5)

with

sol_res = simulate_scheil_solidification(dbf, comps, phases, initial_composition, start_temperature, step_temperature=5, verbose=True)

bocklund commented 3 years ago

One other thing: the if the carbide is a new phase that should precipitate, it should also be present in the list of phases: phases = ['LIQUID','BCC_A2','FCC_A1', 'LAVES']

AndRicci commented 3 years ago

It seems that both compositions of alloy 625 end the simulation with a "Maximum step size reduction exceeded." I attach two scripts with their verbose results. The NbC phase is described in the database as FCC_A1, so I cannot explicitly add it. pycalphad should detect the presence of a second FCC_A1 phase of different composition. Actually when I run an equilibrium calculation setting C=0.1w%, the NbC FCC_A1 phase is correctly accounted for. But when I run a Scheil simulation with the same composition it is not.

625Scheil1.txt 625Scheil2.txt

AndRicci commented 3 years ago

It seems that simulate_scheil_solidification gets confused when there are phases with order/disorder forms. For example, if I simulate a Ni-Si binary alloy with 0.1 Si using the nisi_tok.tdb database available from NIMS, the simulation goes on but the liquid fraction remains 1, until it finally stops with "Maximum step size reduction exceeded". The simulation runs okay if the ordered fcc BETA1 phase is deactivated. By the way, I upgraded to pycalphad 0.9.0.

bocklund commented 3 years ago

@AndRicci, I just merged #16 which fixes the case in the Ni-Si binary where the disordered configurations that become stable when the ordered phase is entered weren't being accounted for. Thanks for reporting this case and for upgrading.

I will also note that stopping with Maximum step size reduction exceeded isn't necessarily a failure mode. There are two ways to that the Scheil simulation can terminate:

  1. When the phase fraction of the liquid becomes lower than what is specified by the stop parameter of simulate_scheil_solidification (defaults to 1e-4).
  2. When the simulation finds steps down in temperature and finds no liquid, it goes back and takes a smaller step. This helps get past peritectic points if the liquid composition at the step above the peritectic reaction doesn't have any liquid phase below the peritectic reaction. Smaller steps will get it closer to the composition where it can continue to make progress. This means that when the simulation reaches an actual end-point, i.e. a eutectic, it will step back and reduce the step size to make sure it can't continue to make progress. At a eutectic, it will eventually reach the maximum step size reduction and exit successfully.

The converged attribute of the SolidificationResult object returned by simulate_scheil_solidificaition programmatically indicates type (1) success. The simulate_scheil_solidification cannot yet distinguish whether a type (2) termination was a successful simulation (exiting at a eutectic) or unsuccessful (exiting for some other reason, such as repeated convergence failures or failing to proceed past a eutectic point), so the converged attribute will be alwaysFalse for this termination mode. With the new minimizer in pycalphad 0.9, we hope to expose the ability to compute invariant reactions directly, so in the future scheil can more robustly handle these cases without relying on step size reduction at all.

bocklund commented 3 years ago

Here's the script and verbose output for the Ni-Si case:

from pycalphad import Database
import pycalphad.variables as v
from scheil import simulate_scheil_solidification

dbf = Database('nisi_Tok.tdb')
comps = ['NI', 'SI', 'VA']
phases = list(dbf.phases.keys())

initial_composition = {v.X('SI'): 0.1}
liquid_phase_name = 'LIQUID'
start_temperature = 1650

sol_res = simulate_scheil_solidification(dbf, comps, phases, initial_composition, start_temperature, step_temperature=5, verbose=True)
building callables... done
T=1650.000, X_SI=0.10, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1645.000, X_SI=0.10, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
T=1640.000, X_SI=0.10, ΔT=5.000, NL: 1.000, NP(LIQUID)=1.000
New phases seen: {'FCC_A1'}. T=1635.000, X_SI=0.10, ΔT=5.000, NL: 0.939, NP(LIQUID)=0.939 NP(FCC_A1)=0.061
T=1630.000, X_SI=0.10, ΔT=5.000, NL: 0.839, NP(LIQUID)=0.893 NP(FCC_A1)=0.107
T=1625.000, X_SI=0.10, ΔT=5.000, NL: 0.754, NP(LIQUID)=0.898 NP(FCC_A1)=0.102
T=1620.000, X_SI=0.11, ΔT=5.000, NL: 0.681, NP(LIQUID)=0.903 NP(FCC_A1)=0.097
T=1615.000, X_SI=0.11, ΔT=5.000, NL: 0.618, NP(LIQUID)=0.907 NP(FCC_A1)=0.093
T=1610.000, X_SI=0.11, ΔT=5.000, NL: 0.563, NP(LIQUID)=0.911 NP(FCC_A1)=0.089
T=1605.000, X_SI=0.12, ΔT=5.000, NL: 0.515, NP(LIQUID)=0.915 NP(FCC_A1)=0.085
T=1600.000, X_SI=0.12, ΔT=5.000, NL: 0.472, NP(LIQUID)=0.918 NP(FCC_A1)=0.082
T=1595.000, X_SI=0.12, ΔT=5.000, NL: 0.435, NP(LIQUID)=0.921 NP(FCC_A1)=0.079
T=1590.000, X_SI=0.12, ΔT=5.000, NL: 0.402, NP(LIQUID)=0.924 NP(FCC_A1)=0.076
T=1585.000, X_SI=0.13, ΔT=5.000, NL: 0.372, NP(LIQUID)=0.926 NP(FCC_A1)=0.074
T=1580.000, X_SI=0.13, ΔT=5.000, NL: 0.346, NP(LIQUID)=0.929 NP(FCC_A1)=0.071
T=1575.000, X_SI=0.13, ΔT=5.000, NL: 0.322, NP(LIQUID)=0.931 NP(FCC_A1)=0.069
T=1570.000, X_SI=0.13, ΔT=5.000, NL: 0.300, NP(LIQUID)=0.933 NP(FCC_A1)=0.067
T=1565.000, X_SI=0.14, ΔT=5.000, NL: 0.281, NP(LIQUID)=0.935 NP(FCC_A1)=0.065
T=1560.000, X_SI=0.14, ΔT=5.000, NL: 0.263, NP(LIQUID)=0.937 NP(FCC_A1)=0.063
T=1555.000, X_SI=0.14, ΔT=5.000, NL: 0.247, NP(LIQUID)=0.938 NP(FCC_A1)=0.062
T=1550.000, X_SI=0.14, ΔT=5.000, NL: 0.232, NP(LIQUID)=0.940 NP(FCC_A1)=0.060
T=1545.000, X_SI=0.15, ΔT=5.000, NL: 0.218, NP(LIQUID)=0.942 NP(FCC_A1)=0.058
T=1540.000, X_SI=0.15, ΔT=5.000, NL: 0.206, NP(LIQUID)=0.943 NP(FCC_A1)=0.057
T=1535.000, X_SI=0.15, ΔT=5.000, NL: 0.195, NP(LIQUID)=0.944 NP(FCC_A1)=0.056
T=1530.000, X_SI=0.15, ΔT=5.000, NL: 0.184, NP(LIQUID)=0.946 NP(FCC_A1)=0.054
T=1525.000, X_SI=0.16, ΔT=5.000, NL: 0.174, NP(LIQUID)=0.947 NP(FCC_A1)=0.053
T=1520.000, X_SI=0.16, ΔT=5.000, NL: 0.165, NP(LIQUID)=0.948 NP(FCC_A1)=0.052
T=1515.000, X_SI=0.16, ΔT=5.000, NL: 0.157, NP(LIQUID)=0.949 NP(FCC_A1)=0.051
T=1510.000, X_SI=0.16, ΔT=5.000, NL: 0.149, NP(LIQUID)=0.950 NP(FCC_A1)=0.050
T=1505.000, X_SI=0.16, ΔT=5.000, NL: 0.142, NP(LIQUID)=0.951 NP(FCC_A1)=0.049
T=1500.000, X_SI=0.17, ΔT=5.000, NL: 0.135, NP(LIQUID)=0.952 NP(FCC_A1)=0.048
T=1495.000, X_SI=0.17, ΔT=5.000, NL: 0.129, NP(LIQUID)=0.953 NP(FCC_A1)=0.047
T=1490.000, X_SI=0.17, ΔT=5.000, NL: 0.123, NP(LIQUID)=0.954 NP(FCC_A1)=0.046
T=1485.000, X_SI=0.17, ΔT=5.000, NL: 0.117, NP(LIQUID)=0.955 NP(FCC_A1)=0.045
T=1480.000, X_SI=0.18, ΔT=5.000, NL: 0.112, NP(LIQUID)=0.956 NP(FCC_A1)=0.044
T=1475.000, X_SI=0.18, ΔT=5.000, NL: 0.107, NP(LIQUID)=0.956 NP(FCC_A1)=0.044
T=1470.000, X_SI=0.18, ΔT=5.000, NL: 0.103, NP(LIQUID)=0.957 NP(FCC_A1)=0.043
T=1465.000, X_SI=0.18, ΔT=5.000, NL: 0.098, NP(LIQUID)=0.958 NP(FCC_A1)=0.042
T=1460.000, X_SI=0.18, ΔT=5.000, NL: 0.094, NP(LIQUID)=0.959 NP(FCC_A1)=0.041
T=1455.000, X_SI=0.18, ΔT=5.000, NL: 0.090, NP(LIQUID)=0.959 NP(FCC_A1)=0.041
T=1450.000, X_SI=0.19, ΔT=5.000, NL: 0.087, NP(LIQUID)=0.960 NP(FCC_A1)=0.040
T=1445.000, X_SI=0.19, ΔT=5.000, NL: 0.083, NP(LIQUID)=0.960 NP(FCC_A1)=0.040
T=1440.000, X_SI=0.19, ΔT=5.000, NL: 0.080, NP(LIQUID)=0.961 NP(FCC_A1)=0.039
T=1435.000, X_SI=0.19, ΔT=5.000, NL: 0.077, NP(LIQUID)=0.961 NP(FCC_A1)=0.039
T=1430.000, X_SI=0.19, ΔT=5.000, NL: 0.074, NP(LIQUID)=0.962 NP(FCC_A1)=0.038
T=1425.000, X_SI=0.20, ΔT=5.000, NL: 0.071, NP(LIQUID)=0.963 NP(FCC_A1)=0.037
T=1420.000, X_SI=0.20, ΔT=5.000, NL: 0.069, NP(LIQUID)=0.963 NP(FCC_A1)=0.037
New phases seen: {'BETA3'}. T=1415.000, X_SI=0.20, ΔT=5.000, NL: 0.066, NP(LIQUID)=0.962 NP(BETA3)=0.003 NP(FCC_A1)=0.034
No liquid phase found at T=1410.000, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=4.167, NL: 0.066, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1410.833, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=3.472, NL: 0.065, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1411.528, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=2.894, NL: 0.065, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1412.106, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=2.411, NL: 0.064, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1412.589, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=2.009, NL: 0.064, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1412.991, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=1.674, NL: 0.064, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1413.326, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=1.395, NL: 0.063, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1413.605, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=1.163, NL: 0.063, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1413.837, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Stepping back and reducing step size.
T=1415.000, X_SI=0.20, ΔT=0.969, NL: 0.062, NP(LIQUID)=0.994 NP(BETA3)=0.003 NP(FCC_A1)=0.003
No liquid phase found at T=1414.031, X_SI=0.20. (Found {'FCC_A1', 'BETA3'}) Maximum step size reduction exceeded. Stopping.

It stopped (successfully) around 1414 K, which corresponds to the eutectic reaction for LIQUID -> FCC_A1 (disordered BETA1) + BETA3.

bocklund commented 3 years ago

With the changes in pycalphad/solve-advance branch that is in the pull request here: https://github.com/pycalphad/pycalphad/pull/373, I re-ran the Inconel 625 example with adaptive sampling enabled in the following script:

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

# Load database
dbf = Database('scratch/mc_ni_v2.034.pycalphad.tdb')
# Define the components and phases
comps = ['NI', 'CR', 'MO', 'NB', 'FE', 'SI', 'MN', 'C', 'VA']
phases = ['LIQUID','BCC_A2','FCC_A1', 'LAVES']
# Create the dictionary of conditions
mass_fracs = {v.W('CR'): 0.21,
              v.W('MO'): 0.08,
              v.W('NB'): 0.035,
              v.W('FE'): 0.04,
              v.W('SI'): 0.5/100,
              v.W('MN'):0.5/100,
              v.W('C'): 0.03/100}
conds = v.get_mole_fractions(mass_fracs, 'NI', dbf)
conds[v.P] = 1e5
conds[v.N] = 1

phase_handles, phasemap = phase_legend(phases)

from scheil import simulate_scheil_solidification

# setup the simulation parameters
initial_composition = v.get_mole_fractions(mass_fracs, 'NI', dbf)
liquid_phase_name = 'LIQUID'
start_temperature = 1650

# sample points for adaptive mode
from pycalphad import Model
from pycalphad.core.calculate import _sample_phase_constitution
from pycalphad.core.utils import point_sample
points_dict = {}
for phase_name in phases:
    mod = Model(dbf, comps, phase_name)
    points_dict[phase_name] = _sample_phase_constitution(mod, point_sample, True, pdens=50)

# perform the simulation
sol_res = simulate_scheil_solidification(dbf, comps, phases, initial_composition, start_temperature, step_temperature=5, verbose=True, eq_kwargs=dict(calc_opts=dict(points=points_dict)))

# plot the result
plt.figure(figsize=(9,6))
for phase_name, amounts in sol_res.cum_phase_amounts.items():
    plt.plot(sol_res.temperatures, amounts, label=phase_name)
plt.plot(sol_res.temperatures, sol_res.fraction_liquid, label='LIQUID')
plt.ylabel('Phase Fraction')
plt.xlabel('Temperature (K)')
plt.title('Alloy 625 Scheil simulation, phase fractions')
plt.legend(loc='best')
plt.savefig('scheil-iss14.png')

scheil-iss14

I verified that there is a miscibility gap where the NbC carbide begins to precipitate as a second FCC_A1 phase. The verbose output is a little misleading because it combines the phase fractions of the two FCC_A1 phases as one unified phase fraction. The simulation is correct, but the combined phase amounts in the printed output and the SolidificationResult object are a little difficult to use in the case of miscibility gaps.

bocklund commented 3 years ago

After #22 was merged, adaptive sampling is now on by default and the above plot can be produced by the following script when on the pycalphad:develop branch:

import matplotlib.pyplot as plt
from pycalphad import Database, equilibrium, variables as v
from scheil import simulate_scheil_solidification

# Load database
dbf = Database('scratch/mc_ni_v2.034.pycalphad.tdb')
# Define the components and phases
comps = ['NI', 'CR', 'MO', 'NB', 'FE', 'SI', 'MN', 'C', 'VA']
phases = ['LIQUID','BCC_A2','FCC_A1', 'LAVES']
# Create the dictionary of conditions
mass_fracs = {v.W('CR'): 0.21,
              v.W('MO'): 0.08,
              v.W('NB'): 0.035,
              v.W('FE'): 0.04,
              v.W('SI'): 0.5/100,
              v.W('MN'):0.5/100,
              v.W('C'): 0.03/100}
conds = v.get_mole_fractions(mass_fracs, 'NI', dbf)
conds[v.P] = 1e5
conds[v.N] = 1

# setup the simulation parameters
initial_composition = v.get_mole_fractions(mass_fracs, 'NI', dbf)
liquid_phase_name = 'LIQUID'
start_temperature = 1650

# perform the simulation
sol_res = simulate_scheil_solidification(dbf, comps, phases, initial_composition, start_temperature, step_temperature=5, verbose=True)

# plot the result
plt.figure(figsize=(9,6))
for phase_name, amounts in sol_res.cum_phase_amounts.items():
    plt.plot(sol_res.temperatures, amounts, label=phase_name)
plt.plot(sol_res.temperatures, sol_res.fraction_liquid, label='LIQUID')
plt.ylabel('Phase Fraction')
plt.xlabel('Temperature (K)')
plt.title('Alloy 625 Scheil simulation, phase fractions')
plt.legend(loc='best')
plt.savefig('scheil-iss14.png')

Closing this issue since it is resolved pending new releases of pycalphad and scheil.