
ID	0036
Class	BUG
Severity	2
Status	ready

Models: —

Found: 2026-10-02, R135 implicit-integrator review of the net-hydrodynamics programme (cage24 case), at e782fc44. Related: FHSIM-0028, MARE-0122.

Evidence

The model is a 24 × 6 cylinder cage Net/NetStructure (1008 states, CastWake="1") with <Integrator Method="BDF"> and no <LinearSolver>. The run aborts before the first step. The same failure occurred with the binaries built before and after the R135 hold change (B2 and B3), so it is not a consequence of the wake work.

Log, at e782fc44 (SUNDIALS 7.5):

[Info]    Jacobian sparsity: 1008 states, NNZ=160776, density=15.823413%, bandwidth=0
[Info]    Linear solver auto-selected: NearTridiagonal (1008 states, NNZ=160776, density=15.823413%,
          bandwidth=0, extra_blocks=0, correction_rank=0)
[Info]    Preconditioner "NearTridiagonal" set up with 1008 states for CVODE SPGMR
[Info]      Method          : cvodes/BDF
[Info]      Linear solver   : NearTridiagonal
[Info]      Preconditioner  : NearTridiagonal
[ERROR][rank 0][.../cvodes/cvodes.c:8110][cvHandleFailure] At t = 0, the setup routine failed in an unrecoverable manner.
SUNDIALS_ERROR: CVode() failed with flag = -6

CVodeGetIntegratorStats at the failure: steps 0, NLS step fails 1, RHS evals 4, initial step 1.5e-8 s. Flag −6 is CV_LSETUP_FAIL, so the SPGMR preconditioner setup callback returned an unrecoverable error on its first call. The error is not reported by FhSim: the log names no cause.

Code path:

src/engine/model/JacobianSparsityAnalyzer.cpp:487-503: SelectLinearSolver returns NearTridiagonal when the density is below 0.5, exactly one SimObject has states (numStateful == 1) and NearTridiagonalSolver::DetectStructure reports detected. Here it reports bandwidth=0, extra_blocks=0, correction_rank=0 for a pattern with 15.8 % density (160 776 non-zeros in 1008 × 1008). That is not the structure the solver is meant for.
src/engine/integrator/sundials/SundialsSolverBuilder.cpp:107-125: for SPGMR, the resolved type NearTridiagonal is used as a preconditioner through CreateLinearSolver, and the log line Preconditioner "NearTridiagonal" set up is written once IsReady() is true. Whether the numeric factorisation at the first gamma succeeds is only known inside the setup callback.

The cause inside NearTridiagonalSolver (detection or factorisation) is not diagnosed here.

Reproduction

At e782fc44, with marine_elements and environment libraries built from the net-hydrodynamics branches (the failure is in the FhSim core, not the wake). The script below writes cage_net.xml and cage_in.xml. It is gen_cage.py --ncirc 24 --nrows 6 --castwake 1 --method BDF --tend 10 --step 0.1 reduced to what the case needs.

python3 gen.py

FhSim cage_in.xml -l log.txt -c 0 # exits 1 with the error above

gen.py:

import math
nc, nr, R, depth, U = 24, 6, 8.0, 12.0, 0.5
bar_l, bar_d = 0.025, 0.003
cw, ch = 2 * math.pi * R / nc, depth / nr
mw, mh = cw / bar_l, ch / bar_l
node = lambda i, j: j * nc + (i % nc)
nn = nc * (nr + 1)
L = ['<Net>\n  <General MaxNumNodes="2147483647" MaxNumPanels="0" MaxNumCables="2147483647"'
     ' MaxElementLength="1.798E+308" MinElementLength="2.225E-308" MaxElementLengthFactor="1.798E+308"'
     ' MaxStepLength="0.1" AdaptationPeriod="0" StepSafetyFactor="50" MaxTension="1E+06" MeanTension="0"'
     ' MinLength="0.1" MassFactor="1" MaxVelocity="0" TwineColour="0.1, 0.3, 0.2, 0.5"/>',
     '  <ExternalNodeMap ' + ' '.join(f'T{i}="{node(i, 0)}" B{i}="{node(i, nr)}"' for i in range(nc)) + '/>',
     '  <NetElements>']
k = 0
for j in range(nr):
    for i in range(nc):
        a, b, c, d = node(i, j), node(i + 1, j), node(i, j + 1), node(i + 1, j + 1)
        u0, u1, v0, v1 = i * mw, (i + 1) * mw, j * mh, (j + 1) * mh
        for conn, m in (((a, b, c), ((u0, v0), (u1, v0), (u0, v1))), ((b, d, c), ((u1, v0), (u1, v1), (u0, v1)))):
            k += 1
            L.append(f'    <El{k} PosAMesh="{m[0][0]:.6f}, {m[0][1]:.6f}" PosBMesh="{m[1][0]:.6f}, {m[1][1]:.6f}"'
                     f' PosCMesh="{m[2][0]:.6f}, {m[2][1]:.6f}" BarD="{bar_d}" BarL="{bar_l}" BarE="2e8" BarRho="1100"'
                     f' RefineFactor="1" TwineColour="0.2, 0.4, 0.2, 0.5" Conn="{conn[0]},{conn[1]},{conn[2]}"/>')
L.append('  </NetElements>\n  <CableElements>')
k = 0
for j in (0, nr):
    for i in range(nc):
        k += 1
        L.append(f'    <El{k} Conn="{node(i, j)},{node(i + 1, j)}" D="{0.03 if j == 0 else 0.02}" L="{cw:.6f}" E="1e10"'
                 ' Rho="2000" MaxTension="1e6" NumSpheres="0" SphereD="" SphereMass="" SpherePos="" SphereMesh=""'
                 ' AddedWeightInWater="0" NumDisks="0" DiskD="" DiskThickness="" DiskMass="" DiskPos="" DiskMesh=""'
                 ' RefineFactor="1"/>')
L.append('  </CableElements>\n</Net>\n')
open('cage_net.xml', 'w').write('\n'.join(L))
fx = -0.5 * 1025 * 0.3 * (2 * R * depth) * U ** 2 * 2 / (2 * nc)
names = ', '.join([f'T{i}' for i in range(nc)] + [f'B{i}' for i in range(nc)])
conn = ' '.join([f'Net0.T{i}Force="{fx:.3f},0,-40"' for i in range(nc)] + [f'Net0.B{i}Force="{fx:.3f},0,40"' for i in range(nc)])
ini = []
for j in range(nr + 1):
    for i in range(nc):
        th = 2 * math.pi * i / nc
        n = node(i, j) + 1
        ini.append(f'Net0.Pos_{n}="{R * math.cos(th):.6f}, {R * math.sin(th):.6f}, {1.0 + j * ch:.6f}" Net0.Vel_{n}="0,0,0"')
open('cage_in.xml', 'w').write(f'''<Contents>
   <OBJECTS>
      <Lib LibName="environment" SimObject="Environment" Name="Env" Waves.Model="Component" Waves.NumWaves="0"
           Bathymetry.Depth="200" Current.Model="Constant" Current.Vector="{U},0,0"/>
      <Lib LibName="marine_elements" SimObject="Net/NetStructure" Name="Net0" File="cage_net.xml"
           NodesInputForce="{names}" NodesOutputPosAndVel="T0" Scale="1" CastWake="1"/>
   </OBJECTS>
   <INTERCONNECTIONS>
      <Connection {conn}/>
   </INTERCONNECTIONS>
   <INITIALIZATION>
      <InitialCondition {' '.join(ini)}/>
   </INITIALIZATION>
   <SIMULATION>
      <Timing TStart="0" TEnd="10.0"/>
      <Integrator Method="BDF">
         <StepControl StepMax="0.1" AbsTol="1e-3" RelTol="1e-3"/>
      </Integrator>
   </SIMULATION>
   <OBSERVERS>
      <FileOutput Select="object.Net0:port.T0Pos"/>
   </OBSERVERS>
</Contents>
''')

Results, one thread, StepMax="0.1", AbsTol=RelTol=1e-3, 10 s:

Configuration	Result
default (no `<LinearSolver>`: SPGMR + auto NearTridiagonal preconditioner)	fails at t = 0, flag −6
`<LinearSolver Type="DENSE"/>`	runs to 10 s (4227 steps)
`<LinearSolver Type="SPGMR" Preconditioner="none"/>`	runs to 10 s (4279 steps)
`<LinearSolver Type="SPGMR" Preconditioner="jacobi"/>`	runs to 10 s
DIRK (ARKode), DENSE	runs to 10 s

Step counts are from the review run with ExternalWakeUpdate at its default. Run times under B3: DENSE 7.0 s, SPGMR none 2.3 s.

Effect

A user who writes <Integrator Method="BDF"> for a net with about 1000 states and no <LinearSolver> gets no run at all. Both cheaper and slower workarounds exist, but nothing in the log points to them. The structural auto-selection also makes the result worse than the unpreconditioned default of FHSIM-0028 and MARE-0122, which at least runs. FHSIM-0028 option (a) (DENSE by default for DenseLU models) would not reach this case, because here the resolved solver is not DenseLU.

Possible fix

(a) Diagnose why DetectStructure accepts this 15.8 %-dense pattern with no bandwidth and no extra blocks, and why the setup fails. Either tighten the detection or fix the factorisation.
(b) If a structural preconditioner setup fails at the first call, fall back to the Jacobi diagonal (as SetupSpgmrPreconditioner already does when prec->IsReady() is false), and log which preconditioner failed and why, instead of aborting with a bare CV_LSETUP_FAIL.
(c) Decide with FHSIM-0028 what the default for a net model should be (DENSE, or SPGMR with a working preconditioner).

Test that would prove it

A Sundials test on a model whose sparsity pattern is dense enough to reach the NearTridiagonal branch with bandwidth=0, extra_blocks=0 (the pattern from the reproduction above, or a reduced one with the same signature, found while diagnosing), Method="BDF", no <LinearSolver>. Assert that the run completes, and that if a preconditioner is rejected, the log names it. It fails today with flag −6 at t = 0.

Risk

Changing the detection can change the auto-selected solver, and therefore the step sequence, of existing models that run with NearTridiagonal today. (b) alone only affects the runs that fail today.