#5198 update definition of hysteresis decay rate

CHANGELOG.md

@@ -23,7 +23,7 @@
 - Fix Bruggeman coefficient computation from BPX porosity and transport efficiency instead of hard-coding, remove redundant values, and add a unit test for verification. ([#5196](https://github.com/pybamm-team/PyBaMM/pull/5196))
 
 ## Breaking changes
-
+- Updates the hysteresis decay rate parameters to a "true" hysteresis decay rate which changes the interpretation of the units of the hysteresis decay rate parameters. ([#5217](https://github.com/pybamm-team/PyBaMM/pull/5217))
 - Changed fundamental variable for all SEI models from thickness to concentration ([#4869](https://github.com/pybamm-team/PyBaMM/pull/4869))
 
 # [v25.8.0](https://github.com/pybamm-team/PyBaMM/tree/v25.8.0) - 2025-08-04
@@ -65,6 +65,7 @@
 ## Breaking changes
 
 - Changed behavior of drive cycle steps in `pybamm.Experiment`s to treat each time point as a discontinuity, consistent with how input interpolants work. This ensures more accurate simulation of drive cycles with rapid changes. ([#5141](https://github.com/pybamm-team/PyBaMM/pull/5141))
+- Makes `A_cc = L_z * L_y * number of layers`, which in turn changes the interpretation of "{domain} electrode capacity [A.h]" variables (and their composite equivalents). The electrode capacity variables now account for all layers, whereas before it was the capacity of a single electrode layer. ([#5138](https://github.com/pybamm-team/PyBaMM/pull/5138))
 - Removed the IREE code from the IDAKLU solver ([#5080](https://github.com/pybamm-team/PyBaMM/pull/5080))
 - Removed support for Python 3.9 ([#5052](https://github.com/pybamm-team/PyBaMM/pull/5052))
 - In OCP hysteresis models, users need to explicitly give the equilibrium, delithiation, and lithiation OCPs when using a hysteresis model. E.g., you must provide all three of "Negative electrode OCP [V]", "Negative electrode delithiation OCP [V]", and "Negative electrode lithiation OCP [V]". ([#4893](https://github.com/pybamm-team/PyBaMM/pull/4893))

docs/source/examples/notebooks/expression_tree/expression-tree.ipynb

@@ -263,7 +263,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pybamm",
+   "display_name": ".env",
    "language": "python",
    "name": "python3"
   },
@@ -277,7 +277,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.12.11"
   },
   "toc": {
    "base_numbering": 1,
@@ -291,11 +291,6 @@
    "toc_position": {},
    "toc_section_display": true,
    "toc_window_display": true
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c"
-   }
   }
  },
  "nbformat": 4,

docs/source/examples/notebooks/models/hysteresis-state-models.ipynb

@@ -9,15 +9,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
       "\n",
-      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.2.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m25.1.1\u001b[0m\n",
+      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m25.1.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m25.2\u001b[0m\n",
       "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
       "Note: you may need to restart the kernel to use updated packages.\n"
      ]
@@ -44,27 +44,25 @@
     "\n",
     "### One-state hysteresis (Axen)\n",
     "In this model, the state variable $h(z,t)$ is both stoichiometry and time-dependent, and is governed by the following ODE:\n",
-    "$$ \\frac{dh(z,t)}{dt} = \\gamma \\left( \\frac{i_{\\text{vol}}}{F \\, c_{\\text{max}} \\, \\epsilon}\\right) (1 - \\text{sgn}(i_{\\text{vol}}) \\, h(z,t)),$$\n",
+    "$$ \\frac{dh(z,t)}{dt} = \\gamma \\left( \\frac{i_{\\text{vol}}}{F \\, c_{\\text{max}} \\, \\epsilon}\\right) \\frac{\\left(1 - \\text{sgn}(i_{\\text{vol}}) \\, h(z,t)\\right)}{2},$$\n",
     "where $i_{\\text{vol}}$ is the volumetric interfacial current density, $F$ is the Faraday constant, $c_{\\text{max}}$ is the maximum lithium concentration in the particles and $\\epsilon$ is the active material volume fraction. In this model, the rate parameter $\\gamma$ is given by\n",
     "$$\\gamma = \\left( \\frac{1-\\text{sgn}(i_{\\text{vol}})}{2}  \\gamma_{\\text{lith}} + \\frac{1+\\text{sgn}(i_{\\text{vol}})}{2} \\gamma_{\\text{delith}}\\right),$$\n",
     "which allows for a different decay rate during delithiation ($\\gamma_{\\text{delith}}$) and lithiation ($\\gamma_{\\text{lith}}$).\n",
     "Note that this implementation is a reformulation of the governing equations in the original paper. Also note that $\\gamma$ is dimensionless.\n",
     "\n",
     "\n",
     "### One-state differential capacity hysteresis (Wycisk)\n",
-    "As in the Axen model, the state variable $h(z,t)$ is both stoichiometry and time-dependent, and is governed by the following ODE:\n",
-    "$$ \\frac{dh(z,t)}{dt} = \\left(\\frac{\\gamma \\cdot i_{\\text{surf}}}{Q_{\\text{cell}}}\\right) \\left(1 - \\text{sgn}(i_{\\text{surf}}) \\, h(z,t)\\right), $$\n",
+    "As in the Axen model, the state variable $h(z,t)$ is both stoichiometry and time-dependent, and is governed by the same ODE\n",
+    "$$ \\frac{dh(z,t)}{dt} = \\gamma \\left( \\frac{i_{\\text{vol}}}{F \\, c_{\\text{max}} \\, \\epsilon}\\right) \\frac{\\left(1 - \\text{sgn}(i_{\\text{vol}}) \\, h(z,t)\\right)}{2},$$\n",
     "\n",
-    "where $\\gamma(z)$ is given by\n",
+    "where now $\\gamma(z)$ is given by\n",
     "\n",
-    "$$ \\gamma(z) = \\Gamma(z) \\cdot \\frac{1}{\\left(C_{\\text{diff}}\\left(z\\right)\\right)^{x}}.$$\n",
+    "$$ \\gamma(z) = \\Gamma(z) \\cdot \\left(\\frac{C_{\\text{diff}}\\left(z\\right)}{C_{\\text{ref,vol}}}\\right)^{-x}.$$\n",
     "\n",
     "Here $C_{\\text{diff}}(z)$ is the differential capacity with respect to potential, expressed as \n",
     "\n",
     "$$ C_{\\text{diff}}(z) = \\frac{dQ}{dU_{\\text{eq}}(z)},$$\n",
-    "where $Q$ is the capacity of the phase of active material experiencing the voltage hysteresis and $U_{\\text{eq}}$ is the average open-circuit potential of the phase of active material experiencing the voltage hysteresis. The remaining parameters are $\\Gamma$ and $x$ which are both fitting parameters that affect the response of the hysteresis state decay when passing charge in either direction.\n",
-    "\n",
-    "Warning: in this implementation, the units of both $\\gamma$ and $\\Gamma$ are unusual. The parameter $\\gamma$ has units $\\text{m}^2$ and $\\Gamma$ has units $\\text{m}^2\\text{V}^{x}\\text{A}^{-x}\\text{h}^{-x}$. Also note that the units of $Q$ are $\\text{A}\\text{h}$, but the time $t$ is expressed in seconds.\n",
+    "where $Q$ is the capacity of the phase of active material experiencing the voltage hysteresis and $U_{\\text{eq}}$ is the average open-circuit potential of the phase of active material experiencing the voltage hysteresis. The remaining parameters are $\\Gamma$ and $x$ which are both fitting parameters that affect the response of the hysteresis state decay when passing charge in either direction. The parameter $C_{\\text{ref,vol}}$ is the reference differential capacity, which in this implementation is hard-coded to 1 $\\text{F}$.\n",
     "\n",
     "### Current Sigmoid (Ai)\n",
     "This approach uses a sigmoid function to switch between the delithiation and lithiation branches of the open-circuit potential, depending on the sign of the current. The original paper gives the following expression for the open-circuit potential:\n",
@@ -88,7 +86,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -129,7 +127,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -166,15 +164,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Wycisk parameters\n",
     "parameters_wycisk = parameters.copy()\n",
     "parameters_wycisk.update(\n",
     "    {\n",
-    "        \"Secondary: Negative particle hysteresis decay rate\": 100,\n",
+    "        \"Secondary: Negative particle hysteresis decay rate\": 0.1,\n",
     "        \"Secondary: Negative particle hysteresis switching factor\": 10,\n",
     "        \"Secondary: Initial hysteresis state in negative electrode\": 1,\n",
     "    },\n",
@@ -208,7 +206,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -234,18 +232,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 17,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/Users/robertwtimms/Documents/PyBaMM/src/pybamm/simulation.py:122: UserWarning: The default solver changed to IDAKLUSolver after the v25.4.0. release. You can swap back to the previous default by using `pybamm.CasadiSolver()` instead.\n",
-      "  self._solver = solver or self._model.default_solver\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "solutions = {}\n",
     "for name, model in models.items():\n",
@@ -264,18 +253,18 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "273b3b86ed5b4a7f8df79459f875414c",
+       "model_id": "5efa147c8cf44c7db96697dfdda4e610",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=3.0033810199725046, step=0.03003381019972504…"
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=3.041045468063589, step=0.03041045468063589)…"
       ]
      },
      "metadata": {},
@@ -284,10 +273,10 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.plotting.quick_plot.QuickPlot at 0x2a640f310>"
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x32a3cf2f0>"
       ]
      },
-     "execution_count": 7,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -323,7 +312,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "env",
+   "display_name": ".env",
    "language": "python",
    "name": "python3"
   },
@@ -337,7 +326,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.12.11"
   }
  },
  "nbformat": 4,

src/pybamm/models/full_battery_models/lithium_ion/util.py

@@ -1,35 +1,46 @@
+import pybamm
+
+
 def check_if_composite(options, electrode):
-    options = options or {}
-    particle_phases = options.get("particle phases", None)
-    if particle_phases is None:
-        return False
+    if not isinstance(options, pybamm.BatteryModelOptions):
+        options = pybamm.BatteryModelOptions(options)
+    domain_options = getattr(options, electrode)
+    particle_phases = domain_options["particle phases"]
     if particle_phases == "2":
         return True
-    if isinstance(particle_phases, tuple) and particle_phases[0] == "2":
-        if electrode == "negative":
-            return True
-    if isinstance(particle_phases, tuple) and particle_phases[1] == "2":
-        if electrode == "positive":
-            return True
-    return False
+    else:
+        return False
 
 
-def _has_hysteresis_composite(options, electrode, phase, hysteresis_options):
-    ocp_index = 0 if electrode == "negative" else 1
-    my_ocp_options = options["open-circuit potential"][ocp_index]
-    if check_if_composite(options, electrode) and isinstance(my_ocp_options, tuple):
+def _has_hysteresis(electrode, options, phase=None):
+    if not isinstance(options, pybamm.BatteryModelOptions):
+        options = pybamm.BatteryModelOptions(options)
+    hysteresis_options = [
+        "current sigmoid",
+        "one-state hysteresis",
+        "one-state differential capacity hysteresis",
+        # Also catch old names
+        "Axen",
+        "Wycisk",
+    ]
+    domain_options = getattr(options, electrode)
+    if isinstance(domain_options["open-circuit potential"], str):
+        return domain_options["open-circuit potential"] in hysteresis_options
+    elif isinstance(domain_options["open-circuit potential"], tuple):
+        ocp_option = domain_options["open-circuit potential"]
         if phase == "primary":
-            if my_ocp_options[0] in hysteresis_options:
-                return True
-            else:
-                return False
+            return ocp_option[0] in hysteresis_options
         elif phase == "secondary":
-            if my_ocp_options[1] in hysteresis_options:
-                return True
-            else:
-                return False
+            return ocp_option[1] in hysteresis_options
         else:
-            raise ValueError("Phase must be either primary or secondary")
+            return any(
+                isinstance(item, str) and item in hysteresis_options
+                for item in ocp_option
+            )
+    else:
+        raise ValueError(
+            f"Invalid open-circuit potential option: {domain_options['open-circuit potential']}"
+        )
 
 
 def _get_lithiation_delithiation(direction, electrode, options, phase=None):
@@ -44,24 +55,4 @@ def _get_lithiation_delithiation(direction, electrode, options, phase=None):
     ):
         return "delithiation"
     else:
-        raise ValueError
-
-
-def _has_hysteresis(electrode, options, phase=None):
-    hysteresis_options = [
-        "current sigmoid",
-        "one-state hysteresis",
-        "one-state differential capacity hysteresis",
-    ]
-    phase = phase or ""
-    if options.get("open-circuit potential") is None:
-        return False
-    if isinstance(options["open-circuit potential"], str):
-        if options["open-circuit potential"] in hysteresis_options:
-            return True
-        else:
-            return False
-    elif isinstance(options["open-circuit potential"], tuple):
-        return _has_hysteresis_composite(options, electrode, phase, hysteresis_options)
-    else:
-        raise ValueError
+        raise ValueError()

src/pybamm/models/submodels/interface/open_circuit_potential/one_state_differential_capacity_hysteresis_ocp.py

@@ -1,7 +1,15 @@
+import warnings
+
 import pybamm
 
 from . import BaseHysteresisOpenCircuitPotential
 
+# Only throw the hysteresis decay rate warning once
+warnings.filterwarnings(
+    "once",
+    message="The definition of the hysteresis decay rate parameter has changed*",
+)
+
 
 class OneStateDifferentialCapacityHysteresisOpenCircuitPotential(
     BaseHysteresisOpenCircuitPotential
@@ -31,6 +39,11 @@ class OneStateDifferentialCapacityHysteresisOpenCircuitPotential(
     def __init__(
         self, param, domain, reaction, options, phase="primary", x_average=False
     ):
+        warnings.warn(
+            "The definition of the hysteresis decay rate parameter has changed in "
+            "PyBaMM v25.10. Please see the CHANGELOG for more details.",
+            stacklevel=2,
+        )
         super().__init__(
             param, domain, reaction, options=options, phase=phase, x_average=x_average
         )
@@ -72,9 +85,11 @@ def set_rhs(self, variables):
         domain, Domain = self.domain_Domain
         phase_name = self.phase_name
 
+        c_max = self.phase_param.c_max
+        eps = self.phase_param.epsilon_s
         if self.x_average is False:
-            i_surf = variables[
-                f"{Domain} electrode {phase_name}interfacial current density [A.m-2]"
+            i_vol = variables[
+                f"{Domain} electrode {phase_name}volumetric interfacial current density [A.m-3]"
             ]
             sto_surf = variables[f"{Domain} {phase_name}particle surface stoichiometry"]
             T = variables[f"{Domain} electrode temperature [K]"]
@@ -83,8 +98,8 @@ def set_rhs(self, variables):
                 f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
             ]
         else:
-            i_surf = variables[
-                f"X-averaged {domain} electrode {phase_name}interfacial current density [A.m-2]"
+            i_vol = variables[
+                f"X-averaged {domain} electrode {phase_name}volumetric interfacial current density [A.m-3]"
             ]
             sto_surf = variables[
                 f"X-averaged {domain} {phase_name}particle surface stoichiometry"
@@ -94,18 +109,16 @@ def set_rhs(self, variables):
             dQdU = variables[
                 f"X-averaged {domain} electrode {phase_name}differential capacity [A.s.V-1]"
             ]
-        # check if composite or not
-        if phase_name != "":
-            Q_cell = variables[f"{Domain} electrode {phase_name}phase capacity [A.h]"]
-        else:
-            Q_cell = variables[f"{Domain} electrode capacity [A.h]"]
+            eps = pybamm.x_average(eps)
 
         Gamma = self.phase_param.hysteresis_decay(sto_surf, T)
         x = self.phase_param.hysteresis_switch
 
-        i_surf_sign = pybamm.sign(i_surf)
-        signed_h = 1 - i_surf_sign * h
-        gamma = Gamma * (1 / dQdU**x)
-        dhdt = gamma * i_surf / Q_cell * signed_h
+        i_vol_sign = pybamm.sign(i_vol)
+        signed_h = (1 - i_vol_sign * h) / 2
+        C_diff = dQdU
+        C_ref_vol = 1  # [F]
+        gamma = Gamma * ((C_diff / C_ref_vol) ** (-x))
+        dhdt = gamma * i_vol / (self.param.F * c_max * eps) * signed_h
 
         self.rhs[h] = dhdt