From 2436806f613627e6313024c88df8fe531922244e Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 21:39:07 +0000 Subject: [PATCH 1/6] Add intial half cell model. --- docs/api/electrochemistry.md | 7 + docs/api/symbolic.md | 12 + .../electrochemistry/vanadium_half_cell_1d.py | 279 +++++++++ src/sumfem/electrochemistry/__init__.py | 6 +- src/sumfem/electrochemistry/electrolyte.py | 8 +- src/sumfem/electrochemistry/half_cell.py | 534 ++++++++++++++++++ src/sumfem/electrochemistry/porous.py | 18 +- src/sumfem/models/porous.py | 20 + src/sumfem/symbolic/terminals.py | 17 +- tests/electrochemistry/half_cell/conftest.py | 434 ++++++++++++++ .../electrochemistry/half_cell/test_models.py | 203 +++++++ .../half_cell/test_reference.py | 353 ++++++++++++ tests/symbolic/test_terminals.py | 13 + tests/test_api_docs.py | 21 +- 14 files changed, 1907 insertions(+), 18 deletions(-) create mode 100644 examples/electrochemistry/vanadium_half_cell_1d.py create mode 100644 src/sumfem/electrochemistry/half_cell.py create mode 100644 tests/electrochemistry/half_cell/conftest.py create mode 100644 tests/electrochemistry/half_cell/test_models.py create mode 100644 tests/electrochemistry/half_cell/test_reference.py diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 0078203..5021441 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -49,3 +49,10 @@ .. autoclass:: MassTransferButlerVolmer :members: ``` + +## Half-cells + +```{eval-rst} +.. autoclass:: RedoxHalfCell + :members: +``` diff --git a/docs/api/symbolic.md b/docs/api/symbolic.md index 3008296..d2cb516 100644 --- a/docs/api/symbolic.md +++ b/docs/api/symbolic.md @@ -121,7 +121,11 @@ .. autoclass:: Interpolation :members: +.. currentmodule:: sumfem.symbolic.transfer + .. autodata:: Transferable + +.. currentmodule:: sumfem.symbolic ``` ## Units and dimensional analysis @@ -140,12 +144,20 @@ .. autoexception:: UnitError +.. currentmodule:: sumfem.symbolic.constants + .. autodata:: BOLTZMANN +.. currentmodule:: sumfem.symbolic.constants + .. autodata:: FARADAY +.. currentmodule:: sumfem.symbolic.constants + .. autodata:: GAS_CONSTANT +.. currentmodule:: sumfem.symbolic + .. autofunction:: boltzmann_constant ``` diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py new file mode 100644 index 0000000..401dff3 --- /dev/null +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -0,0 +1,279 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.13.6 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# # A vanadium half-cell in one dimension +# +# This demo charges, rests and discharges the negative half-cell of a +# vanadium redox flow battery across its porous electrode, as the 1D example +# of [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS). It shows how to: +# +# - Build a half-cell of any redox couple and electrolyte, negative or +# positive, with +# {py:class}`RedoxHalfCell ` +# - Run it through a galvanostatic schedule of charge, rest and discharge, +# with a {py:class}`Simulation ` +# - Follow the voltage and the state of charge of the half-cell +# +# ## The model +# +# The porous electrode holds the solid, of potential $\phi_s$, and in its pores +# of porosity $\varepsilon$ the electrolyte, of potential $\phi_l$, with the +# ions V²⁺, V³⁺, H⁺ and HSO₄⁻ of concentrations $c_k$, and sulfate SO₄²⁻ by +# electroneutrality. At the fibres of specific area $a$, the reaction +# $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ carries the current +# density $i$ from the electrolyte to the solid. Charge is conserved in both +# phases, and each ion is stored, transported and produced: +# +# $$ +# \begin{aligned} +# -\nabla \cdot (\sigma_{\text{eff}} \nabla \phi_s) + a i &= 0, \\ +# \nabla \cdot i_l - a i &= 0, \\ +# \varepsilon \frac{\partial c_k}{\partial t} + \nabla \cdot N_k &= s_k, +# \\ +# \sum_k z_k c_k &= 0, +# \end{aligned} +# $$ (half-cell) +# +# with the effective conductivity of the solid $\sigma_{\text{eff}}$, the +# current density of the electrolyte $i_l$ and the Nernst-Planck fluxes +# $N_k$, whose diffusivities are corrected for the porosity. The reaction +# consumes V²⁺ and produces V³⁺, $s_{\mathrm{V^{2+}}} = -a i / F$ and +# $s_{\mathrm{V^{3+}}} = a i / F$, at the Butler-Volmer current density +# +# $$ +# \begin{aligned} +# i &= F k_0 \left( +# e^{\alpha_a F \eta / (R T)} c_{\mathrm{V^{2+}}} +# - e^{-\alpha_c F \eta / (R T)} c_{\mathrm{V^{3+}}} \right), \\ +# \eta &= \phi_s - \phi_l - E_0, +# \end{aligned} +# $$ (half-cell-kinetics) +# +# with the formal potential $E_0$. The current density $j$ of the cell, +# positive on discharge, enters the solid at the current collector and leaves +# the electrolyte at the membrane, where the electrolyte has the +# concentrations of the reservoir. The voltage of the half-cell is the mean +# solid potential at the current collector less the mean electrolyte +# potential at the membrane. The parameters are those of the RFBniCS +# examples. +# +# We import the relevant packages: + +# + tags=["hide-input"] +from mpi4py import MPI + +import dolfinx.mesh +import matplotlib.pyplot as plt +import numpy as np +import sympy.physics.units as syu + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell + +# - + +# ## The electrode +# +# Across the electrode, from the current collector at $x = 0$ to the +# membrane at $x = W$, the cell is an interval. Its ends are located with a +# tolerance relative to the width, see {ref}`units-length-scale`: + +# + +width, num_cells = 2.5e-3, 250 +dtype = dolfinx.default_real_type +interval = dolfinx.mesh.create_interval( + MPI.COMM_WORLD, num_cells, [0.0, width], dtype=dtype +) +tolerance = 100 * float(np.finfo(dtype).eps) * width +marked = discrete.MeshData.from_locators( + interval, + { + 0: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0, atol=tolerance), + "Membrane": lambda x: np.isclose(x[0], width, atol=tolerance), + } + }, + boundary=True, +) +domain = symbolic.Domain.from_mesh("cell", interval, length_scale=width) +# - + +# ## The half-cell +# +# 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 +# of the reservoir name the ions that are transported, so sulfate, which has +# none, follows from electroneutrality: + +# + +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("H", 1, 9.312e-9 * diffusivity), + Ion("HSO4", -1, 1.33e-9 * diffusivity), + Ion("SO4", -2, 1.065e-9 * diffusivity), +) +reservoir = { + "V2": 300.0 * concentration, + "V3": 300.0 * concentration, + "H": 1000.0 * concentration, + "HSO4": 1000.0 * concentration, +} +# - + +# The porous electrode is a carbon felt, a mat of carbon fibres, with its +# porosity, the conductivity of its fibres, their specific area and diameter, +# see {py:class}`PorousElectrode `. +# The constants of the permeability of its pores and of the inertial drag +# matter only with flow, and this demo has none: + +electrode = PorousElectrode( + "electrode", + domain, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=0.5 * syu.micrometer, + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) + +# {py:class}`RedoxHalfCell ` sums the +# models of {eq}`half-cell` for a couple of these ions, with its standard rate +# constant $k_0$ and formal potential $E_0$. Discharge oxidizes V²⁺ in the +# negative half-cell; the positive half-cell of the battery, with the couple +# VO²⁺/VO₂⁺, reduces its couple on discharge, which reverses the current. The +# concentrations of the reservoir are Dirichlet conditions at the membrane: + +cell = RedoxHalfCell( + domain, + ions, + reservoir, + reduced="V2", + oxidized="V3", + electrode=electrode, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, +) +inlet = cell.inlet_conditions("Membrane") + +# The discretization bounds the quadrature degree of the nonlinear kinetics, +# as RFBniCS does. The initial state is at rest: the concentrations of the +# reservoir, the solid at zero and the electrolyte at minus the formal +# potential: + +disc = discrete.discretize( + [cell.model, *inlet], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +cell.initialize(disc) + +# ## The schedule +# +# The cell is charged, rested and discharged, each with time steps that +# start short, while the cell responds fast, and then grow. Each operation +# gives its mode, the magnitude of the current density in A/m², and its time +# steps, as pairs of a number of steps and their size in seconds. The +# current of the cell is positive on discharge, so the mode gives its sign: + +# + +schedule = ( + ("charge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), + ("rest", 0.0, ((5, 0.1), (10, 1.0), (20, 2.0))), + ("discharge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), +) +sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} +# - + +# A {py:class}`Simulation ` takes the steps with +# backward Euler, checks the units, and solves every time step in +# dimensionless form, see {ref}`units-scaling`. The steps of the schedule +# change in size, so the reference time of the time derivatives is the +# largest of them: a shorter step only makes Newton's absolute tolerance +# stricter, see {doc}`../cooling_plate`: + +# + +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 100}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +# - + +# ## Charge, rest and discharge +# +# Each operation sets the current of the cell, and each part of it the size +# of the time step, before the simulation takes its steps. After every step, +# a callback records the time, the mode, the voltage of the half-cell and its +# state of charge, the share of V²⁺ in the vanadium of the electrode: + +# + +records: list[tuple[float, str, float, float]] = [] +for mode, current, steps in schedule: + disc.set({cell.current: sign[mode] * current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run( + num_steps=count, + callback=lambda t, step, mode=mode: records.append( + (t, mode, cell.voltage(disc), cell.state_of_charge(disc)) + ), + ) +times = np.array([r[0] for r in records]) +modes = np.array([r[1] for r in records]) +voltage = np.array([r[2] for r in records]) +charge = np.array([r[3] for r in records]) +# - + +# 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: + +# + +steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} +print(f"State of charge after charge: {charge[modes == 'charge'][-1]:.6f}") +print(f"State of charge at the end: {charge[-1]:.6f}") +assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0) +rest = charge[modes == "rest"] +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() +# - + +# The voltage and the state of charge through the schedule: + +fig, (top, bottom) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + top.plot(times[modes == mode], 1e3 * voltage[modes == mode], ".", label=mode) + bottom.plot(times[modes == mode], charge[modes == mode], ".", label=mode) +top.set_ylabel("Voltage (mV)") +bottom.set_ylabel("State of charge") +bottom.set_xlabel("Time (s)") +top.legend() +plt.show() diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 7b181b2..51853bb 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -1,4 +1,4 @@ -"""Models of electrochemistry, built on :mod:`sumfem.models`. +r"""Models of electrochemistry, built on :mod:`sumfem.models`. - Electrolytes: :class:`Electrolyte`, a dilute solution of :class:`Ion` and the constants and properties its models share. @@ -9,6 +9,8 @@ - Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer` and :class:`MassTransferButlerVolmer`, on :class:`ElectrodeReaction`, with the Butler-Volmer equation :func:`~sumfem.models.butler_volmer`. +- Half-cells of a redox flow battery: :class:`RedoxHalfCell`, negative or + positive, of any couple and electrolyte. The package imports from the core of :mod:`sumfem`, and never the other way around. @@ -16,6 +18,7 @@ from .charge import ElectrolyteConduction from .electrolyte import Electrolyte, Ion +from .half_cell import RedoxHalfCell from .kinetics import ( ButlerVolmer, ElectrodeReaction, @@ -32,6 +35,7 @@ "Region", "ElectrolyteConduction", "NernstPlanck", + "RedoxHalfCell", "ButlerVolmer", "ElectrodeReaction", "ExchangeCurrentButlerVolmer", diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index dc4f106..116ffd4 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -6,6 +6,7 @@ import types import typing +import sympy import sympy.physics.units as syu import ufl @@ -33,8 +34,9 @@ class Ion: Attributes: name: Name of the ion, e.g. ``"V2"``. charge: Charge number ``z``. - diffusivity: Diffusion coefficient in free electrolyte, in m^2/s. - It is the default value of the parameter the electrolyte creates. + diffusivity: Diffusion coefficient in free electrolyte, in m^2/s or a + quantity with units, e.g. ``1e-9 * syu.m**2 / syu.s``. It is the + 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 @@ -44,7 +46,7 @@ class Ion: name: str charge: float - diffusivity: float + diffusivity: float | sympy.Expr concentration: Field | None = None positive: bool = False diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py new file mode 100644 index 0000000..e15a1bd --- /dev/null +++ b/src/sumfem/electrochemistry/half_cell.py @@ -0,0 +1,534 @@ +"""Half-cells of a redox flow battery. + +:class:`RedoxHalfCell` is a half-cell of any redox couple of dissolved ions, +negative or positive: a porous electrode, optionally beside free channels, +filled with an electrolyte that flows from a reservoir, with charge +conservation in the solid and the electrolyte, transport of every ion, and +Butler-Volmer kinetics of the couple at the fibres. +""" + +from __future__ import annotations + +import dataclasses +import types +import typing + +import sympy +import sympy.physics.units as syu +import ufl + +from sumfem.discrete import Discretization +from sumfem.models import Accumulation, BoundaryFlux, Conduction, Outflow +from sumfem.symbolic import ( + DirichletBC, + Domain, + Field, + Parameter, + Residual, + ResidualSum, +) + +from .charge import ElectrolyteConduction +from .electrolyte import CONCENTRATION, Electrolyte, Ion +from .kinetics import ButlerVolmer, MassTransferButlerVolmer +from .porous import PorousElectrode, Region +from .transport import NernstPlanck + +__all__ = ["RedoxHalfCell"] + + +class RedoxHalfCell: + """A half-cell of a redox flow battery, of a couple of dissolved ions. + + The unknowns are the potentials ``phi_s`` of the solid and ``phi_l`` of + the electrolyte, in V, and the concentration ``c_`` of every ion + with a reservoir concentration, in mol/m^3. The one ion without follows + from electroneutrality. In the porous electrode, the couple + ``reduced <-> oxidized + n e-`` follows Butler-Volmer kinetics, limited + by mass transfer to the fibres when there is flow. In free channels the + electrolyte is not corrected for porosity, and there is no reaction or + solid phase. + + The model, :attr:`model`, is the sum of: + + - ``"solid"``: Ohmic conduction in the solid, + :class:`~sumfem.models.Conduction` with Bruggeman's effective + conductivity of the electrode; + - ``"reaction"``: the kinetics of the couple, + :class:`~sumfem.electrochemistry.ButlerVolmer`, or + :class:`~sumfem.electrochemistry.MassTransferButlerVolmer` with mass + transfer, which moves the current between the phases and consumes and + produces the ions of the couple; + - ``"collector_current"`` and ``"membrane_current"``: the current of the + cell into the solid and out of the electrolyte, + :class:`~sumfem.models.BoundaryFlux`; + - for the electrode's pores and, with a suffix ``"_channels"``, the free + channels: charge conservation in the electrolyte, + :class:`~sumfem.electrochemistry.ElectrolyteConduction`, and for every + transported ion its storage, :class:`~sumfem.models.Accumulation` + weighted with the porosity, and its transport, + :class:`~sumfem.electrochemistry.NernstPlanck`, named + ``"storage_"`` and ``"transport_"``; + - with an outlet, ``"outflow_"``: advection out of the cell, + :class:`~sumfem.models.Outflow`; + - with a proton flux, ``"proton_flux"``: the membrane current carried by + protons, a :class:`~sumfem.models.BoundaryFlux` of the ion ``"H"``. + + The cell is galvanostatic: the parameter ``".current"`` is the + current density through the membrane, in A/m^2, positive on discharge. + On discharge the couple is oxidized in the negative half-cell and + reduced in the positive one, so the current enters the solid at the + current collector of the first, and leaves it at that of the second. + The current collector carries the current density ``current * + area_ratio``, so that the total currents match when the two areas + differ. + + Args: + domain: The domain of the cell. + ions: The ions of the electrolyte, without fields: the cell creates + the field of every ion that has a reservoir concentration. + reservoir: The concentration of each transported ion in the + reservoir, in mol/m^3 or a quantity with units, the default of + the parameter ``".c0_"``. + reduced: Name of the reduced ion of the couple. + oxidized: Name of the oxidized ion of the couple. + electrode: The porous electrode. + oxidation_on_discharge: Whether discharge oxidizes the couple, as in + the negative half-cell. + rate_constant: Standard rate constant of the couple in m/s or with + units, the default of the parameter ``".rate_constant"``. + formal_potential: Formal potential of the couple in V or with + units, the default of the parameter ``".formal_potential"``. + anodic_transfer: Anodic transfer coefficient, the default of the + parameter ``".anodic_transfer"``. + cathodic_transfer: Cathodic transfer coefficient, the default of the + parameter ``".cathodic_transfer"``. + electrons: Number of electrons the couple transfers. + channels: Cell markers of free channels. Empty for none. + velocity: Superficial velocity of the electrolyte, or ``None`` + without flow. + mass_transfer: Limit the kinetics by mass transfer to the fibres, + which needs a velocity. Defaults to whether there is one. + collector: Facet markers of the current collector. + membrane: Facet markers of the membrane. + outlet: Facet markers of an outlet, through which ions leave by + advection. Empty for none. + proton_flux: Carry the membrane current by protons, the ion ``"H"``, + through the membrane. + density: Density of the electrolyte, see + :class:`~sumfem.electrochemistry.Electrolyte`. Pass the + parameter of a flow model to share it. + viscosity: Dynamic viscosity of the electrolyte, given like the + density. + degree: Degree of the Lagrange elements of the unknowns. + name: Prefix of the cell's own parameters. + + Raises: + ValueError: For mass transfer or an outlet without a velocity, or a + reservoir concentration of an ion that the electrolyte does not + have. + """ + + _domain: Domain + """The domain of the cell.""" + + _electrode: PorousElectrode + """The porous electrode.""" + + _regions: tuple[tuple[str, Region], ...] + """The regions filled with electrolyte: the electrode's pores, and the + free channels if any. Each comes with the suffix of its models' names.""" + + _electrolyte: Electrolyte + """The electrolyte.""" + + _reduced: str + """Name of the reduced ion of the couple.""" + + _oxidized: str + """Name of the oxidized ion of the couple.""" + + _phi_s: Field + """Potential of the solid.""" + + _phi_l: Field + """Potential of the electrolyte.""" + + _concentrations: types.MappingProxyType[str, Field] + """Concentration of each transported ion, by ion name.""" + + _inlet: types.MappingProxyType[str, Parameter] + """Reservoir concentration of each transported ion, by ion name.""" + + _current: Parameter + """Current density through the membrane, positive on discharge.""" + + _area_ratio: Parameter + """Ratio of the membrane area to the current collector area.""" + + _formal_potential: Parameter + """Formal potential of the couple.""" + + _reaction: ButlerVolmer + """The reaction in the electrode.""" + + _model: ResidualSum + """The sum of all models of the cell.""" + + _collector: tuple[str, ...] + """Facet markers of the current collector.""" + + _membrane: tuple[str, ...] + """Facet markers of the membrane.""" + + _functionals: types.MappingProxyType[str, ufl.Form] + """The integrals behind the cell's diagnostics, created once so that a + discretization compiles each of them only once.""" + + def __init__( + self, + domain: Domain, + ions: typing.Sequence[Ion], + reservoir: typing.Mapping[str, float | sympy.Expr], + reduced: str, + oxidized: str, + electrode: PorousElectrode, + oxidation_on_discharge: bool, + rate_constant: float | sympy.Expr, + formal_potential: float | sympy.Expr, + anodic_transfer: float = 0.5, + cathodic_transfer: float = 0.5, + electrons: int = 1, + channels: tuple[str, ...] = (), + velocity: ufl.core.expr.Expr | None = None, + mass_transfer: bool | None = None, + collector: tuple[str, ...] = ("CurrentCollector",), + membrane: tuple[str, ...] = ("Membrane",), + outlet: tuple[str, ...] = (), + proton_flux: bool = False, + density: ufl.core.expr.Expr | float = 1399.8, + viscosity: ufl.core.expr.Expr | float = 6.963e-3, + degree: int = 1, + name: str = "cell", + ): + # Mass transfer to the fibres depends on the speed of the flow, so it + # is on by default whenever there is one + mass_transfer = velocity is not None if mass_transfer is None else mass_transfer + if velocity is None and (mass_transfer or outlet): + raise ValueError("Mass transfer and an outlet need a velocity.") + unknown = sorted(set(reservoir) - {ion.name for ion in ions}) + if unknown: + raise ValueError(f"The electrolyte has no ions {unknown}.") + self._domain = domain + self._electrode = electrode + self._reduced, self._oxidized = reduced, oxidized + # The regions filled with electrolyte, with the suffix of their models: + # the pores of the electrode, whose porosity weights storage and + # corrects transport, and free channels, of porosity one + regions = [("", self._electrode.region)] + if channels: + regions.append(("_channels", Region(tuple(channels)))) + self._regions = tuple(regions) + + # The unknowns. A concentration is of the order of a mole per litre, + # its scale for the dimensional analysis and the scaled equations + element = ("Lagrange", degree) + self._phi_s = Field("phi_s", domain, element, unit=syu.volt) + self._phi_l = Field("phi_l", domain, element, unit=syu.volt) + self._concentrations = types.MappingProxyType( + { + ion.name: Field( + f"c_{ion.name}", domain, element, unit=CONCENTRATION, scale=1e3 + ) + for ion in ions + if ion.name in reservoir + } + ) + # The electrolyte, with the field of every transported ion. The ion + # without a reservoir concentration, and so without a field, follows + # from electroneutrality. Its density and viscosity enter the mass + # transfer, and may be shared with a model of the flow + self._electrolyte = Electrolyte( + "electrolyte", + domain, + [ + dataclasses.replace( + ion, concentration=self._concentrations.get(ion.name) + ) + for ion in ions + ], + density=density, + viscosity=viscosity, + ) + # The concentrations of the reservoir: the initial state, and the + # values of the conditions at an inlet, see inlet_conditions + self._inlet = types.MappingProxyType( + { + ion: Parameter( + f"{name}.c0_{ion}", + domain, + default=value, + unit=CONCENTRATION, + doc=f"Reservoir concentration of {ion}", + ) + for ion, value in reservoir.items() + } + ) + # The current density of the cell through the membrane, and the + # ratio of the areas that makes the current through the collector + # equal, see balance_currents + self._current = Parameter( + f"{name}.current", + domain, + default=0.0, + unit=syu.ampere / syu.m**2, + doc="Current density through the membrane, positive on discharge", + ) + self._area_ratio = Parameter( + f"{name}.area_ratio", + domain, + default=1.0, + doc="Membrane area over current collector area", + ) + # The kinetics of the couple, parameters of the cell that the reaction + # shares, and the specific area of the electrode + self._formal_potential = Parameter( + f"{name}.formal_potential", + domain, + default=formal_potential, + unit=syu.volt, + doc=f"Formal potential of {reduced}/{oxidized}", + ) + kinetics: dict[str, typing.Any] = dict( + specific_area=self._electrode.specific_area, + formal_potential=self._formal_potential, + rate_constant=Parameter( + f"{name}.rate_constant", + domain, + default=rate_constant, + unit=syu.m / syu.second, + doc=f"Standard rate constant of {reduced}/{oxidized}", + ), + anodic_transfer=Parameter( + f"{name}.anodic_transfer", + domain, + default=anodic_transfer, + doc=f"Anodic transfer coefficient of {reduced}/{oxidized}", + ), + cathodic_transfer=Parameter( + f"{name}.cathodic_transfer", + domain, + default=cathodic_transfer, + doc=f"Cathodic transfer coefficient of {reduced}/{oxidized}", + ), + electrons=electrons, + markers=self._electrode.markers, + ) + # The reaction in the electrode: Butler-Volmer at the concentrations + # of the pores, or at the surface of the fibres with mass transfer + e = self._electrolyte + if mass_transfer: + assert velocity is not None + self._reaction = MassTransferButlerVolmer( + "reaction", + self._phi_s, + self._phi_l, + e, + reduced, + oxidized, + velocity, + fiber_diameter=self._electrode.fiber_diameter, + **kinetics, + ) + else: + self._reaction = ButlerVolmer( + "reaction", self._phi_s, self._phi_l, e, reduced, oxidized, **kinetics + ) + + # The current into the solid at the collector, and out of the + # electrolyte at the membrane, on discharge. Oxidation on discharge + # releases electrons into the solid, which leave through the + # collector, so the current of the cell enters there; reduction on + # discharge takes them, so the current leaves the solid there + inflow = self._current if oxidation_on_discharge else -self._current + models: list[Residual] = [ + Conduction( + "solid", + self._phi_s, + self._electrode.effective_conductivity(), + markers=self._electrode.markers, + ), + self._reaction, + BoundaryFlux( + "collector_current", self._phi_s, inflow * self._area_ratio, collector + ), + BoundaryFlux("membrane_current", self._phi_l, -inflow, membrane), + ] + # Every region filled with electrolyte conserves charge, and stores + # and transports every ion, with the porosity of the region + for suffix, region in self._regions: + models.append( + ElectrolyteConduction(f"electrolyte{suffix}", self._phi_l, e, region) + ) + for ion, c in self._concentrations.items(): + models.append( + Accumulation( + f"storage_{ion}{suffix}", + c, + weight=region.porosity, + markers=region.markers, + ) + ) + models.append( + NernstPlanck( + f"transport_{ion}{suffix}", + e[ion], + e, + self._phi_l, + region, + velocity, + ) + ) + # The ions that the flow carries out of the cell + if outlet: + assert velocity is not None + for ion, c in self._concentrations.items(): + models.append(Outflow(f"outflow_{ion}", c, velocity, outlet)) + # The membrane passes the current as protons: their flux is the + # current over their charge and Faraday's constant + if proton_flux: + models.append( + BoundaryFlux( + "proton_flux", + self._concentrations["H"], + -inflow / (e["H"].charge * e.F), + membrane, + ) + ) + self._model = ResidualSum(models, name=name) + + # The integrals of the diagnostics: the areas and mean potentials of + # the collector and the membrane, for the voltage, and the amounts of + # the couple in the electrolyte of every region, for the state of + # charge + self._collector = tuple(collector) + self._membrane = tuple(membrane) + ds_collector = domain.ds(*self._collector) + ds_membrane = domain.ds(*self._membrane) + functionals = { + "collector_area": 1.0 * ds_collector, + "membrane_area": 1.0 * ds_membrane, + "collector_potential": self._phi_s * ds_collector, + "membrane_potential": self._phi_l * ds_membrane, + } + for role, ion in (("reduced", reduced), ("oxidized", oxidized)): + c = self._concentrations[ion] + functionals[f"amount_{role}"] = sum( + region.storage * c * domain.dx(*region.markers) + for _, region in self._regions + ) + self._functionals = types.MappingProxyType(functionals) + + @property + def model(self) -> Residual: + """The sum of all models of the cell.""" + return self._model + + @property + def electrode(self) -> PorousElectrode: + """The porous electrode.""" + return self._electrode + + @property + def electrolyte(self) -> Electrolyte: + """The electrolyte.""" + return self._electrolyte + + @property + def reaction(self) -> ButlerVolmer: + """The reaction in the electrode.""" + return self._reaction + + @property + def phi_s(self) -> Field: + """Potential of the solid.""" + return self._phi_s + + @property + def phi_l(self) -> Field: + """Potential of the electrolyte.""" + return self._phi_l + + @property + def concentrations(self) -> typing.Mapping[str, Field]: + """Concentration of each transported ion, by ion name.""" + return self._concentrations + + @property + def inlet(self) -> typing.Mapping[str, Parameter]: + """Reservoir concentration of each transported ion, by ion name.""" + return self._inlet + + @property + def current(self) -> Parameter: + """Current density through the membrane, positive on discharge.""" + return self._current + + @property + def area_ratio(self) -> Parameter: + """Ratio of the membrane area to the current collector area.""" + return self._area_ratio + + def inlet_conditions(self, markers: str | tuple[str, ...]) -> list[DirichletBC]: + """The reservoir concentration of every transported ion on some facets.""" + return [ + DirichletBC(c, markers, self._inlet[ion]) + for ion, c in self._concentrations.items() + ] + + def initialize(self, discretization: Discretization) -> None: + """Set a state at rest, as RFBniCS starts. + + The concentrations start at their reservoir values, the solid at zero + and the electrolyte at minus the formal potential. + """ + values: dict[str | ufl.Coefficient, typing.Any] = { + c: discretization.value(self._inlet[ion]) + for ion, c in self._concentrations.items() + } + values[self._phi_s] = 0.0 + values[self._phi_l] = -discretization.value(self._formal_potential) + discretization.set(values) + + def balance_currents(self, discretization: Discretization) -> None: + """Set the area ratio from the areas of membrane and current collector.""" + assemble = discretization.assemble_scalar + ratio = assemble(self._functionals["membrane_area"]) / assemble( + self._functionals["collector_area"] + ) + discretization.set({self._area_ratio: ratio}) + + def voltage(self, discretization: Discretization) -> float: + """The half-cell voltage in V. + + That is the mean solid potential on the current collector minus the + mean electrolyte potential on the membrane. + """ + assemble = discretization.assemble_scalar + f = self._functionals + return assemble(f["collector_potential"]) / assemble( + f["collector_area"] + ) - assemble(f["membrane_potential"]) / assemble(f["membrane_area"]) + + def state_of_charge(self, discretization: Discretization) -> float: + """The share of the reduced ion in the couple, over the cell.""" + assemble = discretization.assemble_scalar + reduced = assemble(self._functionals["amount_reduced"]) + oxidized = assemble(self._functionals["amount_oxidized"]) + return reduced / max(reduced + oxidized, 1e-30) + + 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._reduced] + c_ox = self._concentrations[self._oxidized] + return c_red / (c_red + c_ox + self._electrolyte.floor) diff --git a/src/sumfem/electrochemistry/porous.py b/src/sumfem/electrochemistry/porous.py index 2fe4b0d..3386866 100644 --- a/src/sumfem/electrochemistry/porous.py +++ b/src/sumfem/electrochemistry/porous.py @@ -4,6 +4,7 @@ import dataclasses +import sympy import sympy.physics.units as syu import ufl @@ -53,7 +54,8 @@ class PorousElectrode: name: Prefix of the parameter names. domain: The domain of the parameters. markers: Cell markers of the electrode. Empty for the whole domain. - porosity: Default porosity, if any. + porosity: Default porosity, if any. Every default with a unit is in + SI units, or a quantity with units, e.g. ``0.5 * syu.micrometer``. conductivity: Default conductivity of the fiber material in S/m, if any. specific_area: Default specific area in 1/m, if any. @@ -90,11 +92,11 @@ def __init__( domain: Domain, markers: tuple[str, ...] = (), porosity: float | None = None, - conductivity: float | None = None, - specific_area: float | None = None, - fiber_diameter: float | None = None, + conductivity: float | sympy.Expr | None = None, + specific_area: float | sympy.Expr | None = None, + fiber_diameter: float | sympy.Expr | None = None, kozeny_carman_constant: float | None = None, - pore_length: float | None = None, + pore_length: float | sympy.Expr | None = None, ): self._markers = tuple(markers) self._porosity = Parameter( @@ -165,11 +167,13 @@ def effective_conductivity(self) -> ufl.core.expr.Expr: return bruggeman(self._conductivity, 1.0 - self._porosity) def permeability(self) -> ufl.core.expr.Expr: - """Kozeny-Carman permeability of the fibers.""" + """Kozeny-Carman permeability of the fibers, see + :func:`~sumfem.models.kozeny_carman`.""" return kozeny_carman( self._fiber_diameter, self._porosity, self._kozeny_carman_constant ) def forchheimer(self) -> ufl.core.expr.Expr: - """Ergun's inertial drag coefficient.""" + """Ergun's inertial drag coefficient, an estimate for fibres, see + :func:`~sumfem.models.ergun_coefficient`.""" return ergun_coefficient(self._porosity, self._pore_length) diff --git a/src/sumfem/models/porous.py b/src/sumfem/models/porous.py index 7435b8f..510a2b7 100644 --- a/src/sumfem/models/porous.py +++ b/src/sumfem/models/porous.py @@ -19,6 +19,11 @@ def bruggeman( ) -> ufl.core.expr.Expr: """Bruggeman's effective property ``value * fraction**exponent``. + It corrects a transport property of a phase for the volume the other + phase takes and for the tortuosity of the paths through it. The exponent + ``3/2`` is that of spherical obstructions, see `Tjaden et al., 2016 + `_. + Args: value: The property of the phase, e.g. a diffusivity or the conductivity of the solid. @@ -36,6 +41,10 @@ def kozeny_carman( ) -> ufl.core.expr.Expr: """Kozeny-Carman permeability ``d^2 eps^3 / (K (1 - eps)^2)`` of a fibrous medium. + An estimate of the permeability of a fibrous porous electrode, such as a + carbon felt, from the diameter of its fibres, with an empirical constant + ``K``. + Args: fiber_diameter: Fiber diameter ``d``. porosity: Porosity ``eps``. @@ -70,6 +79,17 @@ def ergun_coefficient( ) -> ufl.core.expr.Expr: """Ergun's inertial drag coefficient ``1.75 (1 - eps) / (L eps^3)``. + The coefficient ``beta``, in 1/m, of Forchheimer's drag ``rho beta |u| u`` + on the superficial velocity ``u``, from Ergun's relation for randomly + packed beds of spheres of diameter ``L`` (`Ergun and Orning, 1949 + `_). For a fibrous electrode, such + as a carbon felt, it is an estimate only: the model of `Tomadakis and + Robertson, 2005 `_ is made for + random fibre structures, and fits to experiments give the coefficient of + a given electrode, e.g. of the order of Ergun's for the 3D printed + electrodes of `van der Heijden et al., 2023 + `_. + Args: porosity: Porosity ``eps``. length: Characteristic pore length ``L``. diff --git a/src/sumfem/symbolic/terminals.py b/src/sumfem/symbolic/terminals.py index 6d04ac7..2e2c1ae 100644 --- a/src/sumfem/symbolic/terminals.py +++ b/src/sumfem/symbolic/terminals.py @@ -39,7 +39,7 @@ from ..error_handling import ComponentError, check_name, qualified_name from .domains import Domain, domain_of from .elements import ElementSpec, as_element -from .si import Unit, to_si +from .si import Unit, has_units, si_factor, to_si _PIOLA = ( ufl.pullback.contravariant_piola, @@ -442,14 +442,20 @@ class Parameter(ufl.Coefficient, QuantityMixin): named ``"."``. domain: The domain the parameter is defined on. shape: Value shape, ``()`` for a scalar. - default: Default value in ``unit``, used when no other value is set. - A scalar sets every component of a vector or tensor parameter. + default: Default value in ``unit``, used when no other value is set, + or a sympy quantity carrying its own units, e.g. ``0.5 * + syu.millimeter``, converted to ``unit``. A scalar sets every + component of a vector or tensor parameter. unit: Unit of the value, a sympy unit expression. ``None`` means dimensionless. The discretization stores values in SI base units. doc: Description, shown in documentation and archived metadata. The name must be passed positionally: ``ufl.Coefficient.__new__`` reads the first positional argument. + + Raises: + ValueError: If a default with units does not have the dimension of + ``unit``. """ _name: str @@ -467,7 +473,7 @@ def __init__( /, domain: Domain, shape: tuple[int, ...] = (), - default: npt.ArrayLike | None = None, + default: npt.ArrayLike | sympy.Expr | None = None, unit: Unit | None = None, doc: str = "", ): @@ -484,6 +490,9 @@ def __init__( space = ufl.FunctionSpace(domain, element, label=f"{name}#parameter") ufl.Coefficient.__init__(self, space) self._name = name + if has_units(default): + # A quantity with its own units, as a value in the declared unit + default = to_si(default, unit) / si_factor(unit) self._default = ( None if default is None else _validate_value(name, default, shape) ) diff --git a/tests/electrochemistry/half_cell/conftest.py b/tests/electrochemistry/half_cell/conftest.py new file mode 100644 index 0000000..65fc6ea --- /dev/null +++ b/tests/electrochemistry/half_cell/conftest.py @@ -0,0 +1,434 @@ +"""Meshes and the original RFBniCS forms, for comparison with the models. + +The reference residuals are transcribed from the examples of RFBniCS +(https://github.com/Ah-Aghabarari/RFBniCS, MIT licence), with one change: +each unknown is its own function, where RFBniCS splits one blocked +function. The tests pass the functions of a discretization, so both sides +assemble on the same spaces. + +The tests use pytest's importlib import mode, so these are handed to tests +through the ``reference`` fixture instead of being imported. +""" + +import types + +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.mesh +import numpy as np +import pytest +import sympy.physics.units as syu +import ufl + +from sumfem.discrete import MeshData +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.symbolic import Domain + +W = 2.5e-3 +"""Width of the electrode in the examples, in m.""" + +H = 1e-2 +"""Height of the electrode in the 2D example, in m.""" + +PARAMETERS = { + "sigma": 3402.1, + "a": 2.3e4, + "alpha_a": 0.5, + "alpha_c": 0.5, + "k0": 2.5e-7, + "E_formal": -0.36, + "df": 0.5e-6, + "F": 96485.33212, + "R": 8.314462618, + "T": 298.15, + "D_V2": 8.10e-12, + "D_V3": 1.65e-11, + "D_H": 9.312e-9, + "D_HSO4": 1.33e-9, + "D_SO4": 1.065e-9, + "c0_V2": 300.0, + "c0_V3": 300.0, + "c0_H": 1000.0, + "c0_HSO4": 1000.0, + "z_V2": 2.0, + "z_V3": 3.0, + "z_H": 1.0, + "z_HSO4": -1.0, + "z_SO4": -2.0, + "mu": 6.963e-3, + "rho": 1399.8, + "Kc": 5.55, + "epsilon": 0.94, + "u_in": 1.0e-7, + "p_out": 3.0e5, + "Lp": 100.0e-6, +} +"""Default parameters of the RFBniCS examples.""" + +QUADRATURE_DEGREE = 4 +"""The quadrature degree of the examples for degree-one elements.""" + + +def interval(num_cells: int = 20) -> MeshData: + """The 1D example's electrode, with its current collector and membrane.""" + mesh = dolfinx.mesh.create_interval(MPI.COMM_WORLD, num_cells, [0.0, W]) + return MeshData.from_locators( + mesh, + { + 0: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0), + "Membrane": lambda x: np.isclose(x[0], W), + } + }, + boundary=True, + ) + + +def rectangle(nx: int = 5, ny: int = 8) -> MeshData: + """The 2D example's electrode, with flow from the inlet at the bottom.""" + mesh = dolfinx.mesh.create_rectangle( + MPI.COMM_WORLD, + np.array([[0.0, 0.0], [W, H]]), + [nx, ny], + cell_type=dolfinx.mesh.CellType.quadrilateral, + ) + return MeshData.from_locators( + mesh, + { + 1: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0), + "Membrane": lambda x: np.isclose(x[0], W), + "Inlet": lambda x: np.isclose(x[1], 0.0), + "Outlet": lambda x: np.isclose(x[1], H), + } + }, + boundary=True, + ) + + +def box(n: int = 2) -> MeshData: + """A unit cube with an electrode and two channels, as in the 3D example. + + The electrode fills ``x > 1/2``. The inlet and outlet channels fill the + rest, below and above ``y = 1/2``, and open at ``x = 0``. + """ + mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 2 * n, 2 * n, 2 * n) + half = 0.5 + 1e-10 + cells = MeshData.from_locators( + mesh, + { + 3: { + "InletChannel": lambda x: (x[0] <= half) & (x[1] <= half), + "OutletChannel": lambda x: (x[0] <= half) & (x[1] >= 0.5 - 1e-10), + "Electrode": lambda x: x[0] >= 0.5 - 1e-10, + } + }, + ) + sides = { + "Inlet": lambda x: np.isclose(x[0], 0.0) & (x[1] <= half), + "Outlet": lambda x: np.isclose(x[0], 0.0) & (x[1] >= 0.5 - 1e-10), + "CurrentCollector": lambda x: np.isclose(x[2], 0.0) & (x[0] >= 0.5 - 1e-10), + "Membrane": lambda x: np.isclose(x[0], 1.0), + } + + def walls(x): + """The rest of the boundary: markers of facets do not overlap.""" + return ~np.logical_or.reduce([side(x) for side in sides.values()]) + + facets = MeshData.from_locators(mesh, {2: {"Walls": walls, **sides}}, boundary=True) + return MeshData( + mesh, {**cells.tags, **facets.tags}, {**cells.markers, **facets.markers} + ) + + +def measures(marked: MeshData, degree: int | None = QUADRATURE_DEGREE): + """Cell and facet measures of a marked mesh, and a lookup of marker ids.""" + mesh = marked.mesh + tdim = mesh.topology.dim + metadata = None if degree is None else {"quadrature_degree": degree} + dx = ufl.Measure( + "dx", domain=mesh, subdomain_data=marked.tags.get(tdim), metadata=metadata + ) + ds = ufl.Measure( + "ds", domain=mesh, subdomain_data=marked.tags.get(tdim - 1), metadata=metadata + ) + + def ids(dim, *names): + return tuple(marked.markers[dim][name][0] for name in names) + + return dx, ds, ids + + +def _constants(mesh, values): + return { + name: dolfinx.fem.Constant(mesh, np.float64(value)) + for name, value in values.items() + } + + +def electrochemistry( + marked: MeshData, + u, + u_old, + v, + dt, + current, + *, + velocity=None, + regions: bool = False, + area_ratio: float = 1.0, + membrane_flux: bool = False, + outlet: bool = False, +): + """The residual of the RFBniCS examples, in one form. + + Without a velocity it is the 1D example's, with a velocity the 2D + example's, and with regions the 3D example's, where the electrode and + channel cells differ. + + Args: + marked: The mesh. + u: Functions of phi_s, phi_l, c_V2, c_V3, c_H and c_HSO4. + u_old: Functions of the previous concentrations of V2, V3, H, HSO4. + v: Test functions of the six unknowns. + dt: Time step, a number or a constant. + current: Membrane current density, signed, a number or a constant. + velocity: The flow velocity, if any. + regions: Whether the mesh has electrode and channel cells. + area_ratio: Membrane area over current collector area. + membrane_flux: Whether protons carry the membrane current. + outlet: Whether ions leave through an outlet. + """ + mesh = marked.mesh + c = _constants(mesh, PARAMETERS) + tdim = mesh.topology.dim + dx, ds, ids = measures(marked) + if regions: + dxE = dx(ids(tdim, "Electrode")) + dxAll = dxE + dx(ids(tdim, "InletChannel", "OutletChannel")) + Q0 = dolfinx.fem.functionspace(mesh, ("DG", 0)) + mask_E = dolfinx.fem.Function(Q0) + mask_E.x.array[marked.entities(tdim, "Electrode")] = 1.0 + mask_F = dolfinx.fem.Function(Q0) + mask_F.x.array[:] = 1.0 - mask_E.x.array + else: + dxE = dxAll = dx + mask_E = mask_F = None + ds_cc = ds(ids(tdim - 1, "CurrentCollector")) + ds_mem = ds(ids(tdim - 1, "Membrane")) + + phi_s, phi_l, c_V2, c_V3, c_H, c_HSO4 = u + c_V2_old, c_V3_old, c_H_old, c_HSO4_old = u_old + v_s, v_l, v_c_V2, v_c_V3, v_c_H, v_c_HSO4 = v + _F, _R, _T, _eps = c["F"], c["R"], c["T"], c["epsilon"] + dt_c = ( + dt + if isinstance(dt, dolfinx.fem.Constant) + else dolfinx.fem.Constant(mesh, np.float64(dt)) + ) + _i_mem = ( + current + if isinstance(current, dolfinx.fem.Constant) + else dolfinx.fem.Constant(mesh, np.float64(current)) + ) + _i_cc = _i_mem * area_ratio + + # Electroneutrality + c_SO4 = ( + -(c["z_V2"] * c_V2 + c["z_V3"] * c_V3 + c["z_H"] * c_H + c["z_HSO4"] * c_HSO4) + / c["z_SO4"] + ) + eps_c = dolfinx.fem.Constant(mesh, 1e-14) + c_V2_eff = 0.5 * (c_V2 + ufl.sqrt(c_V2 * c_V2 + eps_c * eps_c)) + c_V3_eff = 0.5 * (c_V3 + ufl.sqrt(c_V3 * c_V3 + eps_c * eps_c)) + + delta_phi = phi_s - phi_l - c["E_formal"] + A = c["k0"] * ufl.exp(c["alpha_a"] * _F * delta_phi / (_R * _T)) + B = c["k0"] * ufl.exp(-c["alpha_c"] * _F * delta_phi / (_R * _T)) + if velocity is None: + c_V2_s, c_V3_s = c_V2_eff, c_V3_eff + else: + eps_speed = dolfinx.fem.Constant(mesh, 1e-16) + u_mag = ufl.sqrt(ufl.dot(velocity, velocity) + eps_speed**2) + Re = c["rho"] * u_mag * c["df"] / c["mu"] + k_mR = 7.0 * c["D_V2"] / c["df"] * (Re**0.4) + k_mO = 7.0 * c["D_V3"] / c["df"] * (Re**0.4) + den = 1.0 + A / (k_mR + eps_c) + B / (k_mO + eps_c) + c_V3_s = ( + (1.0 + A / (k_mR + eps_c)) * c_V3_eff + (A / (k_mO + eps_c)) * c_V2_eff + ) / den + c_V2_s = ( + (1.0 + B / (k_mO + eps_c)) * c_V2_eff + (B / (k_mR + eps_c)) * c_V3_eff + ) / den + i_F = _F * (A * c_V2_s - B * c_V3_s) + if mask_E is not None: + i_F = mask_E * i_F + f = c["a"] * i_F + + def effective(D): + if mask_E is None: + return D * _eps**1.5 + return mask_E * (D * _eps**1.5) + mask_F * D + + D_V2, D_V3, D_H, D_HSO4, D_SO4 = ( + effective(c[f"D_{ion}"]) for ion in ("V2", "V3", "H", "HSO4", "SO4") + ) + storage = _eps if mask_E is None else mask_E * _eps + mask_F * 1.0 + + kappa = (_F**2 / (_R * _T)) * ( + c["z_V2"] ** 2 * D_V2 * c_V2_eff + + c["z_V3"] ** 2 * D_V3 * c_V3_eff + + c["z_H"] ** 2 * D_H * c_H + + c["z_HSO4"] ** 2 * D_HSO4 * c_HSO4 + + c["z_SO4"] ** 2 * D_SO4 * c_SO4 + ) + i_l = -kappa * ufl.grad(phi_l) - _F * ( + c["z_V2"] * D_V2 * ufl.grad(c_V2) + + c["z_V3"] * D_V3 * ufl.grad(c_V3) + + c["z_H"] * D_H * ufl.grad(c_H) + + c["z_HSO4"] * D_HSO4 * ufl.grad(c_HSO4) + + c["z_SO4"] * D_SO4 * ufl.grad(c_SO4) + ) + sigma = (1.0 - _eps) ** 1.5 * c["sigma"] + residual = ufl.inner(sigma * ufl.grad(phi_s), ufl.grad(v_s)) * dxE + f * v_s * dxE + residual += ufl.inner(-i_l, ufl.grad(v_l)) * dxAll - ufl.inner(f, v_l) * dxE + + def flux(D, z, conc): + N = -D * ufl.grad(conc) - D * z * _F * conc / (_R * _T) * ufl.grad(phi_l) + return N if velocity is None else N + velocity * conc + + N_V2 = flux(D_V2, c["z_V2"], c_V2) + N_V3 = flux(D_V3, c["z_V3"], c_V3) + N_H = flux(D_H, c["z_H"], c_H) + N_HSO4 = flux(D_HSO4, c["z_HSO4"], c_HSO4) + residual += ( + storage * (c_V2 - c_V2_old) / dt_c * v_c_V2 * dxAll + - ufl.inner(N_V2, ufl.grad(v_c_V2)) * dxAll + + (1.0 / _F) * f * v_c_V2 * dxE + ) + residual += ( + storage * (c_V3 - c_V3_old) / dt_c * v_c_V3 * dxAll + - ufl.inner(N_V3, ufl.grad(v_c_V3)) * dxAll + - (1.0 / _F) * f * v_c_V3 * dxE + ) + residual += ( + storage * (c_H - c_H_old) / dt_c * v_c_H * dxAll + - ufl.inner(N_H, ufl.grad(v_c_H)) * dxAll + ) + residual += ( + storage * (c_HSO4 - c_HSO4_old) / dt_c * v_c_HSO4 * dxAll + - ufl.inner(N_HSO4, ufl.grad(v_c_HSO4)) * dxAll + ) + residual -= _i_cc * v_s * ds_cc + residual += _i_mem * v_l * ds_mem + if membrane_flux: + residual += _i_mem / (c["z_H"] * _F) * v_c_H * ds_mem + if outlet: + n = ufl.FacetNormal(mesh) + ds_out = ds(ids(tdim - 1, "Outlet")) + for conc, test in ( + (c_V2, v_c_V2), + (c_V3, v_c_V3), + (c_H, v_c_H), + (c_HSO4, v_c_HSO4), + ): + residual += ufl.dot(velocity, n) * conc * test * ds_out + return residual + + +def _darcy_coefficient(mesh): + """The Darcy coefficient, permeability over viscosity, of the 2D example.""" + c = _constants(mesh, PARAMETERS) + return ( + (c["df"] ** 2 / (c["Kc"] * c["mu"])) + * (c["epsilon"] ** 3) + / ((1.0 - c["epsilon"]) ** 2) + ) + + +def darcy(marked: MeshData, p, q): + """The pressure equation of the 2D example, with the inflow.""" + mesh = marked.mesh + dx, ds, ids = measures(marked) + u_in = dolfinx.fem.Constant(mesh, np.float64(PARAMETERS["u_in"])) + a = ufl.inner(_darcy_coefficient(mesh) * ufl.grad(p), ufl.grad(q)) * dx + return a - u_in * q * ds(ids(mesh.topology.dim - 1, "Inlet")) + + +def darcy_velocity(marked: MeshData, p): + """The Darcy velocity of the 2D example.""" + return -_darcy_coefficient(marked.mesh) * ufl.grad(p) + + +DIFFUSIVITY = syu.m**2 / syu.second +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("H", 1.0, 9.312e-9 * DIFFUSIVITY), + Ion("HSO4", -1.0, 1.33e-9 * DIFFUSIVITY), + Ion("SO4", -2.0, 1.065e-9 * DIFFUSIVITY), +) +"""The ions of the examples, sulfate by electroneutrality.""" + +RESERVOIR = { + "V2": 300.0 * CONCENTRATION, + "V3": 300.0 * CONCENTRATION, + "H": 1000.0 * CONCENTRATION, + "HSO4": 1000.0 * CONCENTRATION, +} +"""The concentrations of the reservoir of the examples.""" + + +def electrode(domain: Domain, markers: tuple[str, ...] = ()) -> PorousElectrode: + """The porous electrode of the examples.""" + return PorousElectrode( + "electrode", + domain, + markers, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=0.5 * syu.micrometer, + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, + ) + + +def half_cell(domain: Domain, porous=None, **kwargs) -> RedoxHalfCell: + """The negative half-cell of the examples, V2+ oxidized on discharge.""" + return RedoxHalfCell( + domain, + IONS, + RESERVOIR, + "V2", + "V3", + electrode(domain) if porous is None else porous, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + **kwargs, + ) + + +@pytest.fixture +def reference(): + return types.SimpleNamespace( + PARAMETERS=PARAMETERS, + QUADRATURE_DEGREE=QUADRATURE_DEGREE, + W=W, + H=H, + interval=interval, + rectangle=rectangle, + box=box, + measures=measures, + IONS=IONS, + RESERVOIR=RESERVOIR, + electrode=electrode, + half_cell=half_cell, + electrochemistry=electrochemistry, + darcy=darcy, + darcy_velocity=darcy_velocity, + ) diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py new file mode 100644 index 0000000..929ff78 --- /dev/null +++ b/tests/electrochemistry/half_cell/test_models.py @@ -0,0 +1,203 @@ +"""The RFBniCS models: units, structure, schedules and post-processing.""" + +import dolfinx.fem +import numpy as np +import pytest +import sympy.physics.units as syu + +from sumfem.discrete import DomainMaps, discretize +from sumfem.electrochemistry import RedoxHalfCell +from sumfem.models import ( + BoundaryFlux, + BrinkmanForchheimer, + Darcy, + GradDiv, + NavierStokes, +) +from sumfem.post_processing import evaluate_over_line, function_range +from sumfem.solvers import BDF +from sumfem.symbolic import Domain, Field, Parameter, build_system + + +def fluid(domain): + return ( + Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3, + unit=syu.pascal * syu.second, + ), + Parameter( + "electrolyte.density", domain, default=1399.8, unit=syu.kg / syu.m**3 + ), + ) + + +def assert_consistent(model): + """Units are consistent with and without time stepping, and when factorized.""" + build_system(model, check_units=True) + build_system(model, substitutions=BDF(1).substitutions(model), check_units=True) + scaled = build_system(model, factorize_equations=True) + assert all(scale > 0 for scale in scaled.equation_scales.values()) + + +def test_half_cells_have_consistent_units(reference): + line = Domain("cell", "interval") + cell = reference.half_cell(line) + assert [f.name for f in cell.model.equations] == [ + "phi_s", + "phi_l", + "c_V2", + "c_V3", + "c_H", + "c_HSO4", + ] + assert_consistent(cell.model) + + rect = Domain("cell", "quadrilateral") + viscosity, density = fluid(rect) + electrode = reference.electrode(rect) + p = Field("p", rect, ("Lagrange", 1), unit=syu.pascal) + inflow = Parameter("u_in", rect, default=1e-7, unit=syu.m / syu.second) + darcy = Darcy("darcy", p, electrode.permeability(), viscosity) + assert_consistent(darcy + BoundaryFlux("inflow", p, inflow, ("Inlet",))) + cell = reference.half_cell( + rect, + electrode, + velocity=darcy.superficial_velocity, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + assert_consistent(cell.model) + # The flow is data for the electrochemistry + assert list(build_system(cell.model).known_coefficients) == ["p"] + + box = Domain("cell", "tetrahedron") + viscosity, density = fluid(box) + electrode = reference.electrode(box, ("Electrode",)) + channels = ("InletChannel", "OutletChannel") + u = Field("u", box, ("Lagrange", 2, (3,)), unit=syu.m / syu.second, scale=1e-7) + p = Field("p", box, ("Lagrange", 1), unit=syu.pascal) + flow = BrinkmanForchheimer( + "porous_flow", + u, + p, + electrode.porosity, + electrode.permeability(), + electrode.forchheimer(), + viscosity, + density, + markers=("Electrode",), + ) + NavierStokes("channel_flow", u, p, viscosity, density, markers=channels) + flow += GradDiv("grad_div", u, 1e-3 * viscosity) + assert_consistent(flow) + cell = reference.half_cell( + box, + electrode, + channels=channels, + velocity=u, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + assert_consistent(cell.model) + names = [m.name for m in cell.model.components] + assert "transport_V2_channels" in names and "storage_H_channels" in names + assert "electrolyte_channels" in names + + +def test_shared_parameters_are_used_by_flow_and_kinetics(reference): + rect = Domain("cell", "quadrilateral") + viscosity, density = fluid(rect) + p = Field("p", rect, ("Lagrange", 1), unit=syu.pascal) + electrode = reference.electrode(rect) + darcy = Darcy("darcy", p, electrode.permeability(), viscosity) + cell = reference.half_cell( + rect, + electrode, + velocity=darcy.superficial_velocity, + density=density, + viscosity=viscosity, + ) + parameters = build_system(cell.model + darcy).parameters + assert parameters["electrolyte.viscosity"] is viscosity + assert parameters["electrode.porosity"] is electrode.porosity + assert "darcy.viscosity" not in parameters + + +def test_half_cell_errors(reference): + domain = Domain("cell", "interval") + with pytest.raises(ValueError, match="need a velocity"): + reference.half_cell(domain, mass_transfer=True) + with pytest.raises(ValueError, match="need a velocity"): + reference.half_cell(domain, outlet=("Outlet",)) + + +def test_the_couple_is_oxidized_or_reduced_on_discharge(reference): + """On discharge the current enters the solid at the current collector of a + half-cell that oxidizes its couple, and leaves it at that of one that + reduces it: the equations of the potentials, with only the current in + them, are opposite.""" + marked = reference.interval(num_cells=4) + residuals = [] + for oxidation in (True, False): + domain = Domain.from_mesh("cell", marked.mesh) + cell = RedoxHalfCell( + domain, + reference.IONS, + reference.RESERVOIR, + "V2", + "V3", + reference.electrode(domain), + oxidation_on_discharge=oxidation, + rate_constant=0.0, + formal_potential=0.0, + ) + scheme = BDF(1) + system = build_system( + cell.model, substitutions=scheme.substitutions(cell.model) + ) + disc = discretize([system], DomainMaps({domain: marked})) + scheme.set_step(1, 1.0, disc) + # Without gradients, rate or change, only the currents are left in + # the equations of the potentials + disc.set({cell.current: 0.5}) + problem = disc[system] + residuals.append( + [dolfinx.fem.assemble_vector(form).array for form in problem.residual[:2]] + ) + for oxidizing, reducing in zip(*residuals, strict=True): + assert np.abs(oxidizing).max() > 0 + np.testing.assert_allclose(oxidizing, -reducing) + + +def test_the_reservoir_names_ions_of_the_electrolyte(reference): + domain = Domain("cell", "interval") + with pytest.raises(ValueError, match="no ions"): + RedoxHalfCell( + domain, + reference.IONS, + {**reference.RESERVOIR, "Cl": 1.0}, + "V2", + "V3", + reference.electrode(domain), + oxidation_on_discharge=True, + rate_constant=2.5e-7, + formal_potential=-0.36, + ) + + +def test_post_processing(reference): + marked = reference.interval(num_cells=10) + domain = Domain.from_mesh("cell", marked.mesh) + c = Field("c", domain, ("Lagrange", 1)) + disc = discretize([c], DomainMaps({domain: marked})) + disc.set({c: lambda x: 1.0 + x[0] / reference.W}) + function = disc.function(c) + assert function_range(function) == pytest.approx((1.0, 2.0)) + points, values = evaluate_over_line(function, [0.0], [reference.W], 5) + np.testing.assert_allclose(values[:, 0], 1.0 + points[:, 0] / reference.W) + assert values.shape[1] == 1 diff --git a/tests/electrochemistry/half_cell/test_reference.py b/tests/electrochemistry/half_cell/test_reference.py new file mode 100644 index 0000000..8e1ab0c --- /dev/null +++ b/tests/electrochemistry/half_cell/test_reference.py @@ -0,0 +1,353 @@ +"""The models reproduce the forms and solutions of the RFBniCS examples.""" + +import basix.ufl +import dolfinx.fem +import dolfinx.fem.petsc +import numpy as np +import pytest +import sympy.physics.units as syu +import ufl +import ufl.algorithms + +from sumfem.discrete import DomainMaps, discretize +from sumfem.models import BoundaryFlux, Darcy +from sumfem.solvers import BDF, Simulation +from sumfem.symbolic import Clock, Domain, Field, Parameter, build_system + +IONS = ("V2", "V3", "H", "HSO4") + + +def fluid(domain): + """Viscosity and density of the examples, shared by flow and kinetics.""" + return ( + Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3, + unit=syu.pascal * syu.second, + ), + Parameter( + "electrolyte.density", domain, default=1399.8, unit=syu.kg / syu.m**3 + ), + ) + + +def unknowns(cell): + return [cell.phi_s, cell.phi_l, *(cell.concentrations[ion] for ion in IONS)] + + +def set_state(disc, cell, reference): + """A smooth state away from equilibrium, with a different previous step.""" + W = reference.W + disc.set( + { + cell.phi_s: lambda x: 0.01 * np.sin(np.pi * x[0] / W), + cell.phi_l: lambda x: 0.36 + 0.005 * x[0] / W, + cell.concentrations["V2"]: lambda x: ( + 300 * (1 + 0.1 * np.sin(2 * np.pi * x[0] / W)) + ), + cell.concentrations["V3"]: lambda x: ( + 300 * (1 - 0.1 * np.cos(np.pi * x[0] / W)) + ), + cell.concentrations["H"]: lambda x: 1000 * (1 + 0.05 * x[0] / W), + cell.concentrations["HSO4"]: lambda x: 1000 * (1 - 0.02 * x[0] / W), + } + ) + for ion in IONS: + field = cell.concentrations[ion] + previous = disc.function(field.previous(1)) + previous.x.array[:] = 0.99 * disc.function(field).x.array + + +def compare(disc, cell, problem, reference, dt, current, **options): + """Assemble both residuals, block by block, and compare them.""" + fields = unknowns(cell) + spaces = [disc.function_space(f) for f in fields] + tests = ufl.TestFunctions(ufl.MixedFunctionSpace(*spaces)) + marked = disc.maps.mesh_data(cell.phi_s.domain) + residual = reference.electrochemistry( + marked, + [disc.function(f) for f in fields], + [disc.function(cell.concentrations[ion].previous(1)) for ion in IONS], + tests, + dt, + current, + **options, + ) + expected = ufl.extract_blocks(residual) + for field, ours, theirs in zip(fields, problem.residual, expected, strict=True): + a = dolfinx.fem.assemble_vector(ours).array + b = dolfinx.fem.assemble_vector(dolfinx.fem.form(theirs)).array + scale = np.abs(b).max() + assert scale > 0, field.name + np.testing.assert_allclose( + a, b, rtol=1e-9, atol=1e-11 * scale, err_msg=field.name + ) + + +def backward_euler(cell): + """The scheme, and the system of one backward Euler step of the cell's model.""" + scheme = BDF(1) + return scheme, build_system( + cell.model, substitutions=scheme.substitutions(cell.model) + ) + + +def transient(disc, scheme, system, cell, dt, current): + """Bind the step's system, and set the step and current.""" + problem = disc[system] + scheme.set_step(1, dt, disc) + disc.set({cell.current: current}) + return problem + + +def test_1d_residual_matches_rfbnics(reference): + marked = reference.interval() + domain = Domain.from_mesh("cell", marked.mesh) + cell = reference.half_cell(domain) + scheme, system = backward_euler(cell) + disc = discretize( + [system], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + set_state(disc, cell, reference) + problem = transient(disc, scheme, system, cell, dt=0.1, current=-0.5) + compare(disc, cell, problem, reference, dt=0.1, current=-0.5) + + +def test_2d_residual_matches_rfbnics(reference): + marked = reference.rectangle() + domain = Domain.from_mesh("cell", marked.mesh) + viscosity, density = fluid(domain) + electrode = reference.electrode(domain) + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal) + flow = Darcy("darcy", p, electrode.permeability(), viscosity) + cell = reference.half_cell( + domain, + electrode, + velocity=flow.superficial_velocity, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + scheme, system = backward_euler(cell) + disc = discretize( + [system, p], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + set_state(disc, cell, reference) + # A pressure drop of 1 kPa from inlet to outlet + disc.set({p: lambda x: 3e5 + 1e3 * (1 - x[1] / reference.H)}) + problem = transient(disc, scheme, system, cell, dt=0.5, current=0.5) + velocity = reference.darcy_velocity(marked, disc.function(p)) + compare( + disc, + cell, + problem, + reference, + dt=0.5, + current=0.5, + velocity=velocity, + membrane_flux=True, + outlet=True, + ) + + +def test_2d_pressure_matches_rfbnics(reference): + marked = reference.rectangle() + domain = Domain.from_mesh("cell", marked.mesh) + viscosity, _ = fluid(domain) + electrode = reference.electrode(domain) + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal) + inflow = Parameter("u_in", domain, default=1e-7, unit=syu.m / syu.second) + model = Darcy("darcy", p, electrode.permeability(), viscosity) + BoundaryFlux( + "inflow", p, inflow, ("Inlet",) + ) + disc = discretize( + [model], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + disc.set({p: lambda x: 3e5 + 1e3 * x[0] * x[1] / (reference.W * reference.H)}) + problem = disc[model] + q = ufl.TestFunction(disc.function_space(p)) + expected = reference.darcy(marked, disc.function(p), q) + a = dolfinx.fem.assemble_vector(problem.residual).array + b = dolfinx.fem.assemble_vector(dolfinx.fem.form(expected)).array + np.testing.assert_allclose(a, b, rtol=1e-12, atol=1e-12 * np.abs(b).max()) + + +def test_3d_residual_matches_rfbnics(reference): + marked = reference.box() + domain = Domain.from_mesh("cell", marked.mesh) + viscosity, density = fluid(domain) + electrode = reference.electrode(domain, ("Electrode",)) + u = Field("u", domain, ("Lagrange", 2, (3,)), unit=syu.m / syu.second) + cell = reference.half_cell( + domain, + electrode, + channels=("InletChannel", "OutletChannel"), + velocity=u, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + scheme, system = backward_euler(cell) + disc = discretize( + [system], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + # The cube is 1 m wide, so scale the state of the 2.5 mm electrode + set_state(disc, cell, type(reference)(**{**vars(reference), "W": 1.0})) + disc.set({u: lambda x: (1e-4 * (1 + x[1]), 2e-5 * x[2], -1e-5 * x[0])}) + cell.balance_currents(disc) + ratio = float(disc.value(cell.area_ratio)) + assert ratio == pytest.approx(2.0) + problem = transient(disc, scheme, system, cell, dt=0.2, current=-0.5) + compare( + disc, + cell, + problem, + reference, + dt=0.2, + current=-0.5, + velocity=disc.function(u), + regions=True, + area_ratio=ratio, + membrane_flux=True, + outlet=True, + ) + + +def rfbnics_1d(marked, reference, schedule): + """The time loop of the 1D example of RFBniCS, on a blocked function.""" + mesh = marked.mesh + el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1) + V = dolfinx.fem.functionspace(mesh, basix.ufl.blocked_element(el, shape=(6,))) + X = dolfinx.fem.Function(V) + X_old = dolfinx.fem.Function(V) + parameters = reference.PARAMETERS + for k, ion in enumerate(IONS, start=2): + X.sub(k).interpolate( + lambda x, c=parameters[f"c0_{ion}"]: np.full(x.shape[1], c) + ) + X.sub(1).interpolate(lambda x: np.full(x.shape[1], -parameters["E_formal"])) + X_old.x.array[:] = X.x.array + u = ufl.split(X) + dt = dolfinx.fem.Constant(mesh, 0.1) + current = dolfinx.fem.Constant(mesh, 0.0) + residual = reference.electrochemistry( + marked, u, ufl.split(X_old)[2:], ufl.TestFunctions(V), dt, current + ) + tdim = mesh.topology.dim + membrane = marked.entities(tdim - 1, "Membrane") + bcs = [ + dolfinx.fem.dirichletbc( + np.float64(parameters[f"c0_{ion}"]), + dolfinx.fem.locate_dofs_topological(V.sub(k), tdim - 1, membrane), + V.sub(k), + ) + for k, ion in enumerate(IONS, start=2) + ] + problem = dolfinx.fem.petsc.NonlinearProblem( + residual, + X, + bcs=bcs, + petsc_options_prefix="rfbnics_", + petsc_options={ + "snes_type": "newtonls", + "snes_linesearch_type": "bt", + "snes_atol": 1e-8, + "snes_rtol": 1e-12, + "snes_stol": 1e-14, + "snes_error_if_not_converged": True, + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + # No potential is fixed, so the Jacobian is singular + "mat_mumps_icntl_24": 1, + "mat_mumps_icntl_25": 0, + }, + ) + for signed_current, steps in schedule: + current.value = signed_current + for count, step in steps: + dt.value = step + for _ in range(count): + problem.solve() + X_old.x.array[:] = X.x.array + return X, V + + +def test_1d_charge_and_discharge_match_rfbnics(reference): + marked = reference.interval(num_cells=30) + # The signed current density, and the number and size of the time steps + schedule = ( + (-0.5, ((3, 0.1), (2, 1.0))), + (0.0, ((2, 1.0),)), + (0.5, ((3, 1.0),)), + ) + X, V = rfbnics_1d(marked, reference, schedule) + + 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}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + cell.initialize(disc) + clock = Clock(domain) + simulation = Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=BDF(1), + petsc_options={"snes_atol": 1e-8, "snes_rtol": 1e-12, "snes_stol": 1e-14}, + ) + voltages, charge = [], [] + + def monitor(t, step): + voltages.append(cell.voltage(disc)) + charge.append(cell.state_of_charge(disc)) + + for signed_current, steps in schedule: + disc.set({cell.current: signed_current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run(num_steps=count, callback=monitor) + assert len(voltages) == 10 + assert float(disc.value(clock.t)) == pytest.approx(0.3 + 2.0 + 2.0 + 3.0) + values = {} + for k, field in enumerate(unknowns(cell)): + space, dofs = V.sub(k).collapse() + expected = dolfinx.fem.Function(space) + expected.x.array[:] = X.x.array[dofs] + ours = disc.function(field) + # Both are degree-one Lagrange on one interval, so order by vertex + order_ours = np.argsort(ours.function_space.tabulate_dof_coordinates()[:, 0]) + order_theirs = np.argsort( + expected.function_space.tabulate_dof_coordinates()[:, 0] + ) + values[field.name] = ours.x.array[order_ours], expected.x.array[order_theirs] + # Only differences of potentials are defined, so measure both potentials + # from the electrolyte potential at the membrane + for name in ("phi_s", "phi_l"): + values[name] = tuple(values[name][i] - values["phi_l"][i][-1] for i in range(2)) + for name, (ours, theirs) in values.items(): + scale = max(np.abs(theirs).max(), 1e-3) + np.testing.assert_allclose( + ours, theirs, rtol=1e-7, atol=1e-9 * scale, err_msg=name + ) + # Charging turns V3 into V2, and discharging turns it back + assert charge[4] > 0.5 + assert charge[-1] < charge[6] + assert charge[6] == pytest.approx(charge[4], rel=1e-6) + assert np.isfinite(voltages).all() diff --git a/tests/symbolic/test_terminals.py b/tests/symbolic/test_terminals.py index 9efe556..bbe7080 100644 --- a/tests/symbolic/test_terminals.py +++ b/tests/symbolic/test_terminals.py @@ -164,6 +164,19 @@ def test_the_default_of_a_parameter_in_si(domain): np.testing.assert_allclose(ratio.default_si, [1.0, 3.0]) +def test_a_default_with_units_is_converted_to_the_unit(domain): + """A quantity becomes a value in the declared unit, so that constants can + carry their units; another dimension is rejected.""" + length = Parameter("length", domain, default=5 * syu.centimeter, unit=syu.mm) + assert length.default == pytest.approx(50.0) + assert length.default_si == pytest.approx(0.05) + diffusivity = syu.m**2 / syu.s + rate = Parameter("D", domain, default=0.36 * syu.cm**2 / syu.hour, unit=diffusivity) + assert rate.default_si == pytest.approx(1e-8) + with pytest.raises(ValueError, match="does not have the dimension"): + Parameter("bad", domain, default=1.0 * syu.m / syu.s, unit=diffusivity) + + def test_scalar_default_sets_every_component(domain): np.testing.assert_allclose( Parameter("g", domain, shape=(2,), default=3.0).default, [3, 3] diff --git a/tests/test_api_docs.py b/tests/test_api_docs.py index 77fb05d..0ed9a82 100644 --- a/tests/test_api_docs.py +++ b/tests/test_api_docs.py @@ -23,6 +23,15 @@ def documented(page: pathlib.Path) -> set[str]: return names +def resolve(name: str) -> object: + """The object a documented full name refers to, or ``None``.""" + module, _, attribute = name.rpartition(".") + try: + return getattr(importlib.import_module(module), attribute, None) + except ImportError: + return None + + @pytest.mark.parametrize( "module, page", [ @@ -39,12 +48,18 @@ def documented(page: pathlib.Path) -> set[str]: ], ) def test_every_public_name_is_on_its_api_page(module, page): - public = importlib.import_module(module).__all__ + imported = importlib.import_module(module) names = documented(API / f"{page}.md") - # A name may be documented under its module or a package that re-exports it + # A name may be documented under its module, a package that re-exports + # it, or the module that defines it, as a constant must be, whose + # docstring autodoc finds only in the source of the module that defines it parts = module.split(".") owners = [".".join(parts[: k + 1]) for k in range(1, len(parts))] + objects = [resolve(name) for name in names] missing = sorted( - n for n in public if not any(f"{owner}.{n}" in names for owner in owners) + n + for n in imported.__all__ + if not any(f"{owner}.{n}" in names for owner in owners) + and not any(o is getattr(imported, n) for o in objects if o is not None) ) assert not missing, f"{missing} of {module} are missing from docs/api/{page}.md" From d0cf7073904cb73b7038648881395b30231d963a Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 21:40:50 +0000 Subject: [PATCH 2/6] Add to toc --- _toc.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_toc.yml b/_toc.yml index 6eefd2f..d3017c5 100644 --- a/_toc.yml +++ b/_toc.yml @@ -24,6 +24,7 @@ parts: - caption: "Electrochemistry" chapters: - file: "examples/electrochemistry/porous_electrode.py" + - file: "examples/electrochemistry/vanadium_half_cell_1d.py" - caption: "Python API" chapters: - file: "docs/api/index" From 357e7e7f698e1a0e0f3436e1b3a91af525f064df Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 22:28:20 +0000 Subject: [PATCH 3/6] Add 2D and 3D redox models --- _config.yml | 7 + _toc.yml | 2 + docs/api/symbolic.md | 2 + .../electrochemistry/vanadium_half_cell_1d.py | 174 ++++++- .../electrochemistry/vanadium_half_cell_2d.py | 421 +++++++++++++++ .../electrochemistry/vanadium_half_cell_3d.py | 488 ++++++++++++++++++ references.bib | 69 +++ src/sumfem/electrochemistry/electrolyte.py | 8 +- src/sumfem/electrochemistry/half_cell.py | 40 +- src/sumfem/electrochemistry/kinetics.py | 6 +- src/sumfem/models/porous.py | 12 +- src/sumfem/symbolic/terminals.py | 50 ++ tests/electrochemistry/half_cell/conftest.py | 3 +- .../electrochemistry/half_cell/test_models.py | 39 ++ tests/electrochemistry/test_electrolyte.py | 9 +- tests/symbolic/test_terminals.py | 28 + 16 files changed, 1321 insertions(+), 37 deletions(-) create mode 100644 examples/electrochemistry/vanadium_half_cell_2d.py create mode 100644 examples/electrochemistry/vanadium_half_cell_3d.py diff --git a/_config.yml b/_config.yml index 5164b49..42ae3e3 100644 --- a/_config.yml +++ b/_config.yml @@ -61,6 +61,13 @@ sphinx: - fmt: py # A tutorial or demo that raises fails the build, instead of a warning nb_execution_raise_on_error: true + # Pages that cite references cited elsewhere too list them as footnotes, + # {footcite} and {footbibliography}. MyST checks footnotes before + # sphinxcontrib-bibtex creates them, and warns about each; the book has + # no MyST footnotes of its own + suppress_warnings: ["myst.footnote"] + # No rule before the footnotes, which end a page under References + myst_footnote_transition: false autodoc_member_order: bysource # The page tables of contents list objects by name, without their module toc_object_entries_show_parents: hide diff --git a/_toc.yml b/_toc.yml index d3017c5..4b1e5d9 100644 --- a/_toc.yml +++ b/_toc.yml @@ -25,6 +25,8 @@ parts: chapters: - file: "examples/electrochemistry/porous_electrode.py" - file: "examples/electrochemistry/vanadium_half_cell_1d.py" + - file: "examples/electrochemistry/vanadium_half_cell_2d.py" + - file: "examples/electrochemistry/vanadium_half_cell_3d.py" - caption: "Python API" chapters: - file: "docs/api/index" diff --git a/docs/api/symbolic.md b/docs/api/symbolic.md index d2cb516..3de36ed 100644 --- a/docs/api/symbolic.md +++ b/docs/api/symbolic.md @@ -44,6 +44,8 @@ :members: .. autofunction:: as_expr + +.. autofunction:: value_at_defaults ``` ## Components diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py index 401dff3..9eddc8d 100644 --- a/examples/electrochemistry/vanadium_half_cell_1d.py +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -16,7 +16,8 @@ # # This demo charges, rests and discharges the negative half-cell of a # vanadium redox flow battery across its porous electrode, as the 1D example -# of [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS). It shows how to: +# of [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS) +# {footcite}`Aghabarari2026`. It shows how to: # # - Build a half-cell of any redox couple and electrolyte, negative or # positive, with @@ -25,7 +26,7 @@ # with a {py:class}`Simulation ` # - Follow the voltage and the state of charge of the half-cell # -# ## The model +# ## The mathematical model # # The porous electrode holds the solid, of potential $\phi_s$, and in its pores # of porosity $\varepsilon$ the electrolyte, of potential $\phi_l$, with the @@ -65,8 +66,10 @@ # the electrolyte at the membrane, where the electrolyte has the # concentrations of the reservoir. The voltage of the half-cell is the mean # solid potential at the current collector less the mean electrolyte -# potential at the membrane. The parameters are those of the RFBniCS -# examples. +# potential at the membrane. The model is the macro-homogeneous +# porous-electrode model of RFBniCS {footcite}`Aghabarari2026` (Section II), +# after {footcite}`NewmanTiedemann1975`, and its parameters those of Table I +# there, several from {footcite}`MunozPerales2022`. # # We import the relevant packages: @@ -109,14 +112,60 @@ domain = symbolic.Domain.from_mesh("cell", interval, length_scale=width) # - +# The porous electrode, which fills the interval, is a carbon felt, a mat of +# carbon fibres, with its porosity, the conductivity of its fibres, their +# specific area and diameter, see +# {py:class}`PorousElectrode `. +# The fibres of commercial carbon felts, papers and cloths share a diameter +# of several micrometres {footcite}`Barzegari2026`; RFBniCS +# {footcite}`Aghabarari2026` uses a diameter far below that range, kept in a +# comment. The diameter and the +# constants of the permeability of the pores and of the inertial drag matter +# only with flow, and this demo has none: + +electrode = PorousElectrode( + "electrode", + domain, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=8.0 * syu.micrometer, # RFBniCS: 0.5 micrometre + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) + + # ## The half-cell # # 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 -# of the reservoir name the ions that are transported, so sulfate, which has -# none, follows from electroneutrality: +# 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, +# +# $$ +# \begin{aligned} +# 2 c_{\mathrm{V^{2+}}} + 3 c_{\mathrm{V^{3+}}} + c_{\mathrm{H^+}} +# - c_{\mathrm{HSO_4^-}} - 2 c_{\mathrm{SO_4^{2-}}} &= 0, +# \end{aligned} +# $$ (electroneutrality) +# +# is solved for its concentration, which is then no unknown, but an +# expression of the others: +# +# $$ +# \begin{aligned} +# c_{\mathrm{SO_4^{2-}}} &= \frac{1}{2} \left( 2 c_{\mathrm{V^{2+}}} +# + 3 c_{\mathrm{V^{3+}}} + c_{\mathrm{H^+}} +# - c_{\mathrm{HSO_4^-}} \right). +# \end{aligned} +# $$ (sulfate) +# +# It enters the conductivity and the diffusion current of the electrolyte +# like the concentrations of the other ions. The half-cell checks that the +# reservoir lacks exactly one ion, and that this ion is charged: # + diffusivity = syu.m**2 / syu.second @@ -136,27 +185,16 @@ } # - -# The porous electrode is a carbon felt, a mat of carbon fibres, with its -# porosity, the conductivity of its fibres, their specific area and diameter, -# see {py:class}`PorousElectrode `. -# The constants of the permeability of its pores and of the inertial drag -# matter only with flow, and this demo has none: - -electrode = PorousElectrode( - "electrode", - domain, - porosity=0.94, - conductivity=3402.1 * syu.siemens / syu.m, - specific_area=2.3e4 / syu.m, - fiber_diameter=0.5 * syu.micrometer, - kozeny_carman_constant=5.55, - pore_length=100.0 * syu.micrometer, -) - # {py:class}`RedoxHalfCell ` sums the -# models of {eq}`half-cell` for a couple of these ions, with its standard rate -# constant $k_0$ and formal potential $E_0$. Discharge oxidizes V²⁺ in the -# negative half-cell; the positive half-cell of the battery, with the couple +# models of {eq}`half-cell` for a redox couple. The reaction +# $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ of {eq}`half-cell-kinetics` +# is given by three arguments: the ion `reduced`, V²⁺, the ion `oxidized`, +# V³⁺, and the number of `electrons` the reaction transfers, which the +# charges of the ions must match. The kinetics take the standard rate +# constant $k_0$ and the formal potential $E_0$ of the couple. +# The other ions take no part in the reaction; they carry current and keep +# the electrolyte neutral. Discharge oxidizes V²⁺ in this negative +# half-cell; the positive half-cell of the battery, with the couple # VO²⁺/VO₂⁺, reduces its couple on discharge, which reverses the current. The # concentrations of the reservoir are Dirichlet conditions at the membrane: @@ -166,6 +204,7 @@ reservoir, reduced="V2", oxidized="V3", + electrons=1, electrode=electrode, oxidation_on_discharge=True, rate_constant=2.5e-7 * syu.m / syu.second, @@ -185,7 +224,7 @@ ) cell.initialize(disc) -# ## The schedule +# ## The charging and discharging schedule # # The cell is charged, rested and discharged, each with time steps that # start short, while the cell responds fast, and then grow. Each operation @@ -202,6 +241,82 @@ sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} # - +# ## Dimensionless groups +# +# The model is semi-discrete, so its units are checked on its system alone, +# see {ref}`units-checking`, with the length of the whole schedule as the +# reference time of the time derivatives. With the width of the electrode, +# the scales of the concentration of V²⁺ and of the solid potential, +# Faraday's constant, the temperature and the diffusivity of V²⁺ as the +# repeating quantities, {py:func}`pi_groups ` +# gives one group for every other quantity of the model: + +# + +duration = sum(count * dt for *_, steps in schedule for count, dt in steps) +semi_discrete = symbolic.build_system(cell.model, check_units=True, time_scale=duration) +electrolyte = cell.electrolyte +repeating = [ + domain, + cell.concentrations["V2"], + cell.phi_s, + electrolyte.F, + electrolyte.T, + electrolyte.diffusivity("V2"), +] +groups = symbolic.pi_groups(semi_discrete, time_scale=duration, repeating=repeating) +for group in groups: + print(f"{group.expression} = {group.value:.3g}") + + +def value_of(name: str) -> float: + """The value of the one group of the quantity of a name.""" + (value,) = [g.value for g in groups if name in str(g.expression)] + return value + + +# - + +# Most groups are ratios of like quantities: of the diffusivities and the +# concentrations of the ions, the porosity and the transfer coefficients. +# Four describe the behaviour of the cell: +# +# - The group of the reference time is the [mass Fourier +# number](https://dimless.org/docs/numbers/mass-fourier/) $D t / W^2$ of +# V²⁺ over the whole schedule. It is far below one: the ions diffuse only +# a small part of the way across the electrode, so the electrode exchanges +# little with the reservoir at the membrane. +# - The groups of the rate constant and the specific area multiply to the +# Damköhler number $a k_0 W^2 / D$ {footcite}`Fogler2020`, the rate of the +# reaction over that of diffusion across the electrode. Its square root is +# the Thiele modulus {footcite}`Thiele1939`, the width of the electrode over +# the depth a reactant supplied from the membrane would reach before it +# reacts. Both are far above one. +# - The group of the current of the cell is the current over the limiting +# current of diffusion across the electrode, $F D c / W$ +# {footcite}`NewmanBalsara2021`. Above one, the cell runs on the vanadium +# stored in its pores, not on what diffuses in from the reservoir. The +# current has no default, so the analysis takes a unit current density, +# and the group is scaled with the current of the schedule. +# - The groups of the conductivity of the solid and of the gas constant +# multiply to $\sigma R T / (F^2 D c)$, the conductivity of the solid over +# that V²⁺ alone would give the electrolyte. +# +# Neither the Damköhler number nor the Thiele modulus is yet on +# [dimless.org](https://dimless.org/docs/numbers/): + +# + +fourier = value_of("t_ref") +damkohler = value_of("rate_constant") * value_of("specific_area") +largest_current = max(current for _, current, _ in schedule) +current_ratio = value_of("cell.current") * largest_current +conductivity_ratio = value_of("electrode.conductivity") * value_of("electrolyte.R") +print(f"Mass Fourier number of the schedule: {fourier:.3g}") +print(f"Damkoehler number: {damkohler:.3g}, Thiele modulus: {np.sqrt(damkohler):.3g}") +print(f"Current over the limiting current of diffusion: {current_ratio:.3g}") +print(f"Conductivity of the solid over that of V2+: {conductivity_ratio:.3g}") +assert fourier < 1e-2 and damkohler > 1e2 and current_ratio > 1 +# - + # A {py:class}`Simulation ` takes the steps with # backward Euler, checks the units, and solves every time step in # dimensionless form, see {ref}`units-scaling`. The steps of the schedule @@ -277,3 +392,8 @@ bottom.set_xlabel("Time (s)") top.legend() plt.show() + +# ## References +# +# ```{footbibliography} +# ``` diff --git a/examples/electrochemistry/vanadium_half_cell_2d.py b/examples/electrochemistry/vanadium_half_cell_2d.py new file mode 100644 index 0000000..95ae2bd --- /dev/null +++ b/examples/electrochemistry/vanadium_half_cell_2d.py @@ -0,0 +1,421 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.13.6 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# # A vanadium half-cell in two dimensions, with flow +# +# This demo extends {doc}`vanadium_half_cell_1d` with the flow of the +# electrolyte through the porous electrode, as the 2D example of +# [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS) {footcite}`Aghabarari2026`. +# It shows how to: +# +# - Solve Darcy flow through the electrode with +# {py:class}`Darcy `, and check it against its exact +# solution +# - Couple the half-cell to the flow, which carries the ions in at the inlet +# and out at the outlet, and limits the kinetics by mass transfer to the +# fibres +# - Name the dimensionless numbers of the flow from the expressions of the +# model +# - Run the half-cell through a charging and discharging schedule, and plot +# its local state of charge +# +# ## The model +# +# The electrode spans the width $W$ from the current collector at $x = 0$ to +# the membrane at $x = W$, and the height $H$ from the inlet at $y = 0$ to +# the outlet at $y = H$. The electrolyte enters at the inlet with the +# superficial velocity $v_{\text{in}}$ and leaves at the outlet, held at the +# pressure $p_{\text{out}}$. In the pores, Darcy's law gives the +# superficial velocity $\mathbf{u}$ from the pressure $p$, with the +# permeability $K$ of the electrode and the viscosity $\mu$ of the +# electrolyte: +# +# $$ +# \begin{aligned} +# \mathbf{u} &= -\frac{K}{\mu} \nabla p, \\ +# \nabla \cdot \mathbf{u} &= 0, \\ +# -\mathbf{u} \cdot \mathbf{n} &= v_{\text{in}} \quad \text{at the inlet}, \\ +# p &= p_{\text{out}} \quad \text{at the outlet}, +# \end{aligned} +# $$ (darcy) +# +# and no flow through the other sides. Darcy's law is the balance of the +# pressure gradient and the Darcy drag in the volume-averaged Navier-Stokes +# equations, whose inertia, viscous diffusion and Forchheimer drag are +# negligible in the electrode {footcite}`Aghabarari2026` (Section II B 5). The +# permeability follows from the porosity and the fibre diameter by the +# Kozeny-Carman relation, see +# {py:func}`kozeny_carman `. +# +# The half-cell is that of {eq}`half-cell` in {doc}`vanadium_half_cell_1d`, +# with three changes. The Nernst-Planck fluxes carry the ions with +# $\mathbf{u}$. The concentrations of the reservoir hold at the inlet +# instead of the membrane, the ions leave with the flow at the outlet, and +# the membrane passes the current of the cell as protons. The reaction takes +# the concentrations at the surface of the fibres, which mass transfer from +# the pores, at the rate $k_m (c - c_s)$, ties to the concentrations $c$ of +# the pores, with a mass-transfer coefficient $k_m$ from the correlation of +# {footcite}`Krishnamurthy2011`, see +# {py:class}`~sumfem.electrochemistry.MassTransferButlerVolmer`. These are the +# boundary conditions and kinetics of {footcite}`Aghabarari2026` (Sections II B 4 +# and II C). +# The flow is steady and does not depend on the electrochemistry, so it is +# solved once, first. +# +# We import the relevant packages: + +# + tags=["hide-input"] +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.mesh +import matplotlib.pyplot as plt +import numpy as np +import pyvista +import sympy.physics.units as syu +import ufl + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.models import BoundaryFlux, Darcy +from sumfem.plotting import gather_grid + +# - + +# ## The electrode +# +# The cell is a rectangle of quadrilaterals, with its sides marked. They are +# located with a tolerance relative to the height, see +# {ref}`units-length-scale`: + +# + +width, height, cells = 2.5e-3, 1e-2, (25, 100) +dtype = dolfinx.default_real_type +rectangle = dolfinx.mesh.create_rectangle( + MPI.COMM_WORLD, + np.array([[0.0, 0.0], [width, height]]), + list(cells), + cell_type=dolfinx.mesh.CellType.quadrilateral, + dtype=dtype, +) +tolerance = 100 * float(np.finfo(dtype).eps) * height +marked = discrete.MeshData.from_locators( + rectangle, + { + 1: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0, atol=tolerance), + "Membrane": lambda x: np.isclose(x[0], width, atol=tolerance), + "Inlet": lambda x: np.isclose(x[1], 0.0, atol=tolerance), + "Outlet": lambda x: np.isclose(x[1], height, atol=tolerance), + } + }, + boundary=True, +) +domain = symbolic.Domain.from_mesh("cell", rectangle, length_scale=width) +# - + +# The porous electrode is the carbon felt of {doc}`vanadium_half_cell_1d`, +# whose fibres of several micrometres {footcite}`Barzegari2026` now matter: their +# diameter sets the permeability and the mass transfer. RFBniCS +# {footcite}`Aghabarari2026` uses a diameter far below that range, kept in a +# comment: + +electrode = PorousElectrode( + "electrode", + domain, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=8.0 * syu.micrometer, # RFBniCS: 0.5 micrometre + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) + +# ## The flow +# +# The density and the viscosity of the electrolyte are parameters that the +# flow and the half-cell share, as are the inlet velocity and the outlet +# pressure. The flow depends on the gradient of the pressure alone, so the +# pressure is the gauge pressure, zero at the outlet: an offset as large as +# the atmospheric pressure would cost the small pressure drop its digits. + +# + +viscosity = symbolic.Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3 * syu.pascal * syu.second, + unit=syu.pascal * syu.second, + doc="Dynamic viscosity of the electrolyte", +) +density = symbolic.Parameter( + "electrolyte.density", + domain, + default=1399.8 * syu.kg / syu.m**3, + unit=syu.kg / syu.m**3, + doc="Density of the electrolyte", +) +inlet_velocity = symbolic.Parameter( + "flow.inlet_velocity", + domain, + default=1e-7 * syu.m / syu.second, + unit=syu.m / syu.second, + doc="Superficial velocity at the inlet", +) +outlet_pressure = symbolic.Parameter( + "flow.outlet_pressure", + domain, + default=0.0 * syu.pascal, + unit=syu.pascal, + doc="Gauge pressure at the outlet", +) +# - + +# The pressure drop $\mu v_{\text{in}} H / K$ of the uniform flow through the +# electrode is the scale of the pressure. +# {py:func}`value_at_defaults ` evaluates +# it, an expression of parameters, at their defaults: + +# + +pressure_drop = symbolic.value_at_defaults( + viscosity * inlet_velocity * height / electrode.permeability() +) +p = symbolic.Field("p", domain, ("Lagrange", 1), unit=syu.pascal, scale=pressure_drop) +# - + +# {py:class}`Darcy ` is the flow of {eq}`darcy`, and a +# {py:class}`BoundaryFlux ` its inflow. Its system +# is checked for units and solved in dimensionless form, see +# {ref}`units-scaling`: + +darcy = Darcy("darcy", p, electrode.permeability(), viscosity) +flow = darcy + BoundaryFlux("inflow", p, inlet_velocity, ("Inlet",)) +flow_system = symbolic.build_system(flow, check_units=True, factorize_equations=True) +outlet = symbolic.DirichletBC(p, "Outlet", outlet_pressure) + +# ## The half-cell +# +# The ions and the reservoir are those of {doc}`vanadium_half_cell_1d`, with +# sulfate from electroneutrality: + +# + +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("H", 1, 9.312e-9 * diffusivity), + Ion("HSO4", -1, 1.33e-9 * diffusivity), + Ion("SO4", -2, 1.065e-9 * diffusivity), +) +reservoir = { + "V2": 300.0 * concentration, + "V3": 300.0 * concentration, + "H": 1000.0 * concentration, + "HSO4": 1000.0 * concentration, +} +# - + +# The half-cell is the negative one of {doc}`vanadium_half_cell_1d`. The +# superficial velocity of the flow carries the ions and drives the mass +# transfer to the fibres, which a velocity switches on. The flow carries the +# ions out at the `outlet`, the membrane passes the current as protons with +# `proton_flux`, and the reservoir concentrations hold at the inlet: + +cell = RedoxHalfCell( + domain, + ions, + reservoir, + reduced="V2", + oxidized="V3", + electrons=1, + electrode=electrode, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + velocity=darcy.superficial_velocity, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, +) +inlet = cell.inlet_conditions("Inlet") + +# The flow and the half-cell are discretized together, so they share their +# parameters, with the quadrature degree of the nonlinear kinetics bounded as +# in {doc}`vanadium_half_cell_1d`. A {py:class}`Solver ` +# solves the flow: + +disc = discrete.discretize( + [flow_system, outlet, cell.model, *inlet], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +solvers.Solver(disc[flow_system], [outlet]).solve() + +# With a uniform permeability, the flow is uniform, and the pressure falls +# linearly from the inlet to the outlet by $\mu v_{\text{in}} H / K$. Linear +# elements hold it exactly, so the pressure drop and the outflow, the flux of +# the superficial velocity through the outlet, match it to round-off: + +# + +n = ufl.FacetNormal(domain) +ds_inlet, ds_outlet = domain.ds("Inlet"), domain.ds("Outlet") +area = disc.assemble_scalar(1.0 * domain.dx()) +v_in = float(disc.value(inlet_velocity)) +drop = disc.assemble_scalar((p - outlet_pressure) * ds_inlet) / width +outflow = disc.assemble_scalar(ufl.dot(darcy.superficial_velocity, n) * ds_outlet) +print(f"Pressure drop: {drop:.6g} Pa, exact {pressure_drop:.6g} Pa") +print(f"Outflow over inflow: {outflow / (v_in * width):.12f}") +assert np.isclose(drop, pressure_drop, rtol=1e-8) +assert np.isclose(outflow, v_in * width, rtol=1e-8) +# - + +# ## Dimensionless numbers of the flow +# +# The flow past the fibres has the [Reynolds +# number](https://dimless.org/docs/numbers/reynolds/) $\rho |\mathbf{u}| d / +# \mu$ of the fibre diameter $d$, which the kinetics compute, and sets the +# [Sherwood number](https://dimless.org/docs/numbers/sherwood/) $k_m d / D$ +# of the mass transfer by the correlation of +# {py:class}`~sumfem.electrochemistry.MassTransferButlerVolmer`. +# Along the electrode, the [mass Péclet +# number](https://dimless.org/docs/numbers/mass-peclet/) $v_{\text{in}} H / D$ +# of V²⁺ compares the transport by the flow with diffusion. Their means over +# the electrode: + +# + +reaction = cell.reaction +reynolds = disc.assemble_scalar(reaction.reynolds_number() * domain.dx()) / area +d_V2 = float(disc.value(cell.electrolyte.diffusivity("V2"))) +d_fibre = float(disc.value(electrode.fiber_diameter)) +(k_V2, _) = reaction.mass_transfer_coefficients() +sherwood = disc.assemble_scalar(k_V2 * domain.dx()) / area * d_fibre / d_V2 +peclet = v_in * height / d_V2 +print(f"Reynolds number of the fibres: {reynolds:.3g}") +print(f"Sherwood number of V2+: {sherwood:.3g}") +print(f"Mass Peclet number of V2+ along the electrode: {peclet:.3g}") +assert reynolds < 1 and peclet > 1 +# - + +# The flow is creeping, so Darcy's law holds, yet it carries the ions along +# the electrode faster than they diffuse. +# +# ## Charge, rest and discharge +# +# The schedule of {doc}`vanadium_half_cell_1d`, run by a +# {py:class}`Simulation ` in parts, one for each +# current and step size. The simulation's model is the half-cell alone: the +# flow does not depend on the electrochemistry and its inputs are fixed, so +# every step takes the pressure solved above as known, and the velocity from +# it, without solving the flow again. The initial state is at rest, and the +# reference time of the time derivatives is the largest step: + +# + +schedule = ( + ("charge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), + ("rest", 0.0, ((5, 0.1), (5, 1.0), (10, 10.0))), + ("discharge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), +) +sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +cell.initialize(disc) +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 100}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +# - + +# After every step, a callback records the time, the mode, the voltage and +# the state of charge. At the end of the charge, the local state of charge, +# the share of V²⁺ in the vanadium at each point, is kept for a plot: + +# + +local = cell.local_state_of_charge() +space = disc.function_space(cell.concentrations["V2"]) +charged = dolfinx.fem.Function(space, name="soc") +records: list[tuple[float, str, float, float]] = [] +for mode, current, steps in schedule: + disc.set({cell.current: sign[mode] * current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run( + num_steps=count, + callback=lambda t, step, mode=mode: records.append( + (t, mode, cell.voltage(disc), cell.state_of_charge(disc)) + ), + ) + if mode == "charge": + charged.interpolate( + dolfinx.fem.Expression(disc[local], space.element.interpolation_points) + ) +times = np.array([r[0] for r in records]) +modes = np.array([r[1] for r in records]) +voltage = np.array([r[2] for r in records]) +charge = np.array([r[3] for r in records]) +# - + +# 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: + +# + +steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} +rest = charge[modes == "rest"] +print(f"State of charge after charge: {charge[modes == 'charge'][-1]:.6f}") +print(f"State of charge after rest: {rest[-1]:.6f}") +print(f"State of charge at the end: {charge[-1]:.6f}") +assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0) +assert np.all(np.diff(rest) < 0) +assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() +assert np.isfinite(voltage).all() +# - + +# The voltage and the state of charge through the schedule: + +fig, (top, bottom) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + top.plot(times[modes == mode], 1e3 * voltage[modes == mode], ".", label=mode) + bottom.plot(times[modes == mode], charge[modes == mode], ".", label=mode) +top.set_ylabel("Voltage (mV)") +bottom.set_ylabel("State of charge") +bottom.set_xlabel("Time (s)") +top.legend() +plt.show() + +# The local state of charge at the end of the charge, gathered from all +# processes with {py:func}`gather_grid `. The +# electrode is drawn in millimetres, widened to show its width: + +grid = gather_grid(charged) +if MPI.COMM_WORLD.rank == 0: + grid.points[:, 0] *= 4e3 + grid.points[:, 1] *= 1e3 + plotter = pyvista.Plotter(window_size=[480, 640]) + plotter.add_mesh(grid, scalars="soc", scalar_bar_args={"title": "SOC"}) + plotter.view_xy() + plotter.show() + +# ## References +# +# ```{footbibliography} +# ``` diff --git a/examples/electrochemistry/vanadium_half_cell_3d.py b/examples/electrochemistry/vanadium_half_cell_3d.py new file mode 100644 index 0000000..ef79e3c --- /dev/null +++ b/examples/electrochemistry/vanadium_half_cell_3d.py @@ -0,0 +1,488 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.13.6 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# # A vanadium half-cell in three dimensions, with flow channels +# +# This demo runs the half-cell of {doc}`vanadium_half_cell_2d` in a +# flow-through electrode with two channels, as the 3D example of +# [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS) {footcite}`Aghabarari2026`. +# It shows how to: +# +# - Build a mesh with physical groups in [Gmsh](https://gmsh.info), in +# memory, and mark it with them +# - Solve the flow of the electrolyte through the electrode and along the +# channels, Brinkman-Forchheimer in the electrode and Navier-Stokes in the +# channels, from a creeping-flow start +# - Couple the half-cell to the flow in the electrode and the channels, with +# the current of the cell balanced between the membrane and the smaller +# current collector +# +# ## The model +# +# The electrode is a thin plate. Two channels lie against its front: the +# electrolyte enters through the inlet channel, crosses the electrode along +# its height and leaves through the outlet channel. The front of the +# electrode that the channels leave free is the current collector, and its +# back the membrane. In the electrode, of porosity $\varepsilon$, the +# superficial velocity $\mathbf{u}$ and the pressure $p$ follow +# Brinkman-Forchheimer flow, with the Darcy drag of the permeability $K$ +# and the inertial drag of the Forchheimer coefficient $\beta$: +# +# $$ +# \begin{aligned} +# \frac{\rho}{\varepsilon^2} (\mathbf{u} \cdot \nabla) \mathbf{u} +# - \nabla \cdot \left( \frac{2 \mu}{\varepsilon} +# \boldsymbol{\epsilon}(\mathbf{u}) \right) + \nabla p +# + \frac{\mu}{K} \mathbf{u} + \rho \beta |\mathbf{u}| \mathbf{u} +# &= 0, \\ +# \nabla \cdot \mathbf{u} &= 0, +# \end{aligned} +# $$ (brinkman-forchheimer) +# +# with the symmetric gradient $\boldsymbol{\epsilon}$, see +# {py:class}`BrinkmanForchheimer `. In +# the channels, $\varepsilon = 1$ and the drag vanishes: the flow is +# Navier-Stokes flow. The velocity is given at the inlet, vanishes on the +# walls, the current collector and the membrane, and the outlet is free of +# stress. The Forchheimer coefficient is Ergun's, an estimate for fibres, +# see {py:func}`ergun_coefficient `. +# +# The half-cell is that of {doc}`vanadium_half_cell_2d`, with the +# electrolyte also in the channels, where it is stored and transported +# without porosity and without reaction. +# +# We import the relevant packages: + +# + tags=["hide-input"] +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.io.gmsh +import gmsh +import matplotlib.pyplot as plt +import numpy as np +import pyvista +import sympy.physics.units as syu +import ufl + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.models import Brinkman, BrinkmanForchheimer, GradDiv, NavierStokes, Stokes +from sumfem.plotting import gather_grid + +# - + +# ## The mesh +# +# The geometry is the flow-through flow field of RfbFoam +# {footcite}`Barzegari2026`, as RFBniCS {footcite}`Aghabarari2026` (Table III) builds +# it, from whose example (MIT licence) the code below is adapted, in +# millimetres: a plate electrode, and two channels with rounded ends that +# open at its front. A channel is a box with a cylinder at each end: + + +# + +def rounded_channel( + x_start: float, x_depth: float, y_center: float, z_start: float, z_end: float +) -> int: + """A channel of unit width along z with rounded ends, open along x.""" + radius = 0.5 + box = gmsh.model.occ.addBox( + x_start, y_center - radius, z_start, x_depth, 2 * radius, z_end - z_start + ) + caps = [ + gmsh.model.occ.addCylinder(x_start, y_center, z, x_depth, 0.0, 0.0, radius) + for z in (z_start, z_end) + ] + fused, _ = gmsh.model.occ.fuse([(3, box)], [(3, cap) for cap in caps]) + return fused[0][1] + + +# - + +# The electrode and the channels are fragmented into one conforming mesh, +# with the electrode and each channel a physical group of cells. The +# boundary surfaces are told apart by where they lie: the open ends of the +# channels are the inlet and the outlet, the free front of the electrode the +# current collector, its back the membrane, and the rest walls: + + +# + +def flow_through_electrode(sizes: tuple[float, float]) -> None: + """The Gmsh model of the cell, meshed with cells of the given sizes.""" + height, width, thickness, depth = 15.0, 17.0, 0.42, 1.0 + z_start, z_end, y_inlet, y_outlet = 1.5, 13.5, 1.0, 16.0 + electrode = gmsh.model.occ.addBox(depth, 0.0, 0.0, thickness, width, height) + channels = [ + rounded_channel(0.0, depth, y, z_start, z_end) for y in (y_inlet, y_outlet) + ] + gmsh.model.occ.synchronize() + _, fragments = gmsh.model.occ.fragment([(3, electrode)], [(3, c) for c in channels]) + gmsh.model.occ.synchronize() + for name, entities in zip( + ("Electrode", "InletChannel", "OutletChannel"), fragments, strict=True + ): + tags = [tag for dim, tag in entities if dim == 3] + gmsh.model.addPhysicalGroup(3, tags, name=name) + surfaces: dict[str, list[int]] = { + "Inlet": [], + "Outlet": [], + "CurrentCollector": [], + "Membrane": [], + "Walls": [], + } + tol = 1e-4 + for dim, tag in gmsh.model.getEntities(2): + cells, _ = gmsh.model.getAdjacencies(dim, tag) + if len(cells) == 2: + # An interface between the electrode and a channel + continue + xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(dim, tag) + if np.isclose(xmin, 0.0, atol=tol) and np.isclose(xmax, 0.0, atol=tol): + y = 0.5 * (ymin + ymax) + near_inlet = abs(y - y_inlet) < abs(y - y_outlet) + surfaces["Inlet" if near_inlet else "Outlet"].append(tag) + elif np.isclose(xmin, depth, atol=tol) and np.isclose(xmax, depth, atol=tol): + surfaces["CurrentCollector"].append(tag) + elif np.isclose(xmin, depth + thickness, atol=tol) and np.isclose( + xmax, depth + thickness, atol=tol + ): + surfaces["Membrane"].append(tag) + else: + surfaces["Walls"].append(tag) + for name, tags in surfaces.items(): + gmsh.model.addPhysicalGroup(2, tags, name=name) + gmsh.option.setNumber("Mesh.MeshSizeMin", sizes[0]) + gmsh.option.setNumber("Mesh.MeshSizeMax", sizes[1]) + gmsh.model.mesh.generate(3) + + +# - + +# The first process builds the model, and +# {py:func}`model_to_mesh ` distributes it +# without a file. The mesh is coarse, a few cells across the electrode, to +# keep the demo short; RFBniCS meshes far finer. Gmsh scales coordinates only +# when it writes a file, so they are scaled from millimetres to metres here. +# {py:meth}`MeshData.from_gmsh ` names +# the markers after the physical groups: + +# + +gmsh.initialize() +gmsh.option.setNumber("General.Terminal", 0) +if MPI.COMM_WORLD.rank == 0: + flow_through_electrode(sizes=(0.2, 0.5)) +data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3) +gmsh.finalize() +data.mesh.geometry.x[:] *= 1e-3 +marked = discrete.MeshData.from_gmsh(data) +thickness = 0.42e-3 +domain = symbolic.Domain.from_mesh("cell", data.mesh, length_scale=thickness) +num_cells = data.mesh.topology.index_map(3).size_global +print(f"Cells: {num_cells}") +# - + +# ## The electrode and the electrolyte +# +# The electrode of {doc}`vanadium_half_cell_2d` fills the cells marked +# `Electrode`, the channels the others. The density and the viscosity of the +# electrolyte are shared by the flow and the half-cell: + +# + +electrode = PorousElectrode( + "electrode", + domain, + markers=("Electrode",), + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=8.0 * syu.micrometer, # RFBniCS: 0.5 micrometre + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) +channels = ("InletChannel", "OutletChannel") +viscosity = symbolic.Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3 * syu.pascal * syu.second, + unit=syu.pascal * syu.second, + doc="Dynamic viscosity of the electrolyte", +) +density = symbolic.Parameter( + "electrolyte.density", + domain, + default=1399.8 * syu.kg / syu.m**3, + unit=syu.kg / syu.m**3, + doc="Density of the electrolyte", +) +inlet_velocity = symbolic.Parameter( + "flow.inlet_velocity", + domain, + default=1e-7 * syu.m / syu.second, + unit=syu.m / syu.second, + doc="Velocity at the inlet", +) +# - + +# ## The flow +# +# The velocity is of quadratic and the pressure of linear elements, the +# Taylor-Hood pair. The inlet velocity scales the velocity, and the pressure +# drop of Darcy flow at it across the height of the electrode, evaluated at +# the defaults of the parameters with +# {py:func}`value_at_defaults `, scales +# the pressure: + +# + +height = 15e-3 +pressure_scale = symbolic.value_at_defaults( + viscosity * inlet_velocity * height / electrode.permeability() +) +u = symbolic.Field( + "u", + domain, + ("Lagrange", 2, (3,)), + unit=syu.m / syu.second, + scale=symbolic.value_at_defaults(inlet_velocity), +) +p = symbolic.Field("p", domain, ("Lagrange", 1), unit=syu.pascal, scale=pressure_scale) +# - + +# The flow of {eq}`brinkman-forchheimer` is a sum of +# {py:class}`BrinkmanForchheimer ` in the +# electrode and {py:class}`NavierStokes ` in the +# channels, with grad-div stabilization, +# {py:class}`GradDiv `, whose weight RFBniCS takes as +# a small part of the viscosity. Its creeping-flow part, +# {py:class}`Brinkman ` and +# {py:class}`Stokes `, is linear, and its solution the +# start of Newton's method for the whole flow. Both are checked for units +# and solved in dimensionless form: + +# + +grad_div = GradDiv("grad_div", u, 1e-3 * viscosity) +creeping = ( + Brinkman( + "electrode_creeping_flow", + u, + p, + electrode.porosity, + electrode.permeability(), + viscosity, + markers=electrode.markers, + ) + + Stokes("channel_creeping_flow", u, p, viscosity, markers=channels) + + grad_div +) +flow = ( + BrinkmanForchheimer( + "electrode_flow", + u, + p, + electrode.porosity, + electrode.permeability(), + electrode.forchheimer(), + viscosity, + density, + markers=electrode.markers, + ) + + NavierStokes("channel_flow", u, p, viscosity, density, markers=channels) + + grad_div +) +flow_systems = [ + symbolic.build_system(model, check_units=True, factorize_equations=True) + for model in (creeping, flow) +] +flow_bcs = [ + symbolic.DirichletBC(u, "Inlet", ufl.as_vector([inlet_velocity, 0.0, 0.0])), + symbolic.DirichletBC(u, ("Walls", "CurrentCollector", "Membrane"), 0.0), +] +# - + +# ## The half-cell +# +# The ions, the reservoir and the half-cell are those of +# {doc}`vanadium_half_cell_2d`. The electrolyte fills the `channels` too. +# The membrane is larger than the current collector, which the channels +# partly cover, so the cell's current density, given at the membrane, is +# scaled at the collector by the ratio of their areas: + +# + +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("H", 1, 9.312e-9 * diffusivity), + Ion("HSO4", -1, 1.33e-9 * diffusivity), + Ion("SO4", -2, 1.065e-9 * diffusivity), +) +reservoir = { + "V2": 300.0 * concentration, + "V3": 300.0 * concentration, + "H": 1000.0 * concentration, + "HSO4": 1000.0 * concentration, +} +cell = RedoxHalfCell( + domain, + ions, + reservoir, + reduced="V2", + oxidized="V3", + electrons=1, + electrode=electrode, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + channels=channels, + velocity=u, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, +) +inlet = cell.inlet_conditions("Inlet") +# - + +# ## Solving the flow +# +# Everything is discretized together, and the creeping flow and then the +# whole flow are solved: + +# + +disc = discrete.discretize( + [*flow_systems, *flow_bcs, cell.model, *inlet], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +for name, system in zip(("Creeping flow", "Flow"), flow_systems, strict=True): + solver = solvers.Solver(disc[system], flow_bcs) + solver.solve() + print(f"{name}: {solver.iterations} Newton iterations") +# - + +# The flow is so slow that its inertia, and the Forchheimer drag, change it +# by less than Newton's tolerance: the creeping flow is already the flow. It +# is incompressible, so what enters at the inlet leaves at the outlet. The +# velocity of Taylor-Hood elements is not divergence-free cell by cell, but +# the constants are among the test functions of the pressure, so it +# conserves mass over the whole cell, and the outflow matches the inflow to +# round-off: + +# + +n = ufl.FacetNormal(domain) +inflow = -disc.assemble_scalar(ufl.dot(u, n) * domain.ds("Inlet")) +outflow = disc.assemble_scalar(ufl.dot(u, n) * domain.ds("Outlet")) +print(f"Outflow over inflow: {outflow / inflow:.8f}") +assert inflow > 0 and np.isclose(outflow, inflow, rtol=1e-8) +# - + +# ## Charge, rest and discharge +# +# The cell starts at rest, with the current balanced between the membrane +# and the collector, and follows the schedule of +# {doc}`vanadium_half_cell_2d`. As there, the simulation's model is the +# half-cell alone, and every step takes the velocity solved above as known: + +# + +cell.balance_currents(disc) +cell.initialize(disc) +schedule = ( + ("charge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), + ("rest", 0.0, ((5, 0.1), (5, 1.0), (10, 10.0))), + ("discharge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), +) +sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 100}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +# - + +# The simulation runs in parts, as in {doc}`vanadium_half_cell_2d`, and keeps +# the local state of charge at the end of the charge: + +# + +local = cell.local_state_of_charge() +space = disc.function_space(cell.concentrations["V2"]) +charged = dolfinx.fem.Function(space, name="soc") +records: list[tuple[float, str, float, float]] = [] +for mode, current, steps in schedule: + disc.set({cell.current: sign[mode] * current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run( + num_steps=count, + callback=lambda t, step, mode=mode: records.append( + (t, mode, cell.voltage(disc), cell.state_of_charge(disc)) + ), + ) + if mode == "charge": + charged.interpolate( + dolfinx.fem.Expression(disc[local], space.element.interpolation_points) + ) +times = np.array([r[0] for r in records]) +modes = np.array([r[1] for r in records]) +voltage = np.array([r[2] for r in records]) +charge = np.array([r[3] for r in records]) +# - + +# Charging raises the state of charge and discharging lowers it, and the +# voltage is lower on charge than on discharge: + +# + +steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} +print(f"State of charge after charge: {charge[modes == 'charge'][-1]:.6f}") +print(f"State of charge at the end: {charge[-1]:.6f}") +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() +# - + +# The voltage and the state of charge through the schedule: + +fig, (top, bottom) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + top.plot(times[modes == mode], 1e3 * voltage[modes == mode], ".", label=mode) + bottom.plot(times[modes == mode], charge[modes == mode], ".", label=mode) +top.set_ylabel("Voltage (mV)") +bottom.set_ylabel("State of charge") +bottom.set_xlabel("Time (s)") +top.legend() +plt.show() + +# The local state of charge at the end of the charge, on the membrane side +# of the electrode, in millimetres: + +grid = gather_grid(charged) +if MPI.COMM_WORLD.rank == 0: + grid.points[:] *= 1e3 + plotter = pyvista.Plotter(window_size=[640, 640]) + plotter.add_mesh(grid, scalars="soc", scalar_bar_args={"title": "SOC"}) + plotter.view_yz(negative=True) + plotter.show() + +# ## References +# +# ```{footbibliography} +# ``` diff --git a/references.bib b/references.bib index 1b06fae..4d7a89e 100644 --- a/references.bib +++ b/references.bib @@ -257,3 +257,72 @@ @book{CarslawJaeger1959 address = {Oxford}, year = {1959}, } + +@article{Barzegari2026, + author = {Barzegari, Mojtaba and de Waal, Martin and Bordignon, Davide and Rosinaite, Adele and Liu, Baichen and Forner-Cuenca, Antoni}, + title = {{RfbFoam}: An open-source framework for porous electrode modeling in redox flow batteries}, + journal = {ChemRxiv}, + year = {2026}, + doi = {10.26434/chemrxiv.10001658/v1}, +} + +@book{NewmanBalsara2021, + author = {Newman, John and Balsara, Nitash P.}, + title = {Electrochemical Systems}, + edition = {4}, + publisher = {Wiley}, + address = {Hoboken}, + year = {2021}, +} + +@book{Fogler2020, + author = {Fogler, H. Scott}, + title = {Elements of Chemical Reaction Engineering}, + edition = {6}, + publisher = {Pearson}, + address = {Boston}, + year = {2020}, +} + +@article{Thiele1939, + author = {Thiele, E. W.}, + title = {Relation between catalytic activity and size of particle}, + journal = {Industrial \& Engineering Chemistry}, + volume = {31}, + number = {7}, + pages = {916--920}, + year = {1939}, + doi = {10.1021/ie50355a027}, +} + +@misc{Aghabarari2026, + author = {Aghabarari, Amirhossein and Dokken, J{\o}rgen S. and Horsch, Martin T. and Valseth, Eirik and Janssen, Mathijs}, + title = {{RFBniCS}: An open-source simulation framework for redox flow batteries}, + year = {2026}, + eprint = {2610.01840}, + archivePrefix = {arXiv}, + primaryClass = {physics.comp-ph}, + url = {https://arxiv.org/abs/2610.01840}, +} + +@article{MunozPerales2022, + author = {Mu{\~n}oz-Perales, Vanesa and Berling, Sabrina and Garc{\'i}a-Quismondo, Enrique and Garc{\'i}a-Salaberri, Pablo A. and Palma, Jes{\'u}s and Vera, Marcos and Ib{\'a}{\~n}ez, Santiago E.}, + title = {Investigating the Effects of Operation Variables on {All-Vanadium} Redox Flow Batteries Through an Advanced Unit-Cell Model}, + journal = {Journal of The Electrochemical Society}, + volume = {169}, + number = {10}, + pages = {100522}, + year = {2022}, + doi = {10.1149/1945-7111/ac975f}, +} + +@article{Krishnamurthy2011, + author = {Krishnamurthy, Deepak and Johansson, Erik O. and Lee, Jin Wook and Kjeang, Erik}, + title = {Computational modeling of microfluidic fuel cells with flow-through porous electrodes}, + journal = {Journal of Power Sources}, + volume = {196}, + number = {23}, + pages = {10019--10031}, + year = {2011}, + doi = {10.1016/j.jpowsour.2011.08.024}, +} diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index 116ffd4..f06b9d0 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -93,7 +93,8 @@ class Electrolyte: Raises: ValueError: If there are ions but none has a field, more than one - ion has none, or two ions share a name. + ion has none, the ion without one has no charge, or two ions + share a name. """ _name: str @@ -145,6 +146,11 @@ def __init__( "At most one ion can follow from electroneutrality, got " f"{[ion.name for ion in balanced]}." ) + if balanced and balanced[0].charge == 0: + raise ValueError( + f"Ion {balanced[0].name!r} has no charge, so electroneutrality " + "cannot give its concentration. Give it a field." + ) if balanced and len(balanced) == len(ions): raise ValueError( "Electroneutrality needs at least one transported ion, got only " diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py index e15a1bd..5e42a80 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -40,6 +40,10 @@ class RedoxHalfCell: """A half-cell of a redox flow battery, of a couple of dissolved ions. + The macro-homogeneous porous-electrode model of RFBniCS, `Aghabarari et + al., 2026 `_, Section II, for any + couple, electrolyte and sign of the half-cell. + The unknowns are the potentials ``phi_s`` of the solid and ``phi_l`` of the electrolyte, in V, and the concentration ``c_`` of every ion with a reservoir concentration, in mol/m^3. The one ion without follows @@ -90,8 +94,11 @@ class RedoxHalfCell: reservoir: The concentration of each transported ion in the reservoir, in mol/m^3 or a quantity with units, the default of the parameter ``".c0_"``. - reduced: Name of the reduced ion of the couple. - oxidized: Name of the oxidized ion of the couple. + reduced: Name of the reduced ion of the couple. With ``oxidized`` + and ``electrons`` it gives the reaction + ``oxidized + n e- <-> reduced``, e.g. V3 + e- <-> V2. + oxidized: Name of the oxidized ion of the couple, of ``electrons`` + more charge than the reduced one. electrode: The porous electrode. oxidation_on_discharge: Whether discharge oxidizes the couple, as in the negative half-cell. @@ -124,9 +131,11 @@ class RedoxHalfCell: name: Prefix of the cell's own parameters. Raises: - ValueError: For mass transfer or an outlet without a velocity, or a + ValueError: For mass transfer or an outlet without a velocity, a reservoir concentration of an ion that the electrolyte does not - have. + have, a reservoir that lacks not exactly one ion, an ion of the + couple without one, a couple that does not conserve charge, or a + balanced ion without charge. """ _domain: Domain @@ -219,6 +228,29 @@ def __init__( unknown = sorted(set(reservoir) - {ion.name for ion in ions}) if unknown: raise ValueError(f"The electrolyte has no ions {unknown}.") + # The one ion without a reservoir concentration has no field, and + # follows from electroneutrality, which closes the equations + balanced = [ion.name for ion in ions if ion.name not in reservoir] + if len(balanced) != 1: + raise ValueError( + "The reservoir must give the concentrations of all ions but one, " + "whose concentration follows from electroneutrality; it lacks " + f"{balanced}." + ) + # The couple is the reaction oxidized + n e- <-> reduced, which must + # conserve charge + charges = {ion.name: ion.charge for ion in ions} + for ion in (reduced, oxidized): + if ion not in reservoir: + raise ValueError( + f"The ion {ion!r} of the couple needs a reservoir concentration." + ) + if charges[oxidized] - charges[reduced] != electrons: + raise ValueError( + f"{oxidized} + {electrons} e- <-> {reduced} does not conserve " + f"charge: the charges are {charges[oxidized]} and " + f"{charges[reduced]}." + ) self._domain = domain self._electrode = electrode self._reduced, self._oxidized = reduced, oxidized diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index cf1400b..ed09660 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -442,7 +442,11 @@ class MassTransferButlerVolmer(ButlerVolmer): ``k_m (c - c_s)`` from the bulk. The mass-transfer coefficient of each ion is ``k_m = Sh D / d`` with the Sherwood number ``Sh = a Re^b`` and the Reynolds number ``Re = rho |u| d / mu`` of the flow past fibers of - diameter ``d``. ``D`` is the diffusivity in free electrolyte. + diameter ``d``. ``D`` is the diffusivity in free electrolyte. The + defaults ``a = 7`` and ``b = 0.4`` are the correlation of `Krishnamurthy + et al., 2011 `_, and the + surface concentrations those of equation 20 of `Aghabarari et al., 2026 + `_. Args: name: Name of the model. diff --git a/src/sumfem/models/porous.py b/src/sumfem/models/porous.py index 510a2b7..45f0e75 100644 --- a/src/sumfem/models/porous.py +++ b/src/sumfem/models/porous.py @@ -22,7 +22,11 @@ def bruggeman( It corrects a transport property of a phase for the volume the other phase takes and for the tortuosity of the paths through it. The exponent ``3/2`` is that of spherical obstructions, see `Tjaden et al., 2016 - `_. + `_, where the tortuosity + scales with the fraction to the power ``-1/2``. It suits macroporous + electrodes; where the tortuosity is measured, the property + ``value * fraction / tortuosity`` is the alternative of `Barzegari et + al., 2026 `_. Args: value: The property of the phase, e.g. a diffusivity or the @@ -43,7 +47,11 @@ def kozeny_carman( An estimate of the permeability of a fibrous porous electrode, such as a carbon felt, from the diameter of its fibres, with an empirical constant - ``K``. + ``K``. The carbon fibres of commercial cloths, papers and felts have + diameters of 7 to 9 micrometres, see `Barzegari et al., 2026 + `_. The estimate is of + limited use for complex media, whose permeability is best fitted to + measured pressure drops. Args: fiber_diameter: Fiber diameter ``d``. diff --git a/src/sumfem/symbolic/terminals.py b/src/sumfem/symbolic/terminals.py index 2e2c1ae..534bddf 100644 --- a/src/sumfem/symbolic/terminals.py +++ b/src/sumfem/symbolic/terminals.py @@ -39,6 +39,7 @@ from ..error_handling import ComponentError, check_name, qualified_name from .domains import Domain, domain_of from .elements import ElementSpec, as_element +from .expressions import evaluate_constant from .si import Unit, has_units, si_factor, to_si _PIOLA = ( @@ -58,6 +59,7 @@ "Parameter", "TimeDerivative", "field_of", + "value_at_defaults", ] # Map from the part number of a field's test and trial functions to the @@ -493,6 +495,9 @@ def __init__( if has_units(default): # A quantity with its own units, as a value in the declared unit default = to_si(default, unit) / si_factor(unit) + elif isinstance(default, sympy.Expr) and default.is_number: + # A zero quantity, e.g. 0 * syu.pascal, which sympy reduces to 0 + default = float(default) self._default = ( None if default is None else _validate_value(name, default, shape) ) @@ -531,6 +536,51 @@ def doc(self) -> str: return self._doc +def value_at_defaults(expr: ufl.core.expr.Expr) -> float | npt.NDArray[np.float64]: + """The value of an expression of parameters at their defaults, in SI units. + + Use it for a value that follows from the parameters of a model before + anything is discretized, e.g. the permeability of an electrode or the + pressure drop of a flow, as the scale of a field. + + Args: + expr: An expression of parameters and numbers. + + Raises: + ComponentError: If the expression depends on anything but parameters + with defaults, e.g. a field or the spatial coordinate. + + Example: + The scale of a pressure from a viscosity, a velocity, a length and a + permeability, all parameters:: + + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal, + scale=value_at_defaults(mu * v * L / K)) + """ + defaults: dict[ufl.Coefficient, ufl.core.expr.Expr] = {} + for c in ufl.algorithms.extract_coefficients(expr): + default = c.default_si if isinstance(c, Parameter) else None + if default is None: + name = getattr(c, "name", str(c)) + raise ComponentError( + f"The expression depends on {name!r}, which is not a parameter " + "with a default." + ) + defaults[c] = ( + ufl.as_ufl(float(default)) + if default.ndim == 0 + else (ufl.as_tensor(default.tolist())) + ) + replaced = ufl.replace(expr, defaults) if defaults else expr + try: + return evaluate_constant(replaced) + except TypeError as error: + raise ComponentError( + "The expression depends on more than parameters, e.g. on the " + "spatial coordinate." + ) from error + + class KnownCoefficient(ufl.Coefficient, QuantityMixin): """Named, spatially varying data that a model reads but does not solve for. diff --git a/tests/electrochemistry/half_cell/conftest.py b/tests/electrochemistry/half_cell/conftest.py index 65fc6ea..76f9e95 100644 --- a/tests/electrochemistry/half_cell/conftest.py +++ b/tests/electrochemistry/half_cell/conftest.py @@ -1,7 +1,8 @@ """Meshes and the original RFBniCS forms, for comparison with the models. The reference residuals are transcribed from the examples of RFBniCS -(https://github.com/Ah-Aghabarari/RFBniCS, MIT licence), with one change: +(https://github.com/Ah-Aghabarari/RFBniCS, MIT licence; Aghabarari et al., +2026, https://arxiv.org/abs/2610.01840), with one change: each unknown is its own function, where RFBniCS splits one blocked function. The tests pass the functions of a discretization, so both sides assemble on the same spaces. diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py index 929ff78..037a646 100644 --- a/tests/electrochemistry/half_cell/test_models.py +++ b/tests/electrochemistry/half_cell/test_models.py @@ -201,3 +201,42 @@ def test_post_processing(reference): points, values = evaluate_over_line(function, [0.0], [reference.W], 5) np.testing.assert_allclose(values[:, 0], 1.0 + points[:, 0] / reference.W) assert values.shape[1] == 1 + + +def test_the_couple_conserves_charge(reference): + """V3 + 2 e- <-> V2 does not conserve charge, nor V3 + e- <-> SO4.""" + domain = Domain("cell", "interval") + with pytest.raises(ValueError, match="does not conserve charge"): + reference.half_cell(domain, electrons=2) + with pytest.raises(ValueError, match="needs a reservoir"): + RedoxHalfCell( + domain, + reference.IONS, + reference.RESERVOIR, + "SO4", + "V3", + reference.electrode(domain), + oxidation_on_discharge=True, + rate_constant=2.5e-7, + formal_potential=-0.36, + ) + + +@pytest.mark.parametrize("missing", [(), ("H", "SO4")], ids=["none", "two"]) +def test_the_reservoir_lacks_exactly_one_ion(reference, missing): + """One ion, and only one, follows from electroneutrality.""" + domain = Domain("cell", "interval") + reservoir = {**reference.RESERVOIR, "SO4": 100.0} + reservoir = {k: v for k, v in reservoir.items() if k not in missing} + with pytest.raises(ValueError, match="all ions but one"): + RedoxHalfCell( + domain, + reference.IONS, + reservoir, + "V2", + "V3", + reference.electrode(domain), + oxidation_on_discharge=True, + rate_constant=2.5e-7, + formal_potential=-0.36, + ) diff --git a/tests/electrochemistry/test_electrolyte.py b/tests/electrochemistry/test_electrolyte.py index 7c9a790..0990f15 100644 --- a/tests/electrochemistry/test_electrolyte.py +++ b/tests/electrochemistry/test_electrolyte.py @@ -13,7 +13,7 @@ NernstPlanck, PorousElectrode, ) -from sumfem.symbolic import Field, build_system +from sumfem.symbolic import Domain, Field, build_system CONCENTRATION = syu.mol / syu.m**3 @@ -152,3 +152,10 @@ def test_a_vanadium_electrode_has_consistent_units(domain): ) system = build_system(model, check_units=True) assert {f.name for f in system.unknowns} >= {"phi_s", "phi_l", "c2", "c3", "h"} + + +def test_the_balanced_ion_needs_a_charge(): + domain = Domain("cell", "interval") + c = Field("c", domain, ("Lagrange", 1)) + with pytest.raises(ValueError, match="has no charge"): + Electrolyte("e", domain, [Ion("A", 1, 1e-9, c), Ion("N", 0, 1e-9)]) diff --git a/tests/symbolic/test_terminals.py b/tests/symbolic/test_terminals.py index bbe7080..ee4cb3e 100644 --- a/tests/symbolic/test_terminals.py +++ b/tests/symbolic/test_terminals.py @@ -11,6 +11,7 @@ KnownCoefficient, Parameter, as_expr, + value_at_defaults, ) from sumfem.symbolic.elements import as_element from sumfem.symbolic.terminals import FieldRole, FieldTerminal, field_of @@ -235,3 +236,30 @@ def test_as_expr(domain): def test_parameter_names_may_be_qualified(domain): assert Parameter("solid.sigma", domain).name == "solid.sigma" + + +def test_a_zero_quantity_is_a_default(domain): + """Sympy reduces a zero quantity to the number zero, without units.""" + parameter = Parameter("p", domain, default=0.0 * syu.pascal, unit=syu.pascal) + assert parameter.default_si == 0.0 + + +def test_an_expression_of_parameters_has_a_value_at_their_defaults(domain): + mu = Parameter("mu", domain, default=2.0) + v = Parameter("v", domain, default=[1.0, 3.0], shape=(2,)) + assert value_at_defaults(mu * ufl.sqrt(ufl.dot(v, v))) == pytest.approx( + 2.0 * np.sqrt(10.0) + ) + np.testing.assert_allclose(value_at_defaults(mu * v), [2.0, 6.0]) + + +def test_a_value_at_defaults_needs_parameters_with_defaults(domain): + free = Parameter("free", domain) + with pytest.raises(ComponentError, match="'free'"): + value_at_defaults(2.0 * free) + u = Field("u", domain, ("Lagrange", 1)) + with pytest.raises(ComponentError, match="'u'"): + value_at_defaults(u + 1.0) + x = ufl.SpatialCoordinate(domain) + with pytest.raises(ComponentError, match="spatial coordinate"): + value_at_defaults(x[0] * Parameter("a", domain, default=1.0)) From 2f1475f9e3e4f26683bec82a24d440ca3c5af75f Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 22:31:46 +0000 Subject: [PATCH 4/6] Add missing import in expose --- src/sumfem/symbolic/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/sumfem/symbolic/__init__.py b/src/sumfem/symbolic/__init__.py index 41af21f..19a13a6 100644 --- a/src/sumfem/symbolic/__init__.py +++ b/src/sumfem/symbolic/__init__.py @@ -48,7 +48,7 @@ from .domains import Domain, Subdomain from .expressions import as_expr from .system import System, build_system -from .terminals import Field, KnownCoefficient, Parameter +from .terminals import Field, KnownCoefficient, Parameter, value_at_defaults from .time import Clock, Dt from .transfer import Assignment, Interpolation, Transfer, Transferable from .ufl_forms import UFLResidual @@ -105,4 +105,5 @@ "equation_dimensions", "equation_scales", "pi_groups", -] + "value_at_defaults", +] \ No newline at end of file From a854cba621a5bad6d6368328004bb6edae24b0ff Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 22:34:48 +0000 Subject: [PATCH 5/6] Add newlien --- src/sumfem/symbolic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sumfem/symbolic/__init__.py b/src/sumfem/symbolic/__init__.py index 19a13a6..2883bb5 100644 --- a/src/sumfem/symbolic/__init__.py +++ b/src/sumfem/symbolic/__init__.py @@ -106,4 +106,4 @@ "equation_scales", "pi_groups", "value_at_defaults", -] \ No newline at end of file +] From 2b59945740e811d0b8edfd246af11d5d0beae7c6 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 23:02:43 +0000 Subject: [PATCH 6/6] Check the discharge-direction test over all processes With 4 cells on 2 processes, each process holds only one of the collector and membrane facets, so one equation's local residual is zero and the test failed in parallel. The residuals now have their ghost contributions added to their owners, the owned entries are compared, and the maximum and the mismatch are reduced over all processes, so every process reaches the same verdict. Co-Authored-By: Claude Opus 5.5 --- .../electrochemistry/half_cell/test_models.py | 27 +++++++++++++++---- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py index 037a646..045f0f9 100644 --- a/tests/electrochemistry/half_cell/test_models.py +++ b/tests/electrochemistry/half_cell/test_models.py @@ -1,6 +1,9 @@ """The RFBniCS models: units, structure, schedules and post-processing.""" +from mpi4py import MPI + import dolfinx.fem +import dolfinx.la import numpy as np import pytest import sympy.physics.units as syu @@ -166,12 +169,26 @@ def test_the_couple_is_oxidized_or_reduced_on_discharge(reference): # the equations of the potentials disc.set({cell.current: 0.5}) problem = disc[system] - residuals.append( - [dolfinx.fem.assemble_vector(form).array for form in problem.residual[:2]] - ) + residuals.append([owned_residual(form) for form in problem.residual[:2]]) + # Each process may hold only one of the boundaries, so the checks are + # over all processes, which all reach the same verdict + comm = marked.mesh.comm for oxidizing, reducing in zip(*residuals, strict=True): - assert np.abs(oxidizing).max() > 0 - np.testing.assert_allclose(oxidizing, -reducing) + largest = comm.allreduce(np.abs(oxidizing).max(initial=0.0), op=MPI.MAX) + mismatch = comm.allreduce( + np.abs(oxidizing + reducing).max(initial=0.0), op=MPI.MAX + ) + assert largest > 0 + assert mismatch <= 1e-12 * largest + + +def owned_residual(form: dolfinx.fem.Form) -> np.ndarray: + """A residual vector with its ghost contributions added to their owners, + restricted to the entries this process owns.""" + vector = dolfinx.fem.assemble_vector(form) + vector.scatter_reverse(dolfinx.la.InsertMode.add) + size = vector.index_map.size_local * vector.block_size + return vector.array[:size].copy() def test_the_reservoir_names_ions_of_the_electrolyte(reference):