Add some more tests, tweak the guess value and add another escape hatch

lyopronto/design_space.py

@@ -112,6 +112,8 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
                 Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance [cm^2-hr-Torr/g]
 
                 Tsub = fsolve(functions.T_sub_solver_FUN, T0, args = (Pch,vial['Av'],vial['Ap'],Kv,Lpr0,Lck,Rp,Tsh))[0] # Sublimation front temperature [degC]
+                # Use previous value as guess for next iteration
+                T0 = Tsub
                 dmdt = functions.sub_rate(vial['Ap'],Rp,Tsub,Pch)   # Total sublimation rate [kg/hr]
                 if dmdt<0:
                     warn(f"At t={t}hr, shelf temperature Tsh={Tsh} is too low for sublimation.")
@@ -151,8 +153,15 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
             ######################################################
 
             T_max[i_Tsh,i_Pch] = np.max(output_saved[:,1])    # Maximum product temperature [degC]
+
             drying_time[i_Tsh,i_Pch] = t    # Total drying time [hr]
             # TODO: consider whether to make this error rather than return NaN
+            if T_max[i_Tsh,i_Pch] > 0: # exceeds melting temperature, not a physically valid solution
+                warn(f"At Tsh={Tsh} and Pch={Pch}, computed temperatures of {T_max[i_Tsh,i_Pch]} exceed melting point of ice: check inputs.")
+                sub_flux_avg[i_Tsh,i_Pch] = np.nan
+                sub_flux_max[i_Tsh,i_Pch] = np.nan
+                sub_flux_end[i_Tsh,i_Pch] = np.nan
+                continue
             if output_saved.shape[0] <= 2:
                 warn(f"At Tsh={Tsh} and Pch={Pch}, drying completed in single timestep: check inputs.")
                 sub_flux_avg[i_Tsh,i_Pch] = np.nan

tests/test_design_space.py

@@ -209,6 +209,79 @@ def test_design_space_highP(self, design_space_1T1P):
             "Output should contain NaNs for infeasible conditions"
         )
 
+    def test_design_space_highT_highP(self, design_space_1T1P):
+        """Test that design space finishes execution (i.e., no infinite loops) when chamber pressure is too high for drying."""
+        vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial = design_space_1T1P
+        Tshelf["setpt"] = [20.0]  # High shelf temperature
+        Pchamber["setpt"] = [150.0]  # Very high pressure due to incorrect input units
+
+        # Expect a warning about infeasible sublimation
+        with pytest.warns(UserWarning, match="sublimation"):
+            output = design_space.dry(
+                vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial
+            )
+
+        # Calculation should complete anyway
+        assert len(output) == 3
+        assert (
+            output[0].shape[0] == 5
+        )  # [T_max, drying_time, sub_flux_avg, sub_flux_max, sub_flux_end]
+        # But should have some NaNs due to infeasibility
+        assert np.any(np.isnan(output[0])), (
+            "Output should contain NaNs for infeasible conditions in shelf T section"
+        )
+        assert np.any(np.isnan(output[1])), (
+            "Output should contain NaNs for infeasible conditions in product T section"
+        )
+        assert not np.any(np.isnan(output[2])), (
+            "Output should contain no NaNs for infeasible conditions in eq cap section"
+        )
+
+    def test_design_space_hotdry(self, design_space_1T1P):
+        """Test that design space solves correctly when chamber pressure barely exceeds
+        vapor pressure of shelf temperature, but exceeds critical product temperature."""
+        vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial = design_space_1T1P
+        Tshelf["setpt"] = [20.0]
+        # High pressure due to incorrect input units, but just barely lower than Tshelf
+        Pchamber["setpt"] = [15.0]  
+
+        with pytest.warns(UserWarning, match="sublimation"):
+            with pytest.warns(UserWarning, match="melting"):
+                output = design_space.dry(
+                    vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial
+                )
+        check_shape(output, Pchamber, Tshelf)
+        assert np.any(np.isnan(output[0])), (
+            "Output should contain NaNs for infeasible conditions in shelf T section"
+        )
+        assert np.any(np.isnan(output[1])), (
+            "Output should contain NaNs for infeasible conditions in product section"
+        )
+        assert not np.any(np.isnan(output[2])), (
+            "Output should contain no NaNs for infeasible conditions in eq cap section"
+        )
+
+    def test_design_space_subzero_eqcap(self, design_space_1T1P):
+        """Test design space with equipment capability leading to subzero sublimation."""
+        vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial = design_space_1T1P
+        Pchamber["setpt"] = [0.001]  # Pch such that a + b*Pch < 0
+
+        with pytest.warns(UserWarning, match="negative"):
+            output = design_space.dry(
+                vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial
+            )
+
+        check_shape(output, Pchamber, Tshelf)
+        assert not np.any(np.isnan(output[0])), (
+            "Output should contain no NaNs for infeasible conditions in shelf T section"
+        )
+        assert not np.any(np.isnan(output[1])), (
+            "Output should contain no NaNs for infeasible conditions in product T section"
+        )
+        assert np.any(np.isnan(output[2])), (
+            "Output should contain NaNs for infeasible conditions in eq cap section"
+        )
+
     def test_design_space_shelf_ramp_down(self, design_space_1T1P):
         """Test design space with ramp-down in shelf temperature."""
         vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial = design_space_1T1P
@@ -224,7 +297,6 @@ def test_design_space_shelf_ramp_down(self, design_space_1T1P):
     def test_design_space_no_initial_sub(self, design_space_1T1P):
         """Test design space with no sublimation at initial shelf temperature."""
         vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial = design_space_1T1P
-        # Set ramp down in shelf temperature
         Tshelf["init"] = -60.0
         Tshelf["setpt"] = [-30.0]
 
@@ -269,17 +341,6 @@ def test_design_space_fast_completion_Tsh(self, design_space_1T1P):
                 )
         check_shape(output, Pchamber, Tshelf)
 
-    def test_design_space_subzero_eqcap(self, design_space_1T1P):
-        """Test design space with equipment capability leading to subzero sublimation."""
-        vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial = design_space_1T1P
-        Pchamber["setpt"] = [0.001]  # Pch such that a + b*Pch < 0
-
-        with pytest.warns(UserWarning, match="negative"):
-            output = design_space.dry(
-                vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial
-            )
-
-        check_shape(output, Pchamber, Tshelf)
 
 
 class TestDesignSpaceComparison: