migration of tritium module from core pathsim

src/pathsim_chem/__init__.py

@@ -11,3 +11,6 @@
     __version__ = "unknown"
 
 __all__ = ["__version__"]
+
+#for direct block import from main package
+from .tritium import *

src/pathsim_chem/tritium/README.md

@@ -0,0 +1,3 @@
+# Tritium Toolbox
+
+This module contains blocks for nuclear fusion modeling. Especially for tritium cycle modeling.

src/pathsim_chem/tritium/__init__.py

@@ -0,0 +1,4 @@
+from .residencetime import *
+from .splitter import *
+from .bubbler import *
+# from .tcap import *

src/pathsim_chem/tritium/bubbler.py

@@ -0,0 +1,218 @@
+#########################################################################################
+##
+##                                      Bubbler Block
+##                                (blocks/fusion/bubbler.py)
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+
+from ..dynsys import DynamicalSystem
+from ...events.schedule import ScheduleList
+
+
+# BLOCK DEFIINITIONS ====================================================================
+
+class Bubbler4(DynamicalSystem):
+    """
+    Tritium bubbling system with sequential vial collection stages.
+
+    This block models a tritium collection system used in fusion reactor blanket 
+    purge gas processing. The system bubbles tritium-containing gas through a series 
+    of liquid-filled vials to capture and concentrate tritium for measurement and 
+    inventory tracking.
+
+
+    Physical Description
+    --------------------
+    The bubbler consists of two parallel processing chains:
+
+    **Soluble Chain (Vials 1-2):** 
+    Tritium already in soluble forms (HTO, HT) flows sequentially through 
+    vials 1 and 2. Each vial has a collection efficiency :math:`\\eta_{vial}`, 
+    representing the fraction of tritium that dissolves into the liquid phase 
+    and is retained.
+
+    **Insoluble Chain (Vials 3-4):** 
+    Tritium in insoluble forms (T₂, organically bound) first undergoes catalytic 
+    conversion to soluble forms with efficiency :math:`\\alpha_{conv}`. The 
+    converted tritium, along with uncaptured soluble tritium from the first chain, 
+    then flows through vials 3 and 4 with the same collection efficiency.
+
+
+    Mathematical Formulation
+    -------------------------
+    The system is governed by the following differential equations for the 
+    vial inventories :math:`x_i`:
+
+    .. math::
+
+        \\frac{dx_1}{dt} &= \\eta_{vial} \\cdot u_{sol}
+
+        \\frac{dx_2}{dt} &= \\eta_{vial} \\cdot (1-\\eta_{vial}) \\cdot u_{sol}
+
+        \\frac{dx_3}{dt} &= \\eta_{vial} \\cdot [\\alpha_{conv} \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot u_{sol}]
+
+        \\frac{dx_4}{dt} &= \\eta_{vial} \\cdot (1-\\eta_{vial}) \\cdot [\\alpha_{conv} \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot u_{sol}]
+
+
+    The sample output represents uncaptured tritium exiting the system:
+
+    .. math::
+
+        y_{sample} = (1-\\alpha_{conv}) \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot [\\alpha_{conv} \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot u_{sol}]
+
+
+    Where:
+        - :math:`u_{sol}` = soluble tritium input flow rate
+        - :math:`u_{insol}` = insoluble tritium input flow rate  
+        - :math:`\\eta_{vial}` = vial collection efficiency
+        - :math:`\\alpha_{conv}` = conversion efficiency from insoluble to soluble
+        - :math:`x_i` = tritium inventory in vial i
+
+    Parameters
+    ----------
+    conversion_efficiency : float 
+       Conversion efficiency from insoluble to soluble forms (:math:`\\alpha_{conv}`), 
+       between 0 and 1.
+    vial_efficiency : float 
+       Collection efficiency of each vial (:math:`\\eta_{vial}`), between 0 and 1.
+    replacement_times : float | list[float] | list[list[float]] 
+        Times at which each vial is replaced with a fresh one. If None, no 
+        replacement events are created. If a single value is provided, it is 
+        used for all vials. If a single list of floats is provided, it will be 
+        used for all vials. If a list of lists is provided, each sublist 
+        corresponds to the replacement times for each vial.
+
+    Notes
+    -----
+    Vial replacement is modeled as instantaneous reset events that set the 
+    corresponding vial inventory to zero, simulating the physical replacement 
+    of a full vial with an empty one.
+    """
+
+    _port_map_out = {
+        "vial1": 0,
+        "vial2": 1,
+        "vial3": 2,
+        "vial4": 3,
+        "sample_out": 4,
+    }
+    _port_map_in = {
+        "sample_in_soluble": 0,
+        "sample_in_insoluble": 1,
+    }
+
+    def __init__(
+        self,
+        conversion_efficiency=0.9,
+        vial_efficiency=0.9,
+        replacement_times=None,
+        ):
+
+        #bubbler parameters
+        self.replacement_times = replacement_times
+        self.vial_efficiency = vial_efficiency
+        self.conversion_efficiency = conversion_efficiency
+
+        #dynamical component, ode rhs
+        def _fn_d(x, u, t):
+
+            #short
+            ve = self.vial_efficiency
+            ce = self.conversion_efficiency
+
+            #unpack inputs
+            sol, ins = u
+
+            #compute vial content change rates
+            dv1 = ve * sol
+            dv2 = dv1 * (1 - ve)
+            dv3 = ve * (ce * ins + (1 - ve)**2 * sol)
+            dv4 = dv3 * (1 - ve)
+
+            return np.array([dv1, dv2, dv3, dv4])
+
+        #algebraic output component
+        def _fn_a(x, u, t):
+
+            #short
+            ve = self.vial_efficiency
+            ce = self.conversion_efficiency
+
+            #unpack inputs
+            sol, ins = u
+
+            sample_out = (1 - ce) * ins + (1 - ve)**2 * (ce * ins + (1 - ve)**2 * sol)
+
+            return np.hstack([x, sample_out])
+
+        #initialization just like `DynamicalSystem` block
+        super().__init__(func_dyn=_fn_d, func_alg=_fn_a, initial_value=np.zeros(4))
+
+        #create internal vial reset events
+        self._create_reset_events()
+
+
+    def _create_reset_event_vial(self, i, reset_times):
+        """Define event action function and return a `ScheduleList` event 
+        per vial `i` that triggers at predefined `reset_times`. 
+        """
+
+        def reset_vial_i(_):
+            #get the full engine state
+            x = self.engine.get()
+            #set index 'i' to zero
+            x[i] = 0.0
+            #set the full engine state 
+            self.engine.set(x)
+
+        return ScheduleList(
+            times_evt=reset_times, 
+            func_act=reset_vial_i
+            )
+
+
+    def _create_reset_events(self):
+        """Create reset events for all vials based on the replacement times.
+
+        Raises
+        ------
+        ValueError : If reset_times is not valid.
+
+        Returns
+        -------
+        events : list[ScheduleList]
+            list of reset events for vials
+        """
+
+        replacement_times = self.replacement_times
+        self.events = []
+
+        # if reset_times is a single list use it for all vials
+        if replacement_times is None:
+            return 
+
+        if isinstance(replacement_times, (int, float)):
+            replacement_times = [replacement_times]
+
+        # if it's a flat list use it for all vials
+        elif isinstance(replacement_times, list) and all(
+            isinstance(t, (int, float)) for t in replacement_times
+            ):
+            replacement_times = [replacement_times] * 4
+
+        elif isinstance(replacement_times, np.ndarray) and replacement_times.ndim == 1:
+            replacement_times = [replacement_times.tolist()] * 4
+
+        elif isinstance(replacement_times, list) and len(replacement_times) != 4:
+            raise ValueError(
+                "replacement_times must be a single value or a list with the same length as the number of vials"
+            )
+
+        #create the internal events
+        self.events = [
+            self._create_reset_event_vial(i, ts) for i, ts in enumerate(replacement_times)
+            ]

