Merge pull request #29 from pathsim/cleanup/glc-remove-solve

Clean up GLC module: drop solve(), slim results, honest mass balance

Author: milanofthe

src/pathsim_chem/tritium/glc.py

@@ -272,106 +272,34 @@ def _process_results(S_ends, params, phys_props, dim_params):
     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_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
+    # Tritium molar flows from the converged solution (mol/s). Each T2 molecule
+    # carries two tritium atoms.
+    n_T_in_liquid = c_T_in * Q_l
+    n_T_in_gas = 2 * P_T2_in * Q_g / (R * T)
+    n_T_out_liquid = c_T_out * Q_l
     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
+    n_T_out_gas = 2 * P_T2_out * Q_g_out / (R * T)
 
-    # 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)
+    n_T_in = n_T_in_liquid + n_T_in_gas
+    n_T_out = n_T_out_liquid + n_T_out_gas
 
-    results = {
-        "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,
+    # Mass-balance residual reported as a diagnostic. For closed-closed (C-C)
+    # boundary conditions it is at numerical-noise level; for open-closed (O-C)
+    # it reflects the dispersive flux neglected by the open inlet condition.
+    return {
         "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,
-        "total_gas_P_inlet [Pa]": P_in,
+        "extraction_efficiency [fraction]": efficiency,
         "total_gas_P_outlet [Pa]": P_out,
         "liquid_vol_flow [m^3/s]": Q_l,
         "gas_vol_flow_outlet [m^3/s]": Q_g_out,
+        "tritium_out_liquid [mol/s]": n_T_out_liquid,
+        "tritium_out_gas [mol/s]": n_T_out_gas,
+        "Total tritium in [mol/s]": n_T_in,
+        "Total tritium out [mol/s]": n_T_out,
+        "mass_balance_residual [mol/s]": n_T_in - n_T_out,
     }
 
-    # Add all calculated parameters to the results dictionary
-    results.update(phys_props)
-    results.update(dim_params)
-
-    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.
-
-    Returns:
-        list: A list containing:
-              - dict: A dictionary containing the simulation results.
-              - 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)
-
-    # 1. Calculate physical, hydrodynamic, and mass transfer properties
-    phys_props = _calculate_properties(params)
-
-    # Pre-solver check for non-physical outlet pressure
-    P_out = params["P_in"] - (
-        phys_props["rho_l"] * (1 - phys_props["epsilon_g"]) * g * params["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."
-        )
-
-    # 2. Calculate dimensionless groups for the ODE system
-    dim_params = _calculate_dimensionless_groups(params, phys_props)
-
-    # 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 = _process_results(bvp.solution(), params, phys_props, dim_params)
-
-    return [results, bvp]
-
 
 class GLC(BVP1D):
     r"""Counter-current bubble column gas-liquid contactor (GLC) for tritium extraction.

tests/tritium/test_glc.py

@@ -122,9 +122,16 @@ def test_solve_matches_reference_oc(self):
                                x_T_ref * _INPUT[0], places=12)
         self.assertAlmostEqual(res["y_T2_outlet_gas"], y_T2_ref, places=12)
 
+        #the open inlet condition neglects the dispersive flux, so O-C does not
+        #close the mass balance exactly; the residual is small but non-trivial
+        #(a few per mille) and is reported honestly rather than hidden
+        rel = abs(res["mass_balance_residual [mol/s]"]) / res["Total tritium in [mol/s]"]
+        self.assertGreater(rel, 1e-6)
+        self.assertLess(rel, 1e-2)
+
 
     def test_results_physical(self):
-        """The dimensional results are physically sensible and mass conserving."""
+        """The dimensional results are physically sensible."""
         blk = GLC(BCs="C-C", **_BASE)
         blk.inputs.update_from_array(np.array(_INPUT))
         blk.update(0.0)
@@ -138,9 +145,10 @@ def test_results_physical(self):
         #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)
+        #closed-closed boundary conditions conserve tritium to numerical noise;
+        #the residual is reported, not redistributed back into the outputs
+        rel = abs(res["mass_balance_residual [mol/s]"]) / res["Total tritium in [mol/s]"]
+        self.assertLess(rel, 1e-6)
 
 
     def test_output_ports(self):
@@ -197,20 +205,6 @@ def _counting(*a, **k):
             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