Publication prep

README.md

@@ -1,5 +1,6 @@
 # LyoPronto.jl
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://LyoHUB.github.io/LyoPronto.jl/dev)
+[![](https://zenodo.org/badge/DOI/10.5281/zenodo.17373491.svg)](http://doi.org/10.5281/zenodo.17373491)
 
 This package is a Julia complement to [LyoPRONTO](https://github.com/LyoHUB/LyoPronto), an open source Python package.
 It has some overlapping functionality with LyoPRONTO, especially simulation of primary drying for conventional lyophilization.

src/pikal_model.jl

@@ -15,30 +15,31 @@ const end_drying_callback = ContinuousCallback(end_cond, terminate!, save_positi
 # -------------------------------------------
 # Incorporate the nonlinear algebraic part in a DAE formulation.
 # This has the advantage that, afterward, temperatures can be cheaply interpolated by builtin solutions
+"""
+    $(SIGNATURES)
 
-@doc raw"""
-    lyo_1d_dae!(du, u, params, t)
+With the Pikal model, compute the mass flow at `u=[hf, Tf]`, time `t`, and conditions `po`.
+Ideally `po` should be a `ParamObjPikal`.
 
-Internal implementation of the Pikal model.
-See [`lyo_1d_dae_f`](@ref) for the wrapped version, which is more fully documented.
+This gets used inside the DAE solve, but for consistency is provided separately so to avoid
+the risk of making unit mistakes.
 """
-function lyo_1d_dae!(du, u, params, t)
+function calc_md_Q(u, po, t)
 
-    if params isa ParamObj
+    if po isa ParamObj
         (; Rp, hf0, csolid, ρsolution,
-        Kshf, Av, Ap, pch, Tsh) = params
+        Kshf, Av, Ap, pch, Tsh) = po
     else
-        Rp, hf0, csolid, ρsolution = params[1]
-        Kshf, Av, Ap = params[2]
-        pch, Tsh = params[3]
+        Rp, hf0, csolid, ρsolution = po[1]
+        Kshf, Av, Ap = po[2]
+        pch, Tsh = po[3]
     end
-    
+
     td = t*u"hr" # Dimensional time
     hf = u[1]*u"cm"
     Tf = u[2]*u"K"
-    if Tf < 0.0u"K"
-        du .= NaN
-        return nothing
+    if Tf < 0.0u"K" # Unphysical temperature needs an escape hatch
+        return NaN*u"kg/s", NaN*u"W"
     end
 
     hd = hf0 - hf
@@ -47,6 +48,31 @@ function lyo_1d_dae!(du, u, params, t)
     Tsub = Tf - Qshf/k_ice/Ap*hf
     delta_p = calc_psub(Tsub)-pch(td)
     dmdt = - Ap*(delta_p)/Rp(hd)
+    return dmdt, Qshf
+end
+
+@doc raw"""
+    lyo_1d_dae!(du, u, params, t)
+
+Internal implementation of the Pikal model.
+See [`lyo_1d_dae_f`](@ref) for the wrapped version, which is more fully documented.
+"""
+function lyo_1d_dae!(du, u, params, t)
+
+    # Need a handful of parameters in this function.
+    if params isa ParamObj
+        (; csolid, ρsolution, Ap) = params
+    else
+        csolid, ρsolution = params[1][3:4]
+        Ap = params[2][3]
+    end
+    # This logic is carried out in a separate function,
+    # so that it can be reused after the fact for computing mass flow.
+    dmdt, Qshf = calc_md_Q(u, params, t)
+    if isnan(dmdt)
+        du .= NaN
+        return nothing
+    end
     Qsub = uconvert(u"W", dmdt*ΔHsub)
 
     dhf_dt = min(0.0u"cm/hr", dmdt/(ρsolution-csolid)/Ap |> u"cm/hr") # Cap dhf_dt at 0: no desublimation
@@ -214,7 +240,6 @@ function dae_Rp!(du, u, p, tn)
 
     Tf = Tf_interp(t)
 
-
     Q = Kshf(pch(t))*Av*(Tsh(t) - Tf)
     Tsub = Tf - Q/Ap/LyoPronto.k_ice * (hf0-hd)
     md = Q/LyoPronto.ΔH

src/public.jl

@@ -7,4 +7,5 @@ public σ
 public e_0, ϵppf, eppf, εppf, ϵpp_gl, epp_gl, εpp_gl
 public ρ_sucrose, k_sucrose
 
-public extract_ts
+public extract_ts
+public calc_md_Q

src/recipes.jl

@@ -389,11 +389,11 @@ end
 @recipe function f(bsp::BarStackPlot)
     x, ys... = bsp.args
     yprev = zeros(eltype(ys[1]), size(ys[1], 1))
-    ymat = hcat(yprev, ys...)
+    ymat = hcat(yprev, reverse(ys)...)
     ys_cum = cumsum(ymat, dims=2)
     ys_cum = ys_cum ./ maximum(ys_cum, dims=2)
-    pal = [:yellow, :orange, :red]
-    for i in axes(ys_cum, 2)[begin:end-1]
+    pal = [:red, :orange, :yellow]
+    for i in axes(ys_cum, 2)[end-1:-1:begin]
         @series begin
             seriestype := :bar
             fillrange := ys_cum[:,i]