From 85d1540b046bea81c63edc253e7feaaf12f99b39 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Fri, 16 Oct 2026 10:34:06 -0400
Subject: [PATCH 1/8] perf: faster model build and IDAKLU set-up for large
 models

- StateVector stores its slices as index ranges instead of a boolean mask.
- CasADi functions and matrices are shared with pybammsolvers instead of copied.
- IDAKLU no longer builds the rhs-only jacobian, which only ScipySolver uses.
- Unstructured meshes gather faces per vertex.
---
 CHANGELOG.md                                  |   6 +
 .../src/pybamm/expression_tree/array.py       |  50 +++++-
 .../pybamm/expression_tree/concatenations.py  |  20 +--
 .../pybamm/expression_tree/input_parameter.py |   2 +-
 .../pybamm/expression_tree/state_vector.py    | 152 ++++++++++++-----
 .../src/pybamm/expression_tree/symbol.py      |   4 +-
 .../src/pybamm/meshes/unstructured_submesh.py | 160 ++++++++++--------
 .../pybamm/src/pybamm/solvers/base_solver.py  |  17 +-
 .../src/pybamm/solvers/idaklu_solver.py       |  35 ++--
 .../src/pybamm/solvers/nonlinear_solver.py    |   6 +-
 .../pybamm/src/pybamm/solvers/scipy_solver.py |   2 +
 .../unit/test_expression_tree/test_array.py   |  30 ++++
 .../test_expression_tree/test_state_vector.py |  87 ++++++----
 .../test_meshes/test_unstructured_submesh.py  |  11 ++
 .../unit/test_solvers/test_base_solver.py     |  26 +++
 packages/pybammsolvers/CMakeLists.txt         |   6 +
 .../src/pybammsolvers/idaklu.cpp              |  51 ++++++
 17 files changed, 479 insertions(+), 186 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6403e6172a..d63e0998e3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,9 +1,14 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Breaking changes
+
+- `StateVector` and `StateVectorDot` raise `ShapeError` for slices that have a step, overlap or are out of order. Out-of-order slices were previously read only up to the last slice's stop. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
+
 ## Deprecated
 
 - `FiniteVolume2D` is deprecated in favour of `FiniteVolumeUnstructured` and raises a `DeprecationWarning` when created. Remaining gaps before removal are tracked in [#5938](https://github.com/pybamm-team/PyBaMM/issues/5938). ([#5939](https://github.com/pybamm-team/PyBaMM/pull/5939))
 - `pybamm.Magnitude` is deprecated and will be removed in 27.1; use `pybamm.Component(symbol, 0)` / `pybamm.Component(symbol, 1)` for the lr/tb components of a `VectorField`. ([#5943](https://github.com/pybamm-team/PyBaMM/pull/5943))
+- The `evaluation_array` argument of `StateVector` and `StateVectorDot` is deprecated and ignored; the entries read are given by the slices. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
 
 ## Features
 
@@ -32,6 +37,7 @@
 
 - Read the Docs pull request builds drop `linkcheck` and the unpublished PDF, install the released solver wheel instead of compiling it, and are skipped when nothing the docs render changed, saving about 14 of 20 minutes. ([#5986](https://github.com/pybamm-team/PyBaMM/pull/5986))
 - Pull request CI no longer installs TeXLive, or the compilers and OpenBLAS that the prebuilt solver wheels make redundant. The latexify test mocks `sympy.preview` and runs everywhere. TeXLive downloads took 12–50 minutes in 6 of 15 runs. ([#5985](https://github.com/pybamm-team/PyBaMM/pull/5985))
+- Faster model build and `IDAKLUSolver` set-up for large (million-state) models. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
 - Experiments with a different model per step no longer compile Jacobians of the state mappers between models, which were unused. ([#5918](https://github.com/pybamm-team/PyBaMM/pull/5918))
 
 # [v26.10.1.1](https://github.com/pybamm-team/PyBaMM/tree/pybamm-v26.10.1.1) - 2026-10-19
diff --git a/packages/pybamm/src/pybamm/expression_tree/array.py b/packages/pybamm/src/pybamm/expression_tree/array.py
index 55d9ba2cd0..c7c548245d 100644
--- a/packages/pybamm/src/pybamm/expression_tree/array.py
+++ b/packages/pybamm/src/pybamm/expression_tree/array.py
@@ -7,7 +7,8 @@
 import numpy as np
 import numpy.typing as npt
 import sympy
-from scipy.sparse import csr_matrix, issparse
+from pybammsolvers import idaklu
+from scipy.sparse import csc_matrix, csr_matrix, issparse
 
 import pybamm
 from pybamm.type_definitions import AuxiliaryDomainType, DomainsType, DomainType
@@ -132,12 +133,12 @@ def _make_entries_string(self, value: tuple | None) -> tuple:
 
     def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
         """See :meth:`pybamm.Symbol._to_casadi()`."""
-        return casadi.MX(self.evaluate(t, y, y_dot, inputs))
+        return casadi.MX(to_casadi_dm(self.evaluate(t, y, y_dot, inputs)))
 
     def _jac(self, variable) -> pybamm.Matrix:
         """See :meth:`pybamm.Symbol._jac()`."""
         # Return zeros of correct size
-        jac = csr_matrix((self.size, variable.evaluation_array.count(True)))
+        jac = csr_matrix((self.size, variable.size))
         return pybamm.Matrix(jac)
 
     def create_copy(
@@ -217,3 +218,46 @@ def meshgrid(
     """
     [x_grid, y_grid] = np.meshgrid(x.entries, y.entries)
     return pybamm.Array(x_grid), pybamm.Array(y_grid)
+
+
+def to_casadi_dm(entries) -> casadi.DM:
+    """
+    Convert a numpy array or scipy sparse matrix to a casadi.DM.
+
+    When pybammsolvers links the Python casadi library, the buffers are copied
+    once in C++; casadi's own conversion copies element by element.
+
+    Parameters
+    ----------
+    entries : numpy.ndarray or scipy.sparse matrix
+        The values, at most two-dimensional; a 1D array becomes a column.
+
+    Returns
+    -------
+    casadi.DM
+        The values, sparse if ``entries`` is sparse.
+    """
+    if not idaklu.shares_python_casadi:
+        return casadi.DM(entries)
+    matrix = casadi.DM()
+    if issparse(entries):
+        csc = csc_matrix(entries)
+        if not csc.has_canonical_format:
+            csc = csc.copy()
+            csc.sum_duplicates()
+        idaklu.fill_casadi_dm_sparse(
+            int(matrix.this),
+            csc.shape[0],
+            csc.shape[1],
+            csc.indptr.astype(np.int64),
+            csc.indices.astype(np.int64),
+            csc.data.astype(np.float64),
+        )
+    else:
+        values = np.asarray(entries, dtype=np.float64)
+        if values.ndim < 2:
+            values = values.reshape(-1, 1)
+        idaklu.fill_casadi_dm_dense(
+            int(matrix.this), values.shape[0], values.shape[1], values.ravel(order="F")
+        )
+    return matrix
diff --git a/packages/pybamm/src/pybamm/expression_tree/concatenations.py b/packages/pybamm/src/pybamm/expression_tree/concatenations.py
index aac75ca33e..06e233b9ff 100644
--- a/packages/pybamm/src/pybamm/expression_tree/concatenations.py
+++ b/packages/pybamm/src/pybamm/expression_tree/concatenations.py
@@ -695,22 +695,18 @@ def simplified_domain_concatenation(
     # Create the DomainConcatenation to read domain and child domain
     concat = DomainConcatenation(children, mesh, copy_this=copy_this)
     # Simplify Concatenation of StateVectors to a single StateVector
-    # The sum of the evalation arrays of the StateVectors must be exactly 1
+    # The StateVectors must cover their combined span exactly once
     if all(isinstance(child, pybamm.StateVector) for child in children):
         sv_children: list[pybamm.StateVector] = children  # type: ignore[assignment]
-        longest_eval_array = len(sv_children[-1]._evaluation_array)
-        eval_arrays = {}
-        for child in sv_children:
-            eval_arrays[child] = np.concatenate(
-                [
-                    child.evaluation_array,
-                    np.zeros(longest_eval_array - len(child.evaluation_array)),
-                ]
-            )
         first_start = sv_children[0].y_slices[0].start
         last_stop = sv_children[-1].y_slices[-1].stop
-        if all(
-            sum(array for array in eval_arrays.values())[first_start:last_stop] == 1
+        ranges = sorted(
+            (y_slice.start or 0, y_slice.stop)
+            for child in sv_children
+            for y_slice in child.y_slices
+        )
+        if pybamm.expression_tree.state_vector.ranges_tile(
+            ranges, first_start, last_stop
         ):
             return pybamm.StateVector(
                 slice(first_start, last_stop), domains=concat._domains
diff --git a/packages/pybamm/src/pybamm/expression_tree/input_parameter.py b/packages/pybamm/src/pybamm/expression_tree/input_parameter.py
index 639568d314..5e2caf08c6 100644
--- a/packages/pybamm/src/pybamm/expression_tree/input_parameter.py
+++ b/packages/pybamm/src/pybamm/expression_tree/input_parameter.py
@@ -88,7 +88,7 @@ def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
 
     def _jac(self, variable: pybamm.StateVector) -> pybamm.Matrix:
         """See :meth:`pybamm.Symbol._jac()`."""
-        n_variable = variable.evaluation_array.count(True)
+        n_variable = variable.size
         nan_vector = self._evaluate_for_shape()
         if isinstance(nan_vector, numbers.Number):
             n_self = 1
diff --git a/packages/pybamm/src/pybamm/expression_tree/state_vector.py b/packages/pybamm/src/pybamm/expression_tree/state_vector.py
index 3c0603a688..bc5269ef93 100644
--- a/packages/pybamm/src/pybamm/expression_tree/state_vector.py
+++ b/packages/pybamm/src/pybamm/expression_tree/state_vector.py
@@ -3,6 +3,8 @@
 #
 from __future__ import annotations
 
+import warnings
+
 import casadi
 import numpy as np
 import numpy.typing as npt
@@ -12,6 +14,30 @@
 from pybamm.type_definitions import AuxiliaryDomainType, DomainsType, DomainType
 
 
+def ranges_tile(ranges, start: int, stop: int) -> bool:
+    """
+    Whether sorted ``(start, stop)`` ranges cover ``[start, stop)`` exactly once.
+
+    Parameters
+    ----------
+    ranges : iterable of tuple of int
+        Half-open index ranges, sorted by start.
+    start, stop : int
+        The interval to cover.
+
+    Returns
+    -------
+    bool
+        True if the ranges are contiguous, non-overlapping and span the interval.
+    """
+    position = start
+    for range_start, range_stop in ranges:
+        if range_start != position:
+            return False
+        position = range_stop
+    return position == stop
+
+
 class StateVectorBase(pybamm.Symbol):
     """
     Node in the expression tree that holds a slice to read from an external vector type.
@@ -33,13 +59,17 @@ class StateVectorBase(pybamm.Symbol):
         (not both). In future, the 'domain' and 'auxiliary_domains' arguments may be
         deprecated.
     evaluation_array : list, optional
-        List of boolean arrays representing slices. Default is None, in which case the
-        evaluation_array is computed from y_slices.
+        Deprecated and ignored: the entries read are always given by ``y_slices``.
+
+    Raises
+    ------
+    :class:`pybamm.ShapeError`
+        If a slice has a step, or the slices are not in increasing order or overlap.
     """
 
-    __slots__ = ("_evaluation_array", "_first_point", "_last_point", "_y_slices")
-    # the evaluation array determines the slices and end points
-    _id_excluded_fields = ("_y_slices", "_first_point", "_last_point")
+    __slots__ = ("_first_point", "_index_ranges", "_last_point", "_y_slices")
+    # The merged index ranges define the entries read, however the slices split them
+    _id_excluded_fields = ("_first_point", "_last_point", "_y_slices")
 
     def __init__(
         self,
@@ -54,6 +84,13 @@ def __init__(
         for y_slice in y_slices:
             if not isinstance(y_slice, slice):
                 raise TypeError("all y_slices must be slice objects")
+        if evaluation_array is not None:
+            warnings.warn(
+                "The 'evaluation_array' argument of StateVector is deprecated and "
+                "ignored; the entries read are given by the y_slices.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
         if name is None:
             if y_slices[0].start is None:
                 name = base_name + f"[0:{y_slice.stop:d}"
@@ -70,7 +107,7 @@ def __init__(
         self._y_slices = y_slices
         self._first_point = y_slices[0].start
         self._last_point = y_slices[-1].stop
-        self._evaluation_array = self._make_evaluation_array(y_slices, evaluation_array)
+        self._index_ranges = self._merge_slices(y_slices)
         super().__init__(
             name=name,
             domain=domain,
@@ -82,12 +119,7 @@ def __init__(
     def _from_json(cls, snippet: dict):
         y_slices = [slice(s["start"], s["stop"], s["step"]) for s in snippet["y_slice"]]
 
-        return cls(
-            *y_slices,
-            name=snippet["name"],
-            domains=snippet["domains"],
-            evaluation_array=snippet["evaluation_array"],
-        )
+        return cls(*y_slices, name=snippet["name"], domains=snippet["domains"])
 
     @property
     def y_slices(self):
@@ -103,35 +135,73 @@ def last_point(self):
 
     @property
     def evaluation_array(self):
-        """Array to use for evaluating."""
-        return self._evaluation_array
+        """Boolean mask, as long as ``last_point``, of the entries of y read."""
+        array = np.zeros(self._last_point, dtype=bool)
+        for start, stop in self._index_ranges:
+            array[start:stop] = True
+        return array.tolist()
 
     @property
     def size(self):
-        return self.evaluation_array.count(True)
+        return sum(stop - start for start, stop in self._index_ranges)
+
+    def _read(self, vector):
+        """The entries of ``vector`` this state vector reads, as a column."""
+        if len(self._index_ranges) == 1:
+            start, stop = self._index_ranges[0]
+            out = vector[start:stop].copy()
+        else:
+            out = np.concatenate(
+                [vector[start:stop] for start, stop in self._index_ranges]
+            )
+        if isinstance(out, np.ndarray) and out.ndim == 1:
+            out = out[:, np.newaxis]
+        return out
+
+    @staticmethod
+    def _merge_slices(y_slices) -> tuple[tuple[int, int], ...]:
+        """
+        The entries read by ``y_slices``, as ``(start, stop)`` ranges with adjacent
+        slices merged.
+
+        Raises
+        ------
+        :class:`pybamm.ShapeError`
+            If a slice has a step, or the slices are not increasing or overlap.
+        """
+        ranges: list[list[int]] = []
+        for y_slice in y_slices:
+            if y_slice.step not in (None, 1):
+                raise pybamm.ShapeError("y_slices must not have a step")
+            start = y_slice.start or 0
+            if ranges and start < ranges[-1][1]:
+                raise pybamm.ShapeError(
+                    "y_slices must be in increasing order and must not overlap"
+                )
+            if ranges and start == ranges[-1][1]:
+                ranges[-1][1] = y_slice.stop
+            elif y_slice.stop > start:
+                ranges.append([start, y_slice.stop])
+        return tuple((start, stop) for start, stop in ranges)
+
+    def __setstate__(self, state):
+        if "_evaluation_array" in state:  # pickled when the mask was stored
+            state = {
+                key: value for key, value in state.items() if key != "_evaluation_array"
+            }
+            state["_index_ranges"] = self._merge_slices(state["_y_slices"])
+        super().__setstate__(state)
 
     def set_evaluation_array(self, y_slices, evaluation_array):
-        """Set evaluation array using slices (deprecated)."""
+        """Set the slices to read (deprecated); ``evaluation_array`` is ignored."""
         pybamm.expression_tree.symbol._warn_mutation(
             "set_evaluation_array",
             "Construct a new StateVector with the desired slices.",
         )
-        object.__setattr__(
-            self,
-            "_evaluation_array",
-            self._make_evaluation_array(y_slices, evaluation_array),
-        )
-
-    @staticmethod
-    def _make_evaluation_array(y_slices, evaluation_array) -> list[bool]:
-        """The evaluation array: ``evaluation_array`` if given, else computed from
-        the slices."""
-        if evaluation_array is not None and pybamm.settings.debug_mode is False:
-            return evaluation_array
-        array = np.zeros(y_slices[-1].stop)
-        for y_slice in y_slices:
-            array[y_slice] = True
-        return [bool(x) for x in array]
+        object.__setattr__(self, "_y_slices", tuple(y_slices))
+        object.__setattr__(self, "_first_point", y_slices[0].start)
+        object.__setattr__(self, "_last_point", y_slices[-1].stop)
+        object.__setattr__(self, "_index_ranges", self._merge_slices(y_slices))
 
     def _jac_diff_vector(self, variable: pybamm.StateVectorBase):
         """
@@ -209,7 +279,6 @@ def create_copy(
             *self.y_slices,
             name=self.name,
             domains=self._domains,
-            evaluation_array=self.evaluation_array,
         )
 
     def _evaluate_for_shape(self):
@@ -298,15 +367,11 @@ def _base_evaluate(
         """See :meth:`pybamm.Symbol._base_evaluate()`."""
         if y is None:
             raise TypeError("StateVector cannot evaluate input 'y=None'")
-        if y.shape[0] < len(self.evaluation_array):
+        if y.shape[0] < self._last_point:
             raise ValueError(
                 "y is too short, so value with slice is smaller than expected"
             )
-
-        out = (y[: len(self._evaluation_array)])[self._evaluation_array]
-        if isinstance(out, np.ndarray) and out.ndim == 1:
-            out = out[:, np.newaxis]
-        return out
+        return self._read(y)
 
     def diff(self, variable: pybamm.Symbol):
         if variable == self:
@@ -316,7 +381,6 @@ def diff(self, variable: pybamm.Symbol):
                 *self._y_slices,
                 name=self.name + "'",
                 domains=self._domains,
-                evaluation_array=self.evaluation_array,
             )
         else:
             return pybamm.Scalar(0)
@@ -390,15 +454,11 @@ def _base_evaluate(
         """See :meth:`pybamm.Symbol._base_evaluate()`."""
         if y_dot is None:
             raise TypeError("StateVectorDot cannot evaluate input 'y_dot=None'")
-        if y_dot.shape[0] < len(self.evaluation_array):
+        if y_dot.shape[0] < self._last_point:
             raise ValueError(
                 "y_dot is too short, so value with slice is smaller than expected"
             )
-
-        out = (y_dot[: len(self._evaluation_array)])[self._evaluation_array]
-        if isinstance(out, np.ndarray) and out.ndim == 1:
-            out = out[:, np.newaxis]
-        return out
+        return self._read(y_dot)
 
     def diff(self, variable: pybamm.Symbol):
         if variable == self:
diff --git a/packages/pybamm/src/pybamm/expression_tree/symbol.py b/packages/pybamm/src/pybamm/expression_tree/symbol.py
index e42bc30ec3..06a1ddd389 100644
--- a/packages/pybamm/src/pybamm/expression_tree/symbol.py
+++ b/packages/pybamm/src/pybamm/expression_tree/symbol.py
@@ -1533,9 +1533,7 @@ def _compute_shape(self):
         except ValueError:
             unpacker = pybamm.SymbolUnpacker(pybamm.StateVector)
             state_vectors_in_node = unpacker.unpack_symbol(self)
-            min_y_size = max(
-                max(len(x._evaluation_array) for x in state_vectors_in_node), 1
-            )
+            min_y_size = max(max(x.last_point for x in state_vectors_in_node), 1)
             # Pick a y that won't cause RuntimeWarnings
             y = np.nan * np.ones((min_y_size, 1))
             evaluated_self = self.evaluate(0, y, y, inputs="shape test")
diff --git a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
index fb562d518e..5f0ffdf66b 100644
--- a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
+++ b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
@@ -107,13 +107,16 @@ def __init__(self, vertices, elements, coord_sys="cartesian", boundary_faces=Non
     # ------------------------------------------------------------------
 
     def _compute_cell_geometry(self):
-        verts = self.vertices[self.elements]  # (n_cells, n_verts, d)
+        # Contiguous (n_cells, d) array per local vertex
+        verts = [
+            self.vertices[self.elements[:, i]] for i in range(self.elements.shape[1])
+        ]
         # Vertex mean is the exact centroid for simplices (triangles, tets);
         # quads and hexes overwrite it with the true area/volume centroid.
-        self.cell_centroids = verts.mean(axis=1)
+        self.cell_centroids = sum(verts) / len(verts)
 
         if self.element_type == ElementType.TRIANGLE:
-            v0, v1, v2 = verts[:, 0], verts[:, 1], verts[:, 2]
+            v0, v1, v2 = verts
             cross = (v1[:, 0] - v0[:, 0]) * (v2[:, 1] - v0[:, 1]) - (
                 v1[:, 1] - v0[:, 1]
             ) * (v2[:, 0] - v0[:, 0])
@@ -121,8 +124,8 @@ def _compute_cell_geometry(self):
         elif self.element_type == ElementType.QUAD:
             # Shoelace formula for arbitrary simple quadrilaterals
             # Vertices ordered: v0, v1, v2, v3 (counterclockwise or clockwise)
-            x = verts[:, :, 0]  # (n_cells, 4)
-            y = verts[:, :, 1]  # (n_cells, 4)
+            x = np.column_stack([v[:, 0] for v in verts])  # (n_cells, 4)
+            y = np.column_stack([v[:, 1] for v in verts])  # (n_cells, 4)
             # shoelace: sum_i (x_i * y_{i+1} - x_{i+1} * y_i)
             x_next = np.roll(x, -1, axis=1)
             y_next = np.roll(y, -1, axis=1)
@@ -135,7 +138,7 @@ def _compute_cell_geometry(self):
             cy = np.sum((y + y_next) * cross, axis=1) / (6.0 * signed_area)
             self.cell_centroids = np.column_stack([cx, cy])
         elif self.element_type == ElementType.TETRAHEDRON:
-            v0, v1, v2, v3 = verts[:, 0], verts[:, 1], verts[:, 2], verts[:, 3]
+            v0, v1, v2, v3 = verts
             d1 = v1 - v0
             d2 = v2 - v0
             d3 = v3 - v0
@@ -146,20 +149,18 @@ def _compute_cell_geometry(self):
             )
             self.cell_volumes = np.abs(det) / 6.0
         elif self.element_type == ElementType.HEXAHEDRON:
-            # 5-tet decomposition, exact only for planar-faced hexes (warped
-            # faces are rejected in _compute_face_geometry). Tet volumes are
-            # summed signed so a degenerate or inverted cell cannot cancel
-            # into a plausible positive volume.
-            t = verts[:, self._HEX_TETS]  # (n_cells, 5, 4, d)
-            d1 = t[:, :, 1] - t[:, :, 0]
-            d2 = t[:, :, 2] - t[:, :, 0]
-            d3 = t[:, :, 3] - t[:, :, 0]
-            tet_vols = np.einsum("ijk,ijk->ij", d1, np.cross(d2, d3)) / 6.0
-            vol = tet_vols.sum(axis=1)
+            # 5-tet split, exact for planar-faced hexes; signed sum so an
+            # inverted cell cannot cancel into a plausible positive volume
+            vol = np.zeros(len(self.elements))
+            moment = np.zeros((len(self.elements), 3))
+            for tet in self._HEX_TETS:
+                p0, p1, p2, p3 = (verts[i] for i in tet)
+                tet_vol = _dot(p1 - p0, _cross(p2 - p0, p3 - p0)) / 6.0
+                vol += tet_vol
+                # Volume-weighted tet centroids: exact for planar-faced hexes
+                # such as frusta, where the vertex mean is not
+                moment += tet_vol[:, None] * (p0 + p1 + p2 + p3) / 4.0
             self.cell_volumes = np.abs(vol)
-            # Volume-weighted tet centroids: exact for planar-faced hexes
-            # such as frusta, where the vertex mean is not
-            moment = np.einsum("ij,ijk->ik", tet_vols, t.mean(axis=2))
             safe_vol = np.where(vol == 0.0, 1.0, vol)
             nonzero = (vol != 0.0)[:, None]
             self.cell_centroids = np.where(
@@ -178,7 +179,6 @@ def _build_face_connectivity(self):
             n_verts_per_face = self.dimension
 
         elems = self.elements
-        n_cells = len(elems)
 
         # Build all faces at once using local face definitions
         if self.element_type == ElementType.QUAD:
@@ -193,18 +193,10 @@ def _build_face_connectivity(self):
         n_fpc = len(local)
 
         all_faces = elems[:, local].reshape(-1, n_verts_per_face)
-        cell_ids = np.repeat(np.arange(n_cells), n_fpc)
-
-        # Canonical keys: sort vertex indices within each face
-        sorted_faces = np.sort(all_faces, axis=1)
 
-        # Find unique faces and which are shared (internal) vs single (boundary)
-        _, inverse, counts = np.unique(
-            sorted_faces, axis=0, return_inverse=True, return_counts=True
-        )
-        # numpy 2.0.0 (the declared floor) returns inverse with shape (n, 1);
-        # 2.0.1+ returns (n,). Flatten so masks below stay 1-D everywhere.
-        inverse = inverse.reshape(-1)
+        # Canonical keys: sort vertex indices within each face, then group
+        # duplicates; shared faces are internal, single faces boundary
+        order, counts = _group_rows(np.sort(all_faces, axis=1))
 
         # A manifold mesh shares each face between at most two cells. A count of
         # three or more means overlapping or non-conforming elements; such faces
@@ -216,28 +208,18 @@ def _build_face_connectivity(self):
                 "elements in the input mesh."
             )
 
-        is_internal = counts[inverse] == 2
-        is_boundary = counts[inverse] == 1
-
-        # For internal faces, we need owner/neighbor pairs.
-        # Group by unique face index; first occurrence is owner, second is neighbor.
-        internal_mask = is_internal
-        int_inv = inverse[internal_mask]
-        int_cells = cell_ids[internal_mask]
-        int_faces_raw = all_faces[internal_mask]
+        # The grouping is stable, so each internal face's two rows are adjacent
+        # with the lower-numbered (owner) cell first
+        internal_rows = order[counts == 2]
+        owner_rows = internal_rows[0::2]
+        internal_owner = owner_rows // n_fpc
+        internal_neighbor = internal_rows[1::2] // n_fpc
+        internal_face_verts = all_faces[owner_rows]
 
-        # Sort by unique-face-id to pair them up: [owner0, neighbor0, owner1, neighbor1, ...]
-        order = np.argsort(int_inv, kind="stable")
-        int_cells_sorted = int_cells[order]
-        int_faces_sorted = int_faces_raw[order]
-
-        internal_owner = int_cells_sorted[0::2]
-        internal_neighbor = int_cells_sorted[1::2]
-        internal_face_verts = int_faces_sorted[0::2]
-
-        # Boundary faces
-        bnd_face_verts = all_faces[is_boundary]
-        bnd_owners = cell_ids[is_boundary]
+        # Boundary faces keep their original (cell-major) order
+        boundary_rows = np.sort(order[counts == 1])
+        bnd_face_verts = all_faces[boundary_rows]
+        bnd_owners = boundary_rows // n_fpc
 
         n_internal = len(internal_owner)
         n_boundary = len(bnd_owners)
@@ -272,32 +254,31 @@ def _build_face_connectivity(self):
     # ------------------------------------------------------------------
 
     def _compute_face_geometry(self):
-        face_verts = self.vertices[self.faces]
-
-        self.face_centroids = face_verts.mean(axis=1)
+        # Contiguous per-vertex gathers; strided views of an (n, k, d) gather
+        # make every later vector operation slower
+        face_verts = [
+            self.vertices[self.faces[:, i]] for i in range(self.faces.shape[1])
+        ]
+        self.face_centroids = sum(face_verts) / len(face_verts)
 
         if self.dimension == 2:
-            v0, v1 = face_verts[:, 0], face_verts[:, 1]
+            v0, v1 = face_verts
             edge = v1 - v0
             self.face_areas = np.linalg.norm(edge, axis=1)
             normals = np.column_stack([edge[:, 1], -edge[:, 0]])
         elif self.element_type == ElementType.HEXAHEDRON:
             # Face = quad: 4 vertices. Area via cross product of diagonals.
-            v0 = face_verts[:, 0]
-            v1 = face_verts[:, 1]
-            v2 = face_verts[:, 2]
-            v3 = face_verts[:, 3]
+            v0, v1, v2, v3 = face_verts
 
             # Quad areas, normals, and the 5-tet cell volumes are only
             # well-defined when all four vertices of a face are coplanar,
             # so warped (non-planar) faces are rejected outright.
-            plane_normal = np.cross(v1 - v0, v2 - v0)
-            plane_norm = np.linalg.norm(plane_normal, axis=1)
+            plane_normal = _cross(v1 - v0, v2 - v0)
+            plane_norm = np.sqrt(_dot(plane_normal, plane_normal))
             safe_norm = np.where(plane_norm < 1e-30, 1.0, plane_norm)
-            offset = np.abs(
-                np.einsum("ij,ij->i", v3 - v0, plane_normal / safe_norm[:, None])
-            )
-            diag_len = np.linalg.norm(v2 - v0, axis=1)
+            n_hat = plane_normal / safe_norm[:, None]
+            offset = np.abs(_dot(v3 - v0, n_hat))
+            diag_len = np.sqrt(_dot(v2 - v0, v2 - v0))
             warped = offset > 1e-8 * np.maximum(diag_len, 1e-30)
             if warped.any():
                 raise pybamm.GeometryError(
@@ -309,17 +290,16 @@ def _compute_face_geometry(self):
 
             diag1 = v2 - v0
             diag2 = v3 - v1
-            cross = np.cross(diag1, diag2)
-            self.face_areas = 0.5 * np.linalg.norm(cross, axis=1)
+            cross = _cross(diag1, diag2)
+            self.face_areas = 0.5 * np.sqrt(_dot(cross, cross))
             normals = cross
 
             # Area-weighted centroid over the (0,1,2)/(0,2,3) triangle split:
             # exact for planar non-parallelogram quads (e.g. trapezoids),
             # where the vertex mean is not. Signed areas taken along the face
             # normal keep non-convex planar quads exact too.
-            n_hat = plane_normal / safe_norm[:, None]
-            area1 = 0.5 * np.einsum("ij,ij->i", np.cross(v1 - v0, v2 - v0), n_hat)
-            area2 = 0.5 * np.einsum("ij,ij->i", np.cross(v2 - v0, v3 - v0), n_hat)
+            area1 = 0.5 * _dot(plane_normal, n_hat)
+            area2 = 0.5 * _dot(_cross(v2 - v0, v3 - v0), n_hat)
             total = area1 + area2
             safe_total = np.where(np.abs(total) < 1e-30, 1.0, total)
             weighted = (
@@ -329,7 +309,7 @@ def _compute_face_geometry(self):
             self.face_centroids[nonzero] = weighted[nonzero]
         else:
             # Face = triangle: 3 vertices
-            v0, v1, v2 = face_verts[:, 0], face_verts[:, 1], face_verts[:, 2]
+            v0, v1, v2 = face_verts
             cross = np.cross(v1 - v0, v2 - v0)
             self.face_areas = 0.5 * np.linalg.norm(cross, axis=1)
             normals = cross
@@ -1531,3 +1511,39 @@ def _hex_to_tet(x_edges, y_edges, z_edges):
     elements = hex_verts[:, kuhn_tets].reshape(-1, 4)
 
     return nodes, elements
+
+
+def _cross(a, b):
+    """Row-wise cross product of two ``(n, 3)`` arrays, faster than ``np.cross``."""
+    a0, a1, a2 = a[:, 0], a[:, 1], a[:, 2]
+    b0, b1, b2 = b[:, 0], b[:, 1], b[:, 2]
+    return np.column_stack([a1 * b2 - a2 * b1, a2 * b0 - a0 * b2, a0 * b1 - a1 * b0])
+
+
+def _dot(a, b):
+    """Row-wise dot product of two ``(n, d)`` arrays."""
+    return np.einsum("ij,ij->i", a, b)
+
+
+def _group_rows(rows):
+    """
+    Stably sort the rows of a 2D array so equal rows are adjacent, and count them.
+
+    Parameters
+    ----------
+    rows : numpy.ndarray
+        2D array whose rows are compared.
+
+    Returns
+    -------
+    order : numpy.ndarray
+        Row indices in lexicographic order, ties kept in input order.
+    counts : numpy.ndarray
+        Multiplicity of the row at each position of ``order``.
+    """
+    order = np.lexsort(rows.T[::-1])
+    sorted_rows = rows[order]
+    starts_group = np.ones(len(rows), dtype=bool)
+    starts_group[1:] = (sorted_rows[1:] != sorted_rows[:-1]).any(axis=1)
+    group = np.cumsum(starts_group) - 1
+    return order, np.bincount(group)[group]
diff --git a/packages/pybamm/src/pybamm/solvers/base_solver.py b/packages/pybamm/src/pybamm/solvers/base_solver.py
index 1fc94719df..ea704bf78c 100644
--- a/packages/pybamm/src/pybamm/solvers/base_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/base_solver.py
@@ -63,6 +63,9 @@ class BaseSolver:
         Default is False.
     """
 
+    #: Whether the solver uses the jacobian of the rhs on its own
+    uses_rhs_jacobian = False
+
     def __init__(
         self,
         method=None,
@@ -279,7 +282,10 @@ def set_up(
         # Process rhs, algebraic, residual and event expressions
         # and wrap in callables
         rhs, jac_rhs, jacp_rhs, jac_rhs_action = process(
-            model.concatenated_rhs, "RHS", vars_for_processing
+            model.concatenated_rhs,
+            "RHS",
+            vars_for_processing,
+            use_jacobian=model.use_jacobian and self.uses_rhs_jacobian,
         )
 
         algebraic, jac_algebraic, jacp_algebraic, jac_algebraic_action = process(
@@ -351,7 +357,7 @@ def set_up(
                 if model.is_standard_form_dae:
                     explicit_rhs = rhs_eval
                 else:
-                    M_ode = casadi.DM(
+                    M_ode = pybamm.expression_tree.array.to_casadi_dm(
                         model.mass_matrix.entries[: model.len_rhs, : model.len_rhs]
                     )
                     explicit_rhs = casadi.solve(M_ode, rhs_eval)
@@ -2257,7 +2263,12 @@ def stack_inputs(inputs_dict: dict) -> casadi.DM:
     :class:`casadi.DM`
         Column vector of the input values, with zero rows for an empty dict.
     """
-    return casadi.vertcat(*inputs_dict.values())
+    if not inputs_dict:
+        return casadi.DM.zeros(0, 1)
+    stacked = np.concatenate(
+        [np.asarray(value, dtype=float).reshape(-1) for value in inputs_dict.values()]
+    )
+    return pybamm.expression_tree.array.to_casadi_dm(stacked)
 
 
 def process(
diff --git a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
index d8260b3755..b4fae0bd12 100644
--- a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
@@ -74,6 +74,19 @@ def _sensitivity_scales(inputs_dict: dict, sensitivity_names: list[str]) -> np.n
 # integrator has located one or more root function zeros.
 _IDA_ROOT_RETURN = 2
 
+
+def to_idaklu_function(fn: casadi.Function):
+    """
+    Hand a casadi.Function to pybammsolvers.
+
+    When pybammsolvers links the Python casadi package's library, the C++ side
+    shares the function directly; otherwise it is serialised and rebuilt.
+    """
+    if idaklu.shares_python_casadi:
+        return idaklu.shared_function(int(fn.this))
+    return idaklu.generate_function(fn.serialize())
+
+
 # Function entries staged in ``_setup`` while building the C++ solver
 # group; popped after construction (C++ owns them via shared_ptr).
 _SETUP_FCN_KEYS = (
@@ -533,10 +546,9 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
                 p_casadi[name] = casadi.MX.sym(name, value.shape[0])
         p_casadi_stacked = casadi.vertcat(*[p for p in p_casadi.values()])
 
-        jac_times_cjmass_expression = (
-            model.jac_rhs_algebraic_eval(t_casadi, y_casadi, p_casadi_stacked)
-            - cj_casadi * mass_matrix
-        )
+        jac_times_cjmass_expression = model.jac_rhs_algebraic_eval(
+            t_casadi, y_casadi, p_casadi_stacked
+        ) - cj_casadi * pybamm.expression_tree.array.to_casadi_dm(mass_matrix)
         if self._options["jacobian"] == "dense":
             jac_times_cjmass_expression = casadi.densify(jac_times_cjmass_expression)
 
@@ -569,7 +581,9 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
                 np.zeros(model.len_alg),
             )
         else:
-            mass_action_expression = casadi.densify(casadi.DM(mass_matrix) @ v_casadi)
+            mass_action_expression = casadi.densify(
+                pybamm.expression_tree.array.to_casadi_dm(mass_matrix) @ v_casadi
+            )
 
         # also need the action of the mass matrix on a vector
         mass_action = casadi.Function(
@@ -672,9 +686,6 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
             compiled = aot_compile(list(fns.values()))
             fns = dict(zip(fns.keys(), compiled, strict=True))
 
-        def to_idaklu(fn):
-            return idaklu.generate_function(fn.serialize())
-
         self._setup = {
             "number_of_states": len(y0),
             "inputs": len(stacked_inputs),
@@ -702,20 +713,20 @@ def to_idaklu(fn):
         }
 
         for name in _SETUP_FCN_KEYS:
-            self._setup[name] = to_idaklu(fns[name])
+            self._setup[name] = to_idaklu_function(fns[name])
 
         # The idaklu-bound Function isn't callable from Python; keep the
         # original for the closest_event_idx lookup in _post_process_solution.
         self._setup["rootfn_casadi"] = fns["rootfn"]
 
         for key in self.output_variables:
-            self._setup["var_idaklu_fcns"].append(to_idaklu(fns[f"var:{key}"]))
+            self._setup["var_idaklu_fcns"].append(to_idaklu_function(fns[f"var:{key}"]))
             if has_sens:
                 self._setup["dvar_dy_idaklu_fcns"].append(
-                    to_idaklu(fns[f"dvar_dy:{key}"])
+                    to_idaklu_function(fns[f"dvar_dy:{key}"])
                 )
                 self._setup["dvar_dp_idaklu_fcns"].append(
-                    to_idaklu(fns[f"dvar_dp:{key}"])
+                    to_idaklu_function(fns[f"dvar_dp:{key}"])
                 )
 
         self._setup["solver"] = idaklu.create_casadi_solver_group(
diff --git a/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py b/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py
index a38535d521..0bf64cfb17 100644
--- a/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py
@@ -8,7 +8,7 @@
 import pybamm
 from pybamm.codegen.compilation import aot_compile
 from pybamm.solvers.base_solver import flatten_inputs
-from pybamm.solvers.idaklu_solver import _check_klu_pivot_tolerance
+from pybamm.solvers.idaklu_solver import _check_klu_pivot_tolerance, to_idaklu_function
 
 _DEFAULT_OPTIONS = {
     "compile": False,
@@ -186,8 +186,8 @@ def _set_up_root_solver(self, model, inputs_dict, t_eval):
 
     def _build_newton_solver(self, res_fn, jac_fn, len_alg):
         _setup = idaklu.StandaloneNewtonSolver(
-            residual=idaklu.generate_function(res_fn.serialize()),
-            jacobian=idaklu.generate_function(jac_fn.serialize()),
+            residual=to_idaklu_function(res_fn),
+            jacobian=to_idaklu_function(jac_fn),
             atol=np.full(len_alg, float(self.atol)).tolist(),
             rtol=float(self.rtol),
             step_tol=float(self.step_tol),
diff --git a/packages/pybamm/src/pybamm/solvers/scipy_solver.py b/packages/pybamm/src/pybamm/solvers/scipy_solver.py
index e18914dfce..f61ffec848 100644
--- a/packages/pybamm/src/pybamm/solvers/scipy_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/scipy_solver.py
@@ -38,6 +38,8 @@ class ScipySolver(pybamm.BaseSolver):
         for details.
     """
 
+    uses_rhs_jacobian = True
+
     def __init__(
         self,
         method="BDF",
diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_array.py b/packages/pybamm/tests/unit/test_expression_tree/test_array.py
index 59c03d0ace..7f680e1c79 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_array.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_array.py
@@ -3,8 +3,11 @@
 #
 
 
+import casadi
 import numpy as np
+import pytest
 import sympy
+from scipy.sparse import csc_matrix
 
 import pybamm
 
@@ -58,3 +61,30 @@ def test_to_from_json(self, mocker):
 
         # json to array conversion
         assert pybamm.Array._from_json(created_json) == arr
+
+    def test_to_casadi_dm_sparse(self):
+        # Duplicate, unsorted entries must be summed as scipy does
+        matrix = csc_matrix(
+            ([1.0, 2.0, 3.0, 4.0], [2, 2, 0, 1], [0, 0, 3, 3, 4]), shape=(3, 4)
+        )
+        matrix.has_canonical_format = False
+        converted = pybamm.expression_tree.array.to_casadi_dm(matrix)
+        np.testing.assert_array_equal(converted.full(), matrix.toarray())
+        assert converted.nnz() == 3
+        # The input is left untouched
+        assert not matrix.has_canonical_format
+
+    @pytest.mark.parametrize(
+        "entries",
+        [
+            np.array([1.0, 0.0, -2.5]),
+            np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]),
+            np.array(7.0),
+        ],
+    )
+    def test_to_casadi_dm_dense(self, entries):
+        converted = pybamm.expression_tree.array.to_casadi_dm(entries)
+        expected = casadi.DM(entries)
+        assert converted.shape == expected.shape
+        assert converted.is_dense()
+        np.testing.assert_array_equal(converted.full(), expected.full())
diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py b/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py
index b5ebc753ee..57dca2fb49 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py
@@ -51,53 +51,78 @@ def test_name(self):
         )
         assert sv.name == "y[0:10,20:30,...,60:70]"
 
-    def test_pass_evaluation_array(self):
-        # Turn off debug mode for this test
-        original_debug_mode = pybamm.settings.debug_mode
-        pybamm.settings.debug_mode = False
-        # Test that evaluation array gets passed down (doesn't have to be the correct
-        # array for this test)
-        array = np.array([1, 2, 3, 4, 5])
-        sv = pybamm.StateVector(slice(0, 10), evaluation_array=array)
-        np.testing.assert_array_equal(sv.evaluation_array, array)
-        # Turn debug mode back to what is was before
-        pybamm.settings.debug_mode = original_debug_mode
+    def test_evaluation_array_argument_deprecated(self):
+        with pytest.warns(DeprecationWarning, match=r"evaluation_array"):
+            sv = pybamm.StateVector(slice(0, 3), evaluation_array=[True, False, True])
+        # The slices define the entries read
+        assert sv.evaluation_array == [True, True, True]
 
     def test_failure(self):
         with pytest.raises(TypeError, match=r"all y_slices must be slice objects"):
             pybamm.StateVector(slice(0, 10), 1)
 
-    def test_to_from_json(self, mocker):
-        original_debug_mode = pybamm.settings.debug_mode
-        pybamm.settings.debug_mode = False
-
-        array = np.array([1, 2, 3, 4, 5])
-        sv = pybamm.StateVector(slice(0, 10), evaluation_array=array)
-
+    def test_to_from_json(self):
+        sv = pybamm.StateVector(slice(0, 3))
         json_dict = {
-            "name": "y[0:10]",
+            "name": "y[0:3]",
             "domains": {
                 "primary": [],
                 "secondary": [],
                 "tertiary": [],
                 "quaternary": [],
             },
-            "y_slice": [
-                {
-                    "start": 0,
-                    "stop": 10,
-                    "step": None,
-                }
-            ],
-            "evaluation_array": [1, 2, 3, 4, 5],
+            "y_slice": [{"start": 0, "stop": 3, "step": None}],
+            "evaluation_array": [True, True, True],
         }
-
         assert sv.to_json() == json_dict
-
         assert pybamm.StateVector._from_json(json_dict) == sv
 
-        # Turn debug mode back to what is was before
-        pybamm.settings.debug_mode = original_debug_mode
+    def test_several_slices(self):
+        sv = pybamm.StateVector(slice(0, 2), slice(4, 6))
+        assert sv.evaluation_array == [True, True, False, False, True, True]
+        assert sv.size == 4
+        np.testing.assert_array_equal(
+            sv.evaluate(y=np.arange(8.0)), [[0.0], [1.0], [4.0], [5.0]]
+        )
+
+    def test_identity_ignores_how_slices_are_split(self):
+        assert pybamm.StateVector(slice(53)) == pybamm.StateVector(slice(0, 53))
+        assert pybamm.StateVector(
+            slice(0, 2), slice(2, 5), name="y"
+        ) == pybamm.StateVector(slice(0, 5), name="y")
+        assert pybamm.StateVector(slice(0, 2)) != pybamm.StateVector(slice(0, 3))
+
+    @pytest.mark.parametrize(
+        "y_slices",
+        [
+            (slice(5, 8), slice(0, 3)),
+            (slice(0, 3), slice(2, 4)),
+            (slice(0, 6, 2),),
+        ],
+    )
+    def test_invalid_slices(self, y_slices):
+        with pytest.raises(pybamm.ShapeError, match=r"y_slices"):
+            pybamm.StateVector(*y_slices)
+
+    def test_unpickle_evaluation_array_state(self):
+        # Pickles from before the index ranges were stored hold the full mask
+        sv = pybamm.StateVector(slice(0, 2), slice(4, 6))
+        state = sv.__getstate__()
+        del state["_index_ranges"]
+        state["_evaluation_array"] = sv.evaluation_array
+        restored = pybamm.StateVector.__new__(pybamm.StateVector)
+        restored.__setstate__(state)
+        assert restored == sv
+        np.testing.assert_array_equal(
+            restored.evaluate(y=np.arange(6.0)), [[0.0], [1.0], [4.0], [5.0]]
+        )
+
+    def test_ranges_tile(self):
+        ranges_tile = pybamm.expression_tree.state_vector.ranges_tile
+        assert ranges_tile([(0, 2), (2, 5)], 0, 5)
+        assert not ranges_tile([(0, 2), (3, 5)], 0, 5)
+        assert not ranges_tile([(0, 3), (2, 5)], 0, 5)
+        assert not ranges_tile([(0, 2)], 0, 5)
 
 
 class TestStateVectorDot:
diff --git a/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py b/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py
index 17d892597f..520bf6e6c4 100644
--- a/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py
+++ b/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py
@@ -7,6 +7,7 @@
     UnstructuredMeshGenerator,
     UnstructuredSubMesh,
     UserSuppliedUnstructuredMesh,
+    _group_rows,
     _hex_grid,
     _hex_to_tet,
     _quad_to_tri,
@@ -503,6 +504,16 @@ def test_generator_parse_lims_with_string_var(self):
         assert len(spatial_vars) == 2
         assert len(spatial_lims) == 2
 
+    def test_group_rows_matches_unique(self):
+        rng = np.random.default_rng(0)
+        rows = np.sort(rng.integers(0, 6, size=(200, 3)), axis=1)
+        order, counts = _group_rows(rows)
+        _, inverse, unique_counts = np.unique(
+            rows, axis=0, return_inverse=True, return_counts=True
+        )
+        np.testing.assert_array_equal(order, np.lexsort(rows.T[::-1]))
+        np.testing.assert_array_equal(counts, unique_counts[inverse[order]])
+
 
 # ======================================================================
 # TestUnstructuredMeshGenerator
diff --git a/packages/pybamm/tests/unit/test_solvers/test_base_solver.py b/packages/pybamm/tests/unit/test_solvers/test_base_solver.py
index 56b9f73a93..bc154caea8 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_base_solver.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_base_solver.py
@@ -674,6 +674,32 @@ def test_check_restart_sensitivities_allows_full_state(self):
     def test_check_restart_sensitivities_allows_empty_solution(self):
         pybamm.BaseSolver._check_restart_sensitivities(pybamm.EmptySolution(), ["p"])
 
+    def test_stack_inputs(self):
+        inputs = {
+            "scalar": 2.0,
+            "vector": np.array([1.0, 3.0]),
+            "column": np.array([[4.0], [5.0]]),
+            "casadi": casadi.DM([6.0, 7.0]),
+        }
+        stacked = pybamm.solvers.base_solver.stack_inputs(inputs)
+        expected = casadi.vertcat(*inputs.values())
+        assert stacked.shape == expected.shape
+        np.testing.assert_array_equal(stacked.full(), expected.full())
+        assert pybamm.solvers.base_solver.stack_inputs({}).shape == (0, 1)
+
+    def test_rhs_jacobian_only_for_solvers_that_use_it(self):
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model)
+        sim.build()
+        idaklu_model = sim.built_model.new_copy()
+        pybamm.IDAKLUSolver().set_up(idaklu_model)
+        assert idaklu_model.jac_rhs_eval is None
+        with pytest.warns(DeprecationWarning):
+            scipy_solver = pybamm.ScipySolver()
+        scipy_model = sim.built_model.new_copy()
+        scipy_solver.set_up(scipy_model)
+        assert scipy_model.jac_rhs_eval is not None
+
 
 def _independent_decays():
     """u' = -a u and v' = -b v, with u(0) = v(0) = 1."""
diff --git a/packages/pybammsolvers/CMakeLists.txt b/packages/pybammsolvers/CMakeLists.txt
index af58ccc45b..5f59d15922 100644
--- a/packages/pybammsolvers/CMakeLists.txt
+++ b/packages/pybammsolvers/CMakeLists.txt
@@ -235,6 +235,12 @@ else ()
     endif()
 endif ()
 
+if (USE_PYTHON_CASADI)
+    # The extension links the Python casadi package's library, so it can share
+    # casadi::Function objects with Python instead of serialising them
+    target_compile_definitions(idaklu PRIVATE PYBAMM_SHARED_CASADI=1)
+endif()
+
 # Match idaklu's libstdc++ std::string ABI to the linked libcasadi (GNU only;
 # no-op for MSVC and macOS/libc++). CasADi's PyPI casadi-config.cmake forces
 # `_GLIBCXX_USE_CXX11_ABI=0` on consumers: correct for the x86_64 wheel (legacy
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
index 6f9d03d054..d9ece0875d 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
@@ -27,6 +27,40 @@ casadi::Function generate_casadi_function(const std::string &data)
 
 namespace py = pybind11;
 
+#ifdef PYBAMM_SHARED_CASADI
+// address is the casadi::Function behind a Python casadi.Function (int(f.this)),
+// from the same casadi library this extension links; the copy shares its node
+casadi::Function shared_casadi_function(std::uintptr_t address)
+{
+  return *reinterpret_cast<const casadi::Function *>(address);
+}
+
+// Fill the casadi::DM behind a Python casadi.DM (int(dm.this)) with a CSC
+// matrix, copying the buffers once instead of element by element through SWIG
+void fill_casadi_dm_sparse(
+  std::uintptr_t address, casadi_int nrow, casadi_int ncol,
+  const py::array_t<casadi_int, py::array::c_style | py::array::forcecast> &indptr,
+  const py::array_t<casadi_int, py::array::c_style | py::array::forcecast> &indices,
+  const py::array_t<double, py::array::c_style | py::array::forcecast> &data)
+{
+  std::vector<casadi_int> colind(indptr.data(), indptr.data() + indptr.size());
+  std::vector<casadi_int> row(indices.data(), indices.data() + indices.size());
+  std::vector<double> nonzeros(data.data(), data.data() + data.size());
+  *reinterpret_cast<casadi::DM *>(address) =
+    casadi::DM(casadi::Sparsity(nrow, ncol, colind, row), nonzeros);
+}
+
+// As fill_casadi_dm_sparse, for a dense matrix given in column-major order
+void fill_casadi_dm_dense(
+  std::uintptr_t address, casadi_int nrow, casadi_int ncol,
+  const py::array_t<double, py::array::c_style | py::array::forcecast> &data)
+{
+  std::vector<double> nonzeros(data.data(), data.data() + data.size());
+  *reinterpret_cast<casadi::DM *>(address) =
+    casadi::DM(casadi::Sparsity::dense(nrow, ncol), nonzeros);
+}
+#endif
+
 IDAKLUSolverGroup *create_casadi_solver_group(
   int number_of_states,
   int number_of_parameters,
@@ -202,6 +236,23 @@ PYBIND11_MODULE(idaklu, m)
     py::arg("string"),
     py::return_value_policy::take_ownership);
 
+#ifdef PYBAMM_SHARED_CASADI
+  m.attr("shares_python_casadi") = true;
+  m.def("shared_function", &shared_casadi_function,
+    "Share a Python casadi.Function, given int(f.this), without serialising it",
+    py::arg("address"),
+    py::return_value_policy::take_ownership);
+  m.def("fill_casadi_dm_sparse", &fill_casadi_dm_sparse,
+    "Fill a Python casadi.DM, given int(dm.this), with a CSC matrix",
+    py::arg("address"), py::arg("nrow"), py::arg("ncol"),
+    py::arg("indptr"), py::arg("indices"), py::arg("data"));
+  m.def("fill_casadi_dm_dense", &fill_casadi_dm_dense,
+    "Fill a Python casadi.DM, given int(dm.this), with a column-major dense matrix",
+    py::arg("address"), py::arg("nrow"), py::arg("ncol"), py::arg("data"));
+#else
+  m.attr("shares_python_casadi") = false;
+#endif
+
   m.def("sundials_error_message", &sundials_error_message,
     "Get a human-readable message for a SUNDIALS error code",
     py::arg("flag"),

From 4a1d7dca30fd7188eb9846312d62f125ecc5c33c Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Tue, 27 Oct 2026 16:50:48 -0400
Subject: [PATCH 2/8] to_casadi_dm: read CSR matrices directly; shorten
 changelog entry

---
 CHANGELOG.md                                  |  2 +-
 .../src/pybamm/expression_tree/array.py       | 20 +++++----
 .../unit/test_expression_tree/test_array.py   | 16 +++++--
 .../src/pybammsolvers/idaklu.cpp              | 43 ++++++++++++++++---
 4 files changed, 60 insertions(+), 21 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index d63e0998e3..c1325d9514 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -37,7 +37,7 @@
 
 - Read the Docs pull request builds drop `linkcheck` and the unpublished PDF, install the released solver wheel instead of compiling it, and are skipped when nothing the docs render changed, saving about 14 of 20 minutes. ([#5986](https://github.com/pybamm-team/PyBaMM/pull/5986))
 - Pull request CI no longer installs TeXLive, or the compilers and OpenBLAS that the prebuilt solver wheels make redundant. The latexify test mocks `sympy.preview` and runs everywhere. TeXLive downloads took 12–50 minutes in 6 of 15 runs. ([#5985](https://github.com/pybamm-team/PyBaMM/pull/5985))
-- Faster model build and `IDAKLUSolver` set-up for large (million-state) models. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
+- Faster model build and `IDAKLUSolver` set-up for large models. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
 - Experiments with a different model per step no longer compile Jacobians of the state mappers between models, which were unused. ([#5918](https://github.com/pybamm-team/PyBaMM/pull/5918))
 
 # [v26.10.1.1](https://github.com/pybamm-team/PyBaMM/tree/pybamm-v26.10.1.1) - 2026-10-19
diff --git a/packages/pybamm/src/pybamm/expression_tree/array.py b/packages/pybamm/src/pybamm/expression_tree/array.py
index c7c548245d..d6f3147a63 100644
--- a/packages/pybamm/src/pybamm/expression_tree/array.py
+++ b/packages/pybamm/src/pybamm/expression_tree/array.py
@@ -241,17 +241,19 @@ def to_casadi_dm(entries) -> casadi.DM:
         return casadi.DM(entries)
     matrix = casadi.DM()
     if issparse(entries):
-        csc = csc_matrix(entries)
-        if not csc.has_canonical_format:
-            csc = csc.copy()
-            csc.sum_duplicates()
+        compressed_rows = entries.format == "csr"
+        compressed = entries if compressed_rows else csc_matrix(entries)
+        if not compressed.has_canonical_format:
+            compressed = compressed.copy()
+            compressed.sum_duplicates()
         idaklu.fill_casadi_dm_sparse(
             int(matrix.this),
-            csc.shape[0],
-            csc.shape[1],
-            csc.indptr.astype(np.int64),
-            csc.indices.astype(np.int64),
-            csc.data.astype(np.float64),
+            compressed.shape[0],
+            compressed.shape[1],
+            compressed.indptr.astype(np.int64),
+            compressed.indices.astype(np.int64),
+            compressed.data.astype(np.float64),
+            compressed_rows=compressed_rows,
         )
     else:
         values = np.asarray(entries, dtype=np.float64)
diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_array.py b/packages/pybamm/tests/unit/test_expression_tree/test_array.py
index 7f680e1c79..8f84c15c1c 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_array.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_array.py
@@ -7,7 +7,7 @@
 import numpy as np
 import pytest
 import sympy
-from scipy.sparse import csc_matrix
+from scipy.sparse import coo_matrix, csc_matrix, csr_matrix
 
 import pybamm
 
@@ -62,10 +62,13 @@ def test_to_from_json(self, mocker):
         # json to array conversion
         assert pybamm.Array._from_json(created_json) == arr
 
-    def test_to_casadi_dm_sparse(self):
+    @pytest.mark.parametrize(
+        "matrix_format, shape", [(csc_matrix, (3, 4)), (csr_matrix, (4, 3))]
+    )
+    def test_to_casadi_dm_sparse(self, matrix_format, shape):
         # Duplicate, unsorted entries must be summed as scipy does
-        matrix = csc_matrix(
-            ([1.0, 2.0, 3.0, 4.0], [2, 2, 0, 1], [0, 0, 3, 3, 4]), shape=(3, 4)
+        matrix = matrix_format(
+            ([1.0, 2.0, 3.0, 4.0], [2, 2, 0, 1], [0, 0, 3, 3, 4]), shape=shape
         )
         matrix.has_canonical_format = False
         converted = pybamm.expression_tree.array.to_casadi_dm(matrix)
@@ -74,6 +77,11 @@ def test_to_casadi_dm_sparse(self):
         # The input is left untouched
         assert not matrix.has_canonical_format
 
+    def test_to_casadi_dm_other_sparse_format(self):
+        matrix = coo_matrix(([1.0, 2.0], ([0, 2], [1, 0])), shape=(3, 2))
+        converted = pybamm.expression_tree.array.to_casadi_dm(matrix)
+        np.testing.assert_array_equal(converted.full(), matrix.toarray())
+
     @pytest.mark.parametrize(
         "entries",
         [
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
index d9ece0875d..3ad43d6a57 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
@@ -35,17 +35,45 @@ casadi::Function shared_casadi_function(std::uintptr_t address)
   return *reinterpret_cast<const casadi::Function *>(address);
 }
 
-// Fill the casadi::DM behind a Python casadi.DM (int(dm.this)) with a CSC
+// Fill the casadi::DM behind a Python casadi.DM (int(dm.this)) with a CSC or CSR
 // matrix, copying the buffers once instead of element by element through SWIG
 void fill_casadi_dm_sparse(
   std::uintptr_t address, casadi_int nrow, casadi_int ncol,
   const py::array_t<casadi_int, py::array::c_style | py::array::forcecast> &indptr,
   const py::array_t<casadi_int, py::array::c_style | py::array::forcecast> &indices,
-  const py::array_t<double, py::array::c_style | py::array::forcecast> &data)
+  const py::array_t<double, py::array::c_style | py::array::forcecast> &data,
+  bool compressed_rows)
 {
-  std::vector<casadi_int> colind(indptr.data(), indptr.data() + indptr.size());
-  std::vector<casadi_int> row(indices.data(), indices.data() + indices.size());
-  std::vector<double> nonzeros(data.data(), data.data() + data.size());
+  const casadi_int *pointers = indptr.data();
+  const casadi_int *index = indices.data();
+  const double *values = data.data();
+  const auto nnz = static_cast<std::size_t>(data.size());
+  std::vector<casadi_int> colind, row;
+  std::vector<double> nonzeros;
+  if (compressed_rows) {
+    // Transpose by counting sort; rows come out sorted within each column
+    colind.assign(ncol + 1, 0);
+    for (std::size_t k = 0; k < nnz; k++) {
+      colind[index[k] + 1]++;
+    }
+    for (casadi_int j = 0; j < ncol; j++) {
+      colind[j + 1] += colind[j];
+    }
+    std::vector<casadi_int> next(colind.begin(), colind.end() - 1);
+    row.resize(nnz);
+    nonzeros.resize(nnz);
+    for (casadi_int i = 0; i < nrow; i++) {
+      for (casadi_int k = pointers[i]; k < pointers[i + 1]; k++) {
+        const casadi_int position = next[index[k]]++;
+        row[position] = i;
+        nonzeros[position] = values[k];
+      }
+    }
+  } else {
+    colind.assign(pointers, pointers + indptr.size());
+    row.assign(index, index + nnz);
+    nonzeros.assign(values, values + nnz);
+  }
   *reinterpret_cast<casadi::DM *>(address) =
     casadi::DM(casadi::Sparsity(nrow, ncol, colind, row), nonzeros);
 }
@@ -243,9 +271,10 @@ PYBIND11_MODULE(idaklu, m)
     py::arg("address"),
     py::return_value_policy::take_ownership);
   m.def("fill_casadi_dm_sparse", &fill_casadi_dm_sparse,
-    "Fill a Python casadi.DM, given int(dm.this), with a CSC matrix",
+    "Fill a Python casadi.DM, given int(dm.this), with a CSC or CSR matrix",
     py::arg("address"), py::arg("nrow"), py::arg("ncol"),
-    py::arg("indptr"), py::arg("indices"), py::arg("data"));
+    py::arg("indptr"), py::arg("indices"), py::arg("data"),
+    py::arg("compressed_rows") = false);
   m.def("fill_casadi_dm_dense", &fill_casadi_dm_dense,
     "Fill a Python casadi.DM, given int(dm.this), with a column-major dense matrix",
     py::arg("address"), py::arg("nrow"), py::arg("ncol"), py::arg("data"));

From 4324f3d9b66a87ff228f81b764599dd6bc5d4aa2 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Wed, 28 Oct 2026 10:04:36 -0400
Subject: [PATCH 3/8] Changelog: PR number

---
 CHANGELOG.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c1325d9514..b8195b5bc7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,13 +2,13 @@
 
 ## Breaking changes
 
-- `StateVector` and `StateVectorDot` raise `ShapeError` for slices that have a step, overlap or are out of order. Out-of-order slices were previously read only up to the last slice's stop. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
+- `StateVector` and `StateVectorDot` raise `ShapeError` for stepped, overlapping or out-of-order slices. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
 
 ## Deprecated
 
 - `FiniteVolume2D` is deprecated in favour of `FiniteVolumeUnstructured` and raises a `DeprecationWarning` when created. Remaining gaps before removal are tracked in [#5938](https://github.com/pybamm-team/PyBaMM/issues/5938). ([#5939](https://github.com/pybamm-team/PyBaMM/pull/5939))
 - `pybamm.Magnitude` is deprecated and will be removed in 27.1; use `pybamm.Component(symbol, 0)` / `pybamm.Component(symbol, 1)` for the lr/tb components of a `VectorField`. ([#5943](https://github.com/pybamm-team/PyBaMM/pull/5943))
-- The `evaluation_array` argument of `StateVector` and `StateVectorDot` is deprecated and ignored; the entries read are given by the slices. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
+- The `evaluation_array` argument of `StateVector` and `StateVectorDot` is deprecated and ignored. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
 
 ## Features
 
@@ -37,7 +37,7 @@
 
 - Read the Docs pull request builds drop `linkcheck` and the unpublished PDF, install the released solver wheel instead of compiling it, and are skipped when nothing the docs render changed, saving about 14 of 20 minutes. ([#5986](https://github.com/pybamm-team/PyBaMM/pull/5986))
 - Pull request CI no longer installs TeXLive, or the compilers and OpenBLAS that the prebuilt solver wheels make redundant. The latexify test mocks `sympy.preview` and runs everywhere. TeXLive downloads took 12–50 minutes in 6 of 15 runs. ([#5985](https://github.com/pybamm-team/PyBaMM/pull/5985))
-- Faster model build and `IDAKLUSolver` set-up for large models. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX))
+- Faster model build and `IDAKLUSolver` set-up for large models. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
 - Experiments with a different model per step no longer compile Jacobians of the state mappers between models, which were unused. ([#5918](https://github.com/pybamm-team/PyBaMM/pull/5918))
 
 # [v26.10.1.1](https://github.com/pybamm-team/PyBaMM/tree/pybamm-v26.10.1.1) - 2026-10-19

From aedbfed5512b8b0cc037095dc3f2c8acb6795884 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Wed, 28 Oct 2026 10:09:34 -0400
Subject: [PATCH 4/8] Changelog: migration notes

---
 CHANGELOG.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b8195b5bc7..9f9372e0e4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,13 +2,13 @@
 
 ## Breaking changes
 
-- `StateVector` and `StateVectorDot` raise `ShapeError` for stepped, overlapping or out-of-order slices. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
+- `StateVector` and `StateVectorDot` raise `ShapeError` for stepped, overlapping or out-of-order slices; pass sorted, disjoint slices. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
 
 ## Deprecated
 
 - `FiniteVolume2D` is deprecated in favour of `FiniteVolumeUnstructured` and raises a `DeprecationWarning` when created. Remaining gaps before removal are tracked in [#5938](https://github.com/pybamm-team/PyBaMM/issues/5938). ([#5939](https://github.com/pybamm-team/PyBaMM/pull/5939))
 - `pybamm.Magnitude` is deprecated and will be removed in 27.1; use `pybamm.Component(symbol, 0)` / `pybamm.Component(symbol, 1)` for the lr/tb components of a `VectorField`. ([#5943](https://github.com/pybamm-team/PyBaMM/pull/5943))
-- The `evaluation_array` argument of `StateVector` and `StateVectorDot` is deprecated and ignored. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
+- The `evaluation_array` argument of `StateVector` and `StateVectorDot` is deprecated and ignored; drop the argument. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
 
 ## Features
 

From 317dc5c8afefc9499868269660443ff5116df33d Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Wed, 28 Oct 2026 10:18:02 -0400
Subject: [PATCH 5/8] Cast Python casadi addresses in one place

---
 .../pybammsolvers/src/pybammsolvers/idaklu.cpp     | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
index 3ad43d6a57..6dbb72ed48 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
@@ -28,11 +28,19 @@ casadi::Function generate_casadi_function(const std::string &data)
 namespace py = pybind11;
 
 #ifdef PYBAMM_SHARED_CASADI
+// The C++ object behind a Python casadi object, given int(obj.this); Python only
+// exposes it as an integer, so the cast cannot be avoided
+template <class T>
+T *python_casadi_object(std::uintptr_t address)
+{
+  return reinterpret_cast<T *>(address);  // NOLINT(performance-no-int-to-ptr)
+}
+
 // address is the casadi::Function behind a Python casadi.Function (int(f.this)),
 // from the same casadi library this extension links; the copy shares its node
 casadi::Function shared_casadi_function(std::uintptr_t address)
 {
-  return *reinterpret_cast<const casadi::Function *>(address);
+  return *python_casadi_object<const casadi::Function>(address);
 }
 
 // Fill the casadi::DM behind a Python casadi.DM (int(dm.this)) with a CSC or CSR
@@ -74,7 +82,7 @@ void fill_casadi_dm_sparse(
     row.assign(index, index + nnz);
     nonzeros.assign(values, values + nnz);
   }
-  *reinterpret_cast<casadi::DM *>(address) =
+  *python_casadi_object<casadi::DM>(address) =
     casadi::DM(casadi::Sparsity(nrow, ncol, colind, row), nonzeros);
 }
 
@@ -84,7 +92,7 @@ void fill_casadi_dm_dense(
   const py::array_t<double, py::array::c_style | py::array::forcecast> &data)
 {
   std::vector<double> nonzeros(data.data(), data.data() + data.size());
-  *reinterpret_cast<casadi::DM *>(address) =
+  *python_casadi_object<casadi::DM>(address) =
     casadi::DM(casadi::Sparsity::dense(nrow, ncol), nonzeros);
 }
 #endif

From a3311d9711aced98b4c19d6b98f4750623115192 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Wed, 28 Oct 2026 11:50:04 -0400
Subject: [PATCH 6/8] to_casadi_dm: canonicalise sparse input for casadi's own
 conversion too; type-check the SWIG address

---
 .../pybamm/src/pybamm/expression_tree/array.py | 18 ++++++++++++++----
 .../pybamm/src/pybamm/solvers/idaklu_solver.py |  2 ++
 .../unit/test_expression_tree/test_array.py    | 15 ++++++++++++---
 3 files changed, 28 insertions(+), 7 deletions(-)

diff --git a/packages/pybamm/src/pybamm/expression_tree/array.py b/packages/pybamm/src/pybamm/expression_tree/array.py
index d6f3147a63..a1cf46fd08 100644
--- a/packages/pybamm/src/pybamm/expression_tree/array.py
+++ b/packages/pybamm/src/pybamm/expression_tree/array.py
@@ -237,16 +237,19 @@ def to_casadi_dm(entries) -> casadi.DM:
     casadi.DM
         The values, sparse if ``entries`` is sparse.
     """
-    if not idaklu.shares_python_casadi:
-        return casadi.DM(entries)
-    matrix = casadi.DM()
     if issparse(entries):
         compressed_rows = entries.format == "csr"
         compressed = entries if compressed_rows else csc_matrix(entries)
+        # casadi aborts the process on unsorted or duplicate entries
         if not compressed.has_canonical_format:
             compressed = compressed.copy()
             compressed.sum_duplicates()
+        if not idaklu.shares_python_casadi:
+            return casadi.DM(csc_matrix(compressed))
+        matrix = casadi.DM()
         idaklu.fill_casadi_dm_sparse(
+            # SWIG's address attribute is not in casadi's stubs
+            # pyrefly: ignore[missing-attribute]
             int(matrix.this),
             compressed.shape[0],
             compressed.shape[1],
@@ -256,10 +259,17 @@ def to_casadi_dm(entries) -> casadi.DM:
             compressed_rows=compressed_rows,
         )
     else:
+        if not idaklu.shares_python_casadi:
+            return casadi.DM(entries)
+        matrix = casadi.DM()
         values = np.asarray(entries, dtype=np.float64)
         if values.ndim < 2:
             values = values.reshape(-1, 1)
         idaklu.fill_casadi_dm_dense(
-            int(matrix.this), values.shape[0], values.shape[1], values.ravel(order="F")
+            # pyrefly: ignore[missing-attribute]
+            int(matrix.this),
+            values.shape[0],
+            values.shape[1],
+            values.ravel(order="F"),
         )
     return matrix
diff --git a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
index b4fae0bd12..43a4927335 100644
--- a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
@@ -83,6 +83,8 @@ def to_idaklu_function(fn: casadi.Function):
     shares the function directly; otherwise it is serialised and rebuilt.
     """
     if idaklu.shares_python_casadi:
+        # SWIG's address attribute is not in casadi's stubs
+        # pyrefly: ignore[missing-attribute]
         return idaklu.shared_function(int(fn.this))
     return idaklu.generate_function(fn.serialize())
 
diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_array.py b/packages/pybamm/tests/unit/test_expression_tree/test_array.py
index 8f84c15c1c..97f9281f5d 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_array.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_array.py
@@ -12,6 +12,15 @@
 import pybamm
 
 
+@pytest.fixture(params=[True, False], ids=["shared", "copied"])
+def casadi_sharing(request, monkeypatch):
+    """Run with and without the shared casadi library, where it is available."""
+    idaklu = pybamm.expression_tree.array.idaklu
+    if request.param and not idaklu.shares_python_casadi:
+        pytest.skip("pybammsolvers does not share the Python casadi library")
+    monkeypatch.setattr(idaklu, "shares_python_casadi", request.param)
+
+
 class TestArray:
     def test_name(self):
         arr = pybamm.Array(np.array([1, 2, 3]))
@@ -65,7 +74,7 @@ def test_to_from_json(self, mocker):
     @pytest.mark.parametrize(
         "matrix_format, shape", [(csc_matrix, (3, 4)), (csr_matrix, (4, 3))]
     )
-    def test_to_casadi_dm_sparse(self, matrix_format, shape):
+    def test_to_casadi_dm_sparse(self, matrix_format, shape, casadi_sharing):
         # Duplicate, unsorted entries must be summed as scipy does
         matrix = matrix_format(
             ([1.0, 2.0, 3.0, 4.0], [2, 2, 0, 1], [0, 0, 3, 3, 4]), shape=shape
@@ -77,7 +86,7 @@ def test_to_casadi_dm_sparse(self, matrix_format, shape):
         # The input is left untouched
         assert not matrix.has_canonical_format
 
-    def test_to_casadi_dm_other_sparse_format(self):
+    def test_to_casadi_dm_other_sparse_format(self, casadi_sharing):
         matrix = coo_matrix(([1.0, 2.0], ([0, 2], [1, 0])), shape=(3, 2))
         converted = pybamm.expression_tree.array.to_casadi_dm(matrix)
         np.testing.assert_array_equal(converted.full(), matrix.toarray())
@@ -90,7 +99,7 @@ def test_to_casadi_dm_other_sparse_format(self):
             np.array(7.0),
         ],
     )
-    def test_to_casadi_dm_dense(self, entries):
+    def test_to_casadi_dm_dense(self, entries, casadi_sharing):
         converted = pybamm.expression_tree.array.to_casadi_dm(entries)
         expected = casadi.DM(entries)
         assert converted.shape == expected.shape

From 02cfb2fecd6ca3d79af081439ac1764480687df0 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Thu, 29 Oct 2026 14:00:58 -0400
Subject: [PATCH 7/8] Address review: empty/negative StateVector slices,
 rhs-jacobian opt-out default, naming and imports

---
 CHANGELOG.md                                  |  2 +-
 .../pybamm/expression_tree/state_vector.py    | 12 ++++++++----
 .../src/pybamm/meshes/unstructured_submesh.py |  3 +--
 .../pybamm/src/pybamm/solvers/base_solver.py  |  7 ++++---
 .../src/pybamm/solvers/casadi_solver.py       |  2 ++
 .../pybamm/src/pybamm/solvers/dummy_solver.py |  2 ++
 .../src/pybamm/solvers/idaklu_solver.py       | 19 ++++++++++++-------
 .../pybamm/src/pybamm/solvers/jax_solver.py   |  2 ++
 .../src/pybamm/solvers/nonlinear_solver.py    | 11 ++++++++---
 .../pybamm/src/pybamm/solvers/scipy_solver.py |  2 --
 .../test_expression_tree/test_state_vector.py |  7 +++++++
 11 files changed, 47 insertions(+), 22 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9f9372e0e4..097a75db84 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@
 
 ## Breaking changes
 
-- `StateVector` and `StateVectorDot` raise `ShapeError` for stepped, overlapping or out-of-order slices; pass sorted, disjoint slices. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
+- `StateVector` and `StateVectorDot` raise `ShapeError` for stepped, negative, overlapping or out-of-order slices; pass sorted, disjoint, non-negative slices. ([#5994](https://github.com/pybamm-team/PyBaMM/pull/5994))
 
 ## Deprecated
 
diff --git a/packages/pybamm/src/pybamm/expression_tree/state_vector.py b/packages/pybamm/src/pybamm/expression_tree/state_vector.py
index bc5269ef93..d082d9a360 100644
--- a/packages/pybamm/src/pybamm/expression_tree/state_vector.py
+++ b/packages/pybamm/src/pybamm/expression_tree/state_vector.py
@@ -64,7 +64,8 @@ class StateVectorBase(pybamm.Symbol):
     Raises
     ------
     :class:`pybamm.ShapeError`
-        If a slice has a step, or the slices are not in increasing order or overlap.
+        If a slice has a step or a negative bound, or the slices are not in
+        increasing order or overlap.
     """
 
     __slots__ = ("_first_point", "_index_ranges", "_last_point", "_y_slices")
@@ -147,8 +148,8 @@ def size(self):
 
     def _read(self, vector):
         """The entries of ``vector`` this state vector reads, as a column."""
-        if len(self._index_ranges) == 1:
-            start, stop = self._index_ranges[0]
+        if len(self._index_ranges) <= 1:
+            start, stop = self._index_ranges[0] if self._index_ranges else (0, 0)
             out = vector[start:stop].copy()
         else:
             out = np.concatenate(
@@ -167,13 +168,16 @@ def _merge_slices(y_slices) -> tuple[tuple[int, int], ...]:
         Raises
         ------
         :class:`pybamm.ShapeError`
-            If a slice has a step, or the slices are not increasing or overlap.
+            If a slice has a step or a negative bound, or the slices are not
+            increasing or overlap.
         """
         ranges: list[list[int]] = []
         for y_slice in y_slices:
             if y_slice.step not in (None, 1):
                 raise pybamm.ShapeError("y_slices must not have a step")
             start = y_slice.start or 0
+            if start < 0 or y_slice.stop < 0:
+                raise pybamm.ShapeError("y_slices must not have negative bounds")
             if ranges and start < ranges[-1][1]:
                 raise pybamm.ShapeError(
                     "y_slices must be in increasing order and must not overlap"
diff --git a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
index 5f0ffdf66b..b3b742637c 100644
--- a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
+++ b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
@@ -254,8 +254,7 @@ def _build_face_connectivity(self):
     # ------------------------------------------------------------------
 
     def _compute_face_geometry(self):
-        # Contiguous per-vertex gathers; strided views of an (n, k, d) gather
-        # make every later vector operation slower
+        # Contiguous (n_faces, d) array per local vertex
         face_verts = [
             self.vertices[self.faces[:, i]] for i in range(self.faces.shape[1])
         ]
diff --git a/packages/pybamm/src/pybamm/solvers/base_solver.py b/packages/pybamm/src/pybamm/solvers/base_solver.py
index ea704bf78c..4428f19b41 100644
--- a/packages/pybamm/src/pybamm/solvers/base_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/base_solver.py
@@ -14,6 +14,7 @@
 
 import pybamm
 from pybamm import ParameterValues
+from pybamm.expression_tree.array import to_casadi_dm
 from pybamm.expression_tree.binary_operators import _Heaviside
 from pybamm.expression_tree.input_parameter import DUMMY_INPUT_PARAMETER_VALUE
 from pybamm.models.base_model import ModelSolutionObservability
@@ -64,7 +65,7 @@ class BaseSolver:
     """
 
     #: Whether the solver uses the jacobian of the rhs on its own
-    uses_rhs_jacobian = False
+    uses_rhs_jacobian = True
 
     def __init__(
         self,
@@ -357,7 +358,7 @@ def set_up(
                 if model.is_standard_form_dae:
                     explicit_rhs = rhs_eval
                 else:
-                    M_ode = pybamm.expression_tree.array.to_casadi_dm(
+                    M_ode = to_casadi_dm(
                         model.mass_matrix.entries[: model.len_rhs, : model.len_rhs]
                     )
                     explicit_rhs = casadi.solve(M_ode, rhs_eval)
@@ -2268,7 +2269,7 @@ def stack_inputs(inputs_dict: dict) -> casadi.DM:
     stacked = np.concatenate(
         [np.asarray(value, dtype=float).reshape(-1) for value in inputs_dict.values()]
     )
-    return pybamm.expression_tree.array.to_casadi_dm(stacked)
+    return to_casadi_dm(stacked)
 
 
 def process(
diff --git a/packages/pybamm/src/pybamm/solvers/casadi_solver.py b/packages/pybamm/src/pybamm/solvers/casadi_solver.py
index 40576d0b06..a19af37cf0 100644
--- a/packages/pybamm/src/pybamm/solvers/casadi_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/casadi_solver.py
@@ -81,6 +81,8 @@ class CasadiSolver(pybamm.BaseSolver):
         None leaves the number of integrators unbound. Default is 100.
     """
 
+    uses_rhs_jacobian = False
+
     def __init__(
         self,
         mode="safe",
diff --git a/packages/pybamm/src/pybamm/solvers/dummy_solver.py b/packages/pybamm/src/pybamm/solvers/dummy_solver.py
index ac8fc9aee9..907a54c3de 100644
--- a/packages/pybamm/src/pybamm/solvers/dummy_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/dummy_solver.py
@@ -9,6 +9,8 @@
 class DummySolver(pybamm.BaseSolver):
     """Dummy solver class for empty models."""
 
+    uses_rhs_jacobian = False
+
     def __init__(self):
         super().__init__()
         self.name = "Dummy solver"
diff --git a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
index 43a4927335..dc714abb76 100644
--- a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
@@ -13,6 +13,7 @@
 
 import pybamm
 from pybamm.codegen.compilation import aot_compile
+from pybamm.expression_tree.array import to_casadi_dm
 from pybamm.solvers.base_solver import flatten_inputs, stack_inputs
 from pybamm.solvers.observation import OutputAssembly
 
@@ -75,7 +76,7 @@ def _sensitivity_scales(inputs_dict: dict, sensitivity_names: list[str]) -> np.n
 _IDA_ROOT_RETURN = 2
 
 
-def to_idaklu_function(fn: casadi.Function):
+def casadi_to_idaklu_function(fn: casadi.Function):
     """
     Hand a casadi.Function to pybammsolvers.
 
@@ -307,6 +308,8 @@ class IDAKLUSolver(pybamm.BaseSolver):
 
     """
 
+    uses_rhs_jacobian = False
+
     class StateID(IntEnum):
         ALGEBRAIC = 0
         DIFFERENTIAL = 1
@@ -550,7 +553,7 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
 
         jac_times_cjmass_expression = model.jac_rhs_algebraic_eval(
             t_casadi, y_casadi, p_casadi_stacked
-        ) - cj_casadi * pybamm.expression_tree.array.to_casadi_dm(mass_matrix)
+        ) - cj_casadi * to_casadi_dm(mass_matrix)
         if self._options["jacobian"] == "dense":
             jac_times_cjmass_expression = casadi.densify(jac_times_cjmass_expression)
 
@@ -584,7 +587,7 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
             )
         else:
             mass_action_expression = casadi.densify(
-                pybamm.expression_tree.array.to_casadi_dm(mass_matrix) @ v_casadi
+                to_casadi_dm(mass_matrix) @ v_casadi
             )
 
         # also need the action of the mass matrix on a vector
@@ -715,20 +718,22 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
         }
 
         for name in _SETUP_FCN_KEYS:
-            self._setup[name] = to_idaklu_function(fns[name])
+            self._setup[name] = casadi_to_idaklu_function(fns[name])
 
         # The idaklu-bound Function isn't callable from Python; keep the
         # original for the closest_event_idx lookup in _post_process_solution.
         self._setup["rootfn_casadi"] = fns["rootfn"]
 
         for key in self.output_variables:
-            self._setup["var_idaklu_fcns"].append(to_idaklu_function(fns[f"var:{key}"]))
+            self._setup["var_idaklu_fcns"].append(
+                casadi_to_idaklu_function(fns[f"var:{key}"])
+            )
             if has_sens:
                 self._setup["dvar_dy_idaklu_fcns"].append(
-                    to_idaklu_function(fns[f"dvar_dy:{key}"])
+                    casadi_to_idaklu_function(fns[f"dvar_dy:{key}"])
                 )
                 self._setup["dvar_dp_idaklu_fcns"].append(
-                    to_idaklu_function(fns[f"dvar_dp:{key}"])
+                    casadi_to_idaklu_function(fns[f"dvar_dp:{key}"])
                 )
 
         self._setup["solver"] = idaklu.create_casadi_solver_group(
diff --git a/packages/pybamm/src/pybamm/solvers/jax_solver.py b/packages/pybamm/src/pybamm/solvers/jax_solver.py
index 9704a4aa02..8a96401656 100644
--- a/packages/pybamm/src/pybamm/solvers/jax_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/jax_solver.py
@@ -52,6 +52,8 @@ class JaxSolver(pybamm.BaseSolver):
         for details.
     """
 
+    uses_rhs_jacobian = False
+
     def __init__(
         self,
         method="BDF",
diff --git a/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py b/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py
index 0bf64cfb17..abc7bbd576 100644
--- a/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/nonlinear_solver.py
@@ -8,7 +8,10 @@
 import pybamm
 from pybamm.codegen.compilation import aot_compile
 from pybamm.solvers.base_solver import flatten_inputs
-from pybamm.solvers.idaklu_solver import _check_klu_pivot_tolerance, to_idaklu_function
+from pybamm.solvers.idaklu_solver import (
+    _check_klu_pivot_tolerance,
+    casadi_to_idaklu_function,
+)
 
 _DEFAULT_OPTIONS = {
     "compile": False,
@@ -71,6 +74,8 @@ class NonlinearSolver(pybamm.BaseSolver):
         partial pivoting tolerance.
     """
 
+    uses_rhs_jacobian = False
+
     def __init__(
         self,
         atol=1e-6,
@@ -186,8 +191,8 @@ def _set_up_root_solver(self, model, inputs_dict, t_eval):
 
     def _build_newton_solver(self, res_fn, jac_fn, len_alg):
         _setup = idaklu.StandaloneNewtonSolver(
-            residual=to_idaklu_function(res_fn),
-            jacobian=to_idaklu_function(jac_fn),
+            residual=casadi_to_idaklu_function(res_fn),
+            jacobian=casadi_to_idaklu_function(jac_fn),
             atol=np.full(len_alg, float(self.atol)).tolist(),
             rtol=float(self.rtol),
             step_tol=float(self.step_tol),
diff --git a/packages/pybamm/src/pybamm/solvers/scipy_solver.py b/packages/pybamm/src/pybamm/solvers/scipy_solver.py
index f61ffec848..e18914dfce 100644
--- a/packages/pybamm/src/pybamm/solvers/scipy_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/scipy_solver.py
@@ -38,8 +38,6 @@ class ScipySolver(pybamm.BaseSolver):
         for details.
     """
 
-    uses_rhs_jacobian = True
-
     def __init__(
         self,
         method="BDF",
diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py b/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py
index 57dca2fb49..87fd5d8045 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_state_vector.py
@@ -98,12 +98,19 @@ def test_identity_ignores_how_slices_are_split(self):
             (slice(5, 8), slice(0, 3)),
             (slice(0, 3), slice(2, 4)),
             (slice(0, 6, 2),),
+            (slice(-2, 10),),
+            (slice(0, -1),),
         ],
     )
     def test_invalid_slices(self, y_slices):
         with pytest.raises(pybamm.ShapeError, match=r"y_slices"):
             pybamm.StateVector(*y_slices)
 
+    def test_empty_slice(self):
+        sv = pybamm.StateVector(slice(0, 0))
+        assert sv.size == 0
+        assert sv.evaluate(y=np.arange(10.0)).shape == (0, 1)
+
     def test_unpickle_evaluation_array_state(self):
         # Pickles from before the index ranges were stored hold the full mask
         sv = pybamm.StateVector(slice(0, 2), slice(4, 6))

From 8e2701bae8d0d03ae2ab98adbeb73af2d9ff34b3 Mon Sep 17 00:00:00 2001
From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com>
Date: Thu, 29 Oct 2026 14:05:21 -0400
Subject: [PATCH 8/8] Use np.mean for unstructured cell and face centroids

---
 packages/pybamm/src/pybamm/meshes/unstructured_submesh.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
index b3b742637c..a9514f59ce 100644
--- a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
+++ b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
@@ -113,7 +113,7 @@ def _compute_cell_geometry(self):
         ]
         # Vertex mean is the exact centroid for simplices (triangles, tets);
         # quads and hexes overwrite it with the true area/volume centroid.
-        self.cell_centroids = sum(verts) / len(verts)
+        self.cell_centroids = np.mean(verts, axis=0)
 
         if self.element_type == ElementType.TRIANGLE:
             v0, v1, v2 = verts
@@ -258,7 +258,7 @@ def _compute_face_geometry(self):
         face_verts = [
             self.vertices[self.faces[:, i]] for i in range(self.faces.shape[1])
         ]
-        self.face_centroids = sum(face_verts) / len(face_verts)
+        self.face_centroids = np.mean(face_verts, axis=0)
 
         if self.dimension == 2:
             v0, v1 = face_verts