FIX: Raise ValueError if there is no initial liquid in solidification simulations

scheil/simulate.py

@@ -198,7 +198,7 @@ def simulate_scheil_solidification(dbf, comps, phases, composition,
                 print('done')
 
     converged = False
-    phases_seen = {liquid_phase_name, ''}
+    phases_seen = {''}
     liquid_comp = composition
     wks = Workspace(dbf, comps, phases, calc_opts=eq_kwargs.get("calc_opts"))
     last_converged_wks = None
@@ -230,6 +230,10 @@ def simulate_scheil_solidification(dbf, comps, phases, composition,
                 if verbose:
                     print(f'(Convergence failure) ', end='')
             if T_STEP_ORIG / step_temperature > MAXIMUM_STEP_SIZE_REDUCTION:
+                if liquid_phase_name not in phases_seen:
+                    if len(found_ph) == 0:
+                        raise ValueError(f"No equilibrium calculations converged near the starting temperature ({start_temperature} K). The starting temperature must be at or above the liquidus.")
+                    raise ValueError(f"No stable liquid phase ({liquid_phase_name}) found at the starting temperature ({start_temperature} K); found phases {found_ph}. The starting temperature must be at or above the liquidus.")
                 # Only found solid phases and the step size has already been reduced. Stop running without converging.
                 if verbose:
                     print('Maximum step size reduction exceeded. Stopping.')
@@ -431,6 +435,13 @@ def simulate_equilibrium_solidification(dbf, comps, phases, composition,
             temperatures.append(current_T)
             current_T -= step_temperature
         else:
+            if len(temperatures) == 0:
+                # Liquid was not stable at the starting temperature, so there is no liquid-bearing
+                # upper bound for the binary search and no solidification path to compute.
+                found_ph = stable_phases - {''}
+                if len(found_ph) == 0:
+                    raise ValueError(f"No equilibrium calculations converged at the starting temperature ({start_temperature} K). The starting temperature must be at or above the liquidus.")
+                raise ValueError(f"No stable liquid phase ({liquid_phase_name}) found at the starting temperature ({start_temperature} K); found phases {found_ph}. The starting temperature must be at or above the liquidus.")
             # binary search to find the solidus
             T_high = current_T + step_temperature  # High temperature, liquid
             T_low = current_T  # Low temperature, solids only

tests/fixtures.py

@@ -0,0 +1,30 @@
+import os
+import pytest
+from pycalphad import Database
+
+DATABASE_DIR = os.path.join(os.path.dirname(__file__), 'databases')
+
+
+@pytest.fixture(scope="session")
+def load_database(request):
+    """
+    Helper fixture to load a database (parameterized by the value of `request`).
+    """
+    db = Database(os.path.join(DATABASE_DIR, request.param))
+    def _load_database():
+        return db
+    return _load_database
+
+
+def select_database(path):
+    """
+    Decorator to facilitate safe, fast loading of database objects. Use as
+
+    ```
+    @select_database("filename.tdb")  # matches a file in the tests/databases directory
+    def test_name_of_my_test(load_database):
+        dbf = load_database()  # equivalent to `dbf = Database("filename.tdb")`
+        # ... implement test below
+    ```
+    """
+    return pytest.mark.parametrize("load_database", [path], indirect=True)

tests/test_chen2009.py

@@ -5,17 +5,18 @@
 from pycalphad import Database, variables as v
 from scheil import simulate_scheil_solidification
 import pytest
+from fixtures import select_database, load_database, DATABASE_DIR
 
-DB_CHEN = Database(os.path.join(os.path.dirname(__file__), 'sl_chen.tdb'))
 
-
-def test_binary_A_B():
+@select_database("sl_chen.tdb")
+def test_binary_A_B(load_database):
     """Tests for the Scheil properties of the A-B binary peritectic system"""
+    dbf = load_database()
     comp = {v.X('B'): 0.5}
     start = 800  # Kelvin
     step = 5.0
 
-    sol_res = simulate_scheil_solidification(DB_CHEN, ['A', 'B'], ['ALPHA', 'BETA', 'LIQUID'], comp, start, step_temperature=step, stop=1e-8)
+    sol_res = simulate_scheil_solidification(dbf, ['A', 'B'], ['ALPHA', 'BETA', 'LIQUID'], comp, start, step_temperature=step, stop=1e-8)
     print(f"Converged to stopping criteria: {sol_res.converged}")
 
     phase_amnts = sol_res.phase_amounts
@@ -53,7 +54,7 @@ def test_binary_A_B():
 )
 def test_binary_A_B_with_ordering(phases):
     """Tests for the Scheil properties of the A-B binary peritectic system with an ordered phase ORD_BETA"""
-    with open(os.path.join(os.path.dirname(__file__), 'sl_chen.tdb')) as fp:
+    with open(os.path.join(DATABASE_DIR, 'sl_chen.tdb')) as fp:
         tdb_str = fp.read()
     # Append an ordered phase, 'ORD_BETA' with no ordering parameters
     ordered_phase_str = """
@@ -98,13 +99,15 @@ def test_binary_A_B_with_ordering(phases):
     assert np.isclose(sol_res.temperatures[idx_first_beta], 647, atol=step * 1.1)
 
 
-def test_binary_A_C():
+@select_database("sl_chen.tdb")
+def test_binary_A_C(load_database):
     """Tests for the Scheil properties of the A-C binary eutectic system"""
+    dbf = load_database()
     comp = {v.X('C'): 0.5}
     start = 900  # Kelvin
     step = 5.0
 
-    sol_res = simulate_scheil_solidification(DB_CHEN, ['A', 'C'], ['ALPHA', 'BETA', 'LIQUID'], comp, start, step_temperature=step, stop=1e-8)
+    sol_res = simulate_scheil_solidification(dbf, ['A', 'C'], ['ALPHA', 'BETA', 'LIQUID'], comp, start, step_temperature=step, stop=1e-8)
     print(f"Converged to stopping criteria: {sol_res.converged}")
 
     phase_amnts = sol_res.phase_amounts
@@ -126,19 +129,21 @@ def test_binary_A_C():
         assert np.isclose(sol_res.temperatures[-1], 820, atol=step * 1.1)
 
 
-def test_ternary_A_B_C():
+@select_database("sl_chen.tdb")
+def test_ternary_A_B_C(load_database):
     """Tests for the Scheil properties of the A-B-C ternary system
 
     There is a transition point, t, where the solidification changes from the
     eutectic monovariant to peritetic-like. Most of these tests are focused
     around this point.
     """
     t_temp = 683.7  # Eutectic-like to peritectic-like transition temperature, in Kelvin
+    dbf = load_database()
     comp = {v.X('B'): 0.25, v.X('C'): 0.1}
     start = 800  # Kelvin
     step = 5
 
-    sol_res = simulate_scheil_solidification(DB_CHEN, ['A', 'B', 'C'], ['ALPHA', 'BETA', 'LIQUID'], comp, start, step_temperature=step, stop=1e-8)
+    sol_res = simulate_scheil_solidification(dbf, ['A', 'B', 'C'], ['ALPHA', 'BETA', 'LIQUID'], comp, start, step_temperature=step, stop=1e-8)
     print(f"Converged to stopping criteria: {sol_res.converged}")
 
     phase_amnts = sol_res.phase_amounts

tests/test_scheil_solidification.py

@@ -1,18 +1,36 @@
-import os
 import json
 import numpy as np
-from pycalphad import Database, variables as v
+import pytest
+from pycalphad import variables as v
 from scheil import simulate_scheil_solidification, simulate_equilibrium_solidification, SolidificationResult
+from fixtures import select_database, load_database
 
+@select_database("sl_chen.tdb")
+def test_scheil_solidification_raises_without_initial_liquid(load_database):
+    """Scheil solidification should raise a clear error if no liquid is stable at the start temperature."""
+    dbf = load_database()
+    # For X(B)=0.5 the A-B system is fully solid below 600 K, so 500 K has no stable liquid
+    with pytest.raises(ValueError, match="[Ll]iquid"):
+        simulate_scheil_solidification(dbf, ['A', 'B'], ['ALPHA', 'BETA', 'LIQUID'], {v.X('B'): 0.5}, 500.0, step_temperature=5.0)
 
-def test_scheil_solidification_result_properties():
+
+@select_database("sl_chen.tdb")
+def test_equilibrium_solidification_raises_without_initial_liquid(load_database):
+    """Equilibrium solidification should raise a clear error if no liquid is stable at the start temperature."""
+    dbf = load_database()
+    with pytest.raises(ValueError, match="[Ll]iquid"):
+        simulate_equilibrium_solidification(dbf, ['A', 'B'], ['ALPHA', 'BETA', 'LIQUID'], {v.X('B'): 0.5}, 500.0, step_temperature=5.0)
+
+
+@select_database("alzn_mey.tdb")
+def test_scheil_solidification_result_properties(load_database):
     """Test that SolidificationResult objects produced by Scheil solidification have the required properties."""
     # Required properties are
     # 1. the shape of the output arrays are matching
     # 2. the final fraction of solid is 1.0 and fraction of liquid is 0.0
     # 3. the sum total of the final (solid) phase amounts is 1.0
 
-    dbf = Database(os.path.join(os.path.dirname(__file__), 'alzn_mey.tdb'))
+    dbf = load_database()
     comps = ['AL', 'ZN', 'VA']
     phases = sorted(dbf.phases.keys())
 
@@ -60,9 +78,10 @@ def test_scheil_solidification_result_properties():
 
 
 
-def test_scheil_solidification_custom_properties():
+@select_database("alzn_mey.tdb")
+def test_scheil_solidification_custom_properties(load_database):
     """Test that SolidificationResult objects produced by Scheil solidification calculations with custom outputs have the right shape."""
-    dbf = Database(os.path.join(os.path.dirname(__file__), 'alzn_mey.tdb'))
+    dbf = load_database()
     comps = ['AL', 'ZN', 'VA']
     phases = sorted(dbf.phases.keys())
 
@@ -99,21 +118,22 @@ def test_scheil_solidification_custom_properties():
     assert set(sol_res.output.keys()) == set(rnd_trip_sol_res.output.keys())
     for ky in sol_res.output.keys():
         np.testing.assert_almost_equal(rnd_trip_sol_res.output[ky], sol_res.output[ky])
-    
+
     # Test to_dataframe doesn't raise
     sol_res.to_dataframe(include_zero_phases=True)
     df = sol_res.to_dataframe(include_zero_phases=False)
     assert "CPM" in df.columns
     assert "HM" in df.columns
 
 
-def test_equilibrium_solidification_result_properties():
+@select_database("alzn_mey.tdb")
+def test_equilibrium_solidification_result_properties(load_database):
     """Test that SolidificationResult objects produced by equilibrium have the required properties."""
     # Required properties are that the shape of the output arrays are matching
     # NOTE: final phase amounts are not tested because they are not guaranteed
     # to be 0.0 or 1.0 in the same way as in the Scheil simulations.
 
-    dbf = Database(os.path.join(os.path.dirname(__file__), 'alzn_mey.tdb'))
+    dbf = load_database()
     comps = ['AL', 'ZN', 'VA']
     phases = sorted(dbf.phases.keys())
 
@@ -156,9 +176,10 @@ def test_equilibrium_solidification_result_properties():
     sol_res.to_dataframe(include_zero_phases=False)
 
 
-def test_equilibrium_solidification_custom_properties():
+@select_database("alzn_mey.tdb")
+def test_equilibrium_solidification_custom_properties(load_database):
     """Test that SolidificationResult objects produced by equilibrium solidification calculations with custom outputs have the right shape."""
-    dbf = Database(os.path.join(os.path.dirname(__file__), 'alzn_mey.tdb'))
+    dbf = load_database()
     comps = ['AL', 'ZN', 'VA']
     phases = sorted(dbf.phases.keys())