Refactor main script

lyopronto/__init__.py

@@ -5,15 +5,18 @@
 # it under the terms of the GNU General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
-  
+
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
-  
+
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+# ----------------------
+# Import submodules
+
 from . import constant
 from . import freezing
 from . import calc_knownRp
@@ -23,4 +26,32 @@
 from . import opt_Pch
 from . import opt_Tsh
 from . import functions
-from . import plot_styling
+from . import plot_styling
+
+from .high_level import (
+    execute_simulation, 
+    save_inputs_legacy,
+    save_inputs,
+    read_inputs,
+    save_csv,
+    generate_visualizations,
+)
+
+__all__ = [
+    "constant",
+    "freezing",
+    "calc_knownRp",
+    "calc_unknownRp",
+    "design_space",
+    "opt_Pch_Tsh",
+    "opt_Pch",
+    "opt_Tsh",
+    "functions",
+    "plot_styling",
+    "execute_simulation", 
+    "save_inputs_legacy",
+    "save_inputs",
+    "read_inputs", 
+    "save_csv",
+    "generate_visualizations",
+]

lyopronto/calc_knownRp.py

@@ -65,7 +65,7 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt):
         "setpoint(s). Drying cannot proceed.")
         return np.array([[0.0, Tsh_t(0), Tsh_t(0), Tsh_t(0), Pch_t(0), 0.0, 0.0]])
 
-    config = (vial, product, ht, Pch_t, Tsh_t, dt, Lpr0)
+    inputs = (vial, product, ht, Pch_t, Tsh_t, dt, Lpr0)
 
     Lck0 = [0.0]
     T0 = Tsh_t(0)
@@ -104,7 +104,7 @@ def finish(t, L):
     if sol.t[-1] == max_t:# and Lpr0 > sol.y[0, -1]:
         warn("Maximum simulation time (specified by Pchamber and Tshelf) reached before drying completion.")
 
-    output = functions.fill_output(sol, config)
+    output = functions.fill_output(sol, inputs)
 
     return output    
 

lyopronto/design_space.py

@@ -22,7 +22,39 @@
 
 ################# Primary drying at fixed set points ###############
 
+# TODO: document this properly
 def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
+    """Compute quantities necessary for constructing a graphical design space. 
+
+    Args:
+        vial (dict): _description_
+        product (dict): _description_
+        ht (dict): _description_
+        Pchamber (dict): _description_
+        Tshelf (dict): _description_
+        dt (float): _description_
+        eq_cap (dict): _description_
+        nVial (int): _description_
+
+    Returns:
+        ndarray: table of results for shelf isotherms
+        ndarray: table of results for product isotherms
+        ndarray: table of results for equipment capability curve
+
+    The first two returns have 5 rows corresponding to:
+        - Maximum product temperature in degC
+        - Primary drying time in hr
+        - Average sublimation flux in kg/hr/m^2
+        - Maximum/minimum sublimation flux in kg/hr/m^2
+        - Sublimation flux at the end of primary drying in kg/hr/m^2
+    The third return has 3 rows corresponding to the first three of that list.
+
+    With nT setpoints in Tshelf['setpt'] and nP setpoints in Pchamber['setpt'], the returned arrays have the following shapes:
+        - Shelf isotherms: (5, nT, nP) array
+        - Product isotherms: (5, 2) array (for the lowest and highest Pchamber setpoints)
+        - Equipment capability curve: (3, nP) array
+
+    """
 
     T_max = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])
     drying_time = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])

lyopronto/freezing.py

@@ -124,19 +124,21 @@ def freeze(vial,product,h_freezing,Tshelf,dt):
 
     t_last = ts
     Tpr0 = Tpr    # degC
+    Tsh0 = Tsh
     while(t<t_tr[-1]):
 
         i = np.argmax(t_tr>t) # Get first index where time trigger exceeds current time
         if not(i == i_prev):
             i_prev = i
             t_last = t_tr[i-1]
             Tpr0 = Tpr
+            Tsh0 = Tsh
 
         # Evaluate shelf temperature at current time point 
         Tsh = Tshr(t)    # degC
 
         # Product temperature
-        Tpr = functions.lumped_cap_Tpr_ice(t-t_last,Tpr0,V_frozen,h_freezing,vial['Av'],Tsh,Tsh_tr[i-1],r[i])
+        Tpr = functions.lumped_cap_Tpr_ice(t-t_last,Tpr0,V_frozen,h_freezing,vial['Av'],Tsh,Tsh0,r[i])
         # Update record as functions of the cycle time
         freezing_output_saved = np.append(freezing_output_saved, [[t, Tsh, Tpr]],axis=0)
 

lyopronto/functions.py

@@ -328,28 +328,26 @@ def crystallization_time_FUN(V,h,Av,Tf,Tn,Tsh_func,t0):
     hA = h*constant.hr_To_s * Av*constant.cm_To_m**2 # heat transfer coefficient in J/K/hr
     # t = rhoV*(Hf-Cp*(Tf-Tn))/hA/(Tf-Tsh) # time: g*(J/g- J/g/K*K)/(J/m^2/K/hr*m^2*K) = hr
     lhs = rhoV*(Hf-Cp*(Tf-Tn))/hA
-    def integrand(t):
-        return Tf - Tsh_func(t+t0)
-    def resid(t):
-        integral, _ = quad(integrand, 0, t)
+    def integrand(dt):
+        return Tf - Tsh_func(t0+dt)
+    def resid(delta_t):
+        integral, _ = quad(integrand, 0, delta_t)
         return integral - lhs
-    t = brentq(resid, t0, t0+100.0)
-
-
-    return t
+    delta_t = brentq(resid, 0, 100.0)
+    return delta_t
 
 ##
 
 
 
 ################################################################
-def calc_step(t, Lck, config):
+def calc_step(t, Lck, inputs):
     """Calculate the full set of system states at a given time step from ODE solution states.
 
     Args:
         t (float): The current time in hours.
         Lck (float): The cake thickness in cm.
-        config (tuple): A tuple containing the configuration parameters.
+        inputs (tuple): A tuple containing the inputs parameters.
 
     Returns:
         (np.ndarray): The full set of system states at the given time step:
@@ -361,7 +359,7 @@ def calc_step(t, Lck, config):
             5. Sublimation flux [kg/hr/m²],
             6. Drying percent [%]
     """
-    vial, product, ht, Pch_t, Tsh_t, dt, Lpr0 = config
+    vial, product, ht, Pch_t, Tsh_t, dt, Lpr0 = inputs
     Tsh = Tsh_t(t)
     Pch = Pch_t(t)
     Kv = Kv_FUN(ht['KC'],ht['KP'],ht['KD'],Pch)  # Vial heat transfer coefficient in cal/s/K/cm^2
@@ -379,12 +377,12 @@ def calc_step(t, Lck, config):
     col = np.array([t, Tsub, Tbot, Tsh, Pch*constant.Torr_to_mTorr, dmdt/(vial['Ap']*constant.cm_To_m**2), dry_percent])
     return col
 
-def fill_output(sol, config):
+def fill_output(sol, inputs):
     """Fill the output array with the results from the ODE solver.
 
     Args:
         sol (ODESolution): The solution object returned by the ODE solver.
-        config (tuple): A tuple containing the configuration parameters.
+        inputs (tuple): A tuple containing the input parameters.
 
     Returns:
         (np.ndarray): The output array filled with the results from the ODE solver.
@@ -393,11 +391,11 @@ def fill_output(sol, config):
     of points is impractical. Instead, we calculate at the the ODE solver points, and 
     interpolate elsewhere.
     """
-    dt = config[5]
+    dt = inputs[5]
 
     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, config)
+        interp_points[i,:] = calc_step(t, y, inputs)
     # out_t = np.arange(0, sol.t[-1], dt)   
     if dt is None:
         return interp_points