From 875b5c6c43da550b50df23ee22f052996bea9391 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Fri, 12 Jun 2026 10:19:38 -0400
Subject: [PATCH 1/6] Fix RegulariseSqrtAndPower corrupting state-independent
 bases (#5600)

* Fix RegulariseSqrtAndPower corrupting state-independent bases

RegulariseSqrtAndPower replaced every Sqrt/Power node with RegPower,
including bases independent of the state such as pybamm.Parameter and
pybamm.InputParameter rate constants. With the default scale of 1 this
corrupts values below the regularisation tolerance (e.g. a ~1e-9 rate
constant raised to 0.5 evaluated ~470x too small), causing solver
failures or silently wrong overpotentials in user-supplied
exchange-current density functions.

Now only regularise bases that match a registered scale or depend on
the state (variable/state vector/time). State-dependent normalised
bases like (c_e / c_e_ref) ** (1 - alpha) used by built-in sets
(ORegan2022, Chayambuka2022) keep their default-scale regularisation,
and manually-created reg_power nodes (MSMR j0_j) are untouched.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Add PR number to changelog entry

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
(cherry picked from commit 5946b8cfce30d6df9764aea6274a60d2412f2d3d)
---
 CHANGELOG.md                                  |   4 +
 .../expression_tree/operations/regularise.py  |  33 ++--
 .../test_operations/test_regularise.py        | 174 ++++++++++++++++--
 3 files changed, 188 insertions(+), 23 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a080f18316..3e93ab3adb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Bug fixes
+
+- `RegulariseSqrtAndPower` no longer regularises state-independent bases, fixing corrupted small rate constants in exchange-current density functions. ([#5600](https://github.com/pybamm-team/PyBaMM/pull/5600))
+
 # [v26.6.1.1](https://github.com/pybamm-team/PyBaMM/tree/v26.6.1.1) - 2026-06-11
 
 ## Bug fixes
diff --git a/src/pybamm/expression_tree/operations/regularise.py b/src/pybamm/expression_tree/operations/regularise.py
index 98bc7b239c..aeac2621df 100644
--- a/src/pybamm/expression_tree/operations/regularise.py
+++ b/src/pybamm/expression_tree/operations/regularise.py
@@ -11,9 +11,15 @@ class RegulariseSqrtAndPower:
     """
     Callable that replaces Sqrt and Power nodes with RegPower nodes.
 
-    All Sqrt and Power nodes are replaced with RegPower. If the base of the
-    operation matches a symbol in the scales map, that scale is used;
-    otherwise the default scale (None) is used.
+    A Sqrt or Power node is replaced with RegPower when the base of the
+    operation matches a symbol in the scales map (using that scale) or when
+    the base is state-dependent, i.e. contains a variable, state vector or
+    time (using the default scale of 1). Bases that are independent of the
+    state (e.g. ``pybamm.Parameter`` or ``pybamm.InputParameter`` rate
+    constants) are left unchanged: they have no derivative to regularise, and
+    RegPower with the default scale of 1 corrupts values smaller than the
+    regularisation tolerance delta (e.g. a rate constant of ~1e-9 raised to
+    the power 0.5 evaluates ~470x too small).
 
     Parameters
     ----------
@@ -100,20 +106,25 @@ def _process(self, sym, resolved_scales):
 
         new_children = [self._process(child, resolved_scales) for child in sym.children]
 
-        if isinstance(sym, pybamm.Sqrt):
-            child = new_children[0]
-            scale = self._get_scale(child, resolved_scales)
-            return pybamm.RegPower(child, 0.5, scale=scale)
-
-        if isinstance(sym, pybamm.Power):
-            base, exponent = new_children
+        if isinstance(sym, pybamm.Sqrt | pybamm.Power):
+            base = new_children[0]
+            exponent = 0.5 if isinstance(sym, pybamm.Sqrt) else new_children[1]
             scale = self._get_scale(base, resolved_scales)
-            return pybamm.RegPower(base, exponent, scale=scale)
+            if scale is not None or self._depends_on_state(base):
+                return pybamm.RegPower(base, exponent, scale=scale)
 
         if any(n is not o for n, o in zip(new_children, sym.children, strict=True)):
             return sym.create_copy(new_children=new_children)
         return sym
 
+    @staticmethod
+    def _depends_on_state(expr):
+        """Whether the expression depends on the state (or time)."""
+        return any(
+            isinstance(node, pybamm.VariableBase | pybamm.StateVectorBase | pybamm.Time)
+            for node in expr.pre_order()
+        )
+
     def _get_scale(self, expr, resolved_scales):
         """Get scale for an expression, defaulting to None.
 
diff --git a/tests/unit/test_expression_tree/test_operations/test_regularise.py b/tests/unit/test_expression_tree/test_operations/test_regularise.py
index c17ea9555e..12b20ae568 100644
--- a/tests/unit/test_expression_tree/test_operations/test_regularise.py
+++ b/tests/unit/test_expression_tree/test_operations/test_regularise.py
@@ -2,6 +2,8 @@
 # Tests for the RegulariseSqrtAndPower class
 #
 
+import pytest
+
 import pybamm
 
 
@@ -40,16 +42,16 @@ def test_basic_power_replacement(self):
         expr = c_e**0.5
         result = regulariser(expr, inputs=inputs)
 
-        # The original c_e**0.5 Power node should be replaced
-        # (result will have different Power nodes from reg_power formula)
-        assert result != expr
+        # Matched base: RegPower with the registered scale
+        assert isinstance(result, pybamm.RegPower)
+        assert result.children[2] == pybamm.Scalar(1000.0)
 
-    def test_scale_default_to_one(self):
-        """Test that unrecognized expressions get scale=None."""
+    def test_state_dependent_base_default_scale(self):
+        """State-dependent base with no registered scale: regularised with
+        the default scale of 1."""
         c_e = pybamm.Variable("c_e")
         c_s = pybamm.Variable("c_s")
 
-        # Only c_e has a scale, c_s should get scale=None
         inputs = {"c_e": c_e, "c_s": c_s}
         regulariser = pybamm.RegulariseSqrtAndPower(
             {c_e: pybamm.Scalar(1000.0)},
@@ -59,14 +61,162 @@ def test_scale_default_to_one(self):
         expr = pybamm.sqrt(c_s)
         result = regulariser(expr, inputs=inputs)
 
-        # Should be replaced with RegPower (no Sqrt)
-        has_sqrt = any(isinstance(n, pybamm.Sqrt) for n in result.pre_order())
-        assert not has_sqrt
-        # Check it's a RegPower with scale=1 (default)
         assert isinstance(result, pybamm.RegPower)
-        # Scale is the third child
         assert result.children[2] == pybamm.Scalar(1)
 
+    def test_normalised_state_dependent_base_still_regularised(self):
+        """Normalised base (ORegan2022-style (c_e / c_e_ref) ** (1 - alpha))
+        matches no scale but is still regularised with the default scale."""
+        c_e = pybamm.Variable("c_e")
+        inputs = {"c_e": c_e}
+        regulariser = pybamm.RegulariseSqrtAndPower(
+            {c_e: pybamm.Scalar(1000.0)},
+            inputs=inputs,
+        )
+
+        c_e_ref = pybamm.Parameter("Reference concentration")
+        expr = (c_e / c_e_ref) ** 0.208
+        result = regulariser(expr, inputs=inputs)
+
+        assert isinstance(result, pybamm.RegPower)
+        assert result.children[2] == pybamm.Scalar(1)
+
+    def test_parameter_power_not_corrupted(self):
+        """State-independent Parameter base is not regularised: RegPower with
+        the default scale would corrupt a small rate constant."""
+        c_e = pybamm.Variable("c_e")
+        inputs = {"c_e": c_e}
+        regulariser = pybamm.RegulariseSqrtAndPower(
+            {c_e: pybamm.Scalar(1000.0)},
+            inputs=inputs,
+        )
+
+        rate_constant = pybamm.Parameter("Rate constant")
+        expr = rate_constant**0.5
+        result = regulariser(expr, inputs=inputs)
+
+        assert isinstance(result, pybamm.Power)
+        assert not any(isinstance(n, pybamm.RegPower) for n in result.pre_order())
+
+        parameter_values = pybamm.ParameterValues({"Rate constant": 5e-9})
+        value = parameter_values.process_symbol(result).evaluate()
+        assert value == pytest.approx(5e-9**0.5, rel=1e-12)
+
+    def test_input_parameter_power_not_corrupted(self):
+        """State-independent InputParameter base is not regularised."""
+        c_e = pybamm.Variable("c_e")
+        inputs = {"c_e": c_e}
+        regulariser = pybamm.RegulariseSqrtAndPower(
+            {c_e: pybamm.Scalar(1000.0)},
+            inputs=inputs,
+        )
+
+        rate_constant = pybamm.InputParameter("Rate constant")
+        expr = rate_constant**0.5
+        result = regulariser(expr, inputs=inputs)
+
+        assert isinstance(result, pybamm.Power)
+        value = result.evaluate(inputs={"Rate constant": 5e-9})
+        assert value == pytest.approx(5e-9**0.5, rel=1e-12)
+
+    def test_constant_base_state_dependent_exponent(self):
+        """Constant base, state-dependent exponent: outer Power kept, sqrt in
+        the exponent regularised."""
+        c_e = pybamm.Variable("c_e")
+        inputs = {"c_e": c_e}
+        regulariser = pybamm.RegulariseSqrtAndPower(
+            {c_e: pybamm.Scalar(1000.0)},
+            inputs=inputs,
+        )
+
+        base = pybamm.Parameter("Rate constant")
+        expr = base ** pybamm.sqrt(c_e)
+        result = regulariser(expr, inputs=inputs)
+
+        assert isinstance(result, pybamm.Power)
+        assert isinstance(result.children[0], pybamm.Parameter)
+        assert any(isinstance(n, pybamm.RegPower) for n in result.pre_order())
+        assert not any(isinstance(n, pybamm.Sqrt) for n in result.pre_order())
+
+    def test_existing_reg_power_preserved(self):
+        """A manual RegPower node (e.g. MSMR j0_j) is a Function, not a
+        Sqrt/Power, so it is left untouched, even when nested."""
+        c_e = pybamm.Variable("c_e")
+        c_e_ref = pybamm.Parameter("c_e_ref")
+        aj = pybamm.Parameter("aj")
+        inputs = {"c_e": c_e}
+        regulariser = pybamm.RegulariseSqrtAndPower(
+            {c_e: pybamm.Scalar(1000.0)},
+            inputs=inputs,
+        )
+
+        expr = pybamm.reg_power(c_e / c_e_ref, 1 - aj)
+        result = regulariser(expr, inputs=inputs)
+        assert result is expr
+
+        # Nested: reg_power preserved, sibling Power regularised
+        nested = pybamm.reg_power(c_e / c_e_ref, 1 - aj) * c_e**0.5
+        result = regulariser(nested, inputs=inputs)
+        n_reg_power = sum(
+            1 for n in result.pre_order() if isinstance(n, pybamm.RegPower)
+        )
+        assert n_reg_power == 2
+        assert not any(isinstance(n, pybamm.Power) for n in result.pre_order())
+        assert not any(isinstance(n, pybamm.Sqrt) for n in result.pre_order())
+
+    @pytest.mark.parametrize("set_name", ["Chen2020", "ORegan2022", "Ecker2015"])
+    def test_j0_pipeline_produces_regpower(self, set_name):
+        """Built-in j0 pipeline regularises all concentration powers,
+        including ORegan2022's unmatched normalised bases."""
+        param = pybamm.LithiumIonParameters()
+        c_e = pybamm.Variable("c_e")
+        c_s_surf = pybamm.Variable("c_s_surf")
+        T = pybamm.Variable("T")
+        j0 = param.n.prim.j0(c_e, c_s_surf, T)
+
+        parameter_values = pybamm.ParameterValues(set_name)
+        processed = parameter_values.process_symbol(j0)
+
+        n_reg_power = sum(
+            1 for n in processed.pre_order() if isinstance(n, pybamm.RegPower)
+        )
+        n_plain_power = sum(
+            1 for n in processed.pre_order() if isinstance(n, pybamm.Power)
+        )
+        assert n_reg_power == 3, f"{set_name}: expected 3 RegPower, got {n_reg_power}"
+        assert n_plain_power == 0, (
+            f"{set_name}: expected no plain Power nodes, got {n_plain_power}"
+        )
+
+    def test_function_parameter_post_processor_small_parameter(self):
+        """Full j0-style path: a FunctionParameter raising a Parameter to a
+        fractional power must evaluate exactly."""
+
+        def exchange_current(c_e):
+            rate_constant = pybamm.Parameter("Rate constant")
+            return rate_constant**0.5 * c_e**0.5
+
+        c_e = pybamm.Variable("c_e")
+        inputs = {"Electrolyte concentration [mol.m-3]": c_e}
+        regulariser = pybamm.RegulariseSqrtAndPower(
+            {c_e: pybamm.Scalar(1000.0)},
+            inputs=inputs,
+        )
+        function_parameter = pybamm.FunctionParameter(
+            "Exchange-current density [A.m-2]",
+            {"Electrolyte concentration [mol.m-3]": pybamm.Scalar(1000.0)},
+            post_processor=regulariser,
+        )
+        parameter_values = pybamm.ParameterValues(
+            {
+                "Exchange-current density [A.m-2]": exchange_current,
+                "Rate constant": 5e-9,
+            }
+        )
+        value = parameter_values.process_symbol(function_parameter).evaluate()
+        expected = 5e-9**0.5 * 1000.0**0.5
+        assert value == pytest.approx(expected, rel=1e-6)
+
     def test_exact_match_only(self):
         """Test that only exact matches are used for scales."""
         c_s = pybamm.Variable("c_s")
@@ -85,7 +235,7 @@ def test_exact_match_only(self):
         has_sqrt = any(isinstance(n, pybamm.Sqrt) for n in result1.pre_order())
         assert not has_sqrt
 
-        # sqrt(c_s / c_s_max) should also be replaced
+        # sqrt(c_s / c_s_max) should also be replaced (state-dependent base)
         expr2 = pybamm.sqrt(c_s / c_s_max)
         result2 = regulariser(expr2, inputs=inputs)
         has_sqrt = any(isinstance(n, pybamm.Sqrt) for n in result2.pre_order())

From 1cd0d3189762f426c4c67fc37dd455cc390c6550 Mon Sep 17 00:00:00 2001
From: Brady Planden <55357039+BradyPlanden@users.noreply.github.com>
Date: Mon, 15 Jun 2026 23:06:06 +0100
Subject: [PATCH 2/6] Harden serialisation round-trip tests + #5598 (#5599)

* Harden serialisation round-trip tests; fix param-rebuild altitude and MSMR option gaps

Follow-up to the tuple-options fix, addressing issues found while hardening
the serialisation test suite.

Rebuild param at the right altitude. The previous fix rebuilt a loaded
lithium-ion model's `param` with an `isinstance` branch inside
`load_custom_model`. Replace it with a model-owned `_rebuild_param` hook
(no-op on `BaseModel`, overridden by the lithium-ion base model, which now
also uses it in `__init__`). `load_custom_model` calls the hook, so the
options setter stays free of param concerns and the fix lives where the
model owns its derived state (#5598).

Fix two MSMR option-validation gaps surfaced by option fuzzing. Both let an
invalid option set crash later in parameter construction with "invalid
literal for int() with base 10: 'none'" instead of raising a clean
OptionError:
- the all-or-nothing MSMR check missed MSMR inside a per-electrode tuple
  (`("MSMR", "single")`) because it compared the whole value to "MSMR";
- an all-MSMR model left on the default "number of MSMR reactions" ("none")
  passed validation, since the registry only lists that placeholder.

Harden the serialisation test suite:
- generative fuzzing of valid option dicts through to_json/from_json, with
  the construction guard narrowed to (OptionError, NotImplementedError) so
  real bugs surface instead of being silently rejected;
- deterministic tuple-option and int-dimensionality round-trip coverage;
- discretised round-trip now asserts structural identity, solution
  equivalence and event preservation, plus a non-default tuple option.

Co-authored-by: Alec Bills <alexanderallenbills@gmail.com>

* Address review: setter-owned param rebuild, per-electrode MSMR count, symptom assertion

- Move `_rebuild_param` into the `BaseBatteryModel.options` setter so the
  options-derived `param` is kept in sync automatically. Drop the now-redundant
  explicit calls in `lithium_ion.BaseModel.__init__` and `Serialise` deserialise.
  Repurpose the obsolete `test_rebuild_param_syncs_param_with_reassigned_options`
  into `test_options_setter_rebuilds_param`, asserting the new invariant.

- Fix an MSMR regression: validate `number of MSMR reactions` per electrode so a
  mixed cell (MSMR in one electrode, conventional with a "none" count in the
  other) is accepted again, while the genuine all-MSMR "none" default is still
  rejected early. Remove the now-dead `_is_positive_integer_count`; add
  `test_msmr_mixed_electrode_accepted`.

- Strengthen `_assert_param_matches_options` with a symptom-level check: a
  two-phase negative electrode must name its primary active-material fraction
  "Primary: ...", guarding a stale/no-op rebuild beyond the dict comparison.

* fix: MSMR option validation dependent on electrode domains

---------

Co-authored-by: Alec Bills <alexanderallenbills@gmail.com>
(cherry picked from commit 95db6907a151745487acb031ca39481d1f51d510)
---
 CHANGELOG.md                                  |   2 +
 .../expression_tree/operations/serialise.py   |   7 +-
 src/pybamm/models/base_model.py               |   6 +
 .../full_battery_models/base_battery_model.py |  54 ++-
 .../lithium_ion/base_lithium_ion_model.py     |   6 +-
 .../test_base_battery_model.py                |  54 +++
 .../test_serialisation/test_round_trip.py     | 340 ++++++++++++++++++
 7 files changed, 458 insertions(+), 11 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3e93ab3adb..1c9c12711f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,8 @@
 ## Bug fixes
 
 - `RegulariseSqrtAndPower` no longer regularises state-independent bases, fixing corrupted small rate constants in exchange-current density functions. ([#5600](https://github.com/pybamm-team/PyBaMM/pull/5600))
+- Fixed `Serialise.load_custom_model` leaving a loaded lithium-ion model's cached `param` built from default options instead of the restored ones, so a model loaded with non-default options (e.g. composite `"particle phases"`) had parameters inconsistent with its options. ([#5599](https://github.com/pybamm-team/PyBaMM/pull/5599))
+- Fixed `BatteryModelOptions` letting two invalid MSMR option sets pass validation and then crash during parameter construction with `invalid literal for int() with base 10: 'none'`; both now raise a clear `OptionError`. The all-or-nothing MSMR check now detects MSMR inside a per-electrode tuple (e.g. `"open-circuit potential": ("MSMR", "single")`), and an MSMR model must set `"number of MSMR reactions"` to a positive integer rather than the default `"none"`. ([#5599](https://github.com/pybamm-team/PyBaMM/pull/5599))
 
 # [v26.6.1.1](https://github.com/pybamm-team/PyBaMM/tree/v26.6.1.1) - 2026-06-11
 
diff --git a/src/pybamm/expression_tree/operations/serialise.py b/src/pybamm/expression_tree/operations/serialise.py
index b11dac109c..038f33e672 100644
--- a/src/pybamm/expression_tree/operations/serialise.py
+++ b/src/pybamm/expression_tree/operations/serialise.py
@@ -1386,10 +1386,9 @@ def load_custom_model(filename: str | dict) -> pybamm.BaseModel:
         model = base_cls()
         model.name = model_data["name"]
         model.schema_version = schema_version
-        # Restore options so round-trip serialisation produces an equivalent
-        # model. A JSON round-trip turns tuple-valued options (e.g. "particle
-        # phases": ("2", "1")) into lists, which pybamm's options validation
-        # rejects, so convert lists back to tuples before assigning.
+        # JSON turns tuple-valued options into lists (rejected by validation);
+        # convert back. Assigning options rebuilds the options-derived param via
+        # the setter, so it matches the restored options.
         opts = model_data.get("options", {})
         if opts is not None:
             model.options = Serialise._convert_options(opts)
diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index 333ec48e4f..1a3bcd0513 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -641,6 +641,12 @@ def param(self):
     def param(self, values):
         self._param = values
 
+    def _rebuild_param(self):
+        """Rebuild ``param`` from ``self.options``; no-op unless ``param`` is
+        options-derived. Called by the ``options`` setter whenever options are
+        (re)assigned, so subclasses keep ``param`` in sync automatically.
+        """
+
     @property
     def options(self):
         """Returns the model options dictionary that allows customization of the model's behavior."""
diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py
index 8bf74ba715..162a7f6beb 100644
--- a/src/pybamm/models/full_battery_models/base_battery_model.py
+++ b/src/pybamm/models/full_battery_models/base_battery_model.py
@@ -18,6 +18,24 @@ def represents_positive_integer(s):
         return val > 0
 
 
+def _is_msmr(value):
+    """True if an option requests MSMR, including inside a per-electrode tuple."""
+    if isinstance(value, (tuple, list)):
+        return "MSMR" in value
+    return value == "MSMR"
+
+
+def _per_electrode(value, index):
+    """Resolve a possibly per-electrode option to one electrode.
+
+    Per-electrode tuples are ``(negative, positive)``; a scalar applies to
+    both. ``index`` is 0 for the negative electrode, 1 for the positive.
+    """
+    if isinstance(value, (tuple, list)):
+        return value[index]
+    return value
+
+
 def _rename_option(options_dict, option_name, old_name, new_name):
     if option_name not in options_dict:
         return
@@ -621,13 +639,11 @@ def __init__(self, extra_options):
                         f"Option '{name}' not recognised. Best matches are {options.get_best_matches(name)}"
                     )
 
-        # If any of "open-circuit potential", "particle" or "intercalation kinetics" is
-        # "MSMR" then all of them must be "MSMR".
-        # Note: this check is currently performed on full cells, but is loosened for
-        # half-cells where you must pass a tuple of options to only set MSMR models in
-        # the working electrode
+        # All-or-nothing on full cells: if any of OCP/particle/intercalation
+        # kinetics requests MSMR (incl. inside a per-electrode tuple), all must.
+        # Half-cells are loosened -- a tuple sets MSMR in the working electrode.
         msmr_check_list = [
-            options[opt] == "MSMR"
+            _is_msmr(options[opt])
             for opt in ["open-circuit potential", "particle", "intercalation kinetics"]
         ]
         if (
@@ -640,6 +656,29 @@ def __init__(self, extra_options):
                 "'intercalation kinetics' is 'MSMR' then all of them must be 'MSMR'"
             )
 
+        # Validate per electrode so a mixed full cell (MSMR in one electrode,
+        # conventional in the other) is accepted.
+        if options["working electrode"] == "both":
+            electrode_domains = [("negative", 0), ("positive", 1)]
+        else:
+            electrode_domains = [("positive", 1)]
+        for domain, index in electrode_domains:
+            domain_uses_msmr = any(
+                _per_electrode(options[opt], index) == "MSMR"
+                for opt in [
+                    "open-circuit potential",
+                    "particle",
+                    "intercalation kinetics",
+                ]
+            )
+            domain_count = _per_electrode(options["number of MSMR reactions"], index)
+            if domain_uses_msmr and not represents_positive_integer(domain_count):
+                raise pybamm.OptionError(
+                    "'number of MSMR reactions' must be a positive integer for "
+                    f"the {domain} electrode when it uses 'MSMR' "
+                    f"(got {domain_count!r})"
+                )
+
         # If "SEI film resistance" is "distributed" then "total interfacial current
         # density as a state" must be "true"
         if options["SEI film resistance"] == "distributed":
@@ -1209,6 +1248,9 @@ def options(self, extra_options):
                 f"must use surface formulation to solve {self!s} with hydrolysis"
             )
         self._options = options
+        # rebuild whenever options are (re)assigned.
+        # No-op unless the subclass overrides ``_rebuild_param``.
+        self._rebuild_param()
 
     def set_standard_output_variables(self):
         # Time
diff --git a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index 3f48532fda..207b49115f 100644
--- a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -29,13 +29,17 @@ class BaseModel(pybamm.BaseBatteryModel):
 
     def __init__(self, options=None, name="Unnamed lithium-ion model", build=False):
         super().__init__(options, name)
-        self.param = pybamm.LithiumIonParameters(self.options)
 
         self.set_standard_output_variables()
 
         # Li models should calculate esoh by default
         self._calc_esoh = True
 
+    def _rebuild_param(self):
+        # param is options-derived, so rebuild it whenever options are
+        # (re)assigned: at construction and on deserialisation.
+        self.param = pybamm.LithiumIonParameters(self.options)
+
     def set_submodels(self, build):
         self.set_external_circuit_submodel()
         self.set_porosity_submodel()
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index cf25799524..cd14f8b638 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -444,6 +444,60 @@ def test_options(self):
                     "number of MSMR reactions": "1.5",
                 }
             )
+        # MSMR inside a per-electrode tuple must be rejected cleanly, not slip
+        # through to crash later as int('none') in parameter construction.
+        with pytest.raises(pybamm.OptionError, match=r"MSMR"):
+            pybamm.BaseBatteryModel({"open-circuit potential": ("MSMR", "single")})
+        with pytest.raises(pybamm.OptionError, match=r"MSMR"):
+            pybamm.BaseBatteryModel({"particle": ("MSMR", "Fickian diffusion")})
+        # MSMR with the default "none" reaction count must also be rejected cleanly.
+        with pytest.raises(pybamm.OptionError, match=r"MSMR"):
+            pybamm.BaseBatteryModel(
+                {
+                    "open-circuit potential": "MSMR",
+                    "particle": "MSMR",
+                    "intercalation kinetics": "MSMR",
+                }
+            )
+
+    def test_msmr_mixed_electrode_accepted(self):
+        # Verify non-MSMR electrode's "none" reaction count is legitimate
+        options = {
+            "open-circuit potential": ("MSMR", "single"),
+            "particle": ("MSMR", "Fickian diffusion"),
+            "intercalation kinetics": ("MSMR", "symmetric Butler-Volmer"),
+            "number of MSMR reactions": ("3", "none"),
+        }
+        model = pybamm.BaseBatteryModel(options)
+        assert model.options["number of MSMR reactions"] == ("3", "none")
+
+    def test_msmr_half_cell_does_not_validate_counter_electrode(self):
+        # On a half cell the negative electrode is lithium metal, not a porous
+        # MSMR electrode, so a scalar "MSMR" option must not force a reaction
+        # count on it.
+        options = {
+            "working electrode": "positive",
+            "open-circuit potential": "MSMR",
+            "particle": "MSMR",
+            "intercalation kinetics": "MSMR",
+            "number of MSMR reactions": ("none", "6"),
+        }
+        model = pybamm.BaseBatteryModel(options)
+        assert model.options["number of MSMR reactions"] == ("none", "6")
+
+    def test_msmr_half_cell_still_validates_working_electrode(self):
+        # The working electrode is still validated: a "none" count
+        # there must be rejected cleanly rather than slipping through.
+        with pytest.raises(pybamm.OptionError, match=r"positive electrode"):
+            pybamm.BaseBatteryModel(
+                {
+                    "working electrode": "positive",
+                    "open-circuit potential": "MSMR",
+                    "particle": "MSMR",
+                    "intercalation kinetics": "MSMR",
+                    "number of MSMR reactions": ("6", "none"),
+                }
+            )
 
     def test_build_twice(self):
         model = pybamm.lithium_ion.SPM()  # need to pick a model to set vars and build
diff --git a/tests/unit/test_serialisation/test_round_trip.py b/tests/unit/test_serialisation/test_round_trip.py
index 99dc8ffbcc..26dbc42fb5 100644
--- a/tests/unit/test_serialisation/test_round_trip.py
+++ b/tests/unit/test_serialisation/test_round_trip.py
@@ -11,12 +11,17 @@
 
 from __future__ import annotations
 
+import ast
+import inspect
 import json
+import textwrap
 import warnings
 from datetime import datetime
 
 import numpy as np
 import pytest
+from hypothesis import given, reject, settings
+from hypothesis import strategies as st
 
 import pybamm
 from pybamm.expression_tree.operations.serialise import (
@@ -24,6 +29,7 @@
     convert_symbol_from_json,
     convert_symbol_to_json,
 )
+from pybamm.models.full_battery_models.base_battery_model import BatteryModelOptions
 from tests.unit.test_serialisation._helpers import (
     _experiments_equal,
     _solver_init_args_equal,
@@ -500,3 +506,337 @@ def test_preserves_structure(self, symbol):
             f"{type(symbol).__name__} structural id changed across round-trip: "
             f"{symbol!r} -> {restored!r}"
         )
+
+
+# Generative option round-trip fuzzing. JSON has no tuple type, so a tuple option
+# (e.g. "particle phases": ("2", "1")) deserialises as a list, which validation
+# rejects. Fuzzes the build-free to_json/from_json path; the deterministic
+# fixtures below cover to_config/from_config (which builds the model).
+
+# Option registry, sourced from pybamm so the strategy never hardcodes values.
+_POSSIBLE_OPTIONS = BatteryModelOptions({}).possible_options
+
+_STRING_VALUES = {
+    key: [v for v in values if isinstance(v, str)]
+    for key, values in _POSSIBLE_OPTIONS.items()
+}
+
+# Drawable keys. "dimensionality" (the only int option) is excluded: a non-zero
+# value needs companion geometry options the random strategy can't assemble; its
+# int round-trip is covered by test_int_option_survives_json_round_trip_as_int.
+_STRING_OPTION_KEYS = sorted(
+    key for key, values in _STRING_VALUES.items() if key != "dimensionality" and values
+)
+
+# Options accepting a 2-tuple of strings; mirrors the list in BatteryModelOptions
+# (test_tuple_capable_keys_match_validation_source guards drift).
+_TUPLE_CAPABLE_KEYS = frozenset(
+    {
+        "diffusivity",
+        "exchange-current density",
+        "intercalation kinetics",
+        "interface utilisation",
+        "lithium plating",
+        "loss of active material",
+        "number of MSMR reactions",
+        "open-circuit potential",
+        "particle",
+        "particle mechanics",
+        "particle phases",
+        "particle size",
+        "SEI",
+        "SEI on cracks",
+        "stress-induced diffusion",
+    }
+)
+
+# SPM matches the original regression; DFN has a broader option surface.
+_OPTION_MODEL_CLASSES = [pybamm.lithium_ion.SPM, pybamm.lithium_ion.DFN]
+
+
+def _assert_param_matches_options(model):
+    """Assert a loaded lithium-ion model's ``param`` matches its restored options.
+
+    Guards #5598; lead-acid ``param`` is options-independent, so skip it.
+    """
+    if isinstance(model, pybamm.lithium_ion.BaseModel):
+        assert dict(model.param.options) == dict(model.options), (
+            f"param built from {dict(model.param.options)!r} != restored options "
+            f"{dict(model.options)!r}; load_custom_model must rebuild param."
+        )
+        # Symptom-level guard: a two-phase negative electrode names its primary
+        # active-material fraction "Primary: ..."; a param left built from
+        # single-phase options would not.
+        if model.options.negative["particle phases"] == "2":
+            epsilon_s_name = model.param.n.prim.epsilon_s.name
+            assert epsilon_s_name.startswith("Primary:"), (
+                f"negative primary active material fraction is {epsilon_s_name!r}; "
+                "expected a 'Primary:' prefix for a two-phase electrode -- param "
+                "was not rebuilt from the restored composite options."
+            )
+
+
+def test_options_setter_rebuilds_param():
+    """Assigning ``options`` post-construction rebuilds the options-derived
+    ``param`` via the setter, keeping it in sync without an explicit hook call.
+    """
+    # A complete, valid two-phase options set (as a load path would restore).
+    full_opts = dict(pybamm.lithium_ion.SPM({"particle phases": ("2", "1")}).options)
+
+    model = pybamm.lithium_ion.SPM(build=False)  # default: single particle phase
+    model.options = full_opts
+    # The setter rebuilds param, so it tracks the reassigned options.
+    assert dict(model.param.options) == dict(model.options)
+
+
+@st.composite
+def valid_option_dicts(draw):
+    """Draw a small registry-valid option dict, biased to exercise tuples."""
+    keys = draw(
+        st.lists(
+            st.sampled_from(_STRING_OPTION_KEYS),
+            min_size=1,
+            max_size=3,
+            unique=True,
+        )
+    )
+    options = {}
+    for key in keys:
+        choices = _STRING_VALUES[key]
+        if key in _TUPLE_CAPABLE_KEYS and draw(st.booleans()):
+            options[key] = (
+                draw(st.sampled_from(choices)),
+                draw(st.sampled_from(choices)),
+            )
+        else:
+            options[key] = draw(st.sampled_from(choices))
+    return options
+
+
+def test_tuple_capable_keys_match_validation_source():
+    """``_TUPLE_CAPABLE_KEYS`` must match the tuple-capable list in
+    BatteryModelOptions -- AST-extracted so it can't silently drift.
+    """
+    tree = ast.parse(textwrap.dedent(inspect.getsource(BatteryModelOptions)))
+    candidates = [
+        {
+            elt.value
+            for elt in node.elts
+            if isinstance(elt, ast.Constant) and isinstance(elt.value, str)
+        }
+        for node in ast.walk(tree)
+        if isinstance(node, ast.List)
+    ]
+    # The tuple-capable list is the one literal containing these marker keys.
+    matches = [c for c in candidates if {"particle phases", "SEI on cracks"} <= c]
+    assert len(matches) == 1, (
+        "Could not uniquely locate the tuple-capable option list in "
+        "BatteryModelOptions source; update the extraction in this test."
+    )
+    assert matches[0] == set(_TUPLE_CAPABLE_KEYS), (
+        "Tuple-capable options in pybamm have changed. Update "
+        "_TUPLE_CAPABLE_KEYS so option fuzzing keeps generating tuple values "
+        f"for every such key. Source: {sorted(matches[0])}, "
+        f"test constant: {sorted(_TUPLE_CAPABLE_KEYS)}."
+    )
+    assert _TUPLE_CAPABLE_KEYS <= set(_POSSIBLE_OPTIONS), (
+        "_TUPLE_CAPABLE_KEYS contains keys absent from the option registry."
+    )
+
+
+class TestModelOptionsRoundTrip:
+    """Any valid option dict must survive a JSON round-trip with exact types."""
+
+    @pytest.mark.parametrize(
+        "model_cls", _OPTION_MODEL_CLASSES, ids=lambda c: c.__name__
+    )
+    @settings(max_examples=100, deadline=None)
+    @given(options=valid_option_dicts())
+    def test_options_survive_json_round_trip(self, model_cls, options):
+        try:
+            # Construction is the precondition, not the SUT; pass a copy
+            # (pybamm mutates the dict).
+            model = model_cls(dict(options), build=False)
+        except (pybamm.OptionError, NotImplementedError):
+            # Invalid/unsupported option set -> discard. Any other exception is a
+            # real bug that must surface (a broad except here once hid an MSMR
+            # int('none') crash).
+            reject()
+
+        expected = dict(model.options)
+
+        via_json = pybamm.BaseModel.from_json(
+            json.loads(json.dumps(model.to_json(), default=Serialise._json_encoder))
+        )
+        assert _values_equal(expected, dict(via_json.options)), (
+            f"{model_cls.__name__} options lost or mutated via "
+            f"to_json/from_json for input {options!r}: "
+            f"{expected!r} -> {dict(via_json.options)!r}"
+        )
+        # Options surviving isn't enough -- the derived param must match too.
+        _assert_param_matches_options(via_json)
+
+
+def test_int_option_survives_json_round_trip_as_int():
+    """``dimensionality`` (the only int option) must round-trip as an int, not be
+    coerced to/from bool (the #5495 shape, ``True == 1``). Bool-like options are
+    strings "true"/"false" (pybamm rejects Python bools), already fuzzed; a
+    non-zero value needs companion geometry options to construct.
+    """
+    options = {
+        "dimensionality": 1,
+        "current collector": "potential pair",
+        "cell geometry": "pouch",
+    }
+    model = pybamm.lithium_ion.DFN(dict(options), build=False)
+
+    via_json = pybamm.BaseModel.from_json(
+        json.loads(json.dumps(model.to_json(), default=Serialise._json_encoder))
+    )
+    restored = via_json.options["dimensionality"]
+    assert restored == 1
+    assert isinstance(restored, int) and not isinstance(restored, bool)
+
+
+# Buildable tuple options that must round-trip through both load paths --
+# deterministic, always-run coverage of the bug shape.
+_TUPLE_OPTION_FIXTURES = [
+    (pybamm.lithium_ion.SPM, "particle phases", ("2", "1")),
+    (pybamm.lithium_ion.SPM, "SEI", ("none", "constant")),
+    (pybamm.lithium_ion.DFN, "particle size", ("single", "distribution")),
+    (pybamm.lithium_ion.DFN, "particle", ("Fickian diffusion", "uniform profile")),
+]
+
+
+class TestTupleOptionRoundTripRegression:
+    @pytest.mark.parametrize(
+        "model_cls,key,value",
+        _TUPLE_OPTION_FIXTURES,
+        ids=["SPM-particle phases", "SPM-SEI", "DFN-particle size", "DFN-particle"],
+    )
+    def test_tuple_option_survives_both_paths(self, model_cls, key, value):
+        model = model_cls({key: value}, build=False)
+
+        via_json = pybamm.BaseModel.from_json(
+            json.loads(json.dumps(model.to_json(), default=Serialise._json_encoder))
+        )
+        assert via_json.options[key] == value
+        assert isinstance(via_json.options[key], tuple)
+        # Composite "particle phases" is the case (param stuck single-phase).
+        _assert_param_matches_options(via_json)
+
+        via_config = pybamm.BaseModel.from_config(
+            json.loads(json.dumps(model.to_config()))
+        )
+        assert via_config.options[key] == value
+        assert isinstance(via_config.options[key], tuple)
+
+
+# Discretised-model round-trip (save_model/load_model): previously only checked
+# "load and solve". These verify structural identity and solution equivalence
+# (mesh transitively: same equations + solution => same mesh).
+
+_DISCRETISED_MODEL_CLASSES = [pybamm.lithium_ion.SPM, pybamm.lithium_ion.DFN]
+
+
+def _discretised_model(model_cls, options=None):
+    """Build, parameterise and discretise a model; return it and its mesh
+    (the mesh is needed by ``serialise_model`` for named variables).
+    """
+    model = model_cls(options)
+    geometry = model.default_geometry
+    param = model.default_parameter_values
+    param.process_model(model)
+    param.process_geometry(geometry)
+    mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
+    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    disc.process_model(model)
+    return model, mesh
+
+
+def _round_trip_discretised(model, mesh):
+    """Serialise a discretised model and load it back, via dict (no temp file)."""
+    serialised = Serialise().serialise_model(model, mesh)
+    return Serialise().load_model(json.loads(json.dumps(serialised)))
+
+
+@pytest.mark.parametrize(
+    "model_cls", _DISCRETISED_MODEL_CLASSES, ids=lambda c: c.__name__
+)
+class TestDiscretisedModelRoundTrip:
+    def test_derived_state_matches_original(self, model_cls):
+        """Concatenated equations, mass matrix, bounds and options must match."""
+        model, mesh = _discretised_model(model_cls)
+        loaded = _round_trip_discretised(model, mesh)
+
+        for attr in (
+            "concatenated_rhs",
+            "concatenated_algebraic",
+            "concatenated_initial_conditions",
+            "mass_matrix",
+        ):
+            original = getattr(model, attr)
+            restored = getattr(loaded, attr)
+            assert original.id == restored.id, (
+                f"{model_cls.__name__}.{attr} changed across the discretised "
+                f"round-trip: {original.id} != {restored.id}"
+            )
+
+        assert all(
+            np.array_equal(a, b)
+            for a, b in zip(model.bounds, loaded.bounds, strict=True)
+        ), f"{model_cls.__name__} bounds changed across the discretised round-trip"
+        assert dict(model.options) == dict(loaded.options)
+
+        # Events are serialised too -- verify them structurally rather than
+        # relying on the solution check.
+        assert len(model.events) == len(loaded.events)
+        for original_event, restored_event in zip(
+            model.events, loaded.events, strict=True
+        ):
+            assert original_event.name == restored_event.name
+            assert original_event.event_type == restored_event.event_type
+            assert original_event.expression.id == restored_event.expression.id, (
+                f"{model_cls.__name__} event {original_event.name!r} changed "
+                "across the discretised round-trip"
+            )
+
+    def test_solution_matches_original(self, model_cls):
+        """The reloaded model must solve to the same answer as the original."""
+        model, mesh = _discretised_model(model_cls)
+        # Serialise before solving so any solve-time caching cannot leak in.
+        loaded = _round_trip_discretised(model, mesh)
+
+        t_eval = [0, 3600]
+        original_solution = model.default_solver.solve(model, t_eval)
+        loaded_solution = loaded.default_solver.solve(loaded, t_eval)
+
+        assert np.allclose(
+            original_solution.y, loaded_solution.y, rtol=1e-6, atol=1e-8
+        ), (
+            f"{model_cls.__name__} state trajectory diverged after a "
+            "discretised round-trip"
+        )
+        assert np.allclose(
+            original_solution["Voltage [V]"].entries,
+            loaded_solution["Voltage [V]"].entries,
+            rtol=1e-6,
+            atol=1e-8,
+        ), f"{model_cls.__name__} Voltage [V] diverged after a discretised round-trip"
+
+
+def test_nondefault_tuple_option_survives_discretised_round_trip():
+    """A tuple-valued option must survive the discretised save/load path too.
+
+    The parametrised tests above only use default options; this guards
+    ``load_model``'s list->tuple conversion. ``particle`` is used because it
+    discretises with the default parameter set (composite ``particle phases``
+    would need a phase-prefixed one).
+    """
+    options = {"particle": ("Fickian diffusion", "uniform profile")}
+    model, mesh = _discretised_model(pybamm.lithium_ion.DFN, options)
+    loaded = _round_trip_discretised(model, mesh)
+
+    assert loaded.options["particle"] == ("Fickian diffusion", "uniform profile")
+    assert isinstance(loaded.options["particle"], tuple)
+    _assert_param_matches_options(loaded)

From 9ae0ad504d90aed41040930227e16d331b6ddc6d Mon Sep 17 00:00:00 2001
From: Alec Bills <48105066+aabills@users.noreply.github.com>
Date: Tue, 16 Jun 2026 02:19:37 -0700
Subject: [PATCH 3/6] Fix flaky half-cell case in serialisation round-trip
 param check (#5608)

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 tests/unit/test_serialisation/test_round_trip.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/tests/unit/test_serialisation/test_round_trip.py b/tests/unit/test_serialisation/test_round_trip.py
index 26dbc42fb5..d8f1d58fca 100644
--- a/tests/unit/test_serialisation/test_round_trip.py
+++ b/tests/unit/test_serialisation/test_round_trip.py
@@ -566,8 +566,13 @@ def _assert_param_matches_options(model):
         )
         # Symptom-level guard: a two-phase negative electrode names its primary
         # active-material fraction "Primary: ..."; a param left built from
-        # single-phase options would not.
-        if model.options.negative["particle phases"] == "2":
+        # single-phase options would not. Only meaningful on a full cell -- a
+        # half cell's negative electrode is lithium metal with no porous
+        # particle, so ``param.n.prim`` has no ``epsilon_s``.
+        if (
+            model.options["working electrode"] == "both"
+            and model.options.negative["particle phases"] == "2"
+        ):
             epsilon_s_name = model.param.n.prim.epsilon_s.name
             assert epsilon_s_name.startswith("Primary:"), (
                 f"negative primary active material fraction is {epsilon_s_name!r}; "

From 2f13e4397ff6cc0e34817417e5dcab197ad52f6e Mon Sep 17 00:00:00 2001
From: Alec Bills <48105066+aabills@users.noreply.github.com>
Date: Tue, 16 Jun 2026 03:31:38 -0700
Subject: [PATCH 4/6] Preserve custom-model discretisation recipe across
 serialise/load (#5609)

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Brady Planden <55357039+BradyPlanden@users.noreply.github.com>
---
 .../expression_tree/operations/serialise.py   | 149 +++++++++++++-----
 src/pybamm/models/base_model.py               |   8 +-
 .../test_serialisation/test_serialisation.py  | 106 +++++++++++++
 3 files changed, 221 insertions(+), 42 deletions(-)

diff --git a/src/pybamm/expression_tree/operations/serialise.py b/src/pybamm/expression_tree/operations/serialise.py
index 038f33e672..3c83580d3a 100644
--- a/src/pybamm/expression_tree/operations/serialise.py
+++ b/src/pybamm/expression_tree/operations/serialise.py
@@ -534,6 +534,15 @@ def serialise_custom_model(model: pybamm.BaseModel, compress: bool = False) -> d
             },
         }
 
+        # Capture the discretisation recipe (geometry / var_pts / submesh /
+        # spatial methods). These live on subclass ``default_*`` properties and
+        # are otherwise lost when the model is reloaded in an environment where
+        # the defining package can't be imported (the base-class fallback yields
+        # empty defaults), leaving the model undiscretisable.
+        discretisation = Serialise._serialise_discretisation(model)
+        if discretisation:
+            model_content["discretisation"] = discretisation
+
         SCHEMA_VERSION = "1.1"
         model_json = {
             "schema_version": SCHEMA_VERSION,
@@ -628,6 +637,35 @@ def save_custom_model(
         except Exception as e:
             raise ValueError(f"Failed to save custom model: {e}") from e
 
+    @staticmethod
+    def _encode_geometry_value(value):
+        """Recursively convert any :class:`pybamm.Symbol` in a geometry value to
+        its JSON form, descending through nested dicts/lists.
+
+        Geometry limits can nest symbols arbitrarily deep (e.g. current-collector
+        ``tabs -> negative -> z_centre`` is a :class:`pybamm.Parameter`), which a
+        single-level pass leaves as raw, non-JSON-serialisable objects.
+        """
+        if isinstance(value, pybamm.Symbol):
+            return convert_symbol_to_json(value)
+        if isinstance(value, dict):
+            return {k: Serialise._encode_geometry_value(v) for k, v in value.items()}
+        if isinstance(value, (list, tuple)):
+            return [Serialise._encode_geometry_value(v) for v in value]
+        return value
+
+    @staticmethod
+    def _decode_geometry_value(value):
+        """Inverse of :meth:`_encode_geometry_value`: rebuild :class:`pybamm.Symbol`
+        objects from their JSON form at any nesting depth."""
+        if isinstance(value, dict):
+            if "$type" in value or "type" in value:
+                return convert_symbol_from_json(value)
+            return {k: Serialise._decode_geometry_value(v) for k, v in value.items()}
+        if isinstance(value, list):
+            return [Serialise._decode_geometry_value(v) for v in value]
+        return value
+
     @staticmethod
     def serialise_custom_geometry(geometry: pybamm.Geometry) -> dict:
         """
@@ -654,26 +692,15 @@ def serialise_custom_geometry(geometry: pybamm.Geometry) -> dict:
                     geometry_dict_serialized[domain]["symbol_" + key_str] = (
                         convert_symbol_to_json(key)
                     )
-                    # Serialize the value dict
-                    serialized_value = {}
-                    for k, v in value.items():
-                        if isinstance(v, pybamm.Symbol):
-                            serialized_value[k] = convert_symbol_to_json(v)
-                        else:
-                            serialized_value[k] = v
-                    geometry_dict_serialized[domain][key_str] = serialized_value
+                    geometry_dict_serialized[domain][key_str] = (
+                        Serialise._encode_geometry_value(value)
+                    )
                 elif isinstance(key, str):
-                    # String keys (like 'tabs') - keep as is
-                    if isinstance(value, dict):
-                        serialized_value = {}
-                        for k, v in value.items():
-                            if isinstance(v, pybamm.Symbol):
-                                serialized_value[k] = convert_symbol_to_json(v)
-                            else:
-                                serialized_value[k] = v
-                        geometry_dict_serialized[domain][key] = serialized_value
-                    else:
-                        geometry_dict_serialized[domain][key] = value
+                    # String keys (like 'tabs'), whose values may nest Symbols
+                    # arbitrarily deep (e.g. tabs -> negative -> z_centre).
+                    geometry_dict_serialized[domain][key] = (
+                        Serialise._encode_geometry_value(value)
+                    )
 
         SCHEMA_VERSION = "1.1"
         geometry_json = {
@@ -794,29 +821,71 @@ def load_custom_geometry(filename: str | dict) -> pybamm.Geometry:
                 if key in symbol_keys:
                     # Use the reconstructed SpatialVariable as key
                     spatial_var = symbol_keys[key]
-                    reconstructed_value = {}
-                    for k, v in value.items():
-                        if isinstance(v, dict) and ("$type" in v or "type" in v):
-                            # Reconstruct PyBaMM Symbol using convert_symbol_from_json
-                            reconstructed_value[k] = convert_symbol_from_json(v)
-                        else:
-                            reconstructed_value[k] = v
-                    reconstructed_geometry[domain][spatial_var] = reconstructed_value
+                    reconstructed_geometry[domain][spatial_var] = (
+                        Serialise._decode_geometry_value(value)
+                    )
                 else:
-                    # String key (like 'tabs')
-                    if isinstance(value, dict):
-                        reconstructed_value = {}
-                        for k, v in value.items():
-                            if isinstance(v, dict) and ("$type" in v or "type" in v):
-                                reconstructed_value[k] = convert_symbol_from_json(v)
-                            else:
-                                reconstructed_value[k] = v
-                        reconstructed_geometry[domain][key] = reconstructed_value
-                    else:
-                        reconstructed_geometry[domain][key] = value
+                    # String key (like 'tabs'), values may nest Symbols deep
+                    reconstructed_geometry[domain][key] = (
+                        Serialise._decode_geometry_value(value)
+                    )
 
         return pybamm.Geometry(reconstructed_geometry)
 
+    @staticmethod
+    def _serialise_discretisation(model) -> dict:
+        """Serialise a model's discretisation recipe (geometry, var_pts, submesh
+        types, spatial methods) into a JSON-serialisable dict.
+
+        Returns an empty dict when the model exposes no discretisation defaults
+        (e.g. a plain ``pybamm.BaseModel``), so models without custom defaults
+        round-trip exactly as before.
+        """
+        out: dict = {}
+        geometry = getattr(model, "default_geometry", None) or {}
+        if geometry:
+            out["geometry"] = Serialise.serialise_custom_geometry(
+                pybamm.Geometry(geometry)
+            )
+        var_pts = getattr(model, "default_var_pts", None) or {}
+        if var_pts:
+            out["var_pts"] = Serialise.serialise_var_pts(var_pts)
+        submesh_types = getattr(model, "default_submesh_types", None) or {}
+        if submesh_types:
+            out["submesh_types"] = Serialise.serialise_submesh_types(submesh_types)
+        spatial_methods = getattr(model, "default_spatial_methods", None) or {}
+        if spatial_methods:
+            out["spatial_methods"] = Serialise.serialise_spatial_methods(
+                spatial_methods
+            )
+        return out
+
+    @staticmethod
+    def _restore_discretisation(model, discretisation: dict | None) -> None:
+        """Restore a discretisation recipe (from :meth:`_serialise_discretisation`)
+        onto ``model``, so its ``default_*`` properties return the saved values
+        even when the defining subclass could not be imported on load.
+
+        No-op when ``discretisation`` is falsy (older payloads or models without
+        custom discretisation defaults).
+        """
+        if not discretisation:
+            return
+        if "geometry" in discretisation:
+            model._default_geometry = dict(
+                Serialise.load_custom_geometry(discretisation["geometry"])
+            )
+        if "var_pts" in discretisation:
+            model._default_var_pts = Serialise.load_var_pts(discretisation["var_pts"])
+        if "submesh_types" in discretisation:
+            model._default_submesh_types = Serialise.load_submesh_types(
+                discretisation["submesh_types"]
+            )
+        if "spatial_methods" in discretisation:
+            model._default_spatial_methods = Serialise.load_spatial_methods(
+                discretisation["spatial_methods"]
+            )
+
     @staticmethod
     def serialise_spatial_method_item(method) -> dict:
         """Serialise a spatial method. The class is encoded via the kernel's
@@ -1494,6 +1563,10 @@ def load_custom_model(filename: str | dict) -> pybamm.BaseModel:
         # Restore observable state
         model._solution_observable = False
 
+        # Restore the discretisation recipe onto the (possibly fallback) model so
+        # it stays discretisable even when the defining subclass wasn't importable.
+        Serialise._restore_discretisation(model, model_data.get("discretisation"))
+
         return model
 
     @staticmethod
diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index 1a3bcd0513..6432cc39bd 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -688,22 +688,22 @@ def geometry(self):
     @property
     def default_var_pts(self):
         """Returns a dictionary of the default variable points for the model, which is empty by default."""
-        return {}
+        return getattr(self, "_default_var_pts", None) or {}
 
     @property
     def default_geometry(self):
         """Returns a  dictionary of the default geometry for the model, which is empty by default."""
-        return {}
+        return getattr(self, "_default_geometry", None) or {}
 
     @property
     def default_submesh_types(self):
         """Returns a dictionary of the default submesh types for the model, which is empty by default."""
-        return {}
+        return getattr(self, "_default_submesh_types", None) or {}
 
     @property
     def default_spatial_methods(self):
         """Returns a dictionary of the default spatial methods for the model, which is empty by default."""
-        return {}
+        return getattr(self, "_default_spatial_methods", None) or {}
 
     @property
     def default_solver(self):
diff --git a/tests/unit/test_serialisation/test_serialisation.py b/tests/unit/test_serialisation/test_serialisation.py
index c1f7ba4455..f7a966ab21 100644
--- a/tests/unit/test_serialisation/test_serialisation.py
+++ b/tests/unit/test_serialisation/test_serialisation.py
@@ -844,6 +844,112 @@ def test_load_legacy_payload_without_base_class_mro(self, tmp_path):
 
         assert type(loaded) is pybamm.BaseModel
 
+    def test_load_restores_discretisation_after_base_model_fallback(self, tmp_path):
+        """A custom model's discretisation recipe (geometry / var_pts / submesh
+        types / spatial methods) survives a load that falls back to
+        ``pybamm.BaseModel`` because the defining subclass can't be imported."""
+
+        class _DiscModel(pybamm.BaseModel):
+            @property
+            def default_geometry(self):
+                return {
+                    "positive particle": {
+                        "r_p": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1e-5)}
+                    },
+                    "current collector": {"z": {"position": 1}},
+                }
+
+            @property
+            def default_var_pts(self):
+                return {"r_p": 20}
+
+            @property
+            def default_submesh_types(self):
+                return {
+                    "positive particle": pybamm.Uniform1DSubMesh,
+                    "current collector": pybamm.SubMesh0D,
+                }
+
+            @property
+            def default_spatial_methods(self):
+                return {
+                    "positive particle": pybamm.FiniteVolume(),
+                    "current collector": pybamm.ZeroDimensionalSpatialMethod(),
+                }
+
+        model = _DiscModel(name="DiscModel")
+        a = pybamm.Variable("a")
+        model.rhs = {a: a}
+        model.initial_conditions = {a: pybamm.Scalar(1)}
+        model.variables = {"a": a}
+
+        file_path = tmp_path / "model.json"
+        Serialise.save_custom_model(model, filename=str(file_path))
+
+        # The discretisation recipe is captured at serialise time...
+        with open(file_path) as f:
+            data = json.load(f)
+        assert "discretisation" in data["model"]
+        # ...and force the import-failed fallback to a bare pybamm.BaseModel.
+        data["model"]["base_class"] = "nonexistent_module.DiscModel"
+        data["model"]["base_class_mro"] = ["nonexistent_module.DiscModel"]
+        with open(file_path, "w") as f:
+            json.dump(data, f)
+
+        with pytest.warns(UserWarning, match=r"Falling back to pybamm\.BaseModel"):
+            loaded = Serialise.load_custom_model(str(file_path))
+
+        assert type(loaded) is pybamm.BaseModel
+        # The fallback model is still discretisable: defaults are restored.
+        assert list(loaded.default_geometry.keys()) == [
+            "positive particle",
+            "current collector",
+        ]
+        assert loaded.default_var_pts == {"r_p": 20}
+        assert set(loaded.default_submesh_types) == {
+            "positive particle",
+            "current collector",
+        }
+        assert set(loaded.default_spatial_methods) == {
+            "positive particle",
+            "current collector",
+        }
+
+    def test_plain_custom_model_has_no_discretisation_section(self, tmp_path):
+        """A model with no custom discretisation defaults emits no
+        ``discretisation`` section, so existing payloads are unaffected."""
+        model = pybamm.BaseModel(name="Plain")
+        a = pybamm.Variable("a")
+        model.rhs = {a: -a}
+        model.initial_conditions = {a: pybamm.Scalar(1)}
+        model.variables = {"a": a}
+
+        payload = Serialise.serialise_custom_model(model)
+        assert "discretisation" not in payload["model"]
+        # Round-trips without error.
+        Serialise.load_custom_model(payload)
+
+    def test_custom_geometry_round_trip_with_nested_tab_symbols(self):
+        """Current-collector geometry nests symbols arbitrarily deep
+        (``tabs -> negative -> z_centre`` is a Parameter). The geometry must
+        JSON-round-trip without leaking raw, non-serialisable Symbol objects."""
+        options = {
+            "dimensionality": 1,
+            "current collector": "potential pair",
+            "cell geometry": "pouch",
+        }
+        model = pybamm.lithium_ion.DFN(dict(options), build=False)
+        geometry = pybamm.Geometry(model.default_geometry)
+
+        serialised = Serialise.serialise_custom_geometry(geometry)
+        # The whole payload must be JSON-serialisable (no raw Symbols left).
+        serialised = json.loads(json.dumps(serialised, default=Serialise._json_encoder))
+        loaded = Serialise.load_custom_geometry(serialised)
+
+        z_centre = loaded["current collector"]["tabs"]["negative"]["z_centre"]
+        assert isinstance(z_centre, pybamm.Parameter)
+        assert z_centre.name == "Negative tab centre z-coordinate [m]"
+
     def test_function_parameter_with_diff_variable_serialisation(self):
         x = pybamm.Variable("x")
         diff_var = pybamm.Variable("r")

From 96e07630fe8b4383e78dfa533b819afc7b70b027 Mon Sep 17 00:00:00 2001
From: Robert Timms <43040151+rtimms@users.noreply.github.com>
Date: Mon, 8 Jun 2026 10:44:01 +0100
Subject: [PATCH 5/6] docs: fix two dead links breaking linkcheck CI (#5575)

* docs: fix dead links breaking linkcheck CI

Two notebook links 404 after the docs site reorganised, failing the
linkChecker job and the Read the Docs linkcheck pre_build step on every
branch:

- source/api/simulation.html -> source/api/simulation/simulation.html
- installation/gnu-linux-mac.html#optional-jaxsolver ->
  installation/index.html#optional-jaxsolver

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

* fix: ignore bpxstandard.com in lychee (WAF false positive)

The site is live but its LiteSpeed anti-bot WAF returns 415 to CI/datacenter
IPs, breaking the link checker even though the link is valid.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

* ci: re-trigger builds (transient GitHub 504s)

No code changes. Earlier failures (lychee, test-data downloads, RTD
sphinx-linkcheck) were all transient github.com 504 Gateway Timeouts on
live links, not real breakages.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
(cherry picked from commit cf4ae1e4c123e27aa62c39e912e25f63487502a6)
---
 .lycheeignore                                                 | 4 ++++
 .../notebooks/performance/01-simulation-pipeline.ipynb        | 2 +-
 .../examples/notebooks/solvers/idaklu-jax-interface.ipynb     | 2 +-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/.lycheeignore b/.lycheeignore
index 5b74eceeae..8fa7101cd8 100644
--- a/.lycheeignore
+++ b/.lycheeignore
@@ -16,3 +16,7 @@ file:///home/runner/work/PyBaMM/PyBaMM/docs/source/user_guide/fundamentals/pybam
 
 # Telemetry
 https://us.i.posthog.com
+
+# Live site behind a LiteSpeed anti-bot WAF that returns 415 to CI/datacenter
+# IPs (works fine from browsers and locally) — false positive, not a dead link
+https://bpxstandard.com/
diff --git a/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb b/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb
index ca561750fb..33bca4759c 100644
--- a/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb
+++ b/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb
@@ -48,7 +48,7 @@
     "4. **Solver solve** - The next step is to solve the equations. This is done by the numerical solver, which in most cases is one of the solvers provided by the [Sundials](https://sundials.readthedocs.io/en/latest/#) suite of solvers. The solver takes the CasADi functions generated in the previous step and uses them to solve the equations numerically.\n",
     "5. **Post-processing** - The output of the previous step is the solution to the equations at a set of time points. The final step is to post-process this solution to get the variables of interest. This is done by the [`pybamm.Solution`](https://docs.pybamm.org/en/stable/source/api/solvers/solution.html) and [`pybamm.ParameterValues`](https://docs.pybamm.org/en/stable/source/api/parameters/parameter_values.html#pybamm.ParameterValues) classes, which take the solution and the parameter values and return the variables of interest interpolated at the desired time points.\n",
     "\n",
-    "For every simulation you run with PyBaMM, these steps are always performed in order to generate each solution, even if they are hidden from you by various convenience classes like [`pybamm.Simulation`](https://docs.pybamm.org/en/stable/source/api/simulation.html) . Understanding the pipeline is important because it can help you identify where the performance bottlenecks are in your simulations, and how you can optimise them to get the best performance.\n",
+    "For every simulation you run with PyBaMM, these steps are always performed in order to generate each solution, even if they are hidden from you by various convenience classes like [`pybamm.Simulation`](https://docs.pybamm.org/en/stable/source/api/simulation/simulation.html) . Understanding the pipeline is important because it can help you identify where the performance bottlenecks are in your simulations, and how you can optimise them to get the best performance.\n",
     "\n",
     "Below is a simple example that demonstrates the different parts of the pipeline, using the DFN model. While you probably have seen something similar to this script before, the one difference is that we have split the solver setup (step 3) and solve (step 4) parts of the pipeline into two separate steps. Often these are combined into the first call to the solver, but we have split them here to make it clear that they are separate steps."
    ]
diff --git a/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb b/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb
index c77b1596a2..101081f749 100644
--- a/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb
+++ b/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb
@@ -8,7 +8,7 @@
    "source": [
     "# IDAKLU-JAX interface\n",
     "\n",
-    "The IDAKLU-JAX interface requires that PyBaMM is installed with the [optional JAX solver enabled](https://docs.pybamm.org/en/stable/source/user_guide/installation/gnu-linux-mac.html#optional-jaxsolver) (`pip install pybamm[jax]`) and requires at least Python 3.10.\n",
+    "The IDAKLU-JAX interface requires that PyBaMM is installed with the [optional JAX solver enabled](https://docs.pybamm.org/en/stable/source/user_guide/installation/index.html#optional-jaxsolver) (`pip install pybamm[jax]`) and requires at least Python 3.10.\n",
     "\n",
     "PyBaMM provides two mechanisms to interface battery models with JAX. The first (JaxSolver) implements PyBaMM models directly in native JAX, and as such provides the greatest flexibility. However, these models can be very slow to compile, especially during their initial run, and can require large amounts of memory.\n",
     "\n",

From be4a7aa2bd019b7885cbcf18ad9807b0569bd6c8 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Mon, 15 Jun 2026 15:13:47 -0700
Subject: [PATCH 6/6] Update changelog and CITATION.cff for v26.6.1.2

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md | 3 +++
 CITATION.cff | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1c9c12711f..e477d3c6b7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,10 +1,13 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+# [v26.6.1.2](https://github.com/pybamm-team/PyBaMM/tree/v26.6.1.2) - 2026-06-15
+
 ## Bug fixes
 
 - `RegulariseSqrtAndPower` no longer regularises state-independent bases, fixing corrupted small rate constants in exchange-current density functions. ([#5600](https://github.com/pybamm-team/PyBaMM/pull/5600))
 - Fixed `Serialise.load_custom_model` leaving a loaded lithium-ion model's cached `param` built from default options instead of the restored ones, so a model loaded with non-default options (e.g. composite `"particle phases"`) had parameters inconsistent with its options. ([#5599](https://github.com/pybamm-team/PyBaMM/pull/5599))
 - Fixed `BatteryModelOptions` letting two invalid MSMR option sets pass validation and then crash during parameter construction with `invalid literal for int() with base 10: 'none'`; both now raise a clear `OptionError`. The all-or-nothing MSMR check now detects MSMR inside a per-electrode tuple (e.g. `"open-circuit potential": ("MSMR", "single")`), and an MSMR model must set `"number of MSMR reactions"` to a positive integer rather than the default `"none"`. ([#5599](https://github.com/pybamm-team/PyBaMM/pull/5599))
+- Fixed `Serialise.serialise_custom_model` not recording a custom model's discretisation recipe (`geometry`, `var_pts`, `submesh_types`, `spatial_methods`), so a custom model reloaded via the `BaseModel` fallback (when its defining package is unavailable) was no longer discretisable. ([#5609](https://github.com/pybamm-team/PyBaMM/pull/5609))
 
 # [v26.6.1.1](https://github.com/pybamm-team/PyBaMM/tree/v26.6.1.1) - 2026-06-11
 
diff --git a/CITATION.cff b/CITATION.cff
index cacb1087a9..f5c599bafc 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -24,6 +24,6 @@ keywords:
   - "expression tree"
   - "python"
   - "symbolic differentiation"
-version: "26.6.1.1"
+version: "26.6.1.2"
 repository-code: "https://github.com/pybamm-team/PyBaMM"
 title: "Python Battery Mathematical Modelling (PyBaMM)"