Raise error when temperature below melting temperature

.github/workflows/pytest.yaml

@@ -11,7 +11,7 @@ jobs:
       fail-fast: false
       max-parallel: 100
       matrix:
-        python-version: [3.7, 3.8, 3.9]
+        python-version: [3.8, 3.9]
         pycalphad_develop_version: [true, false]
 
     steps:

scheil/simulate.py

@@ -1,5 +1,6 @@
 import sys
 import numpy as np
+from collections import defaultdict
 from pycalphad import equilibrium, variables as v
 from pycalphad.codegen.callables import build_phase_records
 from pycalphad.core.calculate import _sample_phase_constitution
@@ -107,10 +108,17 @@ def simulate_scheil_solidification(dbf, comps, phases, composition,
     solid_phases = sorted((set(phases) | filtered_disordered_phases) - {liquid_phase_name})
     temp = start_temperature
     independent_comps = sorted([str(comp)[2:] for comp in composition.keys()])
+
     x_liquid = {comp: [composition[v.X(comp)]] for comp in independent_comps}
     fraction_solid = [0.0]
     temperatures = [temp]
     phase_amounts = {ph: [0.0] for ph in solid_phases}
+    x_phases = defaultdict(dict)
+    Y_phases = defaultdict(dict)
+    for ph in solid_phases:
+        Y_phases[ph] = [np.nan]
+        for comp in independent_comps:
+            x_phases[ph][comp] = [np.nan]
 
     if adaptive:
         dof_dict = {phase_name: list(map(len, mod.constituents)) for phase_name, mod in models.items()}
@@ -183,12 +191,22 @@ def simulate_scheil_solidification(dbf, comps, phases, composition,
         for solid_phase in solid_phases:
             if solid_phase not in eq_phases:
                 phase_amounts[solid_phase].append(0.0)
+                Y_phases[solid_phase].append(np.nan)
+                for comp in independent_comps:
+                    x_phases[solid_phase][comp].append(np.nan)
                 continue
             np_tieline = np.nansum(eq.isel(vertex=eq_phases.index(solid_phase))["NP"].values.squeeze())
             found_phase_amounts.append((solid_phase, np_tieline))
             delta_fraction_solid = (1 - current_fraction_solid) * np_tieline
             current_fraction_solid += delta_fraction_solid
             phase_amounts[solid_phase].append(delta_fraction_solid)
+            vertex = sorted(np.nonzero(eq["Phase"].values.squeeze().flat == solid_phase))[0]
+            Y_phases[solid_phase].append(list(eq["Y"].isel(vertex=vertex).squeeze().values))
+            for comp in independent_comps:
+                if np.sum(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)==0:
+                    x_phases[solid_phase][comp].append(np.nan)
+                else:
+                    x_phases[solid_phase][comp].append(float(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values))
         fraction_solid.append(current_fraction_solid)
         temperatures.append(temp)
         NL = 1 - fraction_solid[-1]
@@ -217,11 +235,22 @@ def simulate_scheil_solidification(dbf, comps, phases, composition,
         for solid_phase in solid_phases:
             if solid_phase in eq_phases:
                 amount = np.nansum(eq.isel(vertex=eq_phases.index(solid_phase))["NP"].values.squeeze())
-                phase_amounts[solid_phase].append(float(amount) * (1 - current_fraction_solid))
+                try:
+                    phase_amounts[solid_phase].append(float(amount) * (1 - current_fraction_solid))
+                except UnboundLocalError:
+                    phases_solid = set(eq_phases) - {''}
+                    raise ValueError("The input temperature is too low to support a liquid phase, it only contains "+', '.join(phases_solid))
             else:
                 phase_amounts[solid_phase].append(0.0)
+            vertex = sorted(np.nonzero(eq["Phase"].values.squeeze().flat == solid_phase))[0]
+            Y_phases[solid_phase].append(list(eq["Y"].isel(vertex=vertex).squeeze().values))
+            for comp in independent_comps:
+                if np.sum(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)==0:
+                    x_phases[solid_phase][comp].append(np.nan)
+                else:
+                    x_phases[solid_phase][comp].append(float(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values))
 
-    return SolidificationResult(x_liquid, fraction_solid, temperatures, phase_amounts, converged, "scheil")
+    return SolidificationResult(x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, "scheil")
 
 
 def simulate_equilibrium_solidification(dbf, comps, phases, composition,
@@ -291,6 +320,12 @@ def simulate_equilibrium_solidification(dbf, comps, phases, composition,
     fraction_solid = []
     phase_amounts = {ph: [] for ph in solid_phases}  # instantaneous phase amounts
     cum_phase_amounts = {ph: [] for ph in solid_phases}
+    x_phases = defaultdict(dict)
+    Y_phases = defaultdict(dict)
+    for ph in solid_phases:
+        Y_phases[ph] = []
+        for comp in independent_comps:
+            x_phases[ph][comp] = []
     converged = False
     current_T = start_temperature
     if verbose:
@@ -368,8 +403,16 @@ def simulate_equilibrium_solidification(dbf, comps, phases, composition,
                 phase_amounts[solid_phase].append(amount)
             else:
                 phase_amounts[solid_phase].append(amount - cum_phase_amounts[solid_phase][-2])
+            vertex = sorted(np.nonzero(eq["Phase"].values.squeeze().flat == solid_phase))[0]
+            Y_phases[solid_phase].append(list(eq["Y"].isel(vertex=vertex).squeeze().values))
+            for comp in independent_comps:
+                if np.sum(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)==0:
+                    x_phases[solid_phase][comp].append(np.nan)
+                else:
+                    x_phases[solid_phase][comp].append(float(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values))
             current_fraction_solid += amount
+
         fraction_solid.append(current_fraction_solid)
 
     converged = True if np.isclose(fraction_solid[-1], 1.0) else False
-    return SolidificationResult(x_liquid, fraction_solid, temperatures, phase_amounts, converged, "equilibrium")
+    return SolidificationResult(x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, "equilibrium")

scheil/solidification_result.py

@@ -38,8 +38,10 @@ class SolidificationResult():
 
     """
 
-    def __init__(self, x_liquid, fraction_solid, temperatures, phase_amounts, converged, method):
+    def __init__(self, x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, method):
         self.x_liquid = x_liquid
+        self.x_phases = x_phases
+        self.Y_phases = Y_phases
         self.fraction_solid = fraction_solid
         self.fraction_liquid = (1.0 - np.array(fraction_solid)).tolist()
         self.temperatures = temperatures
@@ -57,6 +59,8 @@ def __repr__(self):
     def to_dict(self):
         d = {
             'x_liquid': self.x_liquid,
+            'x_phases': self.x_phases,
+            'Y_phases': self.Y_phases, 
             'fraction_solid': self.fraction_solid,
             'temperatures': self.temperatures,
             'phase_amounts': self.phase_amounts,
@@ -68,9 +72,11 @@ def to_dict(self):
     @classmethod
     def from_dict(cls, d):
         x_liquid = d['x_liquid']
+        x_phases = d['x_phases']
+        Y_phases = d['Y_phases']
         fraction_solid = d['fraction_solid']
         temperatures = d['temperatures']
         phase_amounts = d['phase_amounts']
         converged = d['converged']
         method = d['method']
-        return cls(x_liquid, fraction_solid, temperatures, phase_amounts, converged, method)
+        return cls(x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, method)

tests/test_scheil_solidification.py

@@ -26,6 +26,9 @@ def test_scheil_solidification_result_properties():
     assert num_temperatures == len(sol_res.fraction_solid)
     assert all([num_temperatures == len(nphase) for nphase in sol_res.phase_amounts.values()])
     assert all([num_temperatures == len(liq_comps) for liq_comps in sol_res.x_liquid.values()])
+    for phase in sol_res.x_phases:
+        assert all([num_temperatures == len(phase_comps) for phase_comps in sol_res.x_phases[phase].values()])
+    assert all([num_temperatures == len(phase_sf) for phase_sf in sol_res.Y_phases.values()])    
     assert all([num_temperatures == len(nphase) for nphase in sol_res.cum_phase_amounts.values()])
 
     # final phase amounts are correct
@@ -48,6 +51,8 @@ def test_scheil_solidification_result_properties():
     assert rnd_trip_sol_res.fraction_liquid == sol_res.fraction_liquid
     assert rnd_trip_sol_res.fraction_solid == sol_res.fraction_solid
     assert rnd_trip_sol_res.x_liquid == sol_res.x_liquid
+    assert rnd_trip_sol_res.x_phases == sol_res.x_phases
+    assert rnd_trip_sol_res.Y_phases == sol_res.Y_phases
     assert rnd_trip_sol_res.cum_phase_amounts == sol_res.cum_phase_amounts
     assert rnd_trip_sol_res.phase_amounts == sol_res.phase_amounts
     assert rnd_trip_sol_res.temperatures == sol_res.temperatures