BradyPlanden · Sep 22, 2026
@@ -18,6 +18,7 @@
 - `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))
+- `model.convert_to_format = "rust"` evaluates a model's rhs, algebraic, Jacobian and event functions through compiled Rust evaluators, for `ScipySolver`, `AlgebraicSolver` and EIS. ([#5960](https://github.com/pybamm-team/PyBaMM/pull/5960))
 
 ## Bug fixes
 

@@ -4,7 +4,7 @@ Expression compiler and evaluator for PyBaMM models, filling the role CasADi pla
 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`.
+Models opt in with `model.convert_to_format = "rust"`.
 
 ## Layout
 

@@ -96,10 +96,10 @@ class BaseModel:
     """
 
     _DEFAULT_CONVERT_TO_FORMAT = "casadi"
-    _VALID_CONVERT_TO_FORMATS = ("python", "casadi", "jax")
+    _VALID_CONVERT_TO_FORMATS = ("python", "casadi", "jax", "rust")
 
-    # Set by discretisation, deserialisation and BaseSolver.set_up. Evaluators are
-    # CasADi functions or pybamm.EvaluatorPython/EvaluatorJax, by convert_to_format.
+    # Set by discretisation, deserialisation and BaseSolver.set_up. Evaluators are CasADi
+    # functions, pybamm.EvaluatorPython/EvaluatorJax or RustEvaluators, by convert_to_format.
     len_rhs: int
     len_alg: int
     schema_version: str
@@ -278,6 +278,7 @@ def convert_to_format(self) -> str | None:
             - "python": convert to Python code for evaluating `evaluate(t, y)` on expressions.
             - "casadi": convert to CasADi expression tree for Jacobian calculation.
             - "jax": convert to JAX expression tree.
+            - "rust": convert to a compiled Rust expression graph.
 
         Default is "casadi". ``None`` is deprecated: it evaluates as "python" and
         raises a ``DeprecationWarning`` when set.
@@ -304,8 +305,8 @@ def convert_to_format(self, value: str | None) -> None:
     @property
     def uses_stacked_inputs(self) -> bool:
         """Whether the converted evaluators take the inputs stacked into one
-        column vector rather than as a dict."""
-        return self._convert_to_format == "casadi"
+        vector (a CasADi column, or a flat array for "rust") rather than as a dict."""
+        return self._convert_to_format in ("casadi", "rust")
 
     @property
     def rhs(self):

@@ -124,13 +124,17 @@ def _build_matrix_problem(self, inputs_dict=None):
         model = self._built_model
         inputs_dict = inputs_dict or {}
 
-        # Convert inputs to casadi format for Jacobian evaluation
+        # Inputs in the format the compiled Jacobian expects for each backend
         if model.convert_to_format == "casadi":
             from casadi import vertcat
 
-            casadi_inputs = vertcat(*inputs_dict.values()) if inputs_dict else []
+            eval_inputs = vertcat(*inputs_dict.values()) if inputs_dict else []
+        elif model.convert_to_format == "rust":
+            from pybamm.solvers.base_solver import stack_inputs
+
+            eval_inputs = stack_inputs(inputs_dict, "rust")
         else:
-            casadi_inputs = inputs_dict
+            eval_inputs = inputs_dict
 
         # Only compile Jacobian/model functions on first call; the compiled
         # functions persist on the model and work with any inputs/y0 values.
@@ -150,8 +154,10 @@ def _build_matrix_problem(self, inputs_dict=None):
         self._output_jacobian = Function(
             "eis_outputs", [y], [jacobian(output_expression, y)]
         )(y0).sparse()
-        J_sparse = model.jac_rhs_algebraic_eval(0, y0, casadi_inputs).sparse()
-        neg_J = -csc_matrix(J_sparse)
+        J = model.jac_rhs_algebraic_eval(0, y0, eval_inputs)
+        if model.convert_to_format == "casadi":
+            J = J.sparse()
+        neg_J = -csc_matrix(J)
 
         # State ordering and matrices are fixed until the model is rebuilt.
         if getattr(self, "_matrix_model", None) is not model:

@@ -107,7 +107,7 @@ def _integrate_single(self, model, t_eval, inputs_dict, y0):
         """
         inputs_dict = inputs_dict or {}
         if model.uses_stacked_inputs:
-            inputs = stack_inputs(inputs_dict)
+            inputs = stack_inputs(inputs_dict, model.convert_to_format)
         else:
             inputs = inputs_dict
 

@@ -17,6 +17,7 @@
 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
+from pybamm.solvers.rust_lowering import rust_input_parameter_widths
 
 
 class BaseSolver:
@@ -415,7 +416,7 @@ def _set_initial_conditions(self, model, time, inputs: list[dict]):
         )
         for ipts in inputs:
             if stacked:
-                inputs_y0_ics = stack_inputs(ipts)
+                inputs_y0_ics = stack_inputs(ipts, model.convert_to_format)
             else:
                 inputs_y0_ics = ipts
 
@@ -518,7 +519,11 @@ def _get_vars_for_processing(model, inputs: dict):
             "model": model,
         }
 
-        if model.convert_to_format != "casadi":
+        if model.convert_to_format == "rust":
+            vars_for_processing.update({"input_names": list(inputs.keys())})
+            return vars_for_processing
+
+        elif model.convert_to_format != "casadi":
             # Create Jacobian from concatenated rhs and algebraic
             y = pybamm.StateVector(slice(0, model.len_rhs_and_alg))
             # set up Jacobian object, for re-use of dict
@@ -1309,13 +1314,20 @@ def _get_discontinuity_start_end_indices(self, model, y0, inputs, t_eval):
             pybamm.logger.verbose("No additional discontinuity events found")
             return [0], [len(t_eval)], t_eval
 
-        stacked_inputs = stack_inputs(inputs) if model.uses_stacked_inputs else None
+        stacked_inputs = (
+            stack_inputs(inputs, model.convert_to_format)
+            if model.uses_stacked_inputs
+            else None
+        )
 
         def evaluate(function):
             if stacked_inputs is None:
                 return np.ravel(function(t=0, y=y0, inputs=inputs))
-            # Much faster than numpy's conversion of a CasADi matrix
-            return np.array(function(0, y0, stacked_inputs).elements())
+            value = function(0, y0, stacked_inputs)
+            if isinstance(value, casadi.DM):
+                # Much faster than numpy's conversion of a CasADi matrix
+                return np.array(value.elements())
+            return np.ravel(value)
 
         # The last float before and the first float after each discontinuity
         stops, restarts = [], []
@@ -1386,7 +1398,7 @@ def _check_event_violation(self, t_eval, model, y0, inputs_dict):
 
         inputs = inputs_dict
         if model.uses_stacked_inputs:
-            inputs = stack_inputs(inputs_dict)
+            inputs = stack_inputs(inputs_dict, model.convert_to_format)
 
         events_eval = np.empty(num_terminate_events)
         for idx, event in enumerate(model.terminate_events_eval):
@@ -2155,6 +2167,8 @@ def get_platform_context(self, system_type: str):
     def _set_up_model_inputs(model: pybamm.BaseModel, inputs: dict):
         """Set up input parameters"""
         inputs = ParameterValues.check_parameter_values(inputs)
+        if model.convert_to_format == "rust":
+            validate_rust_input_widths(model, inputs)
 
         # First, check all the strictly required input parameters that must be included
         # for the model to be solvable.
@@ -2203,6 +2217,19 @@ def _set_up_model_inputs(model: pybamm.BaseModel, inputs: dict):
             if name not in inputs_in_model:
                 inputs_in_model[name] = value
 
+        if model.convert_to_format == "rust":
+            # The dummy is one value, but the Rust backend packs each input at its width
+            widths = rust_input_parameter_widths(model)
+            vector_inputs = sorted(
+                missing for missing in missing_inputs if widths.get(missing, 1) > 1
+            )
+            if vector_inputs:
+                raise pybamm.SolverError(
+                    f"No value provided for vector input(s): {vector_inputs}. "
+                    "With convert_to_format='rust', every vector input needs a value "
+                    "of its expected size, even one that only variables use."
+                )
+
         if solution_is_observable and missing_inputs:
             missing_inputs_str = ", ".join(f"'{n}'" for n in missing_inputs)
             pybamm.logger.warning(
@@ -2218,6 +2245,40 @@ def _set_up_model_inputs(model: pybamm.BaseModel, inputs: dict):
         return ordered_inputs
 
 
+def validate_rust_input_widths(model: pybamm.BaseModel, inputs: dict) -> None:
+    """Reject a given input whose number of values is not its packed width."""
+    widths = rust_input_parameter_widths(model)
+    for name, value in inputs.items():
+        expected = widths.get(name)
+        if value is None or expected is None:
+            continue
+        if np.size(value) != expected:
+            raise pybamm.SolverError(
+                f"Input '{name}' has {np.size(value)} value(s), but the model "
+                f"expects {expected}."
+            )
+
+
+def validate_rust_sensitivity_widths(
+    model: pybamm.BaseModel, calculate_sensitivities: list[str]
+) -> None:
+    """Reject sensitivity requests for vector-width inputs on the Rust backend."""
+    widths = rust_input_parameter_widths(model)
+    wide_params = sorted(
+        name for name in calculate_sensitivities if widths.get(name, 1) > 1
+    )
+    if wide_params:
+        raise pybamm.SolverError(
+            "convert_to_format='rust' sensitivities are not supported "
+            f"for vector-width input parameters: {wide_params}. The Rust "
+            "backend's tangent/JVP machinery assumes one scalar direction "
+            "per named parameter, so seeding a sensitivity direction for "
+            "a vector input would silently compute the wrong derivative. "
+            "Use convert_to_format='casadi' for sensitivities with "
+            "respect to these inputs."
+        )
+
+
 def flatten_inputs(inputs_dict: dict) -> npt.NDArray[np.float64]:
     """Flatten ``{name: value}`` into one 1D array.
 
