From 07dd97297058c95bed67421c61a52b72eec59a2b Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 09:59:27 -0500
Subject: [PATCH 01/30] added GLC block

---
 src/pathsim_chem/tritium/glc.py | 362 ++++++++++++++++++++++++++++++++
 1 file changed, 362 insertions(+)
 create mode 100644 src/pathsim_chem/tritium/glc.py

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
new file mode 100644
index 0000000..641ea85
--- /dev/null
+++ b/src/pathsim_chem/tritium/glc.py
@@ -0,0 +1,362 @@
+# GLC block
+import pathsim
+import numpy as np
+from scipy.integrate import solve_bvp
+from scipy import constants as const
+from scipy.optimize import root_scalar
+
+
+def solve(params):
+    def solve_tritium_extraction(dimensionless_params, y_T2_in, elements):
+        """
+        Solves the BVP for tritium extraction in a bubble column.
+
+        Args:
+            params (dict): A dictionary containing the dimensionless parameters:
+                        Bo_l, phi_l, Bo_g, phi_g, psi, nu.
+            y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-).
+
+        Returns:
+            sol: The solution object from scipy.integrate.solve_bvp.
+        """
+
+        Bo_l = dimensionless_params["Bo_l"]
+        phi_l = dimensionless_params["phi_l"]
+        Bo_g = dimensionless_params["Bo_g"]
+        phi_g = dimensionless_params["phi_g"]
+        psi = dimensionless_params["psi"]
+        nu = dimensionless_params["nu"]
+
+        def ode_system(xi, S):
+            """
+            Defines the system of 4 first-order ODEs.
+            S[0] = x_T  (dimensionless liquid concentration)
+            S[1] = dx_T/d(xi)
+            S[2] = y_T2 (dimensionless gas concentration)
+            S[3] = dy_T2/d(xi)
+
+            x_T = c_T / c_T(L+)
+            xi = z / L (dimensionless position)
+            """
+            x_T, dx_T_dxi, y_T2, dy_T2_dxi = S
+
+            # Dimensionless driving force theta. Eq. 12
+            theta = x_T - np.sqrt(((1 - (psi * xi)) / nu) * y_T2)  # MATCHES PAPER
+
+            # Equation for d(S[0])/d(xi) = d(x_T)/d(xi)
+            dS0_dxi = dx_T_dxi
+
+            # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2
+            dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi)
+
+            # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi)
+            dS2_dxi = dy_T2_dxi
+
+            # Equation for d(S[3])/d(xi) = d^2(y_T2)/d(xi)^2 from eq (11)
+            # Avoid division by zero if (1 - psi * xi) is close to zero at xi=1
+
+            term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi  # Part of Eq. 9.3.3 (fourth line)
+            term2 = phi_g * theta  # Part of Eq. 9.3.3 (fourth line)
+            dS3_dxi = (Bo_g / (1 - psi * xi)) * (
+                term1 - term2
+            )  # Eq. 9.3.3 (fourth line)
+
+            return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi))
+
+        def boundary_conditions(Sa, Sb):
+            """
+            Defines the boundary conditions for the BVP.
+            Sa: solution at xi = 0 (liquid outlet)
+            Sb: solution at xi = 1 (liquid inlet)
+            """
+            # Residuals that should be zero for a valid solution.
+            # Based on equations (16) and (17) in the paper.
+
+            # At xi = 0: dx_T/d(xi) = 0
+            res1 = Sa[1]  # Eq. 10.1
+
+            # At xi = 1: x_T(1) = 1 - (1/Bo_l) * dx_T/d(xi)|_1
+            res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1])  # Eq. 10.2
+
+            # At xi = 0: y_T2(0) = y_T2(0-) + (1/Bo_g) * dy_T2/d(xi)|_0
+            res3 = Sa[2] - y_T2_in - (1 / Bo_g) * Sa[3]  # Eq. 10.3
+
+            # At xi = 1: dy_T2/d(xi) = 0
+            res4 = Sb[3]  # Eq. 10.4
+
+            return np.array([res1, res2, res3, res4])
+
+        # Set up the mesh and an initial guess for the solver.
+        xi = np.linspace(0, 1, elements + 1)
+
+        y_guess = np.zeros((4, xi.size))
+
+        # Run the BVP solver
+        sol = solve_bvp(
+            ode_system, boundary_conditions, xi, y_guess, tol=1e-6, max_nodes=10000
+        )
+
+        return sol
+
+    # Unpack parameters
+    c_T_inlet = params["c_T_inlet"]
+    y_T2_in = params["y_T2_in"]
+    P_0 = params["P_0"]
+
+    L = params["L"]
+    D = params["D"]
+
+    flow_l = params["flow_l"]
+    u_g0 = params["u_g0"]
+
+    g = params["g"]
+    T = params["T"]
+
+    elements = params["elements"]  # Number of mesh elements for solver
+
+    # --- Constants ---
+    g = const.g  # m/s^2, Gravitational acceleration
+    R = const.R  # J/(mol·K), Universal gas constant
+    N_A = const.N_A  # 1/mol, Avogadro's number
+    M_LiPb = 2.875e-25  # Kg/molecule, Lipb molecular mass
+
+    # Calculate empirical correlations
+    ρ_l = 10.45e3 * (1 - 1.61e-4 * T)  # kg/m^3, Liquid (LiPb) density
+    σ_l = 0.52 - 0.11e-3 * T  # N/m, Surface tension, liquid (LiPb) - gas (He) interface
+    μ_l = 1.87e-4 * np.exp(11640 / (R * T))  # Pa.s, Dynamic viscosity of liquid LiPb
+    ν_l = μ_l / ρ_l  # m^2/s, Kinematic viscosity of liquid LiPb
+    D_T = 2.5e-7 * np.exp(
+        -27000 / (R * T)
+    )  # m^2/s, Tritium diffusion coefficient in liquid LiPb
+
+    K_s = 2.32e-8 * np.exp(
+        -1350 / (R * T)
+    )  # atfrac*Pa^0.5, Sievert's constant for tritium in liquid LiPb
+    K_s = K_s * (ρ_l / (M_LiPb * N_A))  # mol/(m^3·Pa^0.5)
+
+    A = np.pi * (D / 2) ** 2  # m^2, Cross-sectional area of the column
+
+    # Calculate the volumetric flow rates
+    Q_l = flow_l / ρ_l  # m^3/s, Volumetric flow rate of liquid phase
+    Q_g = u_g0 * A  # m^3/s, Volumetric flow rate of gas phase
+
+    # Calculate the superficial flow velocities
+    # u_g0 =  Q_g / A  # m/s, superficial gas inlet velocity
+    u_l = Q_l / A  # m/s, superficial liquid inlet velocity
+
+    # Calculate Bond, Galilei, Schmidt and Froude numbers
+    Bn = (g * D**2 * ρ_l) / σ_l  # Bond number
+    Ga = (g * D**3) / ν_l**2  # Galilei number
+    Sc = ν_l / D_T  # Schmidt number
+    Fr = u_g0 / (g * D) ** 0.5  # Froude number
+
+    # Calculate dispersion coefficients
+    E_l = (D * u_g0) / (
+        (13 * Fr) / (1 + 6.5 * (Fr**0.8))
+    )  # m^2/s, Effective axial dispersion coefficient, liquid phase
+    E_g = (
+        0.2 * D**2
+    ) * u_g0  # m^2/s, Effective axial dispersion coefficient, gas phase
+
+    # Calculate gas hold-up (phase fraction) & mass transfer coefficient
+    C = 0.2 * (Bn ** (1 / 8)) * (Ga ** (1 / 12)) * Fr  # C = ε_g / (1 - ε_g)^4
+
+    def solveEqn(ε_g, C):
+        # Define the equation to solve
+        eqn = ε_g / (1 - ε_g) ** 4 - C
+        return eqn
+
+    ε_g_initial_guess = 0.1
+    try:
+        # bracket=[0.0001, 0.9999] tells it to *only* look in this range
+        sol = root_scalar(solveEqn, args=(C,), bracket=[0.00001, 0.99999])
+
+        # print(f"--- Using root_scalar (robust method) ---")
+        # print(f"C value was: {C}")
+        # if sol.converged:
+        #     print(f"Solved gas hold-up (εg): {sol.root:.6f}")
+        #     # Verify it
+        #     verification = sol.root / (1 - sol.root)**4
+        #     print(f"Verification (should equal C): {verification:.6f}")
+        # else:
+        #     print("Solver did not converge.")
+
+    except ValueError as e:
+        print(
+            f"Solver failed. This can happen if C is so large that no solution exists between 0 and 1."
+        )
+        print(f"Error: {e}")
+
+    ε_g = sol.root  # Gas phase fraction
+    ε_l = 1 - ε_g  # Liquid phase fraction
+
+    # Calculate outlet pressure hydrostatically & check non-negative
+    P_outlet = P_0 - (ρ_l * (1 - ε_g) * g * L)
+
+    if P_outlet <= 0:
+        raise ValueError(
+            f"Calculated gas outlet pressure P_outlet must be positive, but got {P_outlet:.2e} Pa. Check P_0, rho_l, g, and L are realistic."
+        )
+
+    # Calculate interfacial area
+    d_b = (
+        26 * (Bn**-0.5) * (Ga**-0.12) * (Fr**-0.12)
+    ) * D  # m, Mean bubble diameter AGREES WITH PAPER
+    a = 6 * ε_g / d_b  # m^-1, Specific interfacial area, assuming spherical bubbles
+
+    # Calculate volumetric mass transfer coefficient, liquid-gas
+    h_l_a = (
+        D_T * (0.6 * Sc**0.5 * Bn**0.62 * Ga**0.31 * ε_g**1.1) / (D**2)
+    )  # Volumetric mass transfer coefficient, liquid-gas
+
+    h_l = h_l_a / a  # Mass transfer coefficient
+
+    # Calculate dimensionless values
+
+    # Hydrostatic pressure ratio (Eq. 8.3) # MATCHES PAPER
+    psi = (ρ_l * g * (1 - ε_g) * L) / P_0
+    # Tritium partial pressure ratio (Eq. 8.5) # MATCHES PAPER
+    nu = ((c_T_inlet / K_s) ** 2) / P_0
+    # Bodenstein number, liquid phase (Eq. 8.9) # MATCHES PAPER
+    Bo_l = u_l * L / (ε_l * E_l)
+    # Transfer units parameter, liquid phase (Eq. 8.11) # MATCHES PAPER
+    phi_l = a * h_l * L / u_l
+    # Bodenstein number, gas phase (Eq. 8.10) # MATCHES PAPER ASSUMING u_g0
+    Bo_g = u_g0 * L / (ε_g * E_g)
+    # Transfer units parameter, gas phase (Eq. 8.12) # MATCHES PAPER
+    phi_g = 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0)
+
+    dimensionless_params = {
+        "Bo_l": Bo_l,
+        "phi_l": phi_l,
+        "Bo_g": Bo_g,
+        "phi_g": phi_g,
+        "psi": psi,
+        "nu": nu,
+    }
+
+    # Solve the model
+    solution = solve_tritium_extraction(dimensionless_params, y_T2_in, elements)
+
+    # --- Results ---
+    if solution.success:
+        # --- Dimensionless Results ---
+        x_T_outlet_dimless = solution.y[0, 0]
+        efficiency = 1 - x_T_outlet_dimless
+        y_T2_outlet_gas = solution.y[2, -1]  # y_T2 at xi=1
+
+        # --- Dimensional Results ---
+        # Liquid concentration at outlet (xi=0)
+        c_T_outlet = x_T_outlet_dimless * c_T_inlet
+
+        # Gas partial pressure at outlet (xi=1)
+        P_outlet = P_0 * (
+            1 - dimensionless_params["psi"]
+        )  # Derived from Eq. 8.4 at xi=1
+        P_T2_out = y_T2_outlet_gas * P_outlet
+
+        # Mass transfer consistency check
+        N_A = const.N_A  # Avogadro's number, 1/mol
+        # Tritium molar flow rate into the column via liquid
+        n_T_in_liquid = c_T_inlet * Q_l * N_A  # Triton/s
+
+        # Tritium molar flow rate out of the column via liquid
+        n_T_out_liquid = c_T_outlet * Q_l * N_A  # Tritons/s
+
+        # Tritium molar flow rate into the column via gas
+        P_T2_in = y_T2_in * P_0  # [Pa]
+        n_T2_in_gas = (P_T2_in * Q_g / (const.R * T)) * N_A  # T2/s
+        n_T_in_gas = n_T2_in_gas * 2  # Triton/s
+
+        # Calculate outlet gas volumetric flow rate (gas expands as pressure drops)
+        Q_g_out = (P_0 * Q_g) / P_outlet
+        # Tritium molar flow rate out of the column via gas
+        n_T2_out_gas = (P_T2_out * Q_g_out / (const.R * T)) * N_A  # T2/s
+        n_T_out_gas = n_T2_out_gas * 2  # Triton/s
+
+        T_in = n_T_in_liquid + n_T_in_gas
+        T_out = n_T_out_liquid + n_T_out_gas
+
+        results = {
+            "extraction_efficiency [%]": efficiency * 100,
+            "c_T_inlet [mol/m^3]": c_T_inlet,
+            "c_T_outlet [mol/m^3]": c_T_outlet,
+            "liquid_vol_flow [m^3/s]": Q_l,
+            "P_T2_inlet_gas [Pa]": P_T2_in,
+            "P_T2_outlet_gas [Pa]": P_T2_out,
+            "total_gas_P_outlet [Pa]": P_outlet,
+            "gas_vol_flow [m^3/s]": Q_g,
+            "tritium_out_liquid [mol/s]": n_T_out_liquid / N_A,
+            "tritium_out_gas [mol/s]": n_T_out_gas / N_A,
+        }
+        return results
+    else:
+        raise RuntimeError("BVP solver did not converge.")
+
+
+class GLC(pathsim.blocks.Function):
+    """
+    Gas Liquid Contactor model block. Inherits from Function block.
+    Inputs: c_T_inlet [mol/m^3], P_T2_inlet [Pa]
+    Outputs: n_T_out_liquid [mol/s], n_T_out_gas [mol/s]
+
+    More details about the model can be found in: https://doi.org/10.13182/FST95-A30485
+
+    Args:
+        P_0: Inlet operating pressure [Pa]
+        L: Column height [m]
+        u_g0: Superficial gas inlet velocity [m/s]
+        Q_l: Liquid volumetric flow rate [m^3/s]
+        D: Column diameter [m]
+        T: Temperature [K]
+        g: Gravitational acceleration [m/s^2], default is 9.81
+    """
+
+    _port_map_in = {
+        "c_T_inlet": 0,
+        "y_T2_in": 1,
+    }
+    _port_map_out = {
+        "c_T_outlet": 0,
+        "P_T2_out_gas": 1,
+        "efficiency": 2,
+    }
+
+    def __init__(
+        self,
+        P_0,
+        L,
+        u_g0,
+        flow_l,
+        D,
+        T,
+        g=const.g,
+        initial_nb_of_elements=20,
+    ):
+        self.params = {
+            "P_0": P_0,
+            "L": L,
+            "u_g0": u_g0,
+            "flow_l": flow_l,
+            "g": g,
+            "D": D,
+            "T": T,
+            "elements": initial_nb_of_elements,
+        }
+        super().__init__(func=self.func)
+
+    def func(self, c_T_inlet, y_T2_inlet):
+        new_params = self.params.copy()
+        new_params["c_T_inlet"] = c_T_inlet
+        new_params["y_T2_in"] = y_T2_inlet
+
+        res = solve(new_params)
+
+        c_T_outlet = res["c_T_outlet [mol/m^3]"]
+        P_T2_outlet = res["P_T2_outlet_gas [Pa]"]
+
+        n_T_out_liquid = res["tritium_out_liquid [mol/s]"]
+        n_T_out_gas = res["tritium_out_gas [mol/s]"]
+        eff = res["extraction_efficiency [%]"]
+
+        return c_T_outlet, P_T2_outlet, eff

