BradyPlanden · Sep 22, 2026
@@ -4,6 +4,7 @@
 
 - `pybamm` now requires `pybammsolvers>=0.12.0`, whose IDAKLU backend can evaluate a Rust `CompiledModel`. Upgrade `pybammsolvers` together with `pybamm`. ([#5961](https://github.com/pybamm-team/PyBaMM/pull/5961))
 - `pybamm` now ships platform wheels with its Rust core compiled in as `pybamm.rust._core`: one `cp310-abi3` wheel per platform for Python 3.10-3.15. Building from source or the sdist requires `cargo` 1.89+; `pip install --only-binary pybamm pybamm` forces a wheel. ([#5958](https://github.com/pybamm-team/PyBaMM/pull/5958))
+- Models now default to `convert_to_format = "rust"` instead of `"casadi"`, so every model is lowered to the Rust compute core unless it opts out. Set `model.convert_to_format = "casadi"` to restore the previous backend. ([#5966](https://github.com/pybamm-team/PyBaMM/pull/5966))
 
 ## Deprecated
 

@@ -163,3 +163,9 @@ uv sync --extra all --group dev --reinstall-package pybamm
 
 An ABI mismatch between components — rebuild both extensions from the same
 tree (see "Rebuilding after changes").
+
+## Supported features
+
+Feature coverage (events, DAEs, forward sensitivities, output variables,
+experiments) and the expression types the backend can lower are documented in
+the [Rust backend user guide](../user_guide/rust_backend.rst).
@@ -183,6 +183,22 @@
    "source": [
     "So in this case the speed-up for using multiple threads to solve 1000 SPM simulations is much less than for the DFN simulations, and above 4 threads no speed-up is observed at all."
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The same option on the other backends\n",
+    "\n",
+    "Given several input sets, `num_threads` means the same thing wherever it appears: *this many input sets solve at once*.\n",
+    "\n",
+    "- `IDAKLUSolver` schedules them with OpenMP, one solver per thread, for both `convert_to_format=\"casadi\"` and `\"rust\"`. Work is handed out dynamically, so a sweep whose sets terminate at different times (a current sweep hitting a voltage cut-off, say) does not leave threads idle waiting on the slowest block.\n",
+    "- `DiffsolSolver` takes the same option and schedules with rayon: `pybamm.DiffsolSolver(options={\"num_threads\": 8})`. Pools are shared process-wide by thread count, so ten solvers asking for 8 threads share one pool of 8 rather than spawning eighty. With a single input set it gives the threads to that set's sparse LU instead, which pays off only for large models (around 50k states and up).\n",
+    "\n",
+    "Either way, leave `num_threads` at its default of 1 if something outside PyBaMM is already parallelising over your sweep — a `ThreadPoolExecutor`, joblib, or a parameter-estimation library — or the two thread counts multiply.\n",
+    "\n",
+    "The timings above are from one machine and one PyBaMM version; treat them as the shape of the curve rather than numbers to reproduce."
+   ]
   }
  ],
  "metadata": {

@@ -23,6 +23,7 @@ maxdepth: 2
 fundamentals/index
 fundamentals/battery_models
 fundamentals/public_api
+rust_backend
 ```
 
 ```{toctree}

@@ -0,0 +1,168 @@
+.. _rust-backend:
+
+Rust Backend
+============
+
+PyBaMM ships a Rust compute core (``pybamm.rust``) that compiles a discretised
+model's expression trees: right-hand sides, algebraic residuals, Jacobians,
+events, and output variables into a form evaluated entirely outside Python.
+Two solvers use it:
+
+- :class:`pybamm.IDAKLUSolver` runs the usual SUNDIALS IDA integrator with
+  every callback evaluated by the Rust core instead of CasADi. It supports
+  forward sensitivities, ``output_variables``, events, and experiments, with
+  the differences listed in `Differences from the CasADi backend`_.
+- :class:`pybamm.DiffsolSolver` integrates with the pure-Rust BDF
+  implementation from the `diffsol <https://github.com/martinjrobins/diffsol>`_
+  crate, with no CasADi or SUNDIALS involvement. It supports forward
+  sensitivities under solver error control, ``output_variables``-only solves,
+  events, and ``t_interp``.
+
+Selecting the backend
+---------------------
+
+Which backend a solver uses follows the model's ``convert_to_format``
+attribute (documented on :class:`pybamm.BaseModel`), which defaults to
+``"rust"``:
+
+.. code-block:: python
+
+    import pybamm
+
+    model = pybamm.lithium_ion.SPM()
+
+    sim = pybamm.Simulation(model)  # the default IDAKLUSolver, on the Rust core
+    sol = sim.solve([0, 3600])
+
+Set ``convert_to_format = "casadi"`` before solving to use the CasADi backend
+instead:
+
+.. code-block:: python
+
+    model = pybamm.lithium_ion.SPM()
+    model.convert_to_format = "casadi"
+
+:class:`pybamm.DiffsolSolver` always uses the Rust backend and converts the
+model itself:
+
+.. code-block:: python
+
+    sim = pybamm.Simulation(model, solver=pybamm.DiffsolSolver())
+    sol = sim.solve([0, 3600])
+
+Differences from the CasADi backend
+-----------------------------------
+
+On the Rust backend, :class:`pybamm.IDAKLUSolver`:
+
+- raises :class:`pybamm.SolverError` for sensitivities with respect to an input
+  parameter that holds several values (``pybamm.InputParameter`` with an
+  ``expected_size`` above 1);
+- runs one solver per thread, so it solves with ``num_solvers`` equal to
+  ``num_threads``, with a warning when they differ.
+
+Set ``convert_to_format = "casadi"`` for either.
+
+Output times with ``DiffsolSolver``
+-----------------------------------
+
+Without ``output_variables``, both solvers store the full state at every
+output time; they differ in which times those are. ``IDAKLUSolver`` uses its
+internal integrator steps as the output grid when ``t_interp`` is omitted,
+whereas diffsol evaluates its error-controlled dense output at a grid of its
+own: ``t_eval`` and the times of any discontinuities, merged with 100 points
+uniform over the span. Pass ``t_interp`` to choose the output times exactly;
+they are then ``t_interp`` merged with ``t_eval`` and any discontinuity times.
+
+Off-grid reads — ``sol["Voltage [V]"](t)`` at a ``t`` between output points —
+interpolate with cubic Hermite, as they do for ``IDAKLUSolver``: the solver
+stores the state time derivatives alongside the solution by default. Pass
+``hermite_interpolation=False`` to drop them and halve trajectory memory, at
+the cost of off-grid reads falling back to linear interpolation. IDAKLU's
+Hermite knots are its internal steps, while diffsol's are its output times, so
+on a coarse ``t_interp`` IDAKLU's off-grid reads remain the more accurate of
+the two.
+
+Supported expression types
+--------------------------
+
+The Rust backend supports the following expression types:
+
+**Leaf nodes:**
+
+- ``Scalar`` — constant scalar values
+- ``Array`` / ``Vector`` / ``Matrix`` — dense arrays and matrices
+- ``SparseMatrix`` — CSR sparse matrices
+- ``StateVector`` / ``StateVectorDot`` — state variable slices
+- ``InputParameter`` — named input parameters
+- ``Time`` — simulation time
+
+**Binary operators:**
+
+- Arithmetic: ``+``, ``-``, ``*``, ``/``, ``**``
+- Matrix: ``@`` (matrix multiplication)
+- Comparison: ``minimum``, ``maximum``
+- Other: ``modulo``, ``hypot``, ``EqualHeaviside``, ``NotEqualHeaviside``
+
+**Unary operators:**
+
+- ``-`` (negation), ``abs``
+- Math functions: ``sqrt``, ``exp``, ``log``, ``sin``, ``cos``, ``tanh``,
+  ``sinh``, ``cosh``, ``arcsinh``, ``arctan``, ``erf``, ``sign``, ``floor``,
+  ``ceiling``
+- Reductions: ``max``, ``min`` (over arrays), differentiated to the
+  argmax/argmin subgradient
+
+**Structural:**
+
+- ``Index`` — array slicing
+- ``Concatenation`` — combining arrays
+- ``Conditional`` — branch selection
+- ``VectorField`` — stacked components (read scalar components back with
+  ``pybamm.Component``)
+
+**Interpolation:**
+
+- 1D ``linear``, ``cubic``, and ``pchip`` interpolation
+- 2D and 3D ``linear`` and ``cubic`` regular-grid interpolation
+
+A symbol the backend cannot convert raises an error naming the symbol at
+solver set-up, so an unsupported model fails loudly rather than
+mis-evaluating.
+
+Adding new expression types
+---------------------------
+
+A symbol reaches the Rust core through its ``_to_rust(graph, rust_symbols)``
+method, which allocates nodes on a ``pybamm.rust.ExprGraph``. Subclasses of
+:class:`pybamm.Function` override ``_rust_evaluate(graph, *children)`` instead,
+which receives the children already converted.
+
+1. **If existing graph methods can express it**, only the Python side
+   changes: compose them in ``_rust_evaluate`` or ``_to_rust``, as
+   ``Arcsinh2`` in ``functions.py`` does.
+
+2. **Otherwise, add a primitive** under ``packages/pybamm-rust``:
+
+   - add a ``Node`` variant in ``pybamm-core/src/node.rs``, the vocabulary of
+     the graph Python builds;
+   - lower it to the instruction tape in ``pybamm-core/src/ir/mod.rs``, adding
+     a ``UnaryOp`` or ``BinaryOp`` variant for an element-wise operation;
+   - define that op's value in ``pybamm-core/src/eval.rs``
+     (``with_unary_op!`` or ``with_binary_op!``), which every interpreter of
+     the tape shares;
+   - give it a derivative in ``tangent.rs`` (forward mode) and ``adjoint.rs``
+     (reverse mode), and an arm in each pass whose ``match`` the compiler then
+     reports as non-exhaustive (``simplify.rs``, ``sparsity.rs``,
+     ``zero_propagate.rs``, ``row_extract.rs``, ...);
+   - expose it as an ``ExprGraph`` method in ``pybamm-python/src/expr.rs``,
+     and declare that method in the hand-written stubs,
+     ``packages/pybamm/src/pybamm/rust/_core.pyi``, in the same change
+     (``tests/unit/test_rust_stubs.py`` checks them against the extension).
+
+3. **Test** the lowering in
+   ``packages/pybamm/tests/unit/test_expression_tree/test_operations/test_convert_to_rust.py``
+   and its values against CasADi in
+   ``packages/pybamm/tests/integration/test_rust_parity.py``. A new primitive
+   also joins the generators in ``pybamm-core/tests/common/cases.rs``, so the
+   Rust property tests exercise its values and derivatives.
@@ -5,7 +5,9 @@ PyBaMM's other backends. Python hands over a discretised model as an expression
 this workspace turns it into flat instruction tapes it can evaluate, differentiate and
 solve.
 
-Models opt in with `model.convert_to_format = "rust"`.
+Models use it by default; set `model.convert_to_format = "casadi"` to opt out. The
+[Rust backend user guide](../../docs/source/user_guide/rust_backend.rst) covers choosing
+the backend and how it differs from CasADi.
 
 ## Layout
 

@@ -483,7 +483,9 @@ def serialise_custom_model(model: pybamm.BaseModel, compress: bool = False) -> d
             "base_class": base_cls_str,
             "base_class_mro": base_class_mro,
             "options": getattr(model, "options", {}),
-            "convert_to_format": getattr(model, "convert_to_format", "casadi"),
+            "convert_to_format": getattr(
+                model, "convert_to_format", pybamm.BaseModel._DEFAULT_CONVERT_TO_FORMAT
+            ),
             "rhs": [
                 (
                     convert_symbol_to_json(variable),

@@ -95,7 +95,7 @@ class BaseModel:
         different submodels in the full concatenated solution vector.
     """
 
-    _DEFAULT_CONVERT_TO_FORMAT = "casadi"
+    _DEFAULT_CONVERT_TO_FORMAT = "rust"
     _VALID_CONVERT_TO_FORMATS = ("python", "casadi", "jax", "rust")
 
     # Set by discretisation, deserialisation and BaseSolver.set_up. Evaluators are CasADi
@@ -280,7 +280,7 @@ def convert_to_format(self) -> str | None:
             - "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
+        Default is "rust". ``None`` is deprecated: it evaluates as "python" and
         raises a ``DeprecationWarning`` when set.
         """
         return self._convert_to_format

@@ -601,10 +601,15 @@ def _check_and_prepare_model_inplace(self, model):
 
         if (
             isinstance(self, pybamm.CasadiSolver | pybamm.CasadiAlgebraicSolver)
-        ) and model.convert_to_format != "casadi":
-            pybamm.logger.warning(
-                f"Converting {model.name} to CasADi for solving with CasADi solver"
+            and model.convert_to_format != "casadi"
+        ):
+            # Converting a model left at the default format is routine
+            log = (
+                pybamm.logger.info
+                if model.convert_to_format == model._DEFAULT_CONVERT_TO_FORMAT
+                else pybamm.logger.warning
             )
+            log(f"Converting {model.name} to CasADi for solving with CasADi solver")
             model.convert_to_format = "casadi"
         if isinstance(self.root_method, pybamm.CasadiAlgebraicSolver):
             if model.convert_to_format == "rust":

@@ -186,8 +186,9 @@ def test_overpotentials(self):
         elif self.operating_condition == "off":
             np.testing.assert_allclose(self.eta_r_av(self.t), 0, rtol=1e-7, atol=1e-6)
             np.testing.assert_allclose(self.eta_e_av(self.t), 0, rtol=1e-12, atol=1e-11)
+            # Analytically zero; round-off of ~4 V potentials varies by backend.
             np.testing.assert_allclose(
-                self.delta_phi_s_av(self.t), 0, atol=2e-14, rtol=1e-16
+                self.delta_phi_s_av(self.t), 0, atol=5e-14, rtol=1e-16
             )
 
     def test_ocps(self):

@@ -15,6 +15,8 @@ def test_leading_order_convergence(self):
         # Create models
         leading_order_model = pybamm.lead_acid.LOQS()
         full_model = pybamm.lead_acid.Full()
+        for model in (leading_order_model, full_model):
+            model.convert_to_format = "casadi"
 
         # Same parameters, same geometry
         parameter_values = full_model.default_parameter_values

@@ -10,6 +10,8 @@ class TestCompareBasicModels:
     def test_compare_full(self):
         basic_full = pybamm.lead_acid.BasicFull()
         full = pybamm.lead_acid.Full()
+        for model in (basic_full, full):
+            model.convert_to_format = "casadi"
 
         parameter_values = pybamm.ParameterValues("Sulzer2019")
         parameter_values["Current function [A]"] = 10

@@ -19,6 +19,7 @@ def test_compare_averages_asymptotics(self):
         param = models[0].default_parameter_values
         param.update({"Current function [A]": 1})
         for model in models:
+            model.convert_to_format = "casadi"
             param.process_model(model)
 
         # set mesh
@@ -61,6 +62,7 @@ def test_compare_outputs_surface_form(self):
             param = models[0].default_parameter_values
             param.update({"Current function [A]": 1})
             for model in models:
+                model.convert_to_format = "casadi"
                 param.process_model(model)
 
             # set mesh

@@ -25,6 +25,7 @@ def test_casadi_solver_mode(self):
         for solver in solvers:
             # define model
             model = pybamm.lithium_ion.SPM()
+            model.convert_to_format = "casadi"
 
             # solve simulation
             sim = pybamm.Simulation(
@@ -57,6 +58,7 @@ def test_casadi_fast_matches_safe_default_ode(self):
         solutions = []
         for mode in ["safe", "fast"]:
             model = pybamm.lithium_ion.SPM()
+            model.convert_to_format = "casadi"
             solver = pybamm.CasadiSolver(mode=mode, atol=1e-6, rtol=1e-6)
             sim = pybamm.Simulation(model, solver=solver)
             solutions.append(sim.solve(t_eval))

@@ -323,6 +323,7 @@ class TestIDAKLUSolverAOTCompilation:
     @staticmethod
     def _build_model():
         model = pybamm.lithium_ion.SPMe()
+        model.convert_to_format = "casadi"
         geometry = model.default_geometry
         param = model.default_parameter_values
         param.process_model(model)

@@ -64,6 +64,7 @@ def test_default_voltage_is_expression(self, model_cls):
     @pytest.mark.parametrize("model_cls", REDUCED_MODELS)
     def test_opt_in_solvable_by_casadi_safe(self, model_cls):
         model = model_cls(options=OPT_IN)
+        model.convert_to_format = "casadi"
         sim = pybamm.Simulation(model, solver=pybamm.CasadiSolver(mode="safe"))
         sol = sim.solve([0, 3600])
         v = sol["Voltage [V]"].entries

@@ -176,6 +176,7 @@ def test_simulation_pickled_with_public_convert_to_format(self, monkeypatch):
         # in each model's instance dict under the public name.
         sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
         sim.build()
+        convert_to_format = sim.built_model.convert_to_format
 
         def public_name_state(self):
             state = self.__dict__.copy()
@@ -189,7 +190,7 @@ def public_name_state(self):
         monkeypatch.undo()
 
         loaded = pickle.loads(data)
-        assert loaded.built_model.convert_to_format == "casadi"
+        assert loaded.built_model.convert_to_format == convert_to_format
         solution = loaded.solve([0, 100])
         np.testing.assert_allclose(
             solution["Voltage [V]"](100.0),

@@ -1463,6 +1463,8 @@ def test_run_experiment_cccv_solvers(self):
         solutions = {}
         for name, solver in solvers.items():
             model = pybamm.lithium_ion.SPM()
+            if isinstance(solver, pybamm.CasadiSolver):
+                model.convert_to_format = "casadi"
             sim = pybamm.Simulation(model, experiment=experiment_2step, solver=solver)
             solution = sim.solve()
             assert solution.t[-1] == pytest.approx(3600 * len(experiment_2step.steps))

@@ -12,6 +12,8 @@ class TestQuickPlot:
     @pytest.mark.parametrize("solver", _solver_args)
     def test_simple_ode_model(self, solver):
         model = pybamm.lithium_ion.BaseModel(name="Simple ODE Model")
+        if isinstance(solver, pybamm.CasadiSolver):
+            model.convert_to_format = "casadi"
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         # Create variables: domain is explicitly empty since these variables are only