@@ -2244,19 +2305,26 @@ def flatten_inputs(inputs_dict: dict) -> npt.NDArray[np.float64]:
     )
 
 
-def stack_inputs(inputs_dict: dict) -> casadi.DM:
+def stack_inputs(
+    inputs_dict: dict, convert_to_format: str = "casadi"
+) -> casadi.DM | npt.NDArray[np.float64]:
     """Stack ``{name: value}`` into one column vector, in dict-key order.
 
     Parameters
     ----------
     inputs_dict : dict
         Input values for one input set.
+    convert_to_format : str, optional
+        The format of the evaluators the vector is for (default is "casadi").
 
     Returns
     -------
-    :class:`casadi.DM`
-        Column vector of the input values, with zero rows for an empty dict.
+    :class:`casadi.DM` or numpy.ndarray
+        The input values, with zero rows for an empty dict: a flat numpy array
+        for ``"rust"``, and a column :class:`casadi.DM` otherwise.
     """
+    if convert_to_format == "rust":
+        return flatten_inputs(inputs_dict)
     return casadi.vertcat(*inputs_dict.values())
 
 
@@ -2282,25 +2350,29 @@ def process(
     -------
     func: :class:`pybamm.EvaluatorPython` or
             :class:`pybamm.EvaluatorJax` or
-            :class:`casadi.Function`
+            :class:`casadi.Function` or
+            :class:`pybamm.solvers.rust_evaluator.RustEvaluator`
         evaluator for the function $f(y, t, p)$ given by `symbol`
 
     jac: :class:`pybamm.EvaluatorPython` or
             :class:`pybamm.EvaluatorJaxJacobian` or
-            :class:`casadi.Function`
+            :class:`casadi.Function` or
+            :class:`pybamm.solvers.rust_evaluator.RustEvaluator`
         evaluator for the Jacobian $\\frac{\\partial f}{\\partial y}$
         of the function given by `symbol`
 
     jacp: :class:`pybamm.EvaluatorPython` or
             :class:`pybamm.EvaluatorJaxSensitivities` or
-            :class:`casadi.Function`
+            :class:`casadi.Function` or
+            :class:`pybamm.solvers.rust_evaluator.RustEvaluator`
         evaluator for the parameter sensitivities
         $\frac{\\partial f}{\\partial p}$
         of the function given by `symbol`
 
     jac_action: :class:`pybamm.EvaluatorPython` or
             :class:`pybamm.EvaluatorJax` or
-            :class:`casadi.Function`
+            :class:`casadi.Function` or
+            :class:`pybamm.solvers.rust_evaluator.RustEvaluator`
         evaluator for product of the Jacobian with a vector $v$,
         i.e. $\\frac{\\partial f}{\\partial y} * v$
     """
@@ -2333,6 +2405,49 @@ def report(string):
             jac = None
             jac_action = None
 
+    elif model.convert_to_format == "rust":
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+        from pybamm.solvers.rust_lowering import rust_graph_with_inputs
+
+        if return_jacp_stacked:
+            raise pybamm.SolverError(
+                "return_jacp_stacked is not supported for convert_to_format='rust'. "
+                "Use the per-parameter tuple returned by jacp instead."
+            )
+
+        report(f"Converting {name} to Rust")
+        input_names = vars_for_processing["input_names"]
+
+        graph = rust_graph_with_inputs(model, input_names)
+        rust_expr = symbol.to_rust(graph, {})
+
+        # Full system width, so rectangular sub-groups differentiate wrt all states
+        cf = graph.compile(rust_expr, name=name, n_states=model.len_rhs_and_alg)
+        func = RustEvaluator(cf, "func")
+
+        if model.calculate_sensitivities:
+            report(
+                f"Calculating sensitivities for {name} with respect "
+                f"to parameters {model.calculate_sensitivities} using Rust"
+            )
+            # The graph's registry numbers the parameters; `input_names` only filters
+            sens_indices = []
+            for p in model.calculate_sensitivities:
+                index = graph.input_index(p)
+                if index is not None and p in input_names:
+                    sens_indices.append(index)
+            jacp = RustEvaluator(cf, "jacp", sens_indices=sens_indices)
+        else:
+            jacp = None
+
+        if use_jacobian:
+            report(f"Calculating jacobian for {name} using Rust")
+            jac = RustEvaluator(cf, "jac")
+            jac_action = RustEvaluator(cf, "jac_action")
+        else:
+            jac = None
+            jac_action = None
+
     elif model.convert_to_format != "casadi":
         y = vars_for_processing["y"]
         jacobian = vars_for_processing["jacobian"]

@@ -1074,7 +1074,7 @@ def _rhs_dot_consistent_initialization(self, y0, model, time, inputs_dict):
         """
         inputs_dict = inputs_dict or {}
         if model.uses_stacked_inputs:
-            input_eval = stack_inputs(inputs_dict)
+            input_eval = stack_inputs(inputs_dict, model.convert_to_format)
         else:
             input_eval = inputs_dict