FIX/TST: equilibrium: Phases are correctly suspended when a subsystem…

… is specified that causes the phase to have zero feasible points. Fixes pycalphadgh-503

pycalphad/core/calculate.py

@@ -173,7 +173,7 @@ def _sample_phase_constitution(model, sampler, fixed_grid, pdens, phase_local_co
                 assert np.max(np.abs(constraint_jac.dot(points.T).T - constraint_rhs)) < 1e-6
             else:
                 # No feasible points; return array of nan to preserve shape
-                return np.full((num_points, points.shape[1]), np.nan)
+                return np.full((num_points, constraint_jac.shape[1]), np.nan)
         else:
             points = np.concatenate((points, sampler(sublattice_dof, pdof=pdens)))
 

pycalphad/core/eqsolver.pyx

@@ -112,6 +112,9 @@ def add_nearly_stable(object composition_sets, object phase_records,
             continue
         phase_record = phase_records[phase_name]
         phase_indices = grid.attrs['phase_indices'][phase_name]
+        if phase_indices.start == phase_indices.stop:
+            # Phase has zero feasible grid points to consider
+            continue
         driving_forces_for_phase = driving_forces[phase_indices.start:phase_indices.stop]
         minimum_df_idx = argmax(&driving_forces_for_phase[0], driving_forces_for_phase.shape[0])
         if driving_forces_for_phase[minimum_df_idx] >= minimum_df:

pycalphad/tests/test_equilibrium.py

@@ -1013,6 +1013,33 @@ def test_eq_charge_ndzro(load_database):
     assert np.allclose(Y_PYRO, [9.99970071e-01, 2.99288042e-05, 3.83395063e-02, 9.61660494e-01, 9.93381787e-01,
                                 6.61821340e-03, 1.00000000e+00, 1.39970285e-03, 9.98600297e-01], rtol=5e-4)
 
+@pytest.mark.solver
+def test_issue_503_charged_infeasible_subsystem():
+    "equilibrium suspends a phase with zero feasible points due to internal constraints"
+    tdb = """
+ ELEMENT /-   ELECTRON_GAS              0.0000E+00  0.0000E+00  0.0000E+00!
+ ELEMENT VA   VACUUM                    0.0000E+00  0.0000E+00  0.0000E+00!
+ ELEMENT O    1/2_MOLE_O2(G)            1.5999E+01  4.3410E+03  1.0252E+02!
+ ELEMENT V    BCC_A2                    5.0941E+01  4.5070E+03  3.0890E+01!
+
+ SPECIES O-2                         O1/-2!
+ SPECIES O2                          O2!
+ SPECIES O3                          O3!
+ SPECIES V+2                         V1/+2!
+ SPECIES V+3                         V1/+3!
+
+ PHASE GAS:G %  1  1.0  !
+    CONSTITUENT GAS:G :O,O2,O3 :  !
+
+ PHASE HALITE %  2 1   1 !
+    CONSTITUENT HALITE :V,V+2,V+3,VA : O-2,VA :  !
+"""
+    dbf = Database(tdb)
+    result = equilibrium(dbf, ['O', 'VA'], ['GAS', 'HALITE'], {v.T: 1000, v.P: 1e5})
+    print(result)
+    assert np.all(np.isclose(result.NP.squeeze(), [1.0, np.nan], equal_nan=True))
+    assert np.all(result.Phase.squeeze() == ["GAS", ""])
+
 @pytest.mark.solver
 @select_database("crtiv_ghosh.tdb")
 def test_ternary_three_phase_dilute(load_database):