Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
milanofthe merged 1 commit into from
Jun 25, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 18 additions & 90 deletions src/pathsim_chem/tritium/glc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
30 changes: 12 additions & 18 deletions tests/tritium/test_glc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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):
Expand Down Expand Up @@ -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
Expand Down
Loading