From f29623d8b5a8b25edf2f2324604af6ccd7089a47 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 09:59:39 -0500
Subject: [PATCH 02/30] added to init

---
 src/pathsim_chem/tritium/__init__.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/pathsim_chem/tritium/__init__.py b/src/pathsim_chem/tritium/__init__.py
index 4cae89e..265bbb7 100644
--- a/src/pathsim_chem/tritium/__init__.py
+++ b/src/pathsim_chem/tritium/__init__.py
@@ -1,4 +1,5 @@
 from .residencetime import *
 from .splitter import *
 from .bubbler import *
-# from .tcap import *
\ No newline at end of file
+from .glc import *
+# from .tcap import *

From da3219a0ae40c552b23f000bfc6fed1b193863eb Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 14:30:17 -0500
Subject: [PATCH 03/30] changed model + outputs

---
 src/pathsim_chem/tritium/glc.py | 497 +++++++++++++++-----------------
 1 file changed, 237 insertions(+), 260 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 641ea85..c72c2b2 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -1,297 +1,272 @@
-# GLC block
-import pathsim
+"""
+Bubble column gas-liquid contactor model solver.
+
+This module solves the coupled, non-linear, second-order ordinary differential
+equations that describe tritium transport in a counter-current bubble column,
+based on the model by C. Malara (1995).
+"""
+
 import numpy as np
 from scipy.integrate import solve_bvp
-from scipy import constants as const
 from scipy.optimize import root_scalar
+import scipy.constants as const
+import pathsim
 
+# --- Physical Constants ---
+g = const.g  # m/s^2, Gravitational acceleration
+R = const.R  # J/(mol·K), Universal gas constant
+N_A = const.N_A  # 1/mol, Avogadro's number
+M_LiPb = 2.875e-25  # Kg/molecule, Lipb molecular mass
 
-def solve(params):
-    def solve_tritium_extraction(dimensionless_params, y_T2_in, elements):
-        """
-        Solves the BVP for tritium extraction in a bubble column.
-
-        Args:
-            params (dict): A dictionary containing the dimensionless parameters:
-                        Bo_l, phi_l, Bo_g, phi_g, psi, nu.
-            y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-).
-
-        Returns:
-            sol: The solution object from scipy.integrate.solve_bvp.
-        """
-
-        Bo_l = dimensionless_params["Bo_l"]
-        phi_l = dimensionless_params["phi_l"]
-        Bo_g = dimensionless_params["Bo_g"]
-        phi_g = dimensionless_params["phi_g"]
-        psi = dimensionless_params["psi"]
-        nu = dimensionless_params["nu"]
-
-        def ode_system(xi, S):
-            """
-            Defines the system of 4 first-order ODEs.
-            S[0] = x_T  (dimensionless liquid concentration)
-            S[1] = dx_T/d(xi)
-            S[2] = y_T2 (dimensionless gas concentration)
-            S[3] = dy_T2/d(xi)
-
-            x_T = c_T / c_T(L+)
-            xi = z / L (dimensionless position)
-            """
-            x_T, dx_T_dxi, y_T2, dy_T2_dxi = S
-
-            # Dimensionless driving force theta. Eq. 12
-            theta = x_T - np.sqrt(((1 - (psi * xi)) / nu) * y_T2)  # MATCHES PAPER
-
-            # Equation for d(S[0])/d(xi) = d(x_T)/d(xi)
-            dS0_dxi = dx_T_dxi
 
-            # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2
-            dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi)
+def _calculate_properties(params):
+    """Calculate temperature-dependent and geometry-dependent properties."""
+    T = params["T"]
+    D = params["D"]
+    Flow_l = params["Flow_l"]
+    Flow_g = params["Flow_g"]
+    P_0 = params["P_0"]
 
-            # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi)
-            dS2_dxi = dy_T2_dxi
+    # --- Fluid Properties (Temperature Dependent) ---
+    # TODO add references
+    rho_l = 10.45e3 * (1 - 1.61e-4 * T)  # kg/m^3, Liquid (LiPb) density
+    sigma_l = 0.52 - 0.11e-3 * T  # N/m, Surface tension, liquid-gas interface
+    mu_l = 1.87e-4 * np.exp(11640 / (R * T))  # Pa.s, Dynamic viscosity of liquid
+    nu_l = mu_l / rho_l  # m^2/s, Kinematic viscosity of liquid
+    D_T = 2.5e-7 * np.exp(-27000 / (R * T))  # m^2/s, Tritium diffusion coeff
+    K_s_at = 2.32e-8 * np.exp(-1350 / (R * T))  # atfrac*Pa^0.5, Sievert's const
+    K_s = K_s_at * (rho_l / (M_LiPb * N_A))  # mol/(m^3·Pa^0.5)
+
+    # --- Flow Properties ---
+    A = np.pi * (D / 2) ** 2  # m^2, Cross-sectional area
+    Q_l = Flow_l / rho_l  # m^3/s, Volumetric liquid flow rate
+    Q_g = (Flow_g * R * T) / P_0  # m^3/s, Volumetric gas flow rate at inlet
+    u_l = Q_l / A  # m/s, Superficial liquid velocity
+    u_g0 = Q_g / A  # m/s, Superficial gas velocity at inlet
+
+    # --- Dimensionless Numbers for Correlations ---
+    Bn = (g * D**2 * rho_l) / sigma_l  # Bond number
+    Ga = (g * D**3) / nu_l**2  # Galilei number
+    Sc = nu_l / D_T  # Schmidt number
+    Fr = u_g0 / (g * D) ** 0.5  # Froude number
 
-            # Equation for d(S[3])/d(xi) = d^2(y_T2)/d(xi)^2 from eq (11)
-            # Avoid division by zero if (1 - psi * xi) is close to zero at xi=1
+    # --- Hydrodynamic and Mass Transfer Parameters ---
+    # Gas hold-up (ε_g) from correlation: C = ε_g / (1 - ε_g)^4
+    C = 0.2 * (Bn ** (1 / 8)) * (Ga ** (1 / 12)) * Fr
 
-            term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi  # Part of Eq. 9.3.3 (fourth line)
-            term2 = phi_g * theta  # Part of Eq. 9.3.3 (fourth line)
-            dS3_dxi = (Bo_g / (1 - psi * xi)) * (
-                term1 - term2
-            )  # Eq. 9.3.3 (fourth line)
+    def _f_holdup(e, C_val):
+        return e / (1 - e) ** 4 - C_val
 
-            return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi))
+    try:
+        sol = root_scalar(_f_holdup, args=(C,), bracket=[1e-12, 1 - 1e-12])
+        epsilon_g = sol.root
+    except Exception as exc:
+        raise RuntimeError("Failed to solve for gas hold-up ε_g") from exc
 
-        def boundary_conditions(Sa, Sb):
-            """
-            Defines the boundary conditions for the BVP.
-            Sa: solution at xi = 0 (liquid outlet)
-            Sb: solution at xi = 1 (liquid inlet)
-            """
-            # Residuals that should be zero for a valid solution.
-            # Based on equations (16) and (17) in the paper.
+    epsilon_l = 1 - epsilon_g  # Liquid phase fraction
 
