IDAKLU's full-system Newton initialisation diverges, and the failure reads as IDA_BAD_K #6000

Description

IDAKLU's full-system Newton for the initial conditions solves with IDA's Jacobian at cj = 1 (AlgebraicICBuilder.inl#L32-L34), then zeroes the differential part of the step (NonlinearSolver.inl#L94). With the differential states held fixed, the algebraic step should solve with ∂g/∂y_a. The coupled solve folds ∂g/∂y_d (∂f/∂y_d − M)⁻¹ ∂f/∂y_a into that matrix, and on the DFN the Newton never converges.

IDACalcIC usually recovers, but its flag is ignored (IDAKLUSolverOpenMP.inl#L897). When it also fails, the first step fails, and GetSolutionFull(t0) runs before that step's flag is checked (#L385-L389), so the error reads IDAGetDky: IDA_BAD_K.

Evidence

DFN (Chen2020, default mesh and tolerances), sim.solve([0, 3600]), pybamm 26.10.1.1 / pybammsolvers 0.11.1 on arm64:

newton_mode="auto" converges in 5–8 iterations at every current. With #5966's solver log, IDACalcIC reports IDA_LINESEARCH_FAIL at 10 A. From the same guess, the step at cj = 1 diverges, while at cj = 1e6 it matches the algebraic-block Newton, which converges. The Rust backend in #5966 fails the same way.

Motivation

Non-standard-form DAEs take this path even with newton_mode="auto", so they rely on IDACalcIC alone, and when that fails the error points at the wrong call. The robustness sweep for #5966 hit this on the default DFN at 2C.

Possible implementation

1. Keep ∂f/∂y_d and ∂g/∂y_d out of the full-system matrix, either by dropping them or by using a large cj as IDACalcIC does.
2. Evaluate that Jacobian at t_val; it's currently evaluated at t = 0.
3. Check IDACalcIC's flag, and check the first step's flag before GetSolutionFull(t0).
4. Add a test for the default DFN at 10 A with newton_mode="full", and a changelog entry.

Done when: the full-system Newton converges on the default DFN at 1–10 A, and an initialisation failure raises its own error.

import pybamm

parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.update({"Current function [A]": "[input]"})
for mode in ("auto", "full"):
    solver = pybamm.IDAKLUSolver(options={"newton_mode": mode})
    sim = pybamm.Simulation(
        pybamm.lithium_ion.DFN(), parameter_values=parameter_values, solver=solver
    )
    for current in (1.0, 5.0, 7.5, 10.0):
        try:
            sim.solve([0, 3600], inputs={"Current function [A]": current})
            outcome = "ok"
        except pybamm.SolverError as error:
            outcome = str(error)
        print(f"{mode} {current:4} A: {outcome}")