src/pathsim_chem/tritium/residencetime.py

@@ -0,0 +1,123 @@
+#########################################################################################
+##
+##                           Blocks for residence time modeling
+##                            (blocks/fusion/residencetime.py)
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+
+from ..dynsys import DynamicalSystem
+
+
+# BLOCKS ================================================================================
+
+class ResidenceTime(DynamicalSystem):
+    """Chemical process block with residence time model.
+
+    This block implements an internal 1st order linear ode with 
+    multiple inputs, outputs, an internal constant source term 
+    and no direct passthrough.
+
+    The internal ODE with inputs :math:`u_i` :
+
+    .. math::
+
+        \\dot{x} = - x / \\tau + \\mathrm{src} + \\sum_i \\beta_i u_i
+
+
+    And the output equation for every output `i` :
+
+    .. math::
+
+        y_i = \\gamma_i x
+
+
+    Parameters
+    ----------
+    tau : float
+        residence time, inverse natural frequency (eigenvalue)
+    betas: None | list[float] | np.ndarray[float]
+        weights of inputs that are accumulated in state, optional
+    gammas : None | list[float] | np.ndarray[float]
+        weights of states (fractions) for output, optional
+    initial_value : float
+        initial value of state / initial quantity of process
+    source_term : float
+        constant source term / generation term of the process
+    """
+
+    def __init__(self, tau=1, betas=None, gammas=None, initial_value=0, source_term=0):
+
+        #input validation
+        if np.isclose(tau, 0):
+            raise ValueError(f"'tau' must be nonzero but is {tau}")
+
+        #time constant and input/output weights
+        self.tau = tau
+        self.betas = 1 if betas is None else np.array(betas)
+        self.gammas = 1 if gammas is None else np.array(gammas)
+        self.source_term = source_term
+
+        #rhs of residence time ode
+        def _fn_d(x, u, t):
+            return -x/self.tau + self.source_term + sum(self.betas*u)
+
+        #jacobian of rhs wrt x
+        def _jc_d(x, u, t):
+            return -1/self.tau
+
+        #output function of residence time ode
+        def _fn_a(x, u, t):
+            return self.gammas * x
+
+        #initialization just like `DynamicalSystem` block
+        super().__init__(func_dyn=_fn_d, jac_dyn=_jc_d, func_alg=_fn_a, initial_value=initial_value)
+
+
+class Process(ResidenceTime):
+    """Simplified version of the `ResidenceTime` model block
+    with all inputs being summed equally and only the state 
+    and the flux being returned to the output
+
+    This block implements an internal 1st order linear ode with 
+    multiple inputs, outputs and no direct passthrough.
+
+    The internal ODE with inputs :math:`u_i` :
+
+    .. math::
+
+        \\dot{x} = - x / \\tau + \\mathrm{src} + \\sum_i u_i
+
+
+    And the output equations for output `i=0` and `i=1`:
+
+    .. math::
+
+        y_0 = x
+
+    .. math::
+
+        y_1 = x / \\tau
+
+
+    Parameters
+    ----------
+    tau : float
+        residence time, inverse natural frequency (eigenvalue)
+    initial_value : float
+        initial value of state / initial quantity of process
+    source_term : float
+        constant source term / generation term of the process
+    """
+
+    #max number of ports
+    _n_out_max = 2
+
+    #maps for input and output port labels
+    _port_map_out = {"x": 0, "x/tau": 1}
+
+    def __init__(self, tau=1, initial_value=0, source_term=0):
+        super().__init__(tau, 1, [1, 1/tau], initial_value, source_term)