-            # At xi = 0: dx_T/d(xi) = 0
-            res1 = Sa[1]  # Eq. 10.1
+    # Dispersion coefficients
+    E_l = (D * u_g0) / ((13 * Fr) / (1 + 6.5 * (Fr**0.8)))
+    E_g = (0.2 * D**2) * u_g0
 
-            # At xi = 1: x_T(1) = 1 - (1/Bo_l) * dx_T/d(xi)|_1
-            res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1])  # Eq. 10.2
+    # Interfacial area and mass transfer coefficients
+    d_b = (26 * (Bn**-0.5) * (Ga**-0.12) * (Fr**-0.12)) * D
+    a = 6 * epsilon_g / d_b
+    h_l_a = D_T * (0.6 * Sc**0.5 * Bn**0.62 * Ga**0.31 * epsilon_g**1.1) / (D**2)
+    h_l = h_l_a / a  # Mass transfer coefficient
 
-            # At xi = 0: y_T2(0) = y_T2(0-) + (1/Bo_g) * dy_T2/d(xi)|_0
-            res3 = Sa[2] - y_T2_in - (1 / Bo_g) * Sa[3]  # Eq. 10.3
+    return {
+        "rho_l": rho_l,
+        "sigma_l": sigma_l,
+        "mu_l": mu_l,
+        "nu_l": nu_l,
+        "K_s": K_s,
+        "Q_l": Q_l,
+        "Q_g": Q_g,
+        "u_l": u_l,
+        "u_g0": u_g0,
+        "epsilon_g": epsilon_g,
+        "epsilon_l": epsilon_l,
+        "E_l": E_l,
+        "E_g": E_g,
+        "a": a,
+        "h_l": h_l,
+    }
 
-            # At xi = 1: dy_T2/d(xi) = 0
-            res4 = Sb[3]  # Eq. 10.4
 
-            return np.array([res1, res2, res3, res4])
+def _calculate_dimensionless_groups(params, phys_props):
+    """Calculate the dimensionless groups for the ODE system."""
+    # Unpack parameters
+    L, T, P_0, c_T_inlet = params["L"], params["T"], params["P_0"], params["c_T_inlet"]
+    rho_l, K_s, u_l, u_g0 = (
+        phys_props["rho_l"],
+        phys_props["K_s"],
+        phys_props["u_l"],
+        phys_props["u_g0"],
+    )
+    epsilon_g, epsilon_l, E_l, E_g = (
+        phys_props["epsilon_g"],
+        phys_props["epsilon_l"],
+        phys_props["E_l"],
+        phys_props["E_g"],
+    )
+    a, h_l = phys_props["a"], phys_props["h_l"]
+
+    # Calculate dimensionless groups
+    psi = (rho_l * g * epsilon_l * L) / P_0  # Hydrostatic pressure ratio
+    nu = ((c_T_inlet / K_s) ** 2) / P_0  # Equilibrium ratio
+    Bo_l = u_l * L / (epsilon_l * E_l)  # Bodenstein number, liquid
+    phi_l = a * h_l * L / u_l  # Transfer units, liquid (Eq. 8.11)
+    Bo_g = u_g0 * L / (epsilon_g * E_g)  # Bodenstein number, gas
+    phi_g = 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0)
 
-        # Set up the mesh and an initial guess for the solver.
-        xi = np.linspace(0, 1, elements + 1)
+    return {
+        "Bo_l": Bo_l,
+        "phi_l": phi_l,
+        "Bo_g": Bo_g,
+        "phi_g": phi_g,
+        "psi": psi,
+        "nu": nu,
+    }
 
-        y_guess = np.zeros((4, xi.size))
 
-        # Run the BVP solver
-        sol = solve_bvp(
-            ode_system, boundary_conditions, xi, y_guess, tol=1e-6, max_nodes=10000
-        )
+def _solve_bvp_system(dim_params, y_T2_in, BCs, elements):
+    """Sets up and solves the Boundary Value Problem for tritium extraction."""
+    Bo_l, phi_l, Bo_g, phi_g, psi, nu = (
+        dim_params["Bo_l"],
+        dim_params["phi_l"],
+        dim_params["Bo_g"],
+        dim_params["phi_g"],
+        dim_params["psi"],
+        dim_params["nu"],
+    )
 
-        return sol
+    def ode_system(xi, S):
+        """
+        Defines the system of 4 first-order ODEs.
+        S = [x_T, dx_T/d(xi), y_T2, dy_T2/d(xi)]
+        """
+        x_T, dx_T_dxi, y_T2, dy_T2_dxi = S
+        theta = x_T - np.sqrt(np.maximum(0, (1 - psi * xi) * y_T2 / nu))
+
+        dS0_dxi = dx_T_dxi  # d(x_T)/d(xi)
+        dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi)  # d^2(x_T)/d(xi)^2
+        dS2_dxi = dy_T2_dxi  # d(y_T2)/d(xi)
+        term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi
+        term2 = phi_g * theta
+        dS3_dxi = (Bo_g / (1 - psi * xi)) * (term1 - term2)  # d^2(y_T2)/d(xi)^2
+
+        return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi))
+
+    def boundary_conditions(Sa, Sb):
+        """Defines the boundary conditions at xi=0 (Sa) and xi=1 (Sb)."""
+        if BCs == "C-C":  # Closed-Closed
+            res1 = Sa[1]  # dx_T/d(xi) = 0 at xi=0
+            res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1])  # x_T(1) = 1 - ...
+            res3 = Sa[2] - y_T2_in - (1 / Bo_g) * Sa[3]  # y_T2(0) = y_T2_in + ...
+            res4 = Sb[3]  # dy_T2/d(xi) = 0 at xi=1
+        elif BCs == "O-C":  # Open-Closed
+            res1 = Sa[1]  # dx_T/d(xi) = 0 at xi=0
+            res2 = Sb[0] - 1.0  # x_T(1) = 1
+            res3 = Sa[2] - y_T2_in  # y_T2(0) = y_T2_in
+            res4 = Sb[3]  # dy_T2/d(xi) = 0 at xi=1
+        return np.array([res1, res2, res3, res4])
+
+    xi = np.linspace(0, 1, elements + 1)
+    y_guess = np.zeros((4, xi.size))
+    return solve_bvp(
+        ode_system, boundary_conditions, xi, y_guess, tol=1e-5, max_nodes=10000
+    )
+
+
+def _process_results(solution, params, phys_props, dim_params):
+    """Processes the BVP solution to produce dimensional results."""
+    if not solution.success:
+        print("BVP solver failed to converge.")
+        return {"error": "BVP solver failed"}
 
     # Unpack parameters
-    c_T_inlet = params["c_T_inlet"]
+    c_T_inlet, P_0, T = params["c_T_inlet"], params["P_0"], params["T"]
     y_T2_in = params["y_T2_in"]
-    P_0 = params["P_0"]
-
-    L = params["L"]
-    D = params["D"]
-
-    flow_l = params["flow_l"]
-    u_g0 = params["u_g0"]
 
-    g = params["g"]
-    T = params["T"]
+    # Dimensionless results
+    x_T_outlet_dimless = solution.y[0, 0]
+    Q_l, Q_g = phys_props["Q_l"], phys_props["Q_g"]
+    y_T2_outlet_gas = solution.y[2, -1]
+    efficiency = 1 - x_T_outlet_dimless
+
+    # Dimensional results
+    c_T_outlet = x_T_outlet_dimless * c_T_inlet
+    P_outlet = P_0 * (1 - dim_params["psi"])
+    P_T2_out = y_T2_outlet_gas * P_outlet
+    P_T2_in = y_T2_in * P_0
+
+    # Mass balance check
+    n_T_in_liquid = c_T_inlet * Q_l * N_A
+    n_T_out_liquid = c_T_outlet * Q_l * N_A
+    n_T2_in_gas = (P_T2_in * Q_g / (R * T)) * N_A
+    n_T_in_gas = n_T2_in_gas * 2
+    Q_g_out = (P_0 * Q_g) / P_outlet
+    n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A
+    n_T_out_gas = n_T2_out_gas * 2
+
+    results = {
+        "Total tritium in [T/s]": n_T_in_liquid + n_T_in_gas,
+        "Total tritium out [T/s]": n_T_out_liquid + n_T_out_gas,
+        "extraction_efficiency [fraction]": efficiency,
+        "c_T_inlet [mol/m^3]": c_T_inlet,
+        "c_T_outlet [mol/m^3]": c_T_outlet,
+        "liquid_vol_flow [m^3/s]": Q_l,
+        "P_T2_inlet_gas [Pa]": P_T2_in,
+        "P_T2_outlet_gas [Pa]": P_T2_out,
+        "y_T2_outlet_gas": y_T2_outlet_gas,
+        "gas_vol_flow [m^3/s]": Q_g,
+        "total_gas_P_outlet [Pa]": P_outlet,
+    }
 
-    elements = params["elements"]  # Number of mesh elements for solver
-
-    # --- Constants ---
-    g = const.g  # m/s^2, Gravitational acceleration
-    R = const.R  # J/(mol·K), Universal gas constant
-    N_A = const.N_A  # 1/mol, Avogadro's number
-    M_LiPb = 2.875e-25  # Kg/molecule, Lipb molecular mass
-
-    # Calculate empirical correlations
-    ρ_l = 10.45e3 * (1 - 1.61e-4 * T)  # kg/m^3, Liquid (LiPb) density
-    σ_l = 0.52 - 0.11e-3 * T  # N/m, Surface tension, liquid (LiPb) - gas (He) interface
-    μ_l = 1.87e-4 * np.exp(11640 / (R * T))  # Pa.s, Dynamic viscosity of liquid LiPb
-    ν_l = μ_l / ρ_l  # m^2/s, Kinematic viscosity of liquid LiPb
-    D_T = 2.5e-7 * np.exp(
-        -27000 / (R * T)
-    )  # m^2/s, Tritium diffusion coefficient in liquid LiPb
-
-    K_s = 2.32e-8 * np.exp(
-        -1350 / (R * T)
-    )  # atfrac*Pa^0.5, Sievert's constant for tritium in liquid LiPb
-    K_s = K_s * (ρ_l / (M_LiPb * N_A))  # mol/(m^3·Pa^0.5)
-
-    A = np.pi * (D / 2) ** 2  # m^2, Cross-sectional area of the column
-
-    # Calculate the volumetric flow rates
-    Q_l = flow_l / ρ_l  # m^3/s, Volumetric flow rate of liquid phase
-    Q_g = u_g0 * A  # m^3/s, Volumetric flow rate of gas phase
-
-    # Calculate the superficial flow velocities
-    # u_g0 =  Q_g / A  # m/s, superficial gas inlet velocity
-    u_l = Q_l / A  # m/s, superficial liquid inlet velocity
-
-    # Calculate Bond, Galilei, Schmidt and Froude numbers
-    Bn = (g * D**2 * ρ_l) / σ_l  # Bond number
-    Ga = (g * D**3) / ν_l**2  # Galilei number
-    Sc = ν_l / D_T  # Schmidt number
-    Fr = u_g0 / (g * D) ** 0.5  # Froude number
+    # Add all calculated parameters to the results dictionary
+    results.update(phys_props)
+    results.update(dim_params)
 
-    # Calculate dispersion coefficients
-    E_l = (D * u_g0) / (
-        (13 * Fr) / (1 + 6.5 * (Fr**0.8))
-    )  # m^2/s, Effective axial dispersion coefficient, liquid phase
-    E_g = (
-        0.2 * D**2
-    ) * u_g0  # m^2/s, Effective axial dispersion coefficient, gas phase
+    return [results, solution]
 
