
ID	0028
Class	BUG
Severity	2
Status	blocked

Models: —

Found: 2026-09-27, net-hydrodynamics programme (fhsim_marine_elements test S6), diagnosed against e0bbf4ea

Decision needed: Which default the Sundials backend should use when <LinearSolver> is absent and SPGMR would run unpreconditioned on a model that resolves to DenseLU: option (a), (b) or both (see "Possible fix")

Evidence

At e0bbf4ea:

src/engine/integrator/IntegratorConfig.h:201: SundialsSettings::linearSolver defaults to SundialsLinearSolver::Spgmr, so <Integrator Method="BDF"> without a <LinearSolver> element runs SPGMR.
src/engine/integrator/sundials/SundialsSolverBuilder.cpp:200-205: SUNLinSol_SPGMR(stateVector, precondType, 0, ctx). maxl = 0 selects SUNSPGMR_MAXL_DEFAULT, which is 5 (include/sunlinsol/sunlinsol_spgmr.h of the SUNDIALS packages in the Conan cache). No XML option sets maxl: SundialsOptionApplier.cpp passes CVodeSetEpsLin through but nothing reaches SUNLinSol_SPGMRSetMaxl.
SetupSpgmrPreconditioner (SundialsSolverBuilder.cpp:80-170) builds a structural preconditioner only for BlockTridiagonal, Component, NearTridiagonal and BorderedSchur. It builds the Jacobi diagonal only for those types plus Iterative and IterativeGMRES, or when Preconditioner="jacobi" is given. For a model whose ResolveLinearSolverType() is DenseLU, wantsDiagonal is false, and :136-137 returns nullptr without logging anything (the other unpreconditioned paths do log). precondType is then SUN_PREC_NONE.
The summary prints Linear solver : DenseLU and Preconditioner : none (src/engine/integrator/sundials/FhIntegratorSundials.cpp:322-353). This hides that SPGMR is running. It is filed separately as FHSIM-0029.

Mechanism, traced in gdb on the S6 model below (SUNDIALS 7.5, cvLsSolve): on the second Newton iteration of a step, the unpreconditioned GMRES solve with 5 Krylov vectors does not reach its tolerance and returns SUNLS_RES_REDUCED. Because curiter >= 1, cvLsSolve turns this into SUN_NLS_CONV_RECVR (902), a nonlinear convergence failure, and CVODE cuts the step by 0.25. The instrumented run shows the pattern repeating, e.g. around t = 1.10 s:

NFLAG tn=1.095582959 nflag=902   (conv wrms 1.03 > tol 0.2)
NFLAG tn=1.101832959 nflag=0
NFLAG tn=1.108082959 nflag=902
NFLAG tn=1.113791130 nflag=902
NFLAG tn=1.121860196 nflag=902

Reproduction and numbers

The model is fhsim_marine_elements test S6 (tests/NetValidation_ZhanCylinder_Test.cpp): one NetStructure, 720 states, dense analytical Jacobian, 192 panels casting a gridded wake that is rebuilt and blended in during the run. It runs with <Integrator Method="BDF"> and no <LinearSolver>. Log: Method : cvodes/BDF, Jacobian type : dense, Linear solver : DenseLU, Preconditioner : none, hmax = 0.1.

Configuration	Steps	Wall	Step size in the transients
default (SPGMR, no preconditioner, maxl 5)	761	11 s	falls from 0.1 s to 2e-4 s (first blend) and to 3.5e-5 s (second)
default + `Preconditioner="jacobi"`	—	—	still collapses, to 1.35e-4 s
`<LinearSolver Type="DENSE"/><Jacobian Type="dense"/>`	379	1 s	steady at 0.1 s

gdb counted 77 and 36 SUN_NLS_CONV_RECVR failures inside the two blends and none elsewhere. Wakes are not the cause. A load ramp with the wake switched off collapses the step the same way, to 9.7e-5 s. Any smooth transient that makes the first Newton iterate insufficient triggers it.

Two observations were not diagnosed. They are recorded so that option (a) is not taken as a complete cure. With DENSE, a stiffer variant (stiffness k = 1e2) aborts at t = 2.59 s with repeated error-test failures. With DENSE, a 0.1 s blend still dips to 5e-4 … 1.5e-3 s.

Diagnostic artifacts: fhsim_marine_elements scratch branch diag/bdf-step (cb3af26, never merged; instrumentation in NetStructure.cpp).

Effect

A user who writes <Integrator Method="BDF"> for a dense model gets a solver whose Newton iteration fails in smooth transients. The step collapses by two to three orders of magnitude, and the run is about ten times slower than with a direct solver. Nothing in the log says why: it names DenseLU and no preconditioner.

Possible fix

(a) When <LinearSolver> is absent and SPGMR would run unpreconditioned on a model that resolves to DenseLU (i.e. SetupSpgmrPreconditioner would return nullptr), default to the SUNDIALS DENSE direct solver with a dense Jacobian matrix. This changes the step sequence and timing of every existing BDF/Adams run on such a model, so baselines that rely on the default would move.
(b) Give SPGMR a real preconditioner for DenseLU models, e.g. a dense LU factorisation of I - gamma*J, and/or expose maxl (<LinearSolver Maxl="..."> or a SUNLinSol_SPGMRSetMaxl passthrough). The Jacobi diagonal is not enough here (still 1.35e-4 s).
(c) Log the solver actually in use, and warn when SPGMR runs unpreconditioned by default (FHSIM-0029).

Test that would prove it

A Sundials test on a small dense model with a smooth transient (e.g. a mass-spring chain under a C¹ load ramp), Method="BDF", no <LinearSolver>. Assert that the minimum step during the ramp stays within a small factor of the step before it, and that there are few or no nonlinear convergence failures. It fails today and should pass with (a) or an adequate (b).

Risk

(a) changes the default solver, and so the step sequence, of every existing Sundials run on a DenseLU model that names no linear solver. Regression baselines that use it would need regenerating. (b) keeps results closer to today's for models that already converge.