From d158fdbcb472a4b2d37b36aa743d0096cc3b3055 Mon Sep 17 00:00:00 2001
From: bradyplanden <brady@ionworks.com>
Date: Tue, 22 Sep 2026 22:00:54 +0100
Subject: [PATCH] feat(rust): evaluate rust models through compiled evaluators

"rust" is a valid convert_to_format; BaseSolver.process builds
RustEvaluators; ScipySolver, AlgebraicSolver and EIS run on them.
---
 CHANGELOG.md                                  |   1 +
 packages/pybamm-rust/README.md                |   2 +-
 .../pybamm/src/pybamm/models/base_model.py    |  11 +-
 .../src/pybamm/simulation/eis_simulation.py   |  16 +-
 .../src/pybamm/solvers/algebraic_solver.py    |   2 +-
 .../pybamm/src/pybamm/solvers/base_solver.py  | 141 ++++++-
 .../src/pybamm/solvers/idaklu_solver.py       |   2 +-
 .../src/pybamm/solvers/rust_evaluator.py      | 102 +++++
 .../src/pybamm/solvers/rust_lowering.py       |  24 ++
 .../pybamm/src/pybamm/solvers/scipy_solver.py |   4 +-
 packages/pybamm/tests/conftest_rust.py        |   7 +
 .../test_convert_to_format_rust.py            |  31 ++
 .../pybamm/tests/unit/test_eis_simulation.py  |  12 +
 .../tests/unit/test_models/test_base_model.py |   8 +-
 .../test_serialisation/test_serialisation.py  |   4 +-
 .../test_casadi_algebraic_solver.py           |  13 +
 .../unit/test_solvers/test_casadi_solver.py   |  16 +
 .../unit/test_solvers/test_process_rust.py    | 358 ++++++++++++++++++
 .../unit/test_solvers/test_scipy_solver.py    |  47 ++-
 19 files changed, 768 insertions(+), 33 deletions(-)
 create mode 100644 packages/pybamm/src/pybamm/solvers/rust_evaluator.py
 create mode 100644 packages/pybamm/src/pybamm/solvers/rust_lowering.py
 create mode 100644 packages/pybamm/tests/integration/test_convert_to_format_rust.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1f78d8cd32..62e30f742c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -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
 
diff --git a/packages/pybamm-rust/README.md b/packages/pybamm-rust/README.md
index e93f0af43f..3605f861a5 100644
--- a/packages/pybamm-rust/README.md
+++ b/packages/pybamm-rust/README.md
@@ -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
 
diff --git a/packages/pybamm/src/pybamm/models/base_model.py b/packages/pybamm/src/pybamm/models/base_model.py
index 613d2b230b..2fba62ed5c 100644
--- a/packages/pybamm/src/pybamm/models/base_model.py
+++ b/packages/pybamm/src/pybamm/models/base_model.py
@@ -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):
diff --git a/packages/pybamm/src/pybamm/simulation/eis_simulation.py b/packages/pybamm/src/pybamm/simulation/eis_simulation.py
index 373b2ff872..218115fe1d 100644
--- a/packages/pybamm/src/pybamm/simulation/eis_simulation.py
+++ b/packages/pybamm/src/pybamm/simulation/eis_simulation.py
@@ -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:
diff --git a/packages/pybamm/src/pybamm/solvers/algebraic_solver.py b/packages/pybamm/src/pybamm/solvers/algebraic_solver.py
index a89bfaf48c..6764a779b8 100644
--- a/packages/pybamm/src/pybamm/solvers/algebraic_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/algebraic_solver.py
@@ -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
 
diff --git a/packages/pybamm/src/pybamm/solvers/base_solver.py b/packages/pybamm/src/pybamm/solvers/base_solver.py
index 1fc94719df..78bb1896cb 100644
--- a/packages/pybamm/src/pybamm/solvers/base_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/base_solver.py
@@ -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"]
diff --git a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
index d8260b3755..d85258757c 100644
--- a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
@@ -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
 
diff --git a/packages/pybamm/src/pybamm/solvers/rust_evaluator.py b/packages/pybamm/src/pybamm/solvers/rust_evaluator.py
new file mode 100644
index 0000000000..a56448444f
--- /dev/null
+++ b/packages/pybamm/src/pybamm/solvers/rust_evaluator.py
@@ -0,0 +1,102 @@
+"""CasADi-convention ``(t, y, stacked_inputs)`` callables over Rust tapes."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Literal
+
+import numpy as np
+import numpy.typing as npt
+
+import pybamm
+
+if TYPE_CHECKING:
+    from scipy.sparse import csc_matrix
+
+    from pybamm.rust import CompiledFunction, CompiledJacobian
+
+
+def _as_1d(arr: npt.ArrayLike) -> npt.NDArray[np.float64]:
+    return np.ascontiguousarray(np.asarray(arr, dtype=np.float64)).ravel()
+
+
+class RustEvaluator:
+    """Role-specific view over a per-group ``CompiledFunction``.
+
+    Parameters
+    ----------
+    cf : :class:`pybamm.rust.CompiledFunction`
+        The compiled function to evaluate.
+    role : str
+        ``"func"``, ``"jac"`` (scipy CSC), ``"jac_action"`` or ``"jacp"`` (one
+        ``df/dp_i`` column per ``sens_indices`` entry, as CasADi's ``jacp``).
+    sens_indices : list of int, optional
+        Graph input indices of the sensitivity parameters, read by ``"jacp"``.
+    """
+
+    def __init__(
+        self,
+        cf: CompiledFunction,
+        role: Literal["func", "jac", "jac_action", "jacp"],
+        sens_indices: list[int] | None = None,
+    ) -> None:
+        self._cf = cf
+        self._role = role
+        self._sens_indices = sens_indices
+        # Derived lazily on the first jac/jacp call
+        self._jac: CompiledJacobian | None = None
+
+    def __getstate__(self) -> dict:
+        state = self.__dict__.copy()
+        state["_jac"] = None
+        return state
+
+    def _jacobian(self, wrt: Literal["y", "p"]) -> CompiledJacobian:
+        if self._jac is None:
+            self._jac = self._cf.jacobian(wrt)
+        return self._jac
+
+    def __call__(
+        self,
+        t: float,
+        y: npt.ArrayLike,
+        inputs: npt.ArrayLike,
+        v: npt.ArrayLike | None = None,
+    ) -> npt.NDArray[np.float64] | csc_matrix | tuple[npt.NDArray[np.float64], ...]:
+        """Evaluate this evaluator's role at ``(t, y, inputs)``.
+
+        Parameters
+        ----------
+        t : float
+            The time.
+        y : array-like
+            The state vector.
+        inputs : array-like
+            The stacked input values.
+        v : array-like, optional
+            The vector a ``"jac_action"`` multiplies.
+
+        Returns
+        -------
+        numpy.ndarray, scipy.sparse.csc_matrix or tuple of numpy.ndarray
+            A column for ``"func"`` and ``"jac_action"``, the sparse Jacobian for
+            ``"jac"``, and one column per sensitivity parameter for ``"jacp"``.
+
+        Raises
+        ------
+        pybamm.SolverError
+            If a ``"jac_action"`` evaluator gets no ``v``, or a ``"jacp"`` evaluator
+            has no ``sens_indices``.
+        """
+        t, y, p = float(t), _as_1d(y), _as_1d(inputs)
+        if self._role == "func":
+            return np.asarray(self._cf(t, y, p)).reshape(-1, 1)
+        if self._role == "jac":
+            return self._jacobian("y")(t, y, p)
+        if self._role == "jac_action":
+            if v is None:
+                raise pybamm.SolverError("A 'jac_action' evaluator needs the vector v")
+            return np.asarray(self._cf.jvp(t, y, p, _as_1d(v))).reshape(-1, 1)
+        if self._sens_indices is None:
+            raise pybamm.SolverError("A 'jacp' evaluator needs sens_indices")
+        dense = self._jacobian("p")(t, y, p).toarray()
+        return tuple(dense[:, i].reshape(-1, 1) for i in self._sens_indices)
diff --git a/packages/pybamm/src/pybamm/solvers/rust_lowering.py b/packages/pybamm/src/pybamm/solvers/rust_lowering.py
new file mode 100644
index 0000000000..fc5355301c
--- /dev/null
+++ b/packages/pybamm/src/pybamm/solvers/rust_lowering.py
@@ -0,0 +1,24 @@
+"""Rust expression graphs for a discretised model."""
+
+from __future__ import annotations
+
+import pybamm
+
+
+def rust_input_parameter_widths(model: pybamm.BaseModel) -> dict[str, int]:
+    """Map each input parameter name to its packed width for the Rust backend."""
+    return {ip.name: ip._expected_size or 1 for ip in model.input_parameters}
+
+
+def rust_graph_with_inputs(model: pybamm.BaseModel, input_names):
+    """Build an ``ExprGraph`` with ``input_names`` registered in stacking order.
+
+    Registration precedes any lowering so graph input indices match that order.
+    """
+    from pybamm.rust import ExprGraph
+
+    graph = ExprGraph()
+    widths = rust_input_parameter_widths(model)
+    for name in input_names:
+        graph.input_parameter(name, widths.get(name, 1))
+    return graph
diff --git a/packages/pybamm/src/pybamm/solvers/scipy_solver.py b/packages/pybamm/src/pybamm/solvers/scipy_solver.py
index e18914dfce..0e4a3299b3 100644
--- a/packages/pybamm/src/pybamm/solvers/scipy_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/scipy_solver.py
@@ -98,7 +98,7 @@ def _integrate_single(self, model, t_eval, inputs_dict, y0):
         # Save inputs dictionary, and if necessary convert inputs to a stacked vector
         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
 
@@ -134,7 +134,7 @@ def rhs(t, y):
 
             def event_wrapper(event):
                 def event_fn(t, y):
-                    return event(t, y, inputs)
+                    return np.asarray(event(t, y, inputs)).item()
 
                 event_fn.terminal = True
                 return event_fn
diff --git a/packages/pybamm/tests/conftest_rust.py b/packages/pybamm/tests/conftest_rust.py
index e11c600a91..5c189bf344 100644
--- a/packages/pybamm/tests/conftest_rust.py
+++ b/packages/pybamm/tests/conftest_rust.py
@@ -6,6 +6,13 @@
 import pybamm
 
 
+@pytest.fixture
+def rust_backend(monkeypatch):
+    """Make newly constructed models default to the Rust backend."""
+    monkeypatch.setattr(pybamm.BaseModel, "_DEFAULT_CONVERT_TO_FORMAT", "rust")
+    yield
+
+
 def evaluate_with_rust(expr, t=0.0, y=None, y_dot=None, inputs=None):
     """Evaluate a PyBaMM expression using the Rust backend."""
     from pybamm.rust import ExprGraph
diff --git a/packages/pybamm/tests/integration/test_convert_to_format_rust.py b/packages/pybamm/tests/integration/test_convert_to_format_rust.py
new file mode 100644
index 0000000000..05c3a4a144
--- /dev/null
+++ b/packages/pybamm/tests/integration/test_convert_to_format_rust.py
@@ -0,0 +1,31 @@
+import pytest
+
+import pybamm
+
+
+class TestAlgebraicRustParity:
+    def test_pure_algebraic_rust_vs_casadi(self):
+        results = {}
+        for fmt in ("casadi", "rust"):
+            model = pybamm.BaseModel()
+            v = pybamm.Variable("v")
+            model.algebraic = {v: v**2 - 3 * v + 2}
+            model.initial_conditions = {v: 3.0}
+            disc = pybamm.Discretisation()
+            disc.process_model(model)
+            model.convert_to_format = fmt
+            sol = pybamm.AlgebraicSolver().solve(model, [0])
+            results[fmt] = sol.y[0][0]
+        assert results["rust"] == pytest.approx(results["casadi"], rel=1e-8)
+        assert results["rust"] == pytest.approx(2.0, rel=1e-6)
+
+
+class TestRustRootResolution:
+    def test_default_root_method_keeps_the_rust_format(self):
+        from tests.unit.test_solvers.test_process_rust import _toy_dae
+
+        model = _toy_dae("rust")
+        solver = pybamm.IDAKLUSolver(options={"calc_ic": False})
+        solver._check_and_prepare_model_inplace(model)
+        assert model.convert_to_format == "rust"
+        assert isinstance(solver.root_method, pybamm.NonlinearSolver)
diff --git a/packages/pybamm/tests/unit/test_eis_simulation.py b/packages/pybamm/tests/unit/test_eis_simulation.py
index 4080a52e66..81191bbe7e 100644
--- a/packages/pybamm/tests/unit/test_eis_simulation.py
+++ b/packages/pybamm/tests/unit/test_eis_simulation.py
@@ -313,6 +313,18 @@ def test_initial_soc_forwards_inputs(self):
         z_v = eis_sim.solve(frequencies, inputs={"vf_solid": 0.94}, initial_soc="0.5 V")
         assert z_v.impedance.shape == (5,)
 
+    def test_eis_rust_backend_matches_casadi(self):
+        def impedance(fmt):
+            model = pybamm.lithium_ion.SPM(
+                options={"surface form": "differential"}, name="SPM"
+            )
+            model.convert_to_format = fmt
+            eis_sim = pybamm.EISSimulation(model)
+            frequencies = np.logspace(-2, 2, 10)
+            return eis_sim.solve(frequencies).impedance
+
+        np.testing.assert_allclose(impedance("rust"), impedance("casadi"), rtol=1e-6)
+
     def test_high_freq_intercept_matches_contact_resistance(self):
         model = pybamm.lithium_ion.SPM(
             options={"surface form": "differential", "contact resistance": "true"}
diff --git a/packages/pybamm/tests/unit/test_models/test_base_model.py b/packages/pybamm/tests/unit/test_models/test_base_model.py
index afa7b1aa75..844e7008b0 100644
--- a/packages/pybamm/tests/unit/test_models/test_base_model.py
+++ b/packages/pybamm/tests/unit/test_models/test_base_model.py
@@ -16,6 +16,8 @@
 
 import pybamm
 
+pytest_plugins = ["tests.conftest_rust"]
+
 
 class TestBaseModel:
     def test_rhs_set_get(self):
@@ -1792,7 +1794,7 @@ def test_y0_property(self):
 
     def test_convert_to_format_validated(self):
         model = pybamm.BaseModel()
-        for valid in ("python", "casadi", "jax"):
+        for valid in ("python", "casadi", "jax", "rust"):
             model.convert_to_format = valid
             assert model.convert_to_format == valid
         with pytest.raises(
@@ -1813,6 +1815,7 @@ def test_uses_stacked_inputs(self):
         model = pybamm.BaseModel()
         for convert_to_format, stacked in [
             ("casadi", True),
+            ("rust", True),
             ("python", False),
             ("jax", False),
         ]:
@@ -1845,3 +1848,6 @@ def public_name_state(self):
         assert "convert_to_format" not in vars(restored)
         solution = pybamm.IDAKLUSolver().solve(restored, [0, 1])
         np.testing.assert_allclose(solution["u"](1.0), np.exp(-1.0), rtol=1e-3)
+
+    def test_rust_backend_fixture_drives_default(self, rust_backend):
+        assert pybamm.BaseModel().convert_to_format == "rust"
diff --git a/packages/pybamm/tests/unit/test_serialisation/test_serialisation.py b/packages/pybamm/tests/unit/test_serialisation/test_serialisation.py
index f3d2789cad..2fe4c64966 100644
--- a/packages/pybamm/tests/unit/test_serialisation/test_serialisation.py
+++ b/packages/pybamm/tests/unit/test_serialisation/test_serialisation.py
@@ -242,7 +242,7 @@ def test_save_load_model(self, mocker, tmp_path, request):
         newest_solver = newest_model.default_solver
         newest_solver.solve(newest_model, [0, 3600])
 
-    @pytest.mark.parametrize("convert_to_format", ["python", "jax"])
+    @pytest.mark.parametrize("convert_to_format", ["python", "jax", "rust"])
     def test_discretised_model_round_trips_convert_to_format(self, convert_to_format):
         model = pybamm.BaseModel()
         a = pybamm.Variable("a")
@@ -261,7 +261,7 @@ def test_discretised_model_round_trips_convert_to_format(self, convert_to_format
         del data["convert_to_format"]
         assert Serialise().load_model(data).convert_to_format == "casadi"
 
-    @pytest.mark.parametrize("convert_to_format", ["python", "jax"])
+    @pytest.mark.parametrize("convert_to_format", ["python", "jax", "rust"])
     def test_custom_model_round_trips_convert_to_format(self, convert_to_format):
         model = pybamm.BaseModel()
         a = pybamm.Variable("a")
diff --git a/packages/pybamm/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/packages/pybamm/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index 4aa36d92e3..beadb8894a 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -36,6 +36,19 @@ def test_simple_root_find(self):
         solution = solver.solve(model, np.linspace(0, 1, 10))
         np.testing.assert_array_equal(solution.y, -2)
 
+    def test_rust_model_is_converted_to_casadi(self):
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + 2}
+        model.initial_conditions = {var: 2}
+        model.convert_to_format = "rust"
+        pybamm.Discretisation().process_model(model)
+
+        solution = pybamm.CasadiAlgebraicSolver().solve(model, np.linspace(0, 1, 10))
+
+        assert model.convert_to_format == "casadi"
+        np.testing.assert_array_equal(solution.y, -2)
+
     def test_algebraic_root_solver_reuse(self):
         # Create a model with an input parameter and
         # check that the algebraic root solver is reused
diff --git a/packages/pybamm/tests/unit/test_solvers/test_casadi_solver.py b/packages/pybamm/tests/unit/test_solvers/test_casadi_solver.py
index a941074446..08acb90a29 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_casadi_solver.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_casadi_solver.py
@@ -27,6 +27,22 @@ def test_bad_mode(self):
         with pytest.raises(ValueError, match=r"invalid mode"):
             pybamm.CasadiSolver(mode="bad mode")
 
+    def test_rust_model_is_converted_to_casadi(self):
+        model = pybamm.BaseModel()
+        var = pybamm.Variable("var")
+        model.rhs = {var: 0.1 * var}
+        model.initial_conditions = {var: 1}
+        model.convert_to_format = "rust"
+        pybamm.Discretisation().process_model(model)
+
+        t_eval = np.linspace(0, 1, 10)
+        solution = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8).solve(model, t_eval)
+
+        assert model.convert_to_format == "casadi"
+        np.testing.assert_allclose(
+            solution.y.full()[0], np.exp(0.1 * t_eval), rtol=1e-6, atol=1e-8
+        )
+
     def test_model_solver(self):
         # Create model
         model = pybamm.BaseModel()
diff --git a/packages/pybamm/tests/unit/test_solvers/test_process_rust.py b/packages/pybamm/tests/unit/test_solvers/test_process_rust.py
index ec2245809a..48cd2fa081 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_process_rust.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_process_rust.py
@@ -1,8 +1,11 @@
+import logging
+
 import numpy as np
 import pytest
 import scipy.sparse
 
 import pybamm
+from tests import get_discretisation_for_testing
 
 casadi = pytest.importorskip("casadi")
 
@@ -159,6 +162,25 @@ def test_named_outputs_eval_trajectory_hermite_accepts_loose_t_query(self):
             np.testing.assert_allclose(np.asarray(out), np.asarray(expected))
 
 
+class TestDenseMatrixToRust:
+    def test_dense_matrix_matmul_to_rust(self):
+        from pybamm.rust import ExprGraph
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        A = np.arange(6.0).reshape(2, 3)
+        y = pybamm.StateVector(slice(0, 3))
+        expr = pybamm.Matrix(A) @ y
+        graph = ExprGraph()
+        cf = graph.compile(expr.to_rust(graph, {}), name="f", n_states=3)
+        yv = np.array([1.0, 2.0, 3.0])
+        np.testing.assert_allclose(
+            np.asarray(cf(0.0, yv, np.array([]))).ravel(), A @ yv
+        )
+        # Jacobian path: d(A @ y)/dy == A.
+        ev = RustEvaluator(cf, "jac")
+        np.testing.assert_allclose(ev(0.0, yv, np.array([])).toarray(), A)
+
+
 class TestMatMulLhsValidation:
     def test_matmul_non_constant_lhs_raises_clean_error(self):
         from pybamm.rust import ExprGraph
@@ -193,6 +215,328 @@ def test_incompatible_binary_widths_raise_before_evaluation(self):
             graph.compile(invalid, name="invalid", n_states=5)
 
 
+class TestRustEvaluatorWrappers:
+    def _cf(self):
+        from pybamm.rust import ExprGraph
+
+        graph = ExprGraph()
+        a = graph.input_parameter("a")
+        y0 = graph.state_vector(0, 1)
+        y1 = graph.state_vector(1, 2)
+        expr = graph.concat([y0 * a - y1, y0 + y1 * y1])
+        return graph.compile(expr, name="rhs", n_states=2)
+
+    def test_func_returns_column(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        out = RustEvaluator(cf, "func")(0.0, np.array([[1.5], [-2.0]]), np.array([3.0]))
+        assert out.shape == (2, 1)
+        np.testing.assert_allclose(out.ravel(), [6.5, 5.5])
+
+    def test_jac_returns_scipy_csc(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        jac = RustEvaluator(cf, "jac")(0.0, np.array([1.5, -2.0]), np.array([3.0]))
+        assert scipy.sparse.issparse(jac)
+        np.testing.assert_allclose(jac.toarray(), [[3.0, -1.0], [1.0, -4.0]])
+
+    def test_jac_action_matches_jac_matvec(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        y, p, v = np.array([1.5, -2.0]), np.array([3.0]), np.array([0.3, 0.7])
+        jv = RustEvaluator(cf, "jac_action")(0.0, y, p, v)
+        full = RustEvaluator(cf, "jac")(0.0, y, p).toarray() @ v
+        np.testing.assert_allclose(jv.ravel(), full)
+
+    def test_jacp_tuple_per_param(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        out = RustEvaluator(cf, "jacp", sens_indices=[0])(
+            0.0, np.array([1.5, -2.0]), np.array([3.0])
+        )
+        assert isinstance(out, tuple) and len(out) == 1
+        np.testing.assert_allclose(out[0].ravel(), [1.5, 0.0])
+
+    def test_jac_action_without_v_and_jacp_without_indices_raise(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        y, p = np.array([1.5, -2.0]), np.array([3.0])
+        with pytest.raises(pybamm.SolverError, match="needs the vector v"):
+            RustEvaluator(cf, "jac_action")(0.0, y, p)
+        with pytest.raises(pybamm.SolverError, match="needs sens_indices"):
+            RustEvaluator(cf, "jacp")(0.0, y, p)
+
+    def test_jacobian_is_derived_lazily(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        ev = RustEvaluator(cf, "jac")
+        assert ev._jac is None
+        ev(0.0, np.array([1.5, -2.0]), np.array([3.0]))
+        assert ev._jac is not None
+
+    def test_func_accepts_eval_only(self):
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        cf = self._cf()
+        out = RustEvaluator(cf, "func")(0.0, np.array([1.5, -2.0]), np.array([3.0]))
+        assert out.shape == (2, 1)
+        np.testing.assert_allclose(out.ravel(), [6.5, 5.5])
+
+
+def _jacp_as_list(result):
+    """Normalise a jacp result to a list of 1-D arrays (casadi returns a bare DM)."""
+    if isinstance(result, casadi.DM):
+        result = (result,)
+    return [np.asarray(j).ravel() for j in result]
+
+
+def _toy_dae(convert_to_format):
+    model = pybamm.BaseModel()
+    u = pybamm.Variable("u")
+    v = pybamm.Variable("v")
+    a = pybamm.InputParameter("a")
+    model.rhs = {u: -2 * u + a * v}
+    model.algebraic = {v: 2 * u - v}
+    model.initial_conditions = {u: 1.0, v: 2.0}
+    model.events = [pybamm.Event("u-0.05", u - 0.05)]
+    model.variables = {"u": u, "v": v}
+    disc = pybamm.Discretisation()
+    disc.process_model(model)
+    model.convert_to_format = convert_to_format
+    return model
+
+
+class TestProcessRustParity:
+    @pytest.fixture
+    def pair(self):
+        from pybamm.solvers.base_solver import BaseSolver
+
+        inputs = {"a": 0.5}
+        out = {}
+        for fmt in ("casadi", "rust"):
+            model = _toy_dae(fmt)
+            model.calculate_sensitivities = ["a"]
+            vars_ = BaseSolver._get_vars_for_processing(model, inputs)
+            rhs_alg = pybamm.numpy_concatenation(
+                model.concatenated_rhs, model.concatenated_algebraic
+            )
+            out[fmt] = (model, vars_, rhs_alg)
+        return out
+
+    def _stack(self, fmt, inputs):
+        if fmt == "casadi":
+            return casadi.vertcat(*inputs.values())
+        return np.array(list(inputs.values()), dtype=np.float64)
+
+    def test_rhs_algebraic_func_jac_action_jacp(self, pair):
+        from pybamm.solvers.base_solver import process
+
+        t, y, inputs = 0.3, np.array([1.1, 0.4]), {"a": 0.5}
+        results = {}
+        for fmt, (_model, vars_, rhs_alg) in pair.items():
+            func, jac, jacp, jac_action = process(rhs_alg, "rhs_algebraic", vars_)
+            p = self._stack(fmt, inputs)
+            jac_mat = jac(t, y, p)
+            results[fmt] = {
+                "f": np.asarray(func(t, y, p)).ravel(),
+                "jac": np.asarray(
+                    jac_mat.toarray() if hasattr(jac_mat, "toarray") else jac_mat
+                ),
+                "jv": np.asarray(jac_action(t, y, p, np.array([0.2, -0.3]))).ravel(),
+                "jacp": _jacp_as_list(jacp(t, y, p)),
+            }
+        for key in ("f", "jac", "jv"):
+            np.testing.assert_allclose(
+                results["rust"][key], results["casadi"][key], rtol=1e-12, atol=1e-14
+            )
+        for r, c in zip(
+            results["rust"]["jacp"], results["casadi"]["jacp"], strict=True
+        ):
+            np.testing.assert_allclose(r, c, rtol=1e-12, atol=1e-14)
+
+    def test_event_and_ic_eval_only(self, pair):
+        from pybamm.solvers.base_solver import process
+
+        t, y, inputs = 0.0, np.array([1.0, 2.0]), {"a": 0.5}
+        vals = {}
+        for fmt, (model, vars_, _) in pair.items():
+            ev = process(
+                model.events[0].expression, "event_0", vars_, use_jacobian=False
+            )[0]
+            ic = process(
+                model.concatenated_initial_conditions,
+                "initial_conditions",
+                vars_,
+                use_jacobian=False,
+            )[0]
+            p = self._stack(fmt, inputs)
+            vals[fmt] = (
+                float(np.asarray(ev(t, y, p)).item()),
+                np.asarray(ic(t, np.zeros((2, 1)), p)).ravel(),
+            )
+        assert vals["rust"][0] == pytest.approx(vals["casadi"][0], rel=1e-12)
+        np.testing.assert_allclose(vals["rust"][1], vals["casadi"][1], rtol=1e-12)
+
+    def test_rectangular_group_jacobian_matches_casadi(self, pair):
+        from pybamm.solvers.base_solver import process
+
+        t, y = 0.3, np.array([1.1, 0.4])
+        mats = {}
+        for fmt, (model, vars_, _) in pair.items():
+            _, jac, _, jac_action = process(model.concatenated_rhs, "RHS", vars_)
+            assert jac is not None and jac_action is not None
+            p = self._stack(fmt, inputs={"a": 0.5})
+            jm = jac(t, y, p)
+            mats[fmt] = np.asarray(jm.toarray() if hasattr(jm, "toarray") else jm)
+        np.testing.assert_allclose(mats["rust"], mats["casadi"], rtol=1e-12, atol=1e-14)
+
+
+class TestProcessRustDuplication:
+    def test_lazy_jac_not_built_for_uncalled_groups(self):
+        from pybamm.solvers.base_solver import BaseSolver, process
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        model = _toy_dae("rust")
+        model.calculate_sensitivities = []
+        vars_ = BaseSolver._get_vars_for_processing(model, {"a": 0.5})
+        _, jac, _, jac_action = process(model.concatenated_rhs, "RHS", vars_)
+        assert (
+            isinstance(jac, RustEvaluator)
+            and jac._jac is None
+            and jac_action is not None
+        )
+
+
+class TestStackedInputPredicate:
+    def test_stack_inputs_rust_is_ndarray(self):
+        from pybamm.solvers.base_solver import stack_inputs
+
+        out = stack_inputs({"a": 1.0, "b": np.array([2.0, 3.0])}, "rust")
+        np.testing.assert_allclose(out, [1.0, 2.0, 3.0])
+        assert stack_inputs({}, "rust").size == 0
+
+    def test_stack_inputs_defaults_to_a_casadi_column(self):
+        from pybamm.solvers.base_solver import stack_inputs
+
+        out = stack_inputs({"a": 1.0, "b": np.array([2.0, 3.0])})
+        assert isinstance(out, casadi.DM)
+        np.testing.assert_allclose(out.full(), [[1.0], [2.0], [3.0]])
+
+    def test_set_initial_conditions_and_event_check_rust(self):
+        from pybamm.solvers.base_solver import BaseSolver
+
+        class _DaeSolver(BaseSolver):
+            def _run(self, *a, **kw):  # pragma: no cover
+                raise NotImplementedError
+
+        model = _toy_dae("rust")
+        solver = _DaeSolver()
+        solver.set_up(model, inputs=[{"a": 0.5}], ics_only=True)
+        y0 = np.asarray(model.y0_list[0]).ravel()
+        np.testing.assert_allclose(y0, [1.0, 2.0])
+        model.terminate_events_eval = []
+        solver._check_event_violation([0.0], model, y0, {"a": 0.5})
+
+
+class TestRustSolverGuards:
+    def test_casadi_solver_converts_rust_model_with_a_warning(self, caplog):
+        model = _toy_dae("rust")
+        with caplog.at_level(logging.WARNING, logger=pybamm.logger.name):
+            pybamm.CasadiSolver()._check_and_prepare_model_inplace(model)
+        assert model.convert_to_format == "casadi"
+        assert "Converting" in caplog.text
+
+    @pytest.mark.parametrize("width", [1, 3])
+    def test_an_input_of_another_width_is_a_solver_error(self, width):
+        model = pybamm.BaseModel()
+        v = pybamm.Variable("v")
+        b = pybamm.InputParameter("b", expected_size=2)
+        model.algebraic = {v: v - pybamm.Index(b, 0) - pybamm.Index(b, 1)}
+        model.initial_conditions = {v: 0}
+        pybamm.Discretisation().process_model(model)
+        model.convert_to_format = "rust"
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model, [0], inputs={"b": np.array([1.0, 2.0])})
+        np.testing.assert_allclose(solution.y, 3.0, rtol=1e-6, atol=1e-6)
+
+        with pytest.raises(
+            pybamm.SolverError,
+            match=rf"Input 'b' has {width} value\(s\), but the model expects 2",
+        ):
+            solver.solve(model, [0], inputs={"b": np.ones(width)})
+
+    def test_a_vector_input_only_variables_use_must_be_given(self):
+        model = pybamm.BaseModel()
+        v = pybamm.Variable("v")
+        model.algebraic = {v: v - 1}
+        model.initial_conditions = {v: 0}
+        model.variables = {"v plus b": v + pybamm.Parameter("b")}
+        parameter_values = pybamm.ParameterValues(
+            {"b": pybamm.InputParameter("b", expected_size=2)}
+        )
+        parameter_values.process_model(model)
+        pybamm.Discretisation().process_model(model)
+        model.convert_to_format = "rust"
+        solver = pybamm.AlgebraicSolver()
+
+        with pytest.raises(
+            pybamm.SolverError,
+            match=r"No value provided for vector input\(s\): \['b'\]",
+        ):
+            solver.solve(model, [0])
+        solution = solver.solve(model, [0], inputs={"b": np.array([1.0, 2.0])})
+        np.testing.assert_allclose(solution.y, 1.0, rtol=1e-6, atol=1e-6)
+
+    def test_return_jacp_stacked_is_rejected_before_compiling(self, monkeypatch):
+        from pybamm.solvers import rust_lowering
+        from pybamm.solvers.base_solver import BaseSolver, process
+
+        model = _toy_dae("rust")
+        vars_ = BaseSolver._get_vars_for_processing(model, {"a": 0.5})
+
+        def no_lowering(*args, **kwargs):  # pragma: no cover
+            raise AssertionError("lowered before rejecting return_jacp_stacked")
+
+        monkeypatch.setattr(rust_lowering, "rust_graph_with_inputs", no_lowering)
+        with pytest.raises(pybamm.SolverError, match=r"return_jacp_stacked"):
+            process(model.concatenated_rhs, "RHS", vars_, return_jacp_stacked=True)
+
+
+class TestRustAlgebraicMinimize:
+    def test_minimize_converges_on_rust_model(self):
+        # jac_norm's broadcasting must hold for numpy jacobians as well as casadi DM.
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        model.algebraic = {var1: var1 - 3, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+        model.convert_to_format = "rust"
+
+        solution = pybamm.AlgebraicSolver("minimize", tol=1e-8).solve(model)
+        np.testing.assert_allclose(
+            model.get_processed_variable("var1").evaluate(t=None, y=solution.y),
+            3,
+            rtol=1e-7,
+            atol=1e-6,
+        )
+        np.testing.assert_allclose(
+            model.get_processed_variable("var2").evaluate(t=None, y=solution.y),
+            6,
+            rtol=1e-7,
+            atol=1e-6,
+        )
+
+
 class TestVectorWidthInputParameter:
     def test_default_width_is_scalar(self):
         from pybamm.rust import ExprGraph
@@ -359,3 +703,17 @@ def test_rebuild_rejects_unsupported_root(self):
         mm = graph.matmul(graph.state_vector(0, 2), graph.state_vector(0, 2))
         with pytest.raises(NotImplementedError, match=r"MatMul left operand"):
             CompiledFunction._rebuild(graph, mm.id, None, None)
+
+    def test_rust_evaluator_pickle_rederives_jacobian(self):
+        import pickle
+
+        from pybamm.rust import ExprGraph
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        graph = ExprGraph()
+        expr = (3 * pybamm.StateVector(slice(0, 2))).to_rust(graph, {})
+        ev = RustEvaluator(graph.compile(expr, name="f", n_states=2), "jac")
+        y, p = np.array([1.0, 2.0]), np.array([])
+        expected = ev(0.0, y, p).toarray()  # populates the _jac cache
+        ev2 = pickle.loads(pickle.dumps(ev))
+        np.testing.assert_allclose(ev2(0.0, y, p).toarray(), expected)
diff --git a/packages/pybamm/tests/unit/test_solvers/test_scipy_solver.py b/packages/pybamm/tests/unit/test_solvers/test_scipy_solver.py
index 99779c68ac..73ea321937 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_scipy_solver.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_scipy_solver.py
@@ -372,10 +372,53 @@ def test_model_solver_with_event_with_casadi(self):
             )
             np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
-    def test_model_solver_with_inputs_with_casadi(self):
+    def test_rust_evaluators_receive_numpy_inputs(self, monkeypatch):
+        # A CasADi DM would be converted to numpy again on every evaluation
+        from pybamm.solvers.rust_evaluator import RustEvaluator
+
+        input_types = []
+        evaluate = RustEvaluator.__call__
+
+        def recording_call(self, t, y, inputs, v=None):
+            input_types.append(type(inputs))
+            return evaluate(self, t, y, inputs, v)
+
+        monkeypatch.setattr(RustEvaluator, "__call__", recording_call)
+        model = pybamm.BaseModel()
+        model.convert_to_format = "rust"
+        var = pybamm.Variable("var")
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = [pybamm.Event("var=0.5", var - 0.5)]
+        pybamm.Discretisation().process_model(model)
+
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8)
+        solver.solve(model, np.linspace(0, 10, 20), inputs={"rate": 0.1})
+
+        assert input_types
+        assert set(input_types) == {np.ndarray}
+
+    def test_rust_model_with_an_input_dependent_discontinuity(self):
+        # The switch time is evaluated through Rust, with the inputs stacked for it
+        model = pybamm.BaseModel()
+        model.convert_to_format = "rust"
+        var = pybamm.Variable("var")
+        t_switch = pybamm.InputParameter("t_switch")
+        model.rhs = {var: -(1 + (pybamm.t > t_switch)) * var}
+        model.initial_conditions = {var: 1}
+        pybamm.Discretisation().process_model(model)
+
+        solver = pybamm.ScipySolver(rtol=1e-10, atol=1e-10)
+        solution = solver.solve(model, np.linspace(0, 1, 11), inputs={"t_switch": 0.45})
+        t = solution.t
+        expected = np.where(t <= 0.45, np.exp(-t), np.exp(-0.45 - 2 * (t - 0.45)))
+        np.testing.assert_allclose(solution.y[0], expected, rtol=1e-7, atol=1e-9)
+
+    @pytest.mark.parametrize("convert_to_format", ["casadi", "rust"])
+    def test_model_solver_with_inputs_and_event(self, convert_to_format):
         # Create model
         model = pybamm.BaseModel()
-        model.convert_to_format = "casadi"
+        model.convert_to_format = convert_to_format
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: -pybamm.InputParameter("rate") * var}