Resolve https://github.com/scientificcomputing/sumfem/issues/20

examples/electrochemistry/vanadium_half_cell_1d.py

@@ -139,8 +139,8 @@
 #
 # The ions of the electrolyte are given with their charges and diffusivities,
 # in their units, which sumfem converts to SI, see {ref}`units-checking`. The
-# vanadium ions take the positive part of their concentrations in the
-# kinetics and the conductivity, as they may be depleted. The concentrations
+# kinetics take the positive part of the vanadium concentrations, as they may
+# be depleted. The concentrations
 # of the reservoir name the ions that are transported, each with a field and
 # the equations of {eq}`half-cell`. Sulfate has none: electroneutrality,
 # the last equation of {eq}`half-cell`, with the charges of the ions,
@@ -171,8 +171,8 @@
 diffusivity = syu.m**2 / syu.second
 concentration = syu.mol / syu.m**3
 ions = (
-    Ion("V2", 2, 8.10e-12 * diffusivity, positive=True),
-    Ion("V3", 3, 1.65e-11 * diffusivity, positive=True),
+    Ion("V2", 2, 8.10e-12 * diffusivity),
+    Ion("V3", 3, 1.65e-11 * diffusivity),
     Ion("H", 1, 9.312e-9 * diffusivity),
     Ion("HSO4", -1, 1.33e-9 * diffusivity),
     Ion("SO4", -2, 1.065e-9 * diffusivity),

examples/electrochemistry/vanadium_half_cell_2d.py

@@ -212,8 +212,8 @@
 diffusivity = syu.m**2 / syu.second
 concentration = syu.mol / syu.m**3
 ions = (
-    Ion("V2", 2, 8.10e-12 * diffusivity, positive=True),
-    Ion("V3", 3, 1.65e-11 * diffusivity, positive=True),
+    Ion("V2", 2, 8.10e-12 * diffusivity),
+    Ion("V3", 3, 1.65e-11 * diffusivity),
     Ion("H", 1, 9.312e-9 * diffusivity),
     Ion("HSO4", -1, 1.33e-9 * diffusivity),
     Ion("SO4", -2, 1.065e-9 * diffusivity),

examples/electrochemistry/vanadium_half_cell_3d.py

@@ -322,8 +322,8 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None:
 diffusivity = syu.m**2 / syu.second
 concentration = syu.mol / syu.m**3
 ions = (
-    Ion("V2", 2, 8.10e-12 * diffusivity, positive=True),
-    Ion("V3", 3, 1.65e-11 * diffusivity, positive=True),
+    Ion("V2", 2, 8.10e-12 * diffusivity),
+    Ion("V3", 3, 1.65e-11 * diffusivity),
     Ion("H", 1, 9.312e-9 * diffusivity),
     Ion("HSO4", -1, 1.33e-9 * diffusivity),
     Ion("SO4", -2, 1.065e-9 * diffusivity),

src/sumfem/electrochemistry/electrolyte.py

@@ -39,16 +39,12 @@ class Ion:
             default value of the parameter the electrolyte creates.
         concentration: The field of a transported ion. ``None`` for the one
             ion whose concentration follows from electroneutrality.
-        positive: Use a smooth positive part of the concentration in the
-            conductivity and in reaction kinetics. Use it for ions that can
-            be depleted, such as the active species.
     """
 
     name: str
     charge: float
     diffusivity: float | sympy.Expr
     concentration: Field | None = None
-    positive: bool = False
 
     @property
     def field(self) -> Field:
@@ -278,26 +274,30 @@ def __getitem__(self, name: str) -> Ion:
         return self._ions[name]
 
     def positive_part(self, c: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
-        """The smooth positive part ``(c + sqrt(c^2 + floor^2)) / 2``."""
+        """The smooth positive part ``(c + sqrt(c^2 + floor^2)) / 2``.
+
+        Reaction kinetics take it of the concentrations they react, see
+        :class:`~sumfem.electrochemistry.ButlerVolmer`, so that a
+        concentration that undershoots zero neither reverses a reaction nor
+        enters a logarithm.
+
+        Args:
+            c: A concentration.
+        """
         return 0.5 * (c + ufl.sqrt(c * c + self._floor * self._floor))
 
-    def concentration(
-        self, name: str, positive: bool | None = None
-    ) -> ufl.core.expr.Expr:
-        """Concentration of an ion, as used in conductivity and kinetics.
+    def concentration(self, name: str) -> ufl.core.expr.Expr:
+        """Concentration of an ion, with its sign.
 
         That is the field of a transported ion, and the expression given by
-        electroneutrality for the ion without one.
+        electroneutrality for the ion without one. The conductivity uses it
+        as it is, as Nernst-Planck migration does, so the current is the
+        charge the ions carry, see :meth:`current_density`.
 
         Args:
             name: Name of the ion.
-            positive: Take the positive part of the concentration. Defaults
-                to :attr:`Ion.positive` of the ion.
         """
-        concentration = self._concentrations[name]
-        if self._ions[name].positive if positive is None else positive:
-            return self.positive_part(concentration)
-        return concentration
+        return self._concentrations[name]
 
     def diffusivity(self, name: str) -> Parameter:
         """Diffusivity of an ion in free electrolyte, the parameter ``D_<ion>``.
@@ -331,9 +331,7 @@ def diffusion_current(
             correction: See :meth:`conductivity`.
         """
         total = sum(
-            s.charge
-            * self.diffusivity(s.name)
-            * ufl.grad(self.concentration(s.name, positive=False))
+            s.charge * self.diffusivity(s.name) * ufl.grad(self.concentration(s.name))
             for s in self._ions.values()
         )
         return self._F * correction * total
@@ -345,6 +343,12 @@ def current_density(
     ) -> ufl.core.expr.Expr:
         """Current density ``-kappa grad(phi) - F sum_j z_j D_j grad(c_j)``.
 
+        It is the charge the ions carry, ``F sum_j z_j N_j``, with the
+        Nernst-Planck fluxes ``N_j`` of every ion, also the one given by
+        electroneutrality: the conductivity and the migration use the same
+        concentrations, and advection carries no charge, as
+        ``sum_j z_j c_j = 0``.
+
         Args:
             potential: The electrolyte potential.
             correction: See :meth:`conductivity`.

src/sumfem/electrochemistry/half_cell.py

@@ -814,7 +814,7 @@ def minimum_balanced_concentration(self, discretization: Discretization) -> floa
         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)]
+        expression = discretization[e.concentration(e.balanced)]
         balanced.interpolate(
             dolfinx.fem.Expression(expression, space.element.interpolation_points)
         )

src/sumfem/electrochemistry/kinetics.py

@@ -371,10 +371,7 @@ def __init__(
                 )
 
         def concentrations(ions: dict[str, int]) -> dict[ufl.core.expr.Expr, int]:
-            return {
-                electrolyte.concentration(ion, positive=False): nu
-                for ion, nu in ions.items()
-            }
+            return {electrolyte.concentration(ion): nu for ion, nu in ions.items()}
 
         super().__init__(
             name,
@@ -433,8 +430,8 @@ def _bulk(self) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]:
         """Positive parts of the bulk concentrations of the reduced and oxidized ion."""
         e = self._electrolyte
         return (
-            e.concentration(self._reduced.name, positive=True),
-            e.concentration(self._oxidized.name, positive=True),
+            e.positive_part(e.concentration(self._reduced.name)),
+            e.positive_part(e.concentration(self._oxidized.name)),
         )
 
     def rate_constants(self) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]:

tests/electrochemistry/half_cell/conftest.py

@@ -366,8 +366,8 @@ def darcy_velocity(marked: MeshData, p):
 CONCENTRATION = syu.mol / syu.m**3
 
 IONS = (
-    Ion("V2", 2.0, 8.10e-12 * DIFFUSIVITY, positive=True),
-    Ion("V3", 3.0, 1.65e-11 * DIFFUSIVITY, positive=True),
+    Ion("V2", 2.0, 8.10e-12 * DIFFUSIVITY),
+    Ion("V3", 3.0, 1.65e-11 * DIFFUSIVITY),
     Ion("H", 1.0, 9.312e-9 * DIFFUSIVITY),
     Ion("HSO4", -1.0, 1.33e-9 * DIFFUSIVITY),
     Ion("SO4", -2.0, 1.065e-9 * DIFFUSIVITY),

tests/electrochemistry/half_cell/test_models.py

@@ -257,8 +257,8 @@ def test_the_reservoir_lacks_exactly_one_ion(reference, missing):
 DIFFUSIVITY_UNIT = syu.m**2 / syu.second
 
 POSITIVE_IONS = (
-    Ion("VO", 2, 3.9e-10 * DIFFUSIVITY_UNIT, positive=True),
-    Ion("VO2", 1, 3.9e-10 * DIFFUSIVITY_UNIT, positive=True),
+    Ion("VO", 2, 3.9e-10 * DIFFUSIVITY_UNIT),
+    Ion("VO2", 1, 3.9e-10 * DIFFUSIVITY_UNIT),
     Ion("H", 1, 9.312e-9 * DIFFUSIVITY_UNIT),
     Ion("HSO4", -1, 1.33e-9 * DIFFUSIVITY_UNIT),
     Ion("H2O", 0, 2.3e-9 * DIFFUSIVITY_UNIT),

tests/electrochemistry/test_electrochemistry_from_ufl.py

@@ -52,7 +52,7 @@ def test_conduction_and_transport_in_a_porous_region(cell, box, compare):
     electroneutrality, and Nernst-Planck transport of the three, with
     advection of V2+, in the pores of an electrode on the left half.
 
-    The conductivity takes the positive parts of the vanadium ions, and every
+    The conductivity takes the concentrations as they are, and every
     diffusivity Bruggeman's correction of the porosity.
     """
     gdim = ufl.Cell(cell).topological_dimension
@@ -61,19 +61,14 @@ def test_conduction_and_transport_in_a_porous_region(cell, box, compare):
     W = ufl.MixedFunctionSpace(V, V, V, V)
     phi, c2, c3, h = (ufl.Coefficient(V) for _ in range(4))
     v_phi, w2, w3, wh = ufl.TestFunctions(W)
-    F, R, T, floor, eps, D2, D3, DH, DS = (ufl.Constant(mesh) for _ in range(9))
+    F, R, T, eps, D2, D3, DH, DS = (ufl.Constant(mesh) for _ in range(8))
     sulfate = -(2 * c2 + 3 * c3 + h) / -2
     correction = eps**1.5
     conductivity = (
         F**2
         / (R * T)
         * correction
-        * (
-            4 * D2 * positive_part(c2, floor)
-            + 9 * D3 * positive_part(c3, floor)
-            + DH * h
-            + 4 * DS * sulfate
-        )
+        * (4 * D2 * c2 + 9 * D3 * c3 + DH * h + 4 * DS * sulfate)
     )
     diffusion = (
         F
@@ -103,7 +98,7 @@ def flux(c, z, D, advected=False):
     names = {
         DOMAIN: mesh,
         **{n: c for n, c in zip(("phi", "c2", "c3", "h"), (phi, c2, c3, h))},
-        **dict(F=F, R=R, T=T, floor=floor, eps=eps, D2=D2, D3=D3, DH=DH, DS=DS),
+        **dict(F=F, R=R, T=T, eps=eps, D2=D2, D3=D3, DH=DH, DS=DS),
     }
     plain = UFLResidual(
         "plain",
@@ -117,12 +112,12 @@ def flux(c, z, D, advected=False):
     potential = Field("phi", domain, ("Lagrange", 1))
     fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c2", "c3", "h")]
     ions = [
-        Ion("V2", 2, 2.4e-10, fields[0], positive=True),
-        Ion("V3", 3, 2.4e-10, fields[1], positive=True),
+        Ion("V2", 2, 2.4e-10, fields[0]),
+        Ion("V3", 3, 2.4e-10, fields[1]),
         Ion("H", 1, 9.3e-9, fields[2]),
         Ion("SO4", -2, 1.065e-9),
     ]
-    electrolyte = Electrolyte("e", domain, ions, concentration_floor=0.1)
+    electrolyte = Electrolyte("e", domain, ions)
     region = PorousElectrode("el", domain, ("LeftHalf",), porosity=0.8).region
     model = ElectrolyteConduction("charge", potential, electrolyte, region)
     model += NernstPlanck("V2", ions[0], electrolyte, potential, region, velocity=u)
@@ -133,7 +128,6 @@ def flux(c, z, D, advected=False):
         "F": electrolyte.F.default_si,
         "R": electrolyte.R.default_si,
         "T": electrolyte.T.default_si,
-        "floor": 0.1,
         "eps": 0.8,
         "D2": 2.4e-10,
         "D3": 2.4e-10,
@@ -198,8 +192,8 @@ def test_kinetics_of_a_redox_couple(cell, mass_transfer, box, compare):
     potentials = [Field(name, domain, ("Lagrange", 1)) for name in ("phi_s", "phi_l")]
     fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c2", "c3")]
     ions = [
-        Ion("V2", 2, 2.4e-10, fields[0], positive=True),
-        Ion("V3", 3, 2.4e-10, fields[1], positive=True),
+        Ion("V2", 2, 2.4e-10, fields[0]),
+        Ion("V3", 3, 2.4e-10, fields[1]),
     ]
     electrolyte = Electrolyte(
         "e", domain, ions, density=1.3e3, viscosity=5e-3, concentration_floor=0.1
@@ -286,8 +280,8 @@ def test_positive_couple_with_protons_water_and_piston_velocities(cell, box, com
     potentials = [Field(name, domain, ("Lagrange", 1)) for name in ("phi_s", "phi_l")]
     fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c4", "c5", "ch", "cw")]
     ions = [
-        Ion("VO", 2, 3.9e-10, fields[0], positive=True),
-        Ion("VO2", 1, 3.9e-10, fields[1], positive=True),
+        Ion("VO", 2, 3.9e-10, fields[0]),
+        Ion("VO2", 1, 3.9e-10, fields[1]),
         Ion("H", 1, 9.3e-9, fields[2]),
         Ion("H2O", 0, 2.3e-9, fields[3]),
         Ion("SO4", -2, 1.1e-9),

tests/electrochemistry/test_electrolyte.py

@@ -3,6 +3,7 @@
 
 import pytest
 import sympy.physics.units as syu
+import ufl
 
 from sumfem.discrete import DomainMaps, discretize
 from sumfem.electrochemistry import (
@@ -31,8 +32,8 @@ def vanadium(domain):
         for name in ("c2", "c3", "h")
     )
     return [
-        Ion("V2", 2, 2.4e-10, c2, positive=True),
-        Ion("V3", 3, 2.4e-10, c3, positive=True),
+        Ion("V2", 2, 2.4e-10, c2),
+        Ion("V3", 3, 2.4e-10, c3),
         Ion("H", 1, 9.3e-9, h),
         Ion("SO4", -2, 1.065e-9),
     ]
@@ -56,7 +57,7 @@ def uniform(domain, square, values):
     ions = vanadium(domain)
     electrolyte = Electrolyte("e", domain, ions)
     _, phi_l = potentials(domain)
-    # The conduction holds the parameters of every ion, and the floor
+    # The conduction holds the parameters of every ion
     model = ElectrolyteConduction("charge", phi_l, electrolyte)
     for ion in ions[:3]:
         model += NernstPlanck(ion.name, ion, electrolyte, phi_l)
@@ -89,15 +90,64 @@ def test_the_balanced_ion_follows_from_electroneutrality(domain, square):
     assert sulfate == pytest.approx(-charge / electrolyte["SO4"].charge)
 
 
-def test_positive_ions_take_the_positive_part(domain, square):
+def test_concentrations_keep_their_sign(domain, square):
+    """The electrolyte gives a concentration that undershoots zero as it is,
+    and the balanced one from it."""
     electrolyte, disc = uniform(domain, square, {**VALUES, "V2": -1.0})
-    assert mean(disc, domain, electrolyte.concentration("V2")) == pytest.approx(
-        0.0, abs=1e-12
+    assert mean(disc, domain, electrolyte.concentration("V2")) == pytest.approx(-1.0)
+    sulfate = (2 * -1.0 + 3 * VALUES["V3"] + VALUES["H"]) / 2
+    assert mean(disc, domain, electrolyte.concentration("SO4")) == pytest.approx(
+        sulfate
     )
-    unclipped = electrolyte.concentration("V2", positive=False)
-    assert mean(disc, domain, unclipped) == pytest.approx(-1.0)
-    # An ion that is not positive keeps its sign unless asked
-    assert mean(disc, domain, electrolyte.concentration("H")) == pytest.approx(3.0)
+
+
+def test_the_current_is_the_charge_the_ions_carry(domain, square):
+    """The electrolyte current equals ``F sum_j z_j N_j``, with the
+    Nernst-Planck fluxes of every ion, sulfate's from electroneutrality, in
+    a porous region with advection, also where V2+ undershoots zero.
+
+    A regression test: the conductivity once took the positive part of the
+    vanadium concentrations while migration did not, so conduction did not
+    carry the charge the ions moved where a concentration was negative.
+    """
+    ions = vanadium(domain)
+    electrolyte = Electrolyte("e", domain, ions)
+    _, phi_l = potentials(domain)
+    region = Region(porosity=0.6)
+    velocity = ufl.as_vector((1e-3, 2e-3))
+    transport = [
+        NernstPlanck(ion.name, ion, electrolyte, phi_l, region, velocity)
+        for ion in ions[:3]
+    ]
+    model = ElectrolyteConduction("charge", phi_l, electrolyte, region)
+    for summand in transport:
+        model += summand
+    disc = discretize([model], DomainMaps({domain: square()}))
+    profiles = {
+        ions[0].field: lambda x: -1.0 + 3.0 * x[0],
+        ions[1].field: lambda x: 2.0 + x[1],
+        ions[2].field: lambda x: 3.0 + x[0] * x[1],
+        phi_l: lambda x: x[0] ** 2 + 0.5 * x[1],
+    }
+    for field, profile in profiles.items():
+        disc.function(field).interpolate(profile)
+
+    e, sulfate = electrolyte, ions[3]
+    c = e.concentration(sulfate.name)
+    D = region.transport * e.diffusivity(sulfate.name)
+    migration = D * sulfate.charge * e.F * c / (e.R * e.T) * ufl.grad(phi_l)
+    fluxes = [m.flux() for m in transport]
+    fluxes.append(-D * ufl.grad(c) - migration + velocity * c)
+    charges = [e.F * ion.charge * N for ion, N in zip(ions, fluxes)]
+    current = e.current_density(phi_l, region.transport)
+
+    # The charge fluxes of the ions largely cancel, so the round-off is
+    # relative to them, not to the current
+    difference = current - sum(charges)
+    error = mean(disc, domain, ufl.inner(difference, difference))
+    size = sum(mean(disc, domain, ufl.inner(q, q)) for q in charges)
+    assert mean(disc, domain, ufl.inner(current, current)) > 0
+    assert error <= 1e-24 * size
 
 
 @pytest.mark.parametrize(