From 5f4605d1ed8e432c914d41c3297e809188d2d0a4 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Thu, 25 Jun 2026 07:02:57 +0200 Subject: [PATCH 1/2] Rebuild GLC block as native BVP1D subclass --- src/pathsim_chem/tritium/glc.py | 346 ++++++++++++++++++++------------ tests/tritium/test_glc.py | 237 ++++++++++++++++++++++ 2 files changed, 456 insertions(+), 127 deletions(-) create mode 100644 tests/tritium/test_glc.py diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py index cf6b3ff..1f26f9b 100644 --- a/src/pathsim_chem/tritium/glc.py +++ b/src/pathsim_chem/tritium/glc.py @@ -3,15 +3,17 @@ 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). +based on the model by C. Malara (1995). The boundary value problem is solved +with the native :class:`~pathsim.blocks.BVP1D` block, which warm-starts the +mesh between evaluations and skips the re-solve when the boundary data is +unchanged. """ import numpy as np -from scipy.integrate import solve_bvp from scipy.optimize import root_scalar import scipy.constants as const -import pathsim -from pathsim.utils.register import Register + +from pathsim.blocks import BVP1D # --- Physical Constants --- g = const.g # m/s^2, Gravitational acceleration @@ -162,71 +164,79 @@ def _calculate_dimensionless_groups(params, phys_props): } -def _solve_bvp_system(dim_params, y_T2_in, BCs, elements): - """ - Set up and solve the Boundary Value Problem for tritium extraction. +def _ode_system(xi, S, dim): + """First-order ODE system (4 equations) of the dimensionless GLC model. - This function defines the system of ordinary differential equations (ODEs) - and the corresponding boundary conditions, then solves the BVP using - `scipy.integrate.solve_bvp`. + The state is ``S = [x_T, dx_T/d(xi), y_T2, dy_T2/d(xi)]``. The signature + matches what the native :class:`~pathsim.blocks.BVP1D` block expects, with + the spatial mesh ``xi`` of shape ``(m,)`` and the state ``S`` of shape + ``(4, m)``. 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. + xi (numpy.ndarray): Dimensionless axial coordinate, shape ``(m,)``. + S (numpy.ndarray): State, shape ``(4, m)``. + dim (dict): Dimensionless groups (Bo_l, phi_l, Bo_g, phi_g, psi, nu). Returns: - scipy.integrate.OdeSolution: The solution object from `solve_bvp`. + numpy.ndarray: Derivatives ``dS/d(xi)``, shape ``(4, m)``. """ 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"], + dim["Bo_l"], + dim["phi_l"], + dim["Bo_g"], + dim["phi_g"], + dim["psi"], + dim["nu"], ) - 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 - ) + 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 _process_results(solution, params, phys_props, dim_params): +def _boundary_conditions(Sa, Sb, dim, y_T2_in, BCs): + """Boundary-condition residuals at xi=0 (Sa) and xi=1 (Sb). + + Args: + Sa (numpy.ndarray): State at xi=0 (liquid outlet). + Sb (numpy.ndarray): State at xi=1 (liquid inlet). + dim (dict): Dimensionless groups (uses Bo_l, Bo_g). + y_T2_in (float): Inlet mole fraction of T2 in the gas phase. + BCs (str): Boundary-condition type, ``"C-C"`` or ``"O-C"``. + + Returns: + numpy.ndarray: The four boundary-condition residuals. + + Raises: + ValueError: If ``BCs`` is not a recognised type. + """ + Bo_l, Bo_g = dim["Bo_l"], dim["Bo_g"] + + 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 + else: + raise ValueError(f"Unknown boundary condition type: {BCs!r}") + + return np.array([res1, res2, res3, res4]) + + +def _process_results(S_ends, params, phys_props, dim_params): """ Process the BVP solution to produce dimensional results. @@ -235,26 +245,25 @@ def _process_results(solution, params, phys_props, dim_params): mass balance check, and aggregates all results. Args: - solution (scipy.integrate.OdeSolution): The BVP solution object. + S_ends (numpy.ndarray): Solution sampled at the domain endpoints, + shape ``(4, 2)`` with column 0 at xi=0 (liquid outlet) and column 1 + at xi=1 (liquid inlet). 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. + dict: The dimensional results dictionary. """ # 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: - raise RuntimeError("BVP solver failed to converge.") - - # Dimensionless results - x_T_outlet_dimless = solution.y[0, 0] + # Dimensionless results (xi=0 -> column 0, xi=1 -> column 1) + x_T_outlet_dimless = S_ends[0, 0] Q_l, Q_g = phys_props["Q_l"], phys_props["Q_g"] - y_T2_out = solution.y[2, -1] + y_T2_out = S_ends[2, 1] efficiency = 1 - x_T_outlet_dimless # Dimensional results @@ -297,13 +306,16 @@ def _process_results(solution, params, phys_props, dim_params): results.update(phys_props) results.update(dim_params) - return results, solution + return results def solve(params): """ Main solver function for the bubble column model. + Builds the physical properties and dimensionless groups, then solves the + boundary value problem with the native :class:`~pathsim.blocks.BVP1D` block. + Args: params (dict): A dictionary of all input parameters for the model, including operational conditions and geometry. @@ -311,8 +323,12 @@ def solve(params): Returns: list: A list containing: - dict: A dictionary containing the simulation results. - - scipy.integrate.OdeSolution: The raw solution object from the - BVP solver. + - pathsim.blocks.BVP1D: The BVP block, exposing the refined mesh + (``.x``) and the sampled solution (``.solution()``). + + Raises: + ValueError: If the calculated gas outlet pressure is non-positive. + RuntimeError: If the BVP solver fails to converge. """ # Adjust inlet gas concentration to avoid numerical instability at zero y_T2_in = max(params["y_T2_in"], 1e-20) @@ -333,16 +349,31 @@ def solve(params): # 2. Calculate dimensionless groups for the ODE system dim_params = _calculate_dimensionless_groups(params, phys_props) - # 3. Solve the boundary value problem - solution = _solve_bvp_system(dim_params, y_T2_in, params["BCs"], params["elements"]) + # 3. Solve the boundary value problem with the native BVP1D block, sampling + # the two domain endpoints (xi=0 and xi=1) + bvp = BVP1D( + fun=lambda x, y, p, u: _ode_system(x, y, dim_params), + bc=lambda ya, yb, p, u: _boundary_conditions( + ya, yb, dim_params, y_T2_in, params["BCs"] + ), + n=4, + domain=(0.0, 1.0), + n_nodes=params["elements"] + 1, + x_eval=np.array([0.0, 1.0]), + tol=1e-5, + max_nodes=10000, + ) + bvp.update(0.0) + if not bvp.success: + raise RuntimeError("BVP solver failed to converge.") # 4. Process and return the results in a dimensional format - [results, solution] = _process_results(solution, params, phys_props, dim_params) + results = _process_results(bvp.solution(), params, phys_props, dim_params) - return [results, solution] + return [results, bvp] -class GLC(pathsim.blocks.Function): +class GLC(BVP1D): r"""Counter-current bubble column gas-liquid contactor (GLC) for tritium extraction. Solves the coupled, non-linear, second-order boundary value problem that @@ -352,10 +383,22 @@ class GLC(pathsim.blocks.Function): transfer via Sieverts' law, and hydrostatic pressure variation along the column. - The block is intended for steady-state tritium extraction calculations - in fusion blanket systems. At each evaluation it computes - temperature-dependent fluid properties, dimensionless groups, and solves - the BVP using ``scipy.integrate.solve_bvp``. + The block is a thin specialisation of the native + :class:`~pathsim.blocks.BVP1D` block: the constructor seeds the parent with + the Malara right-hand side and boundary conditions, and the four block + inputs supply the per-evaluation boundary data. The hydrodynamic + correlations and dimensionless groups are computed from the current input + inside the collocation callbacks. As for any ``BVP1D``, the solve is + warm-started from the previous mesh and skipped entirely when the input is + unchanged. + + The block output is the dimensionless solution sampled at the two domain + endpoints (``xi=0`` liquid outlet, ``xi=1`` liquid inlet), in the + :class:`~pathsim.blocks.BVP1D` row-major layout + ``[x_T(0), x_T(1), x_T'(0), x_T'(1), y_T2(0), y_T2(1), y_T2'(0), y_T2'(1)]``. + Use :meth:`results` to post-process the current solution into dimensional + quantities (concentrations, extraction efficiency, molar flows, mass + balance). Reference: https://doi.org/10.13182/FST95-A30485 @@ -365,16 +408,6 @@ class GLC(pathsim.blocks.Function): ``y_T2_inlet`` -- T₂ mole fraction in inlet gas [-], ``flow_g`` -- gas mass flow rate [kg/s]. - **Output ports:** - ``c_T_out`` -- dissolved tritium concentration in liquid outlet [mol/m³], - ``y_T2_out`` -- T₂ mole fraction in outlet gas [-], - ``eff`` -- extraction efficiency [-], - ``P_out`` -- total gas outlet pressure [Pa], - ``Q_l`` -- liquid volumetric flow rate [m³/s], - ``Q_g_out`` -- gas volumetric flow rate at outlet [m³/s], - ``n_T_out_liquid`` -- tritium molar flow in liquid outlet [mol/s], - ``n_T_out_gas`` -- tritium molar flow in gas outlet [mol/s]. - Parameters ---------- P_in : float @@ -393,6 +426,11 @@ class GLC(pathsim.blocks.Function): Gravitational acceleration [m/s²]. Default: ``scipy.constants.g``. initial_nb_of_elements : int, optional Number of mesh elements for the initial BVP grid. Default: 20. + tol : float, optional + Solver tolerance for the BVP. Default: 1e-5. + max_nodes : int, optional + Maximum number of mesh nodes allowed during BVP refinement. Default: + 10000. """ input_port_labels = { @@ -401,16 +439,6 @@ class GLC(pathsim.blocks.Function): "y_T2_inlet": 2, "flow_g": 3, } - output_port_labels = { - "c_T_out": 0, - "y_T2_out": 1, - "eff": 2, - "P_out": 3, - "Q_l": 4, - "Q_g_out": 5, - "n_T_out_liquid": 6, - "n_T_out_gas": 7, - } def __init__( self, @@ -421,7 +449,11 @@ def __init__( BCs, g=const.g, initial_nb_of_elements=20, + tol=1e-5, + max_nodes=10000, ): + #fixed operating point and geometry; the four block inputs supply the + #per-evaluation boundary data (c_T_in, flow_l, y_T2_inlet, flow_g) self.params = { "P_in": P_in, "L": L, @@ -431,33 +463,93 @@ def __init__( "elements": initial_nb_of_elements, "BCs": BCs, } - super().__init__(func=self.func) - - def func(self, c_T_in, flow_l, y_T2_inlet, flow_g): - new_params = self.params.copy() - 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 - - res, _ = solve(new_params) - - 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]"] - Q_g_out = res["gas_vol_flow_outlet [m^3/s]"] - n_T_out_liquid = res["tritium_out_liquid [mol/s]"] - n_T_out_gas = res["tritium_out_gas [mol/s]"] - - return ( - c_T_out, - y_T2_out, - eff, - P_total_outlet, - Q_l, - Q_g_out, - n_T_out_liquid, - n_T_out_gas, + self.BCs = BCs + + #cache of (physical properties, dimensionless groups) keyed on the input + #so the hydrodynamic correlations run once per operating point instead of + #on every collocation mesh point + self._u_cache = None + self._phys = None + self._dim = None + + #seed the native BVP block with the Malara model physics; sample the two + #domain endpoints so the dimensionless outlet/inlet states are available + super().__init__( + fun=self._ode, + bc=self._bc, + n=4, + domain=(0.0, 1.0), + n_nodes=initial_nb_of_elements + 1, + x_eval=np.array([0.0, 1.0]), + tol=tol, + max_nodes=max_nodes, ) + + def _physics(self, u): + """Physical properties and dimensionless groups for the input `u`. + + Cached on `u` so the hydrodynamic correlations (and the gas hold-up + root solve) run once per operating point rather than on every mesh + point of the collocation solve. + """ + if self._u_cache is not None and np.array_equal(u, self._u_cache): + return self._phys, self._dim + + c_T_in, flow_l, y_T2_inlet, flow_g = u + p = dict(self.params) + p["c_T_in"] = c_T_in + p["flow_l"] = flow_l + p["flow_g"] = flow_g + p["y_T2_in"] = y_T2_inlet + + phys = _calculate_properties(p) + + #guard against a non-physical (non-positive) gas outlet pressure + P_out = p["P_in"] - phys["rho_l"] * (1 - phys["epsilon_g"]) * g * p["L"] + if P_out <= 0: + raise ValueError( + f"Calculated gas outlet pressure is non-positive ({P_out:.2e} Pa). " + "Check input parameters P_in, L, etc." + ) + + dim = _calculate_dimensionless_groups(p, phys) + + self._u_cache = np.array(u, dtype=float) + self._phys, self._dim = phys, dim + return phys, dim + + def _ode(self, xi, S, p, u): + """Malara right-hand side, with dimensionless groups derived from `u`.""" + _, dim = self._physics(u) + return _ode_system(xi, S, dim) + + def _bc(self, Sa, Sb, p, u): + """Malara boundary conditions, with groups derived from `u`.""" + _, dim = self._physics(u) + y_T2_in = max(u[2], 1e-20) + return _boundary_conditions(Sa, Sb, dim, y_T2_in, self.BCs) + + def results(self): + """Post-process the current BVP solution into dimensional results. + + Returns + ------- + dict, None + Dimensional results (outlet concentration, extraction efficiency, + outlet pressure, volumetric flows, tritium molar flows and the mass + balance) for the most recent input, or ``None`` if no successful + solve has happened yet. + """ + if not self.success: + return None + + u = self.inputs.to_array() + c_T_in, flow_l, y_T2_inlet, flow_g = u + p = dict(self.params) + p["c_T_in"] = c_T_in + p["flow_l"] = flow_l + p["flow_g"] = flow_g + p["y_T2_in"] = y_T2_inlet + + phys, dim = self._physics(u) + return _process_results(self.solution(), p, phys, dim) diff --git a/tests/tritium/test_glc.py b/tests/tritium/test_glc.py new file mode 100644 index 0000000..8d5ef83 --- /dev/null +++ b/tests/tritium/test_glc.py @@ -0,0 +1,237 @@ +######################################################################################## +## +## TESTS FOR +## 'tritium.glc.py' +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from scipy.integrate import solve_bvp + +from pathsim_chem.tritium import GLC +from pathsim_chem.tritium import glc + +from pathsim import Simulation, Connection +from pathsim.blocks import BVP1D, Constant, Scope + + +# HELPERS ============================================================================= + +#fixed operating point and geometry used across the tests +_BASE = dict(P_in=2e5, L=1.0, D=0.1, T=623.0, g=9.81) + +#per-evaluation boundary data (c_T_in, flow_l, y_T2_inlet, flow_g) +_INPUT = (1e-3, 1.0, 0.0, 1e-4) + + +def _reference(BCs): + """Independent gold-standard solve of the Malara model with scipy.solve_bvp. + + Returns the dimensionless outlet liquid concentration ``x_T(0)`` and the + dimensionless outlet gas fraction ``y_T2(1)``. + """ + c_T_in, flow_l, y_T2_inlet, flow_g = _INPUT + p = dict(_BASE) + p.update(elements=20, BCs=BCs, c_T_in=c_T_in, flow_l=flow_l, + flow_g=flow_g, y_T2_in=y_T2_inlet) + + phys = glc._calculate_properties(p) + dim = glc._calculate_dimensionless_groups(p, phys) + y2in = max(p["y_T2_in"], 1e-20) + + Bo_l, phi_l, Bo_g, phi_g, psi, nu = ( + dim["Bo_l"], dim["phi_l"], dim["Bo_g"], dim["phi_g"], dim["psi"], dim["nu"] + ) + + def ode(xi, S): + x_T, dx, y_T2, dy = S + th = x_T - np.sqrt(np.maximum(0, (1 - psi * xi) * y_T2 / nu)) + return np.vstack(( + dx, + Bo_l * (phi_l * th - dx), + dy, + (Bo_g / (1 - psi * xi)) * ((1 + 2 * psi / Bo_g) * dy - phi_g * th), + )) + + def bc(Sa, Sb): + if BCs == "C-C": + return np.array([Sa[1], Sb[0] - (1 - (1 / Bo_l) * Sb[1]), + Sa[2] - y2in - (1 / Bo_g) * Sa[3], Sb[3]]) + return np.array([Sa[1], Sb[0] - 1.0, Sa[2] - y2in, Sb[3]]) + + xi = np.linspace(0, 1, 21) + sol = solve_bvp(ode, bc, xi, np.zeros((4, 21)), tol=1e-5, max_nodes=10000) + return sol.y[0, 0], sol.y[2, -1] + + +# TESTS =============================================================================== + +class TestGLC(unittest.TestCase): + """ + Test the implementation of the 'GLC' bubble-column gas-liquid contactor + block from the fusion toolbox. The block inherits from the native `BVP1D` + block and seeds it with the Malara (1995) tritium-extraction model. + """ + + def test_init(self): + """The block is a BVP1D subclass seeded with the GLC problem dimensions.""" + blk = GLC(BCs="C-C", **_BASE) + self.assertIsInstance(blk, BVP1D) + self.assertEqual(blk.n, 4) + self.assertEqual(blk.domain, (0.0, 1.0)) + #four named input ports supply the boundary data + self.assertEqual(len(blk.inputs), 4) + for label in ("c_T_in", "flow_l", "y_T2_inlet", "flow_g"): + self.assertIn(label, blk.input_port_labels) + #no successful solve yet, so no dimensional results + self.assertIsNone(blk.results()) + + + def test_solve_matches_reference_cc(self): + """The C-C solve matches an independent scipy.solve_bvp reference.""" + blk = GLC(BCs="C-C", **_BASE) + blk.inputs.update_from_array(np.array(_INPUT)) + blk.update(0.0) + + self.assertTrue(blk.success) + x_T_ref, y_T2_ref = _reference("C-C") + res = blk.results() + self.assertAlmostEqual(res["c_T_outlet [mol/m^3]"], + x_T_ref * _INPUT[0], places=12) + self.assertAlmostEqual(res["y_T2_outlet_gas"], y_T2_ref, places=12) + + + def test_solve_matches_reference_oc(self): + """The O-C solve matches an independent scipy.solve_bvp reference.""" + blk = GLC(BCs="O-C", **_BASE) + blk.inputs.update_from_array(np.array(_INPUT)) + blk.update(0.0) + + self.assertTrue(blk.success) + x_T_ref, y_T2_ref = _reference("O-C") + res = blk.results() + self.assertAlmostEqual(res["c_T_outlet [mol/m^3]"], + x_T_ref * _INPUT[0], places=12) + self.assertAlmostEqual(res["y_T2_outlet_gas"], y_T2_ref, places=12) + + + def test_results_physical(self): + """The dimensional results are physically sensible and mass conserving.""" + blk = GLC(BCs="C-C", **_BASE) + blk.inputs.update_from_array(np.array(_INPUT)) + blk.update(0.0) + res = blk.results() + + #extraction efficiency in [0, 1] and outlet below inlet concentration + self.assertGreaterEqual(res["extraction_efficiency [fraction]"], 0.0) + self.assertLessEqual(res["extraction_efficiency [fraction]"], 1.0) + self.assertLess(res["c_T_outlet [mol/m^3]"], _INPUT[0]) + + #hydrostatic head drops the gas pressure below the inlet pressure + self.assertLess(res["total_gas_P_outlet [Pa]"], _BASE["P_in"]) + + #closed tritium mass balance + self.assertAlmostEqual(res["Total tritium in [mol/s]"], + res["Total tritium out [mol/s]"], places=12) + + + def test_output_layout(self): + """The block output is the dimensionless solution at the two endpoints.""" + blk = GLC(BCs="C-C", **_BASE) + blk.inputs.update_from_array(np.array(_INPUT)) + blk.update(0.0) + + #row-major [x_T(0), x_T(1), ...]; output[0] is the dimensionless outlet + #liquid concentration, i.e. c_T_out / c_T_in + res = blk.results() + self.assertAlmostEqual(blk.outputs[0], + res["c_T_outlet [mol/m^3]"] / _INPUT[0], places=12) + + + def test_unknown_bc_raises(self): + """An unknown boundary-condition type is rejected at solve time.""" + blk = GLC(BCs="bogus", **_BASE) + blk.inputs.update_from_array(np.array(_INPUT)) + with self.assertRaises(ValueError): + blk.update(0.0) + + + def test_input_unchanged_skips_resolve(self): + """A re-evaluation with unchanged input does not re-solve the BVP.""" + import pathsim.blocks.bvp as bvpmod + + blk = GLC(BCs="C-C", **_BASE) + blk.inputs.update_from_array(np.array(_INPUT)) + + calls = [0] + orig = bvpmod.solve_bvp + def _counting(*a, **k): + calls[0] += 1 + return orig(*a, **k) + bvpmod.solve_bvp = _counting + try: + blk.update(0.0) #first solve + blk.update(0.0) #unchanged input -> skipped + self.assertEqual(calls[0], 1) + + blk.inputs[3] = 2e-4 #flow_g + blk.update(0.0) #changed input -> solve again + self.assertEqual(calls[0], 2) + finally: + bvpmod.solve_bvp = orig + + + def test_solve_function(self): + """The standalone solve() helper returns results and the BVP block.""" + p = dict(_BASE) + p.update(elements=20, BCs="C-C", c_T_in=_INPUT[0], flow_l=_INPUT[1], + y_T2_in=_INPUT[2], flow_g=_INPUT[3]) + results, bvp = glc.solve(p) + + self.assertIsInstance(bvp, BVP1D) + self.assertTrue(bvp.success) + x_T_ref, _ = _reference("C-C") + self.assertAlmostEqual(results["c_T_outlet [mol/m^3]"], + x_T_ref * _INPUT[0], places=12) + + + def test_non_physical_pressure_raises(self): + """A column tall enough to drive the outlet pressure non-positive fails.""" + #very tall column at low inlet pressure -> hydrostatic head exceeds P_in + blk = GLC(P_in=1e4, L=10.0, D=0.1, T=623.0, BCs="C-C") + blk.inputs.update_from_array(np.array(_INPUT)) + with self.assertRaises(ValueError): + blk.update(0.0) + + + def test_simulation(self): + """The block runs inside a Simulation driven by constant sources.""" + srcs = [Constant(v) for v in _INPUT] + block = GLC(BCs="C-C", **_BASE) + sco = Scope() + + sim = Simulation( + blocks=[*srcs, block, sco], + connections=( + [Connection(srcs[i], block[i]) for i in range(4)] + + [Connection(block[0], sco[0])] + ), + log=False, + ) + sim.run(0.05) + + self.assertTrue(block.success) + res = block.results() + np.testing.assert_allclose( + block.outputs[0], res["c_T_outlet [mol/m^3]"] / _INPUT[0], atol=1e-12 + ) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From 7e42d4b0e44fd87439a1eaaff09d39571d839574 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Thu, 25 Jun 2026 07:11:02 +0200 Subject: [PATCH 2/2] Restore dimensional output ports on GLC via update override --- src/pathsim_chem/tritium/glc.py | 81 ++++++++++++++++++++++++++------- tests/tritium/test_glc.py | 31 +++++++++---- 2 files changed, 88 insertions(+), 24 deletions(-) diff --git a/src/pathsim_chem/tritium/glc.py b/src/pathsim_chem/tritium/glc.py index 1f26f9b..0764775 100644 --- a/src/pathsim_chem/tritium/glc.py +++ b/src/pathsim_chem/tritium/glc.py @@ -383,22 +383,18 @@ class GLC(BVP1D): transfer via Sieverts' law, and hydrostatic pressure variation along the column. - The block is a thin specialisation of the native - :class:`~pathsim.blocks.BVP1D` block: the constructor seeds the parent with - the Malara right-hand side and boundary conditions, and the four block - inputs supply the per-evaluation boundary data. The hydrodynamic - correlations and dimensionless groups are computed from the current input - inside the collocation callbacks. As for any ``BVP1D``, the solve is - warm-started from the previous mesh and skipped entirely when the input is - unchanged. - - The block output is the dimensionless solution sampled at the two domain - endpoints (``xi=0`` liquid outlet, ``xi=1`` liquid inlet), in the - :class:`~pathsim.blocks.BVP1D` row-major layout - ``[x_T(0), x_T(1), x_T'(0), x_T'(1), y_T2(0), y_T2(1), y_T2'(0), y_T2'(1)]``. - Use :meth:`results` to post-process the current solution into dimensional - quantities (concentrations, extraction efficiency, molar flows, mass - balance). + The block is a specialisation of the native :class:`~pathsim.blocks.BVP1D` + block: the constructor seeds the parent with the Malara right-hand side and + boundary conditions, and the four block inputs supply the per-evaluation + boundary data. The hydrodynamic correlations and dimensionless groups are + computed from the current input inside the collocation callbacks. As for any + ``BVP1D``, the solve is warm-started from the previous mesh and skipped + entirely when the input is unchanged. After each solve the dimensionless + endpoint solution is post-processed into the dimensional output ports. + + Use :meth:`results` to retrieve the full dimensional result dictionary, + which additionally contains partial pressures, physical properties and the + dimensionless groups. Reference: https://doi.org/10.13182/FST95-A30485 @@ -408,6 +404,16 @@ class GLC(BVP1D): ``y_T2_inlet`` -- T₂ mole fraction in inlet gas [-], ``flow_g`` -- gas mass flow rate [kg/s]. + **Output ports:** + ``c_T_out`` -- dissolved tritium concentration in liquid outlet [mol/m³], + ``y_T2_out`` -- T₂ mole fraction in outlet gas [-], + ``eff`` -- extraction efficiency [-], + ``P_out`` -- total gas outlet pressure [Pa], + ``Q_l`` -- liquid volumetric flow rate [m³/s], + ``Q_g_out`` -- gas volumetric flow rate at outlet [m³/s], + ``n_T_out_liquid`` -- tritium molar flow in liquid outlet [mol/s], + ``n_T_out_gas`` -- tritium molar flow in gas outlet [mol/s]. + Parameters ---------- P_in : float @@ -439,6 +445,16 @@ class GLC(BVP1D): "y_T2_inlet": 2, "flow_g": 3, } + output_port_labels = { + "c_T_out": 0, + "y_T2_out": 1, + "eff": 2, + "P_out": 3, + "Q_l": 4, + "Q_g_out": 5, + "n_T_out_liquid": 6, + "n_T_out_gas": 7, + } def __init__( self, @@ -529,6 +545,39 @@ def _bc(self, Sa, Sb, p, u): y_T2_in = max(u[2], 1e-20) return _boundary_conditions(Sa, Sb, dim, y_T2_in, self.BCs) + def update(self, t): + """Solve the BVP and expose the dimensional results at the output ports. + + Delegates the solve to :class:`~pathsim.blocks.BVP1D` (warm-started, and + skipped when the input is unchanged), then post-processes the + dimensionless endpoint solution into the eight dimensional output ports. + + Parameters + ---------- + t : float + evaluation time + + Raises + ------ + RuntimeError + if the BVP solve did not converge + """ + super().update(t) + if not self.success: + raise RuntimeError("BVP solver failed to converge.") + + res = self.results() + self.outputs.update_from_array(np.array([ + res["c_T_outlet [mol/m^3]"], + res["y_T2_outlet_gas"], + res["extraction_efficiency [fraction]"], + res["total_gas_P_outlet [Pa]"], + res["liquid_vol_flow [m^3/s]"], + res["gas_vol_flow_outlet [m^3/s]"], + res["tritium_out_liquid [mol/s]"], + res["tritium_out_gas [mol/s]"], + ])) + def results(self): """Post-process the current BVP solution into dimensional results. diff --git a/tests/tritium/test_glc.py b/tests/tritium/test_glc.py index 8d5ef83..6cc7407 100644 --- a/tests/tritium/test_glc.py +++ b/tests/tritium/test_glc.py @@ -87,6 +87,10 @@ def test_init(self): self.assertEqual(len(blk.inputs), 4) for label in ("c_T_in", "flow_l", "y_T2_inlet", "flow_g"): self.assertIn(label, blk.input_port_labels) + #eight named dimensional output ports + for label in ("c_T_out", "y_T2_out", "eff", "P_out", "Q_l", "Q_g_out", + "n_T_out_liquid", "n_T_out_gas"): + self.assertIn(label, blk.output_port_labels) #no successful solve yet, so no dimensional results self.assertIsNone(blk.results()) @@ -139,17 +143,25 @@ def test_results_physical(self): res["Total tritium out [mol/s]"], places=12) - def test_output_layout(self): - """The block output is the dimensionless solution at the two endpoints.""" + def test_output_ports(self): + """The eight dimensional output ports match the results dictionary.""" blk = GLC(BCs="C-C", **_BASE) blk.inputs.update_from_array(np.array(_INPUT)) blk.update(0.0) - #row-major [x_T(0), x_T(1), ...]; output[0] is the dimensionless outlet - #liquid concentration, i.e. c_T_out / c_T_in res = blk.results() - self.assertAlmostEqual(blk.outputs[0], - res["c_T_outlet [mol/m^3]"] / _INPUT[0], places=12) + for label, key in ( + ("c_T_out", "c_T_outlet [mol/m^3]"), + ("y_T2_out", "y_T2_outlet_gas"), + ("eff", "extraction_efficiency [fraction]"), + ("P_out", "total_gas_P_outlet [Pa]"), + ("Q_l", "liquid_vol_flow [m^3/s]"), + ("Q_g_out", "gas_vol_flow_outlet [m^3/s]"), + ("n_T_out_liquid", "tritium_out_liquid [mol/s]"), + ("n_T_out_gas", "tritium_out_gas [mol/s]"), + ): + idx = blk.output_port_labels[label] + self.assertAlmostEqual(blk.outputs[idx], res[key], places=12) def test_unknown_bc_raises(self): @@ -218,7 +230,9 @@ def test_simulation(self): blocks=[*srcs, block, sco], connections=( [Connection(srcs[i], block[i]) for i in range(4)] - + [Connection(block[0], sco[0])] + #connect the named efficiency output port, as a downstream + #consumer would + + [Connection(block["eff"], sco[0])] ), log=False, ) @@ -227,7 +241,8 @@ def test_simulation(self): self.assertTrue(block.success) res = block.results() np.testing.assert_allclose( - block.outputs[0], res["c_T_outlet [mol/m^3]"] / _INPUT[0], atol=1e-12 + block.outputs[block.output_port_labels["eff"]], + res["extraction_efficiency [fraction]"], atol=1e-12 )