[BradyPlanden]

Description

Adds a DAE solver to pybamm-core: a solver module behind a diffsol feature, on by default. It integrates a CompiledModel in-process with diffsol's BDF method and faer's sparse LU, with forward sensitivities, termination events, output-variable trajectories and a rayon batch layer for input sets. It depends on diffsol 0.18.0 from crates.io, faer 0.24 and rayon. Python reaches it in #5965.

• PreparedSolver:
  • It is built once per model and holds everything a solve reuses: the compiled model, the Jacobian and mass sparsity, the mass matrix, and the symbolic LU of the algebraic block.
  • It holds no per-solve mutable state, so it is Send + Sync and one instance serves concurrent solves.
  • Each solve builds its own Equations (diffsol's OdeEquations). Equations owns a ModelEvaluator plus the solve's inputs, initial state and initial sensitivities, and the operator views borrow the evaluator's RefCell<Workspace> per call.
• Requests and outcomes:
  • A solve takes a SolveRequest (stop times, outputs, sensitivities and LU threads) and an InputSet (inputs, y0, and an optional dy0/dp seed for a continued step).
  • It returns a SolveOutcome with the times, states, yp and sensitivities (whose last column is the end state), the start and event states and their sensitivities, the index of the event that ended the solve, and SolverStatistics.
• Sensitivities:
  • They are integrated under error control, with sens_rtol = rtol and parameter scales |p|, or 1 where p is 0 or not finite. On differential states sens_atol = atol * sens_atol_factor; algebraic states keep the state atol. The default factor of 1.0 matches IDAS's IDASensEEtolerances.
  • If integration fails under that control, the solve is retried once from the same initial state with sensitivities excluded from error control, and sens_error_control_relaxed is set. An initial-condition failure is not retried.
• Initial conditions and stop times:
  • Initial conditions come from diffsol's bdf_state.
  • At a discontinuity, which PyBaMM brackets in t_eval with two times one ulp apart, the solver steps across to the far side of the bracket, re-solves y and y' with a line-searched Newton iteration and then the sensitivity system, and resets the step size, as IDAKLU's HandleBreakpoint does. Other stop times on the grid restart where the previous segment ended, and stop times off the grid are integrated through.
  • For a model with algebraic states whose rhs reads t, derivative_wrt_time corrects the stored yp at each segment start.
• Outputs and events:
  • Output variables are evaluated in windows of 128 points, with sensitivities projected as dg/dp + dg/dy·S.
  • A termination event that is already at or below zero at a segment start ends the solve there, naming the most negative event, with a consistent yp on that column and no repeat of the previous column's time.
• Batches and options:
  • solve_batch runs input sets on a rayon pool, with every sparse LU of each set, the algebraic block's included, on one thread to avoid oversubscription. A single solve's LUs take the request's thread count, never faer's global setting, so the solver also runs where a host disables faer's global parallelism. A failing set returns its own error, and results are bitwise equal to a serial loop.
  • SolverOptions exposes 15 of diffsol's options with diffsol's defaults, except max_nonlinear_solver_failures = 100_000.

Each solve builds its OdeSolverProblem afresh. Equations holds a RefCell, so it is not Sync, while solve_batch shares one &PreparedSolver across workers. The time this takes is reported per solve as problem_setup_time_secs.

Limits:

• Linear algebra: faer's sparse LU is the only linear solver, and df/dp is assembled from a dense pattern.
• Events: they terminate the solve; none resets the state.
• Stored yp: only for state trajectories of up to 4096 points, never for output variables.
• Algebraic states: diffsol picks them from the mass matrix diagonal and PyBaMM from zero columns, so they disagree on a zero diagonal entry beside an off-diagonal nonzero.

pybamm-python builds pybamm-core with its default features, so the wheel compiles diffsol, faer and rayon from this PR on, before #5965 uses them.

Review on #5732

• Equations duplicates ModelEvaluator: Equations owns a ModelEvaluator, the type the C ABI uses, and adds only the solve's own state.
• Rc where a RefCell would do: ModelEvaluator holds its workspace in a RefCell, which the operator views borrow through a shared reference. The C ABI holds the evaluator exclusively and uses get_mut, with no borrow-flag check.
• Use diffsol's UnitCallable: Reset and Out are UnitCallable.
• Resets through the unified experiment mode: events terminate the solve. Experiment steps carry gradients by seeding dy0/dp from the previous step (InputSet::with_sens_seed), which feat(rust): pybamm.DiffsolSolver, integrating rust models with diffsol #5965 tests against CasADi.
• Redundant docstring paragraph: the RhsOp docstring is one line.
• A new OdeSolverProblem per solve: this keeps PreparedSolver Send + Sync for solve_batch, as above. The cost is reported as problem_setup_time_secs; per solve it is 0.7 µs for the SPM, 1.2 µs for the SPMe and 4.5 µs for the DFN.

Type of change

New feature, under ## Features in CHANGELOG.md.

Key checklist

• No style issues: pre-commit run, including cargo fmt and cargo clippy
• Rust tests pass (578 with --all-features, 1 ignored), and the MSRV check passes (cargo +1.89 check --locked --all-targets --all-features)
• Tests added: 95 Rust tests.
  • Accuracy:
    • two-state analytic decay to 1e-6;
    • sensitivities against central finite differences.
  • Restarts and initial conditions:
    • restarts across algebraic jumps, with and without sensitivities and events;
    • consistent yp, including time-dependent constraints;
    • diagonal and general mass matrices.
  • Input validation: empty, unsorted or non-finite t_eval, lengths, and tolerances IDAS would reject.
  • Sensitivities:
    • the dy0/dp seed and sensitivity subsets;
    • the 128-point output window;
    • event and end-state sensitivities;
    • sens_atol and its factor;
    • the retry through solve(), which from the shared initial state is bitwise equal to a relaxed solve from scratch, and a retry that fails too (SensRetryFailed).
  • Batches and threading:
    • batches bitwise equal to the serial loop at 1, 2 and 8 threads, with parallel LU, sensitivities and a failing set;
    • PreparedSolver is Send + Sync;
    • batches and a multithreaded LU with faer's global parallelism disabled, in their own test binary.
  • Events: the termination names the event that ended the solve, and an event at a segment start adds no repeated column.
  • Benches:
    • solver_setup: preparing, reusing, and setting up afresh;
    • sens_params: 1, 4 and 16 of 16 sensitivity columns.

Rust execution engine stack: proposal #5667, showcase #5732

Each PR is one commit on the one above it.

#	PR	Adds
0	#5954	Function.create_copy keeps the function's class
0b	#5955	CasADi's RegPower and uniform-grid interpolant lookups match evaluate
1	#5956	pybamm-core: expression graph, typed IR, interpreters, forward-mode AD
2	#5957	pybamm-core: CompiledModel, coloured Jacobians, the C ABI
3	#5958	pybamm.rust bindings, built into the pybamm wheel
4	#5959	Symbol.to_rust
5	#5960	convert_to_format = "rust" for ScipySolver, AlgebraicSolver and EIS
6	#5961	pybammsolvers: IDAKLU over the C ABI
7	#5962	IDAKLUSolver on rust models
8	#5963	variables and sensitivities read in Rust
9	#5964	pybamm-core: diffsol DAE solver ← this PR
10	#5965	pybamm.DiffsolSolver
11	#5966	"rust" as the default convert_to_format

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
⚠️ Please upload report for BASE (stack/08-native-observation@52b69b1). Learn more about missing BASE report.

Additional details and impacted files

@@                      Coverage Diff                       @@
##             stack/08-native-observation    #5964   +/-   ##
==============================================================
  Coverage                               ?   98.20%           
==============================================================
  Files                                  ?      353           
  Lines                                  ?    37029           
  Branches                               ?        0           
==============================================================
  Hits                                   ?    36364           
  Misses                                 ?      665           
  Partials                               ?        0           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[github-actions]

Bencher Report

Project	PyBaMM
Branch	stack/09-rust-core-diffsol
Testbed	bare-metal

⚠️ WARNING: Truncated view!

The full continuous benchmarking report exceeds the maximum length allowed on this platform.

🐰 View full continuous benchmarking report in Bencher