diff --git a/lyopronto/functions.py b/lyopronto/functions.py index 53a76a8..7abe061 100644 --- a/lyopronto/functions.py +++ b/lyopronto/functions.py @@ -400,6 +400,8 @@ def fill_output(sol, inputs): """ dt = inputs[5] + Pch_t, Tsh_t = inputs[3], inputs[4] + interp_points = np.zeros((len(sol.t), 7)) for i,(t, y) in enumerate(zip(sol.t, sol.y[0])): interp_points[i,:] = calc_step(t, y, inputs) @@ -415,4 +417,6 @@ def fill_output(sol, inputs): fullout[i,:] = interp_points[sol.t == t, :] else: fullout[i,:] = interp_func(t) + fullout[i, 3] = Tsh_t(t) # Ensure shelf temp is exact + fullout[i, 4] = Pch_t(t)*constant.Torr_to_mTorr # Ensure chamber pressure is exact return fullout diff --git a/tests/test_calc_knownRp.py b/tests/test_calc_knownRp.py index f081eee..51fc72d 100644 --- a/tests/test_calc_knownRp.py +++ b/tests/test_calc_knownRp.py @@ -2,7 +2,7 @@ import pytest import numpy as np -from lyopronto import calc_knownRp, constant +from lyopronto import calc_knownRp, constant, functions from lyopronto.high_level import execute_simulation from .utils import ( assert_physically_reasonable_output, @@ -274,6 +274,33 @@ def test_high_resistance_product(self, knownRp_standard_setup): # Note: May not take >20 hours depending on other parameters assert_physically_reasonable_output(output) + def test_sharp_corners(self, knownRp_standard_setup): + """Test that solve_ivp interpolation does not smooth out shelf temperature profile.""" + vial, _, ht, __, ___, ____ = knownRp_standard_setup + # Larger resistance is not crucial, but this example shows some smoothing + product = {"R0": 0.5, "A1": 40.0, "A2": 2.0, "cSolid": 0.05} + # More steps -> more smoothing + Tshelf = { + "init": -45.0, + "setpt": [10.0, 40.0], + "dt_setpt": [180.0, 120.0], + "ramp_rate": 1.0, + } + Pchamber = {"setpt": [0.005], "dt_setpt": [1800.0], "ramp_rate": 0.5} + + output = calc_knownRp.dry(vial, product, ht, Pchamber, Tshelf, 0.01) + + assert_complete_drying(output) + assert_physically_reasonable_output(output) + ri = functions.RampInterpolator(Tshelf) + for i, t in enumerate(output[:, 0]): + expected_Tsh = ri(t) + actual_Tsh = output[i, 3] + assert actual_Tsh == pytest.approx(expected_Tsh, abs=1e-4), ( + f"At time {t:.1f} hr, expected shelf temp {expected_Tsh:.2f} °C, " + f"but got {actual_Tsh:.2f} °C" + ) + class TestRegression: """