-    # Calculate gas hold-up (phase fraction) & mass transfer coefficient
-    C = 0.2 * (Bn ** (1 / 8)) * (Ga ** (1 / 12)) * Fr  # C = ε_g / (1 - ε_g)^4
 
-    def solveEqn(ε_g, C):
-        # Define the equation to solve
-        eqn = ε_g / (1 - ε_g) ** 4 - C
-        return eqn
+def solve(params):
+    """
+    Main solver function for the bubble column model.
 
-    ε_g_initial_guess = 0.1
-    try:
-        # bracket=[0.0001, 0.9999] tells it to *only* look in this range
-        sol = root_scalar(solveEqn, args=(C,), bracket=[0.00001, 0.99999])
-
-        # print(f"--- Using root_scalar (robust method) ---")
-        # print(f"C value was: {C}")
-        # if sol.converged:
-        #     print(f"Solved gas hold-up (εg): {sol.root:.6f}")
-        #     # Verify it
-        #     verification = sol.root / (1 - sol.root)**4
-        #     print(f"Verification (should equal C): {verification:.6f}")
-        # else:
-        #     print("Solver did not converge.")
-
-    except ValueError as e:
-        print(
-            f"Solver failed. This can happen if C is so large that no solution exists between 0 and 1."
-        )
-        print(f"Error: {e}")
+    Args:
+        params (dict): A dictionary of all input parameters for the model.
 
-    ε_g = sol.root  # Gas phase fraction
-    ε_l = 1 - ε_g  # Liquid phase fraction
+    Returns:
+        dict: A dictionary containing the simulation results.
+    """
+    # Adjust inlet gas concentration to avoid numerical instability at zero
+    y_T2_in = max(params["y_T2_in"], 1e-20)
 
-    # Calculate outlet pressure hydrostatically & check non-negative
-    P_outlet = P_0 - (ρ_l * (1 - ε_g) * g * L)
+    # 1. Calculate physical, hydrodynamic, and mass transfer properties
+    phys_props = _calculate_properties(params)
 
+    # Pre-solver check for non-physical outlet pressure
+    P_outlet = params["P_0"] - (
+        phys_props["rho_l"] * (1 - phys_props["epsilon_g"]) * g * params["L"]
+    )
     if P_outlet <= 0:
         raise ValueError(
-            f"Calculated gas outlet pressure P_outlet must be positive, but got {P_outlet:.2e} Pa. Check P_0, rho_l, g, and L are realistic."
+            f"Calculated gas outlet pressure is non-positive ({P_outlet:.2e} Pa). "
+            "Check input parameters P_0, L, etc."
         )
 
-    # Calculate interfacial area
-    d_b = (
-        26 * (Bn**-0.5) * (Ga**-0.12) * (Fr**-0.12)
-    ) * D  # m, Mean bubble diameter AGREES WITH PAPER
-    a = 6 * ε_g / d_b  # m^-1, Specific interfacial area, assuming spherical bubbles
+    # 2. Calculate dimensionless groups for the ODE system
+    dim_params = _calculate_dimensionless_groups(params, phys_props)
 
-    # Calculate volumetric mass transfer coefficient, liquid-gas
-    h_l_a = (
-        D_T * (0.6 * Sc**0.5 * Bn**0.62 * Ga**0.31 * ε_g**1.1) / (D**2)
-    )  # Volumetric mass transfer coefficient, liquid-gas
-
-    h_l = h_l_a / a  # Mass transfer coefficient
+    # 3. Solve the boundary value problem
+    solution = _solve_bvp_system(dim_params, y_T2_in, params["BCs"], params["elements"])
 
-    # Calculate dimensionless values
-
-    # Hydrostatic pressure ratio (Eq. 8.3) # MATCHES PAPER
-    psi = (ρ_l * g * (1 - ε_g) * L) / P_0
-    # Tritium partial pressure ratio (Eq. 8.5) # MATCHES PAPER
-    nu = ((c_T_inlet / K_s) ** 2) / P_0
-    # Bodenstein number, liquid phase (Eq. 8.9) # MATCHES PAPER
-    Bo_l = u_l * L / (ε_l * E_l)
-    # Transfer units parameter, liquid phase (Eq. 8.11) # MATCHES PAPER
-    phi_l = a * h_l * L / u_l
-    # Bodenstein number, gas phase (Eq. 8.10) # MATCHES PAPER ASSUMING u_g0
-    Bo_g = u_g0 * L / (ε_g * E_g)
-    # Transfer units parameter, gas phase (Eq. 8.12) # MATCHES PAPER
-    phi_g = 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0)
-
-    dimensionless_params = {
-        "Bo_l": Bo_l,
-        "phi_l": phi_l,
-        "Bo_g": Bo_g,
-        "phi_g": phi_g,
-        "psi": psi,
-        "nu": nu,
-    }
+    # 4. Process and return the results in a dimensional format
+    [results, solution] = _process_results(solution, params, phys_props, dim_params)
 
-    # Solve the model
-    solution = solve_tritium_extraction(dimensionless_params, y_T2_in, elements)
-
-    # --- Results ---
-    if solution.success:
-        # --- Dimensionless Results ---
-        x_T_outlet_dimless = solution.y[0, 0]
-        efficiency = 1 - x_T_outlet_dimless
-        y_T2_outlet_gas = solution.y[2, -1]  # y_T2 at xi=1
-
-        # --- Dimensional Results ---
-        # Liquid concentration at outlet (xi=0)
-        c_T_outlet = x_T_outlet_dimless * c_T_inlet
-
-        # Gas partial pressure at outlet (xi=1)
-        P_outlet = P_0 * (
-            1 - dimensionless_params["psi"]
-        )  # Derived from Eq. 8.4 at xi=1
-        P_T2_out = y_T2_outlet_gas * P_outlet
-
-        # Mass transfer consistency check
-        N_A = const.N_A  # Avogadro's number, 1/mol
-        # Tritium molar flow rate into the column via liquid
-        n_T_in_liquid = c_T_inlet * Q_l * N_A  # Triton/s
-
-        # Tritium molar flow rate out of the column via liquid
-        n_T_out_liquid = c_T_outlet * Q_l * N_A  # Tritons/s
-
-        # Tritium molar flow rate into the column via gas
-        P_T2_in = y_T2_in * P_0  # [Pa]
-        n_T2_in_gas = (P_T2_in * Q_g / (const.R * T)) * N_A  # T2/s
-        n_T_in_gas = n_T2_in_gas * 2  # Triton/s
-
-        # Calculate outlet gas volumetric flow rate (gas expands as pressure drops)
-        Q_g_out = (P_0 * Q_g) / P_outlet
-        # Tritium molar flow rate out of the column via gas
-        n_T2_out_gas = (P_T2_out * Q_g_out / (const.R * T)) * N_A  # T2/s
-        n_T_out_gas = n_T2_out_gas * 2  # Triton/s
-
-        T_in = n_T_in_liquid + n_T_in_gas
-        T_out = n_T_out_liquid + n_T_out_gas
-
-        results = {
-            "extraction_efficiency [%]": efficiency * 100,
-            "c_T_inlet [mol/m^3]": c_T_inlet,
-            "c_T_outlet [mol/m^3]": c_T_outlet,
-            "liquid_vol_flow [m^3/s]": Q_l,
-            "P_T2_inlet_gas [Pa]": P_T2_in,
-            "P_T2_outlet_gas [Pa]": P_T2_out,
-            "total_gas_P_outlet [Pa]": P_outlet,
-            "gas_vol_flow [m^3/s]": Q_g,
-            "tritium_out_liquid [mol/s]": n_T_out_liquid / N_A,
-            "tritium_out_gas [mol/s]": n_T_out_gas / N_A,
-        }
-        return results
-    else:
-        raise RuntimeError("BVP solver did not converge.")
+    return [results, solution]
 
 
 class GLC(pathsim.blocks.Function):
