BradyPlanden · Sep 22, 2026
@@ -17,6 +17,7 @@
 - `pybamm-core` compiles a model into a `CompiledModel`: a DAE with a symbolic Jacobian, per-output sparsity, coloured and lane-batched Jacobian sweeps, fused event and output tapes, and a versioned C ABI (`ffi.rs`) that a drift test pins against pybammsolvers' header. ([#5957](https://github.com/pybamm-team/PyBaMM/pull/5957))
 - `pybamm-core`'s reverse-mode AD values the primitives through the same definitions as the interpreters, and a test holds its partials bitwise to them. ([#5957](https://github.com/pybamm-team/PyBaMM/pull/5957))
 - `pybamm.rust` exposes `ExprGraph`, `CompiledFunction` (one expression, or a named set via `compile_outputs`, with lane-batched trajectory and Hermite evaluation) and `CompiledJacobian`, with hand-written type stubs that `mypy.stubtest` pins against the built extension. ([#5958](https://github.com/pybamm-team/PyBaMM/pull/5958))
+- `Symbol.to_rust` lowers a discretised expression tree into a `pybamm.rust.ExprGraph`, covering the operators, matrix products, concatenations, `Conditional`, the specific functions, `RegPower` and 1D-3D interpolants. Spatial and tensor operators, such as `Norm` and `Component`, raise `TypeError`. ([#5959](https://github.com/pybamm-team/PyBaMM/pull/5959))
 
 ## Bug fixes
 

@@ -1,8 +1,8 @@
 # PyBaMM Rust Core
 
 Expression compiler and evaluator for PyBaMM models, filling the role CasADi plays on
-PyBaMM's other backends. Python hands over an expression DAG; this workspace turns it
-into flat instruction tapes it can evaluate and differentiate.
+PyBaMM's other backends. Python hands over a discretised model as an expression DAG;
+this workspace turns it into flat instruction tapes it can evaluate and differentiate.
 
 Python reaches it through `pybamm.rust`.
 

@@ -3,6 +3,8 @@
 #
 from __future__ import annotations
 
+from typing import cast
+
 import casadi
 import numpy as np
 import numpy.typing as npt
@@ -134,6 +136,45 @@ 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))
 
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        entries = self.evaluate()
+        # A Rust sparse matrix only feeds a matrix product, so a value is dense
+        if issparse(entries):
+            entries = cast(csr_matrix, entries).toarray()
+        entries = np.asarray(entries, dtype=float)
+        if entries.ndim == 2 and entries.shape[1] > 1:
+            return graph.dense_matrix(
+                np.ascontiguousarray(entries).ravel(),
+                entries.shape[0],
+                entries.shape[1],
+            )
+        return graph.array(entries.flatten())
+
+    def _to_rust_matmul_left(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        entries = self.evaluate()
+        if not issparse(entries):
+            return self._to_rust_inner(graph, rust_symbols)
+        # Keyed apart from the dense value of the same symbol
+        key = (self.id, "matmul left")
+        cached = rust_symbols.get(key)
+        if cached is None:
+            # Rows sum in stored order; in column order a stencil's rounding telescopes
+            csr = csr_matrix(entries, copy=True)
+            csr.sum_duplicates()
+            cached = graph.sparse_matrix(
+                csr.indptr.tolist(),
+                csr.indices.tolist(),
+                csr.data.astype(float),
+                csr.shape[0],
+                csr.shape[1],
+            )
+            rust_symbols[key] = cached
+        return cached
+
     def _jac(self, variable) -> pybamm.Matrix:
         """See :meth:`pybamm.Symbol._jac()`."""
         # Return zeros of correct size

@@ -237,6 +237,48 @@ def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
         )
         return self._casadi_evaluate(converted_left, converted_right)
 
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        if isinstance(self, MatrixMultiplication):
+            left = self.left._to_rust_matmul_left(graph, rust_symbols)
+            right = self.right._to_rust_inner(graph, rust_symbols)
+            return graph.matmul(left, right)
+        left, right = self._children_to_rust(graph, rust_symbols)
+        if isinstance(self, Addition):
+            return graph.add(left, right)
+        elif isinstance(self, Subtraction):
+            return graph.sub(left, right)
+        elif isinstance(self, Multiplication):
+            return graph.mul(left, right)
+        elif isinstance(self, Division):
+            return graph.div(left, right)
+        elif isinstance(self, Power):
+            return graph.pow(left, right)
+        elif isinstance(self, Inner):
+            return graph.mul(left, right)
+        elif isinstance(self, Minimum):
+            return graph.minimum(left, right)
+        elif isinstance(self, Maximum):
+            return graph.maximum(left, right)
+        elif isinstance(self, Modulo):
+            return graph.modulo(left, right)
+        elif isinstance(self, Hypot):
+            return graph.hypot(left, right)
+        elif isinstance(self, EqualHeaviside):
+            return graph.equal_heaviside(left, right)
+        elif isinstance(self, NotEqualHeaviside):
+            return graph.not_equal_heaviside(left, right)
+        elif isinstance(self, Equality):
+            return graph.equality(left, right)
+        else:
+            raise TypeError(
+                f"Cannot convert binary operator of type '{type(self).__name__}' "
+                "to the Rust backend. Use a CasADi-backed model (set "
+                "`model.convert_to_format = 'casadi'`) for models using "
+                "this operator."
+            )
+
     def _evaluates_on_edges(self, dimension: str) -> bool | str:
         """See :meth:`pybamm.Symbol._evaluates_on_edges()`."""
         return self.left.evaluates_on_edges(dimension) or self.right.evaluates_on_edges(

@@ -261,6 +261,12 @@ def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
         )
         return casadi.vertcat(*converted_children)
 
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        converted_children = self._children_to_rust(graph, rust_symbols)
+        return graph.concat(converted_children)
+
     def _concatenation_jac(self, children_jacs):
         """See :meth:`pybamm.Concatenation.concatenation_jac()`."""
         children = self._children
@@ -459,6 +465,29 @@ def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
             all_child_vectors.extend(v for _, v in child_vectors)
         return casadi.vertcat(*all_child_vectors)
 
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        converted_children = self._children_to_rust(graph, rust_symbols)
+        all_child_exprs = []
+        for i in range(self.secondary_dimensions_npts):
+            child_exprs = []
+            for child_var, slices in zip(
+                converted_children, self._children_slices, strict=True
+            ):
+                for child_dom, child_slice in slices.items():
+                    child_exprs.append(
+                        (
+                            self._slices[child_dom][i].start,
+                            graph.index(
+                                child_var, child_slice[i].start, child_slice[i].stop
+                            ),
+                        )
+                    )
+            child_exprs.sort(key=itemgetter(0))
+            all_child_exprs.extend(expr for _, expr in child_exprs)
+        return graph.concat(all_child_exprs)
+
     def _concatenation_jac(self, children_jacs):
         """See :meth:`pybamm.Concatenation.concatenation_jac()`."""
         # note that this assumes that the children are in the right order and only have
@@ -556,6 +585,12 @@ def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
         )
         return casadi.vertcat(*converted_children)
 
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        converted_children = self._children_to_rust(graph, rust_symbols)
+        return graph.concat(converted_children)
+
     def _concatenation_new_copy(self, children, perform_simplifications=True):
         """See :meth:`pybamm.Concatenation._concatenation_new_copy()`."""
         return SparseStack(*children)

@@ -195,6 +195,15 @@ def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
         index = casadi.if_else(on_boundary, -1, floored - 1)
         return switch(index, *shared)
 
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        converted_selector = self.selector._to_rust_inner(graph, rust_symbols)
+        converted_branches = [
+            branch._to_rust_inner(graph, rust_symbols) for branch in self.branches
+        ]
+        return graph.conditional(converted_selector, converted_branches)
+
     def to_equation(self):
         if self.print_name is not None:
             return sympy.Symbol(self.print_name)

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from typing import Any
 
 import numpy as np
@@ -125,3 +127,10 @@ def sum_times(self):
     def _unary_evaluate(self, child):
         # return result of evaluating the child, we'll only implement the sum once the model is solved (in pybamm.ProcessedVariable)
         return child
+
+    def _to_rust(
+        self, graph: pybamm.rust.ExprGraph, rust_symbols: dict
+    ) -> pybamm.rust.Expr:
+        # Pass-through: actual summation happens in ProcessedVariable
+        (child,) = self._children_to_rust(graph, rust_symbols)
+        return child