Add balance check for electroneutrality

examples/electrochemistry/vanadium_half_cell_1d.py

@@ -367,7 +367,10 @@ def value_of(name: str) -> float:
 # Charging reduces V³⁺ to V²⁺, so the state of charge grows, and discharging
 # turns it back. At rest only the exchange with the reservoir at the membrane
 # changes the electrode, slowly. The voltage is lower on charge than on
-# discharge, by the overpotentials that drive the reaction either way:
+# discharge, by the overpotentials that drive the reaction either way.
+# Sulfate, which follows from electroneutrality, stays positive, see
+# {py:meth}`minimum_balanced_concentration
+# <sumfem.electrochemistry.RedoxHalfCell.minimum_balanced_concentration>`:
 
 # +
 steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")}
@@ -378,6 +381,7 @@ def value_of(name: str) -> float:
 assert np.abs(np.diff(rest)).max() < np.abs(steps["charge"]).min()
 assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean()
 assert np.isfinite(voltage).all()
+assert cell.minimum_balanced_concentration(disc) > 0
 # -
 
 # The voltage and the state of charge through the schedule:

examples/electrochemistry/vanadium_half_cell_2d.py

@@ -380,7 +380,10 @@
 # Charging raises the state of charge and discharging lowers it. The flow
 # now renews the electrolyte, so at rest it carries the charged electrolyte
 # out and brings the reservoir's in: the state of charge falls back towards
-# that of the reservoir:
+# that of the reservoir. Sulfate, which follows from electroneutrality, stays
+# positive, see
+# {py:meth}`minimum_balanced_concentration
+# <sumfem.electrochemistry.RedoxHalfCell.minimum_balanced_concentration>`:
 
 # +
 steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")}
@@ -392,6 +395,7 @@
 assert np.all(np.diff(rest) < 0)
 assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean()
 assert np.isfinite(voltage).all()
+assert cell.minimum_balanced_concentration(disc) > 0
 # -
 
 # The voltage and the state of charge through the schedule:

examples/electrochemistry/vanadium_half_cell_3d.py

@@ -450,7 +450,10 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None:
 # -
 
 # Charging raises the state of charge and discharging lowers it, and the
-# voltage is lower on charge than on discharge:
+# voltage is lower on charge than on discharge. Sulfate, which follows from
+# electroneutrality, stays positive, see
+# {py:meth}`minimum_balanced_concentration
+# <sumfem.electrochemistry.RedoxHalfCell.minimum_balanced_concentration>`:
 
 # +
 steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")}
@@ -459,6 +462,7 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None:
 assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0)
 assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean()
 assert np.isfinite(voltage).all()
+assert cell.minimum_balanced_concentration(disc) > 0
 # -
 
 # The voltage and the state of charge through the schedule:

src/sumfem/electrochemistry/electrolyte.py

@@ -74,9 +74,14 @@ class Electrolyte:
     an electrolyte share these parameters.
 
     At most one ion may have no field. Its concentration follows from
-    electroneutrality, ``sum_j z_j c_j = 0``. Without ions, the electrolyte
-    has a uniform composition, and only provides the constants, the
-    temperature and the properties of the solution.
+    electroneutrality, ``sum_j z_j c_j = 0``, see :attr:`balanced`. Nothing
+    keeps that concentration positive, so choose the most abundant ion of
+    the supporting electrolyte, e.g. sulfate in a vanadium electrolyte,
+    never an ion of a reaction, and check its minimum, e.g. with
+    :meth:`~sumfem.electrochemistry.RedoxHalfCell.minimum_balanced_concentration`.
+    Without ions, the electrolyte has a uniform composition, and only
+    provides the constants, the temperature and the properties of the
+    solution.
 
     Args:
         name: Prefix of the parameter names.
@@ -231,6 +236,14 @@ def ions(self) -> typing.Mapping[str, Ion]:
         """The ions by name, in the given order, as given."""
         return self._ions
 
+    @property
+    def balanced(self) -> str | None:
+        """Name of the ion that follows from electroneutrality, if any."""
+        for ion in self._ions.values():
+            if ion.concentration is None:
+                return ion.name
+        return None
+
     @property
     def F(self) -> Parameter:
         """Faraday's constant."""

src/sumfem/electrochemistry/half_cell.py

@@ -13,6 +13,10 @@
 import types
 import typing
 
+from mpi4py import MPI
+
+import dolfinx.fem
+import numpy as np
 import sympy
 import sympy.physics.units as syu
 import ufl
@@ -790,6 +794,34 @@ def state_of_charge(self, discretization: Discretization) -> float:
         oxidized = assemble(self._functionals["amount_oxidized"])
         return reduced / max(reduced + oxidized, 1e-30)
 
+    def minimum_balanced_concentration(self, discretization: Discretization) -> float:
+        """The smallest concentration of the ion given by electroneutrality.
+
+        Nothing keeps that concentration positive. A negative minimum means
+        the electrolyte cannot stay neutral with the computed
+        concentrations, e.g. after a time step that is too long, or because
+        the balanced ion is not abundant enough, see
+        :class:`~sumfem.electrochemistry.Electrolyte`.
+
+        Args:
+            discretization: The discretization of the cell.
+
+        Returns:
+            The minimum over the nodes of the concentrations' space, over all
+            processes.
+        """
+        e = self._electrolyte
+        assert e.balanced is not None
+        space = discretization.function_space(next(iter(self._concentrations.values())))
+        balanced = dolfinx.fem.Function(space)
+        expression = discretization[e.concentration(e.balanced, positive=False)]
+        balanced.interpolate(
+            dolfinx.fem.Expression(expression, space.element.interpolation_points)
+        )
+        owned = space.dofmap.index_map.size_local * space.dofmap.index_map_bs
+        local = balanced.x.array[:owned].min(initial=np.inf)
+        return float(space.mesh.comm.allreduce(local, op=MPI.MIN))
+
     def local_state_of_charge(self) -> ufl.core.expr.Expr:
         """The share of the reduced ion in the couple, ``c_red / (c_red + c_ox)``."""
         c_red = self._concentrations[self._couple.reduced]

tests/electrochemistry/half_cell/test_models.py

@@ -401,3 +401,17 @@ def test_two_half_cells_share_the_electrolyte_potential_and_water(reference):
     assert {"negative_c_H", "positive_c_H", "c_VO", "c_V2"} <= set(names)
     assert "negative.rate_constant" in system.parameters
     assert "positive_electrolyte.D_H" in system.parameters
+
+
+def test_the_balanced_ion_is_checked_for_its_minimum(reference):
+    """At rest the vanadium cell holds sulfate at the concentration that
+    neutralizes the others, (2 c_V2 + 3 c_V3 + c_H - c_HSO4) / 2."""
+    marked = reference.interval(num_cells=8)
+    domain = Domain.from_mesh("cell", marked.mesh)
+    cell = reference.half_cell(domain)
+    inlet = cell.inlet_conditions("Membrane")
+    disc = discretize([cell.model, *inlet], DomainMaps({domain: marked}))
+    cell.initialize(disc)
+    assert cell.electrolyte.balanced == "SO4"
+    expected = (2 * 300.0 + 3 * 300.0 + 1000.0 - 1000.0) / 2
+    assert cell.minimum_balanced_concentration(disc) == pytest.approx(expected)