@@ -326,7 +301,7 @@ def __init__(
         self,
         P_0,
         L,
-        u_g0,
+        flow_g,
         flow_l,
         D,
         T,
@@ -336,12 +311,13 @@ def __init__(
         self.params = {
             "P_0": P_0,
             "L": L,
-            "u_g0": u_g0,
-            "flow_l": flow_l,
+            "Flow_l": flow_l,
+            "Flow_g": flow_g,
             "g": g,
             "D": D,
             "T": T,
             "elements": initial_nb_of_elements,
+            "BCs": "C-C",  # hard coded for now
         }
         super().__init__(func=self.func)
 
@@ -353,10 +329,11 @@ def func(self, c_T_inlet, y_T2_inlet):
         res = solve(new_params)
 
         c_T_outlet = res["c_T_outlet [mol/m^3]"]
-        P_T2_outlet = res["P_T2_outlet_gas [Pa]"]
+        y_T2_outlet = res["y_T2_outlet_gas"]
+        P_total_outlet = res["total_gas_P_outlet [Pa]"]
 
         n_T_out_liquid = res["tritium_out_liquid [mol/s]"]
         n_T_out_gas = res["tritium_out_gas [mol/s]"]
         eff = res["extraction_efficiency [%]"]
 
-        return c_T_outlet, P_T2_outlet, eff
+        return c_T_outlet, y_T2_outlet, P_total_outlet, eff

From 3498f5461c71faef4b0754e2964ec1b2b6f4d924 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 14:32:34 -0500
Subject: [PATCH 04/30] change out ports aliases

---
 src/pathsim_chem/tritium/glc.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index c72c2b2..7fe734a 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -293,8 +293,9 @@ class GLC(pathsim.blocks.Function):
     }
     _port_map_out = {
         "c_T_outlet": 0,
-        "P_T2_out_gas": 1,
-        "efficiency": 2,
+        "y_T2_out": 1,
+        "P_out_gas": 2,
+        "efficiency": 3,
     }
 
     def __init__(

From 80db752e8b3642e97d49e65d495c3d28eaebafad Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 14:45:58 -0500
Subject: [PATCH 05/30] fixed type returned by solve

---
 src/pathsim_chem/tritium/glc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 7fe734a..f0d1d76 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -327,7 +327,7 @@ def func(self, c_T_inlet, y_T2_inlet):
         new_params["c_T_inlet"] = c_T_inlet
         new_params["y_T2_in"] = y_T2_inlet
 
-        res = solve(new_params)
+        [res, _] = solve(new_params)
 
         c_T_outlet = res["c_T_outlet [mol/m^3]"]
         y_T2_outlet = res["y_T2_outlet_gas"]

From 3ca7ed4c67952412e966904b5a98a3a42f37e816 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 14:49:02 -0500
Subject: [PATCH 06/30] removed unused variables

---
 src/pathsim_chem/tritium/glc.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index f0d1d76..cac9287 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -333,8 +333,6 @@ def func(self, c_T_inlet, y_T2_inlet):
         y_T2_outlet = res["y_T2_outlet_gas"]
         P_total_outlet = res["total_gas_P_outlet [Pa]"]
 
-        n_T_out_liquid = res["tritium_out_liquid [mol/s]"]
-        n_T_out_gas = res["tritium_out_gas [mol/s]"]
         eff = res["extraction_efficiency [%]"]
 
         return c_T_outlet, y_T2_outlet, P_total_outlet, eff

From 3bed188ee0db28be72a6b7452dd316cd852f447f Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 15:03:38 -0500
Subject: [PATCH 07/30] fraction instead of %

---
 src/pathsim_chem/tritium/glc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index cac9287..96ea097 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -333,6 +333,6 @@ def func(self, c_T_inlet, y_T2_inlet):
         y_T2_outlet = res["y_T2_outlet_gas"]
         P_total_outlet = res["total_gas_P_outlet [Pa]"]
 
-        eff = res["extraction_efficiency [%]"]
+        eff = res["extraction_efficiency [fraction]"]
 
         return c_T_outlet, y_T2_outlet, P_total_outlet, eff

From c957a45b4fcf4e9d2a0232ecd215c372430d4d5b Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 15:04:05 -0500
Subject: [PATCH 08/30] raise error instead

---
 src/pathsim_chem/tritium/glc.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 96ea097..5c4adc1 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -182,8 +182,7 @@ def boundary_conditions(Sa, Sb):
 def _process_results(solution, params, phys_props, dim_params):
     """Processes the BVP solution to produce dimensional results."""
     if not solution.success:
-        print("BVP solver failed to converge.")
-        return {"error": "BVP solver failed"}
+        raise RuntimeError("BVP solver failed to converge.")
 
     # Unpack parameters
     c_T_inlet, P_0, T = params["c_T_inlet"], params["P_0"], params["T"]

From 422e07e765b1568b53eb9739fcc372360258e873 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Thu, 13 Nov 2025 20:56:35 +0000
Subject: [PATCH 09/30] collected physical properties unpacking

---
 src/pathsim_chem/tritium/glc.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 5c4adc1..dd7d09c 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -98,19 +98,20 @@ def _calculate_dimensionless_groups(params, phys_props):
     """Calculate the dimensionless groups for the ODE system."""
     # Unpack parameters
     L, T, P_0, c_T_inlet = params["L"], params["T"], params["P_0"], params["c_T_inlet"]
-    rho_l, K_s, u_l, u_g0 = (
+    
+    # Unpack physical properties
+    rho_l, K_s, u_l, u_g0, epsilon_g, epsilon_l, E_l, E_g, a, h_l = (
         phys_props["rho_l"],
         phys_props["K_s"],
         phys_props["u_l"],
         phys_props["u_g0"],
-    )
-    epsilon_g, epsilon_l, E_l, E_g = (
         phys_props["epsilon_g"],
         phys_props["epsilon_l"],
         phys_props["E_l"],
         phys_props["E_g"],
+        phys_props["a"], 
+        phys_props["h_l"]
     )
-    a, h_l = phys_props["a"], phys_props["h_l"]
 
     # Calculate dimensionless groups
     psi = (rho_l * g * epsilon_l * L) / P_0  # Hydrostatic pressure ratio

From 0d1b27fcbfbe5cfd70d9a2988b0b0c7a08012a8a Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 16:01:21 -0500
Subject: [PATCH 10/30] added n_T_out_gas as output

---
 src/pathsim_chem/tritium/glc.py | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 5c4adc1..adc945d 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -205,13 +205,15 @@ def _process_results(solution, params, phys_props, dim_params):
     n_T_out_liquid = c_T_outlet * Q_l * N_A
     n_T2_in_gas = (P_T2_in * Q_g / (R * T)) * N_A
     n_T_in_gas = n_T2_in_gas * 2
-    Q_g_out = (P_0 * Q_g) / P_outlet
-    n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A
-    n_T_out_gas = n_T2_out_gas * 2
+    # Q_g_out = (P_0 * Q_g) / P_outlet
+    # n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A
+    n_T_out_gas = efficiency * n_T_in_liquid
 
     results = {
         "Total tritium in [T/s]": n_T_in_liquid + n_T_in_gas,
         "Total tritium out [T/s]": n_T_out_liquid + n_T_out_gas,
+        "tritium_out_liquid [T/s]": n_T_out_liquid,
+        "tritium_out_gas [T/s]": n_T_out_gas,
         "extraction_efficiency [fraction]": efficiency,
         "c_T_inlet [mol/m^3]": c_T_inlet,
         "c_T_outlet [mol/m^3]": c_T_outlet,
@@ -227,7 +229,7 @@ def _process_results(solution, params, phys_props, dim_params):
     results.update(phys_props)
     results.update(dim_params)
 
-    return [results, solution]
+    return results, solution
 
 
 def solve(params):
@@ -295,6 +297,7 @@ class GLC(pathsim.blocks.Function):
         "y_T2_out": 1,
         "P_out_gas": 2,
         "efficiency": 3,
+        "n_T_out_gas": 4,
     }
 
     def __init__(
@@ -326,12 +329,13 @@ def func(self, c_T_inlet, y_T2_inlet):
         new_params["c_T_inlet"] = c_T_inlet
         new_params["y_T2_in"] = y_T2_inlet
 
-        [res, _] = solve(new_params)
+        res, _ = solve(new_params)
 
         c_T_outlet = res["c_T_outlet [mol/m^3]"]
         y_T2_outlet = res["y_T2_outlet_gas"]
         P_total_outlet = res["total_gas_P_outlet [Pa]"]
+        n_T_out_gas = res["tritium_out_gas [T/s]"]
 
         eff = res["extraction_efficiency [fraction]"]
 
-        return c_T_outlet, y_T2_outlet, P_total_outlet, eff
+        return c_T_outlet, y_T2_outlet, P_total_outlet, eff, n_T_out_gas

From 6ef07dbf26ed08037f69b0d5644b94fa691e6837 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 16:04:56 -0500
Subject: [PATCH 11/30] everything in mol

---
 src/pathsim_chem/tritium/glc.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index adc945d..462ce0d 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -201,19 +201,19 @@ def _process_results(solution, params, phys_props, dim_params):
     P_T2_in = y_T2_in * P_0
 
     # Mass balance check
-    n_T_in_liquid = c_T_inlet * Q_l * N_A
-    n_T_out_liquid = c_T_outlet * Q_l * N_A
-    n_T2_in_gas = (P_T2_in * Q_g / (R * T)) * N_A
+    n_T_in_liquid = c_T_inlet * Q_l
+    n_T_out_liquid = c_T_outlet * Q_l
+    n_T2_in_gas = P_T2_in * Q_g / (R * T)
     n_T_in_gas = n_T2_in_gas * 2
     # Q_g_out = (P_0 * Q_g) / P_outlet
-    # n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A
+    # n_T2_out_gas = (P_T2_out * Q_g_out / (R * T))
     n_T_out_gas = efficiency * n_T_in_liquid
 
     results = {
-        "Total tritium in [T/s]": n_T_in_liquid + n_T_in_gas,
-        "Total tritium out [T/s]": n_T_out_liquid + n_T_out_gas,
-        "tritium_out_liquid [T/s]": n_T_out_liquid,
-        "tritium_out_gas [T/s]": n_T_out_gas,
+        "Total tritium in [mol/s]": n_T_in_liquid + n_T_in_gas,
+        "Total tritium out [mol/s]": n_T_out_liquid + n_T_out_gas,
+        "tritium_out_liquid [mol/s]": n_T_out_liquid,
+        "tritium_out_gas [mol/s]": n_T_out_gas,
         "extraction_efficiency [fraction]": efficiency,
         "c_T_inlet [mol/m^3]": c_T_inlet,
         "c_T_outlet [mol/m^3]": c_T_outlet,

From a216a0385252a6ea5dd1153ea5166f4f33f83806 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 16:05:52 -0500
Subject: [PATCH 12/30] fixed T/s to mol/s

---
 src/pathsim_chem/tritium/glc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 462ce0d..482a4fc 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -334,7 +334,7 @@ def func(self, c_T_inlet, y_T2_inlet):
         c_T_outlet = res["c_T_outlet [mol/m^3]"]
         y_T2_outlet = res["y_T2_outlet_gas"]
         P_total_outlet = res["total_gas_P_outlet [Pa]"]
-        n_T_out_gas = res["tritium_out_gas [T/s]"]
+        n_T_out_gas = res["tritium_out_gas [mol/s]"]
 
         eff = res["extraction_efficiency [fraction]"]
 

From 53572095c53c889e164f696cb76c850328825f65 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Thu, 13 Nov 2025 21:07:43 +0000
Subject: [PATCH 13/30] Improved function descriptions

---
 src/pathsim_chem/tritium/glc.py | 83 ++++++++++++++++++++++++++++-----
 1 file changed, 71 insertions(+), 12 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index dd7d09c..da5f630 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -20,7 +20,20 @@
 
 
 def _calculate_properties(params):
-    """Calculate temperature-dependent and geometry-dependent properties."""
+    """
+    Calculate temperature-dependent and geometry-dependent physical properties.
+
+    This function computes fluid properties, flow characteristics,
+    dimensionless numbers, and hydrodynamic parameters based on the input
+    geometry and operating conditions.
+
+    Args:
+        params (dict): Dictionary of input parameters including T, D, Flow_l,
+                       Flow_g, and P_0.
+
+    Returns:
+        dict: A dictionary containing the calculated physical properties.
+    """
     T = params["T"]
     D = params["D"]
     Flow_l = params["Flow_l"]
@@ -95,7 +108,19 @@ def _f_holdup(e, C_val):
 
 
 def _calculate_dimensionless_groups(params, phys_props):
-    """Calculate the dimensionless groups for the ODE system."""
+    """
+    Calculate the dimensionless groups for the ODE system.
+
+    Args:
+        params (dict): Dictionary of input parameters including L, T, P_0,
+                       and c_T_inlet.
+        phys_props (dict): Dictionary of physical properties calculated by
+                           _calculate_properties.
+
+    Returns:
+        dict: A dictionary containing the dimensionless groups (Bo_l, phi_l,
+              Bo_g, phi_g, psi, nu).
+    """
     # Unpack parameters
     L, T, P_0, c_T_inlet = params["L"], params["T"], params["P_0"], params["c_T_inlet"]
     
@@ -132,7 +157,22 @@ def _calculate_dimensionless_groups(params, phys_props):
 
 
 def _solve_bvp_system(dim_params, y_T2_in, BCs, elements):
-    """Sets up and solves the Boundary Value Problem for tritium extraction."""
+    """
+    Set up and solve the Boundary Value Problem for tritium extraction.
+
+    This function defines the system of ordinary differential equations (ODEs)
+    and the corresponding boundary conditions, then solves the BVP using
+    `scipy.integrate.solve_bvp`.
+
+    Args:
+        dim_params (dict): Dictionary of dimensionless groups.
+        y_T2_in (float): Inlet mole fraction of T2 in the gas phase.
+        BCs (str): String specifying the boundary conditions (e.g., "C-C").
+        elements (int): Number of elements for the initial mesh.
+
+    Returns:
+        scipy.integrate.OdeSolution: The solution object from `solve_bvp`.
+    """
     Bo_l, phi_l, Bo_g, phi_g, psi, nu = (
         dim_params["Bo_l"],
         dim_params["phi_l"],
@@ -181,7 +221,22 @@ def boundary_conditions(Sa, Sb):
 
 
 def _process_results(solution, params, phys_props, dim_params):
-    """Processes the BVP solution to produce dimensional results."""
+    """
+    Process the BVP solution to produce dimensional results.
+
+    This function converts the dimensionless solution of the BVP back into
+    dimensional quantities, calculates the extraction efficiency, performs a
+    mass balance check, and aggregates all results.
+
+    Args:
+        solution (scipy.integrate.OdeSolution): The BVP solution object.
+        params (dict): The original input parameters.
+        phys_props (dict): The calculated physical properties.
+        dim_params (dict): The calculated dimensionless groups.
+
+    Returns:
+        list: A list containing the results dictionary and the solution object.
+    """
     if not solution.success:
         raise RuntimeError("BVP solver failed to converge.")
 
@@ -216,12 +271,14 @@ def _process_results(solution, params, phys_props, dim_params):
         "extraction_efficiency [fraction]": efficiency,
         "c_T_inlet [mol/m^3]": c_T_inlet,
         "c_T_outlet [mol/m^3]": c_T_outlet,
-        "liquid_vol_flow [m^3/s]": Q_l,
         "P_T2_inlet_gas [Pa]": P_T2_in,
         "P_T2_outlet_gas [Pa]": P_T2_out,
         "y_T2_outlet_gas": y_T2_outlet_gas,
-        "gas_vol_flow [m^3/s]": Q_g,
+        "total_gas_P_inlet [Pa]": P_0,
         "total_gas_P_outlet [Pa]": P_outlet,
+        "liquid_vol_flow [m^3/s]": Q_l,
+        "gas_vol_flow [m^3/s]": Q_g,
+
     }
 
     # Add all calculated parameters to the results dictionary
@@ -236,10 +293,14 @@ def solve(params):
     Main solver function for the bubble column model.
 
     Args:
-        params (dict): A dictionary of all input parameters for the model.
+        params (dict): A dictionary of all input parameters for the model,
+                       including operational conditions and geometry.
 
     Returns:
-        dict: A dictionary containing the simulation results.
+        list: A list containing:
+              - dict: A dictionary containing the simulation results.
+              - scipy.integrate.OdeSolution: The raw solution object from the
+                BVP solver.
     """
     # Adjust inlet gas concentration to avoid numerical instability at zero
     y_T2_in = max(params["y_T2_in"], 1e-20)
@@ -272,16 +333,14 @@ def solve(params):
 class GLC(pathsim.blocks.Function):
     """
     Gas Liquid Contactor model block. Inherits from Function block.
-    Inputs: c_T_inlet [mol/m^3], P_T2_inlet [Pa]
-    Outputs: n_T_out_liquid [mol/s], n_T_out_gas [mol/s]
 
     More details about the model can be found in: https://doi.org/10.13182/FST95-A30485
 
     Args:
         P_0: Inlet operating pressure [Pa]
         L: Column height [m]
-        u_g0: Superficial gas inlet velocity [m/s]
-        Q_l: Liquid volumetric flow rate [m^3/s]
+        flow_g: Gas mass flow rate [kg/s]
+        flow_l: Liquid mass flow rate [kg/s]
         D: Column diameter [m]
         T: Temperature [K]
         g: Gravitational acceleration [m/s^2], default is 9.81

From 859cd76febd31fb6e0414da57371606c441087da Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 16:13:01 -0500
Subject: [PATCH 14/30] update tritium output calculations and adjust for mass
 balance error

---
 src/pathsim_chem/tritium/glc.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 482a4fc..7e4ecac 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -205,9 +205,15 @@ def _process_results(solution, params, phys_props, dim_params):
     n_T_out_liquid = c_T_outlet * Q_l
     n_T2_in_gas = P_T2_in * Q_g / (R * T)
     n_T_in_gas = n_T2_in_gas * 2
-    # Q_g_out = (P_0 * Q_g) / P_outlet
-    # n_T2_out_gas = (P_T2_out * Q_g_out / (R * T))
-    n_T_out_gas = efficiency * n_T_in_liquid
+    # n_T_out_gas = efficiency * n_T_in_liquid
+    Q_g_out = (P_0 * Q_g) / P_outlet
+    n_T2_out_gas = P_T2_out * Q_g_out / (R * T)
+    n_T_out_gas = n_T2_out_gas * 2
+
+    # Adjust for any mass balance error
+    mass_balance_error = (n_T_in_liquid + n_T_in_gas) - (n_T_out_liquid + n_T_out_gas)
+    n_T_out_gas += mass_balance_error * efficiency
+    n_T_out_liquid += mass_balance_error * (1 - efficiency)
 
     results = {
         "Total tritium in [mol/s]": n_T_in_liquid + n_T_in_gas,

From 8a96f1f6177a5a3663888c568410069e0ca25928 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 13 Nov 2025 16:34:46 -0500
Subject: [PATCH 15/30] comments

---
 src/pathsim_chem/tritium/glc.py | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index cba9e09..c47cd7a 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -123,7 +123,7 @@ def _calculate_dimensionless_groups(params, phys_props):
     """
     # Unpack parameters
     L, T, P_0, c_T_inlet = params["L"], params["T"], params["P_0"], params["c_T_inlet"]
-    
+
     # Unpack physical properties
     rho_l, K_s, u_l, u_g0, epsilon_g, epsilon_l, E_l, E_g, a, h_l = (
         phys_props["rho_l"],
@@ -134,8 +134,8 @@ def _calculate_dimensionless_groups(params, phys_props):
         phys_props["epsilon_l"],
         phys_props["E_l"],
         phys_props["E_g"],
-        phys_props["a"], 
-        phys_props["h_l"]
+        phys_props["a"],
+        phys_props["h_l"],
     )
 
     # Calculate dimensionless groups
@@ -257,14 +257,13 @@ def _process_results(solution, params, phys_props, dim_params):
     P_T2_in = y_T2_in * P_0
 
     # Mass balance check
-    n_T_in_liquid = c_T_inlet * Q_l
-    n_T_out_liquid = c_T_outlet * Q_l
-    n_T2_in_gas = P_T2_in * Q_g / (R * T)
-    n_T_in_gas = n_T2_in_gas * 2
-    # n_T_out_gas = efficiency * n_T_in_liquid
-    Q_g_out = (P_0 * Q_g) / P_outlet
-    n_T2_out_gas = P_T2_out * Q_g_out / (R * T)
-    n_T_out_gas = n_T2_out_gas * 2
+    n_T_in_liquid = c_T_inlet * Q_l  # mol/s
+    n_T_out_liquid = c_T_outlet * Q_l  # mol/s
+    n_T2_in_gas = P_T2_in * Q_g / (R * T)  # mol/s
+    n_T_in_gas = n_T2_in_gas * 2  # mol/s
+    Q_g_out = (P_0 * Q_g) / P_outlet  # m3/s
+    n_T2_out_gas = P_T2_out * Q_g_out / (R * T)  # mol/s
+    n_T_out_gas = n_T2_out_gas * 2  # mol/s
 
     # Adjust for any mass balance error
     mass_balance_error = (n_T_in_liquid + n_T_in_gas) - (n_T_out_liquid + n_T_out_gas)
@@ -286,7 +285,6 @@ def _process_results(solution, params, phys_props, dim_params):
         "total_gas_P_outlet [Pa]": P_outlet,
         "liquid_vol_flow [m^3/s]": Q_l,
         "gas_vol_flow [m^3/s]": Q_g,
-
     }
 
     # Add all calculated parameters to the results dictionary

From 0cff1d4466e13575bb1e4d92f19724e1ca326f1c Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Mon, 17 Nov 2025 20:59:13 +0000
Subject: [PATCH 16/30] added Q_l & Q_g_out to GLC block

---
 src/pathsim_chem/tritium/glc.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index c47cd7a..c8d4486 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -284,7 +284,7 @@ def _process_results(solution, params, phys_props, dim_params):
         "total_gas_P_inlet [Pa]": P_0,
         "total_gas_P_outlet [Pa]": P_outlet,
         "liquid_vol_flow [m^3/s]": Q_l,
-        "gas_vol_flow [m^3/s]": Q_g,
+        "gas_vol_flow_outlet [m^3/s]": Q_g,
     }
 
     # Add all calculated parameters to the results dictionary
@@ -359,9 +359,10 @@ class GLC(pathsim.blocks.Function):
     _port_map_out = {
         "c_T_outlet": 0,
         "y_T2_out": 1,
-        "P_out_gas": 2,
-        "efficiency": 3,
-        "n_T_out_gas": 4,
+        "eff": 2,
+        "P_out_gas": 3,
+        "Q_l": 4,
+        "Q_g_out": 5,
     }
 
     def __init__(
@@ -397,9 +398,9 @@ def func(self, c_T_inlet, y_T2_inlet):
 
         c_T_outlet = res["c_T_outlet [mol/m^3]"]
         y_T2_outlet = res["y_T2_outlet_gas"]
-        P_total_outlet = res["total_gas_P_outlet [Pa]"]
-        n_T_out_gas = res["tritium_out_gas [mol/s]"]
-
         eff = res["extraction_efficiency [fraction]"]
+        P_total_outlet = res["total_gas_P_outlet [Pa]"]
+        Q_l = res["liquid_vol_flow [m^3/s]"]
+        Q_g_out = res["gas_vol_flow_outlet [m^3/s]"]
 
-        return c_T_outlet, y_T2_outlet, P_total_outlet, eff, n_T_out_gas
+        return c_T_outlet, y_T2_outlet, eff, P_total_outlet, Q_l, Q_g_out

From f2d6b20856a54f42c37d1234042242d206e72d2a Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Mon, 17 Nov 2025 21:28:01 +0000
Subject: [PATCH 17/30] corrected port map names

---
 src/pathsim_chem/tritium/glc.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index c8d4486..d2f048e 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -353,14 +353,14 @@ class GLC(pathsim.blocks.Function):
     """
 
     _port_map_in = {
-        "c_T_inlet": 0,
+        "c_T_in": 0,
         "y_T2_in": 1,
     }
     _port_map_out = {
-        "c_T_outlet": 0,
+        "c_T_out": 0,
         "y_T2_out": 1,
         "eff": 2,
-        "P_out_gas": 3,
+        "P_out": 3,
         "Q_l": 4,
         "Q_g_out": 5,
     }

From 27e9f7b16d05c490b4cb6090bbb52b75dbe68629 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Mon, 17 Nov 2025 21:49:14 +0000
Subject: [PATCH 18/30] renamed P_in

---
 src/pathsim_chem/tritium/glc.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index d2f048e..aa8991f 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -343,7 +343,7 @@ class GLC(pathsim.blocks.Function):
     More details about the model can be found in: https://doi.org/10.13182/FST95-A30485
 
     Args:
-        P_0: Inlet operating pressure [Pa]
+        P_in: Inlet operating pressure [Pa]
         L: Column height [m]
         flow_g: Gas mass flow rate [kg/s]
         flow_l: Liquid mass flow rate [kg/s]
@@ -367,7 +367,7 @@ class GLC(pathsim.blocks.Function):
 
     def __init__(
         self,
-        P_0,
+        P_in,
         L,
         flow_g,
         flow_l,
@@ -377,7 +377,7 @@ def __init__(
         initial_nb_of_elements=20,
     ):
         self.params = {
-            "P_0": P_0,
+            "P_in": P_in,
             "L": L,
             "Flow_l": flow_l,
             "Flow_g": flow_g,

From b5d528174d6c8751692b4473da3587a275481fd5 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Mon, 17 Nov 2025 21:53:44 +0000
Subject: [PATCH 19/30] Updated parameter name

---
 src/pathsim_chem/tritium/glc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index aa8991f..2e87fed 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -38,7 +38,7 @@ def _calculate_properties(params):
     D = params["D"]
     Flow_l = params["Flow_l"]
     Flow_g = params["Flow_g"]
-    P_0 = params["P_0"]
+    P_0 = params["P_in"]
 
     # --- Fluid Properties (Temperature Dependent) ---
     # TODO add references

From ff9140d11fae7d9a54e5ae65a8b42f41986faf44 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Mon, 17 Nov 2025 22:05:06 +0000
Subject: [PATCH 20/30] Updated P_out names

---
 src/pathsim_chem/tritium/glc.py | 43 ++++++++++++++++++---------------
 1 file changed, 24 insertions(+), 19 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 2e87fed..25fd220 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -38,7 +38,7 @@ def _calculate_properties(params):
     D = params["D"]
     Flow_l = params["Flow_l"]
     Flow_g = params["Flow_g"]
-    P_0 = params["P_in"]
+    P_in = params["P_in"]
 
     # --- Fluid Properties (Temperature Dependent) ---
     # TODO add references
@@ -53,7 +53,7 @@ def _calculate_properties(params):
     # --- Flow Properties ---
     A = np.pi * (D / 2) ** 2  # m^2, Cross-sectional area
     Q_l = Flow_l / rho_l  # m^3/s, Volumetric liquid flow rate
-    Q_g = (Flow_g * R * T) / P_0  # m^3/s, Volumetric gas flow rate at inlet
+    Q_g = (Flow_g * R * T) / P_in  # m^3/s, Volumetric gas flow rate at inlet
     u_l = Q_l / A  # m/s, Superficial liquid velocity
     u_g0 = Q_g / A  # m/s, Superficial gas velocity at inlet
 
@@ -112,7 +112,7 @@ def _calculate_dimensionless_groups(params, phys_props):
     Calculate the dimensionless groups for the ODE system.
 
     Args:
-        params (dict): Dictionary of input parameters including L, T, P_0,
+        params (dict): Dictionary of input parameters including L, T, P_in,
                        and c_T_inlet.
         phys_props (dict): Dictionary of physical properties calculated by
                            _calculate_properties.
@@ -122,7 +122,12 @@ def _calculate_dimensionless_groups(params, phys_props):
               Bo_g, phi_g, psi, nu).
     """
     # Unpack parameters
-    L, T, P_0, c_T_inlet = params["L"], params["T"], params["P_0"], params["c_T_inlet"]
+    L, T, P_in, c_T_inlet = (
+        params["L"],
+        params["T"],
+        params["P_in"],
+        params["c_T_inlet"],
+    )
 
     # Unpack physical properties
     rho_l, K_s, u_l, u_g0, epsilon_g, epsilon_l, E_l, E_g, a, h_l = (
@@ -139,12 +144,12 @@ def _calculate_dimensionless_groups(params, phys_props):
     )
 
     # Calculate dimensionless groups
-    psi = (rho_l * g * epsilon_l * L) / P_0  # Hydrostatic pressure ratio
-    nu = ((c_T_inlet / K_s) ** 2) / P_0  # Equilibrium ratio
+    psi = (rho_l * g * epsilon_l * L) / P_in  # Hydrostatic pressure ratio
+    nu = ((c_T_inlet / K_s) ** 2) / P_in  # Equilibrium ratio
     Bo_l = u_l * L / (epsilon_l * E_l)  # Bodenstein number, liquid
-    phi_l = a * h_l * L / u_l  # Transfer units, liquid (Eq. 8.11)
+    phi_l = a * h_l * L / u_l  # Transfer units, liquid
     Bo_g = u_g0 * L / (epsilon_g * E_g)  # Bodenstein number, gas
-    phi_g = 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0)
+    phi_g = 0.5 * (R * T * c_T_inlet / P_in) * (a * h_l * L / u_g0)
 
     return {
         "Bo_l": Bo_l,
@@ -241,7 +246,7 @@ def _process_results(solution, params, phys_props, dim_params):
         raise RuntimeError("BVP solver failed to converge.")
 
     # Unpack parameters
-    c_T_inlet, P_0, T = params["c_T_inlet"], params["P_0"], params["T"]
+    c_T_inlet, P_in, T = params["c_T_inlet"], params["P_in"], params["T"]
     y_T2_in = params["y_T2_in"]
 
     # Dimensionless results
@@ -252,16 +257,16 @@ def _process_results(solution, params, phys_props, dim_params):
 
     # Dimensional results
     c_T_outlet = x_T_outlet_dimless * c_T_inlet
-    P_outlet = P_0 * (1 - dim_params["psi"])
-    P_T2_out = y_T2_outlet_gas * P_outlet
-    P_T2_in = y_T2_in * P_0
+    P_out = P_in * (1 - dim_params["psi"])
+    P_T2_out = y_T2_outlet_gas * P_out
+    P_T2_in = y_T2_in * P_in
 
     # Mass balance check
     n_T_in_liquid = c_T_inlet * Q_l  # mol/s
     n_T_out_liquid = c_T_outlet * Q_l  # mol/s
     n_T2_in_gas = P_T2_in * Q_g / (R * T)  # mol/s
     n_T_in_gas = n_T2_in_gas * 2  # mol/s
-    Q_g_out = (P_0 * Q_g) / P_outlet  # m3/s
+    Q_g_out = (P_in * Q_g) / P_out  # m3/s
     n_T2_out_gas = P_T2_out * Q_g_out / (R * T)  # mol/s
     n_T_out_gas = n_T2_out_gas * 2  # mol/s
 
@@ -281,8 +286,8 @@ def _process_results(solution, params, phys_props, dim_params):
         "P_T2_inlet_gas [Pa]": P_T2_in,
         "P_T2_outlet_gas [Pa]": P_T2_out,
         "y_T2_outlet_gas": y_T2_outlet_gas,
-        "total_gas_P_inlet [Pa]": P_0,
-        "total_gas_P_outlet [Pa]": P_outlet,
+        "total_gas_P_inlet [Pa]": P_in,
+        "total_gas_P_outlet [Pa]": P_out,
         "liquid_vol_flow [m^3/s]": Q_l,
         "gas_vol_flow_outlet [m^3/s]": Q_g,
     }
@@ -315,13 +320,13 @@ def solve(params):
     phys_props = _calculate_properties(params)
 
     # Pre-solver check for non-physical outlet pressure
-    P_outlet = params["P_0"] - (
+    P_out = params["P_in"] - (
         phys_props["rho_l"] * (1 - phys_props["epsilon_g"]) * g * params["L"]
     )
-    if P_outlet <= 0:
+    if P_out <= 0:
         raise ValueError(
-            f"Calculated gas outlet pressure is non-positive ({P_outlet:.2e} Pa). "
-            "Check input parameters P_0, L, etc."
+            f"Calculated gas outlet pressure is non-positive ({P_out:.2e} Pa). "
+            "Check input parameters P_in, L, etc."
         )
 
     # 2. Calculate dimensionless groups for the ODE system

From 78ee20c4dac88fa1c8b371a08664e93761464a16 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 14:17:37 +0000
Subject: [PATCH 21/30] corrected Q_g_out result

---
 src/pathsim_chem/tritium/glc.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 25fd220..19bbf20 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -252,13 +252,13 @@ def _process_results(solution, params, phys_props, dim_params):
     # Dimensionless results
     x_T_outlet_dimless = solution.y[0, 0]
     Q_l, Q_g = phys_props["Q_l"], phys_props["Q_g"]
-    y_T2_outlet_gas = solution.y[2, -1]
+    y_T2_out = solution.y[2, -1]
     efficiency = 1 - x_T_outlet_dimless
 
     # Dimensional results
     c_T_outlet = x_T_outlet_dimless * c_T_inlet
     P_out = P_in * (1 - dim_params["psi"])
-    P_T2_out = y_T2_outlet_gas * P_out
+    P_T2_out = y_T2_out * P_out
     P_T2_in = y_T2_in * P_in
 
     # Mass balance check
@@ -266,7 +266,7 @@ def _process_results(solution, params, phys_props, dim_params):
     n_T_out_liquid = c_T_outlet * Q_l  # mol/s
     n_T2_in_gas = P_T2_in * Q_g / (R * T)  # mol/s
     n_T_in_gas = n_T2_in_gas * 2  # mol/s
-    Q_g_out = (P_in * Q_g) / P_out  # m3/s
+    Q_g_out = Q_g * (P_in / P_out)  # m3/s
     n_T2_out_gas = P_T2_out * Q_g_out / (R * T)  # mol/s
     n_T_out_gas = n_T2_out_gas * 2  # mol/s
 
@@ -281,15 +281,14 @@ def _process_results(solution, params, phys_props, dim_params):
         "tritium_out_liquid [mol/s]": n_T_out_liquid,
         "tritium_out_gas [mol/s]": n_T_out_gas,
         "extraction_efficiency [fraction]": efficiency,
-        "c_T_inlet [mol/m^3]": c_T_inlet,
         "c_T_outlet [mol/m^3]": c_T_outlet,
         "P_T2_inlet_gas [Pa]": P_T2_in,
         "P_T2_outlet_gas [Pa]": P_T2_out,
-        "y_T2_outlet_gas": y_T2_outlet_gas,
+        "y_T2_outlet_gas": y_T2_out,
         "total_gas_P_inlet [Pa]": P_in,
         "total_gas_P_outlet [Pa]": P_out,
         "liquid_vol_flow [m^3/s]": Q_l,
-        "gas_vol_flow_outlet [m^3/s]": Q_g,
+        "gas_vol_flow_outlet [m^3/s]": Q_g_out,
     }
 
     # Add all calculated parameters to the results dictionary

From 2433cb6c4e8b9d916c9d37d3aaa8fe70a8ffc44f Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 14:30:52 +0000
Subject: [PATCH 22/30] added n_T outputs

---
 src/pathsim_chem/tritium/glc.py | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 19bbf20..57cc0b9 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -367,6 +367,8 @@ class GLC(pathsim.blocks.Function):
         "P_out": 3,
         "Q_l": 4,
         "Q_g_out": 5,
+        "n_T_out_liquid": 6,
+        "n_T_out_gas": 7,
     }
 
     def __init__(
@@ -406,5 +408,16 @@ def func(self, c_T_inlet, y_T2_inlet):
         P_total_outlet = res["total_gas_P_outlet [Pa]"]
         Q_l = res["liquid_vol_flow [m^3/s]"]
         Q_g_out = res["gas_vol_flow_outlet [m^3/s]"]
-
-        return c_T_outlet, y_T2_outlet, eff, P_total_outlet, Q_l, Q_g_out
+        n_T_out_liquid = res["tritium_out_liquid [mol/s]"]
+        n_T_out_gas = res["tritium_out_gas [mol/s]"]
+
+        return (
+            c_T_outlet,
+            y_T2_outlet,
+            eff,
+            P_total_outlet,
+            Q_l,
+            Q_g_out,
+            n_T_out_liquid,
+            n_T_out_gas,
+        )

From c0acda3a0a7029312040f3fd98fac16733d0e38e Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 20:28:06 +0000
Subject: [PATCH 23/30] Added flow_g and flow_l as block inputs

---
 src/pathsim_chem/tritium/glc.py | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 57cc0b9..35ba968 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -358,8 +358,11 @@ class GLC(pathsim.blocks.Function):
 
     _port_map_in = {
         "c_T_in": 0,
-        "y_T2_in": 1,
+        "flow_l":1,
+        "y_T2_in": 2,
+        "flow_g":3,
     }
+
     _port_map_out = {
         "c_T_out": 0,
         "y_T2_out": 1,
@@ -375,30 +378,29 @@ def __init__(
         self,
         P_in,
         L,
-        flow_g,
-        flow_l,
         D,
         T,
         g=const.g,
         initial_nb_of_elements=20,
+        BCs="C-C",
     ):
         self.params = {
             "P_in": P_in,
             "L": L,
-            "Flow_l": flow_l,
-            "Flow_g": flow_g,
             "g": g,
             "D": D,
             "T": T,
             "elements": initial_nb_of_elements,
-            "BCs": "C-C",  # hard coded for now
+            "BCs": BCs,
         }
         super().__init__(func=self.func)
 
-    def func(self, c_T_inlet, y_T2_inlet):
+    def func(self, c_T_inlet, y_T2_inlet, flow_g, flow_l):
         new_params = self.params.copy()
         new_params["c_T_inlet"] = c_T_inlet
         new_params["y_T2_in"] = y_T2_inlet
+        new_params["flow_g"] = flow_g
+        new_params["flow_l"] = flow_l
 
         res, _ = solve(new_params)
 

From 87c2c39951eef1b14439492d785fef190ec50b73 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 20:44:19 +0000
Subject: [PATCH 24/30] renamed Flow_g to flow_g

---
 src/pathsim_chem/tritium/glc.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 35ba968..e3667b7 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -36,8 +36,8 @@ def _calculate_properties(params):
     """
     T = params["T"]
     D = params["D"]
-    Flow_l = params["Flow_l"]
-    Flow_g = params["Flow_g"]
+    flow_l = params["flow_l"]
+    flow_g = params["flow_g"]
     P_in = params["P_in"]
 
     # --- Fluid Properties (Temperature Dependent) ---
@@ -52,8 +52,8 @@ def _calculate_properties(params):
 
     # --- Flow Properties ---
     A = np.pi * (D / 2) ** 2  # m^2, Cross-sectional area
-    Q_l = Flow_l / rho_l  # m^3/s, Volumetric liquid flow rate
-    Q_g = (Flow_g * R * T) / P_in  # m^3/s, Volumetric gas flow rate at inlet
+    Q_l = flow_l / rho_l  # m^3/s, Volumetric liquid flow rate
+    Q_g = (flow_g * R * T) / P_in  # m^3/s, Volumetric gas flow rate at inlet
     u_l = Q_l / A  # m/s, Superficial liquid velocity
     u_g0 = Q_g / A  # m/s, Superficial gas velocity at inlet
 
@@ -349,11 +349,11 @@ class GLC(pathsim.blocks.Function):
     Args:
         P_in: Inlet operating pressure [Pa]
         L: Column height [m]
-        flow_g: Gas mass flow rate [kg/s]
-        flow_l: Liquid mass flow rate [kg/s]
         D: Column diameter [m]
         T: Temperature [K]
         g: Gravitational acceleration [m/s^2], default is 9.81
+        initial_nb_of_elements: Initial number of elements for BVP solver
+        BCs: Boundary conditions type, "C-C" (Closed-Closed) or "O-C" (Open-Closed), default is "C-C"
     """
 
     _port_map_in = {

From 9543edb33f4fd5fc59816b5b6c9a6c8b064855df Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 22:48:36 +0000
Subject: [PATCH 25/30] bug fix

---
 src/pathsim_chem/tritium/glc.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index e3667b7..17bd548 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -380,9 +380,9 @@ def __init__(
         L,
         D,
         T,
+        BCs,
         g=const.g,
         initial_nb_of_elements=20,
-        BCs="C-C",
     ):
         self.params = {
             "P_in": P_in,
@@ -395,12 +395,13 @@ def __init__(
         }
         super().__init__(func=self.func)
 
-    def func(self, c_T_inlet, y_T2_inlet, flow_g, flow_l):
+    def func(self, c_T_in, flow_l, y_T2_inlet, flow_g):
         new_params = self.params.copy()
-        new_params["c_T_inlet"] = c_T_inlet
+        new_params["c_T_in"] = c_T_in
+        new_params["flow_l"] = flow_l
         new_params["y_T2_in"] = y_T2_inlet
         new_params["flow_g"] = flow_g
-        new_params["flow_l"] = flow_l
+
 
         res, _ = solve(new_params)
 

From f55266da7f66ae98b792b6189354984d7137e692 Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 22:54:27 +0000
Subject: [PATCH 26/30] more bug fixes

---
 src/pathsim_chem/tritium/glc.py | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 17bd548..ca7dc52 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -122,11 +122,11 @@ def _calculate_dimensionless_groups(params, phys_props):
               Bo_g, phi_g, psi, nu).
     """
     # Unpack parameters
-    L, T, P_in, c_T_inlet = (
+    L, T, P_in, c_T_in = (
         params["L"],
         params["T"],
         params["P_in"],
-        params["c_T_inlet"],
+        params["c_T_in"],
     )
 
     # Unpack physical properties
@@ -145,11 +145,11 @@ def _calculate_dimensionless_groups(params, phys_props):
 
     # Calculate dimensionless groups
     psi = (rho_l * g * epsilon_l * L) / P_in  # Hydrostatic pressure ratio
-    nu = ((c_T_inlet / K_s) ** 2) / P_in  # Equilibrium ratio
+    nu = ((c_T_in / K_s) ** 2) / P_in  # Equilibrium ratio
     Bo_l = u_l * L / (epsilon_l * E_l)  # Bodenstein number, liquid
     phi_l = a * h_l * L / u_l  # Transfer units, liquid
     Bo_g = u_g0 * L / (epsilon_g * E_g)  # Bodenstein number, gas
-    phi_g = 0.5 * (R * T * c_T_inlet / P_in) * (a * h_l * L / u_g0)
+    phi_g = 0.5 * (R * T * c_T_in / P_in) * (a * h_l * L / u_g0)
 
     return {
         "Bo_l": Bo_l,
@@ -246,7 +246,7 @@ def _process_results(solution, params, phys_props, dim_params):
         raise RuntimeError("BVP solver failed to converge.")
 
     # Unpack parameters
-    c_T_inlet, P_in, T = params["c_T_inlet"], params["P_in"], params["T"]
+    c_T_in, P_in, T = params["c_T_in"], params["P_in"], params["T"]
     y_T2_in = params["y_T2_in"]
 
     # Dimensionless results
@@ -256,14 +256,14 @@ def _process_results(solution, params, phys_props, dim_params):
     efficiency = 1 - x_T_outlet_dimless
 
     # Dimensional results
-    c_T_outlet = x_T_outlet_dimless * c_T_inlet
+    c_T_out = x_T_outlet_dimless * c_T_in
     P_out = P_in * (1 - dim_params["psi"])
     P_T2_out = y_T2_out * P_out
     P_T2_in = y_T2_in * P_in
 
     # Mass balance check
-    n_T_in_liquid = c_T_inlet * Q_l  # mol/s
-    n_T_out_liquid = c_T_outlet * Q_l  # mol/s
+    n_T_in_liquid = c_T_in * Q_l  # mol/s
+    n_T_out_liquid = c_T_out * Q_l  # mol/s
     n_T2_in_gas = P_T2_in * Q_g / (R * T)  # mol/s
     n_T_in_gas = n_T2_in_gas * 2  # mol/s
     Q_g_out = Q_g * (P_in / P_out)  # m3/s
@@ -281,7 +281,7 @@ def _process_results(solution, params, phys_props, dim_params):
         "tritium_out_liquid [mol/s]": n_T_out_liquid,
         "tritium_out_gas [mol/s]": n_T_out_gas,
         "extraction_efficiency [fraction]": efficiency,
-        "c_T_outlet [mol/m^3]": c_T_outlet,
+        "c_T_outlet [mol/m^3]": c_T_out,
         "P_T2_inlet_gas [Pa]": P_T2_in,
         "P_T2_outlet_gas [Pa]": P_T2_out,
         "y_T2_outlet_gas": y_T2_out,
@@ -405,8 +405,8 @@ def func(self, c_T_in, flow_l, y_T2_inlet, flow_g):
 
         res, _ = solve(new_params)
 
-        c_T_outlet = res["c_T_outlet [mol/m^3]"]
-        y_T2_outlet = res["y_T2_outlet_gas"]
+        c_T_out = res["c_T_outlet [mol/m^3]"]
+        y_T2_out = res["y_T2_outlet_gas"]
         eff = res["extraction_efficiency [fraction]"]
         P_total_outlet = res["total_gas_P_outlet [Pa]"]
         Q_l = res["liquid_vol_flow [m^3/s]"]
@@ -415,8 +415,8 @@ def func(self, c_T_in, flow_l, y_T2_inlet, flow_g):
         n_T_out_gas = res["tritium_out_gas [mol/s]"]
 
         return (
-            c_T_outlet,
-            y_T2_outlet,
+            c_T_out,
+            y_T2_out,
             eff,
             P_total_outlet,
             Q_l,

From 1b0b282a271a89347a5e940b94055d54ed2ff05b Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 22:57:51 +0000
Subject: [PATCH 27/30] debug print

---
 src/pathsim_chem/tritium/glc.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index ca7dc52..1a1ea6c 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -328,6 +328,8 @@ def solve(params):
             "Check input parameters P_in, L, etc."
         )
 
+    print(params)
+
     # 2. Calculate dimensionless groups for the ODE system
     dim_params = _calculate_dimensionless_groups(params, phys_props)
 

From b8b8895d84044bb6c57056b5dea1b1356451864e Mon Sep 17 00:00:00 2001
From: rossmacdonald98 <rossmacdonald98@gmail.com>
Date: Tue, 18 Nov 2025 23:39:59 +0000
Subject: [PATCH 28/30] removing debug print

---
 src/pathsim_chem/tritium/glc.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index 1a1ea6c..ca7dc52 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -328,8 +328,6 @@ def solve(params):
             "Check input parameters P_in, L, etc."
         )
 
-    print(params)
-
     # 2. Calculate dimensionless groups for the ODE system
     dim_params = _calculate_dimensionless_groups(params, phys_props)
 

From 64acb257f379cfd4890181bdf7162fcf73f76db7 Mon Sep 17 00:00:00 2001
From: Ross MacDonald <rossmacdonald98@gmail.com>
Date: Wed, 26 Nov 2025 19:58:56 +0000
Subject: [PATCH 29/30] added debug line to BVP failure message

---
 src/pathsim_chem/tritium/glc.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index ca7dc52..ab4f78b 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -242,13 +242,16 @@ def _process_results(solution, params, phys_props, dim_params):
     Returns:
         list: A list containing the results dictionary and the solution object.
     """
-    if not solution.success:
-        raise RuntimeError("BVP solver failed to converge.")
 
     # Unpack parameters
     c_T_in, P_in, T = params["c_T_in"], params["P_in"], params["T"]
     y_T2_in = params["y_T2_in"]
 
+    if not solution.success:
+        print("c_T_in: ",c_T_in)
+        print("y_T2_in :", y_T2_in)
+        raise RuntimeError("BVP solver failed to converge.")
+
     # Dimensionless results
     x_T_outlet_dimless = solution.y[0, 0]
     Q_l, Q_g = phys_props["Q_l"], phys_props["Q_g"]

From 30f595c5dabb8d19f0bf2fe2d1114bd4097b64c6 Mon Sep 17 00:00:00 2001
From: Ross MacDonald <rossmacdonald98@gmail.com>
Date: Wed, 26 Nov 2025 21:10:46 +0000
Subject: [PATCH 30/30] more debug prints

---
 src/pathsim_chem/tritium/glc.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py
index ab4f78b..25df5ff 100644
--- a/src/pathsim_chem/tritium/glc.py
+++ b/src/pathsim_chem/tritium/glc.py
@@ -247,9 +247,14 @@ def _process_results(solution, params, phys_props, dim_params):
     c_T_in, P_in, T = params["c_T_in"], params["P_in"], params["T"]
     y_T2_in = params["y_T2_in"]
 
+    ########### Debugging ##########
+    print("c_T_in: ",c_T_in)
+    print("y_T2_in :", y_T2_in)
+    print("flow_l: ",params["flow_l"])
+    print("flow_g: ",params["flow_g"])
+
+
     if not solution.success:
-        print("c_T_in: ",c_T_in)
-        print("y_T2_in :", y_T2_in)
         raise RuntimeError("BVP solver failed to converge.")
 
     # Dimensionless results