fix(statics): a restart at full stiffness runs on while its residual falls

CMakeLists.txt

@@ -1935,6 +1935,17 @@ if(Python3_Interpreter_FOUND)
     TIMEOUT 900
     ENVIRONMENT "${_cabledyn_python_test_environment}")
 
+  # The finite-EI statics of a clamped lazy wave do not depend on its azimuth: its restart at full
+  # stiffness runs while the residual keeps falling instead of to a fixed count.
+  add_test(NAME hermite_static_azimuth
+           COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_hermite_static_azimuth.py
+                   --exe $<TARGET_FILE:cabledyn>
+                   --output ${CMAKE_BINARY_DIR}/hermite_static_azimuth)
+  set_tests_properties(hermite_static_azimuth PROPERTIES
+    LABELS "python;hermite;static;robustness"
+    TIMEOUT 900
+    ENVIRONMENT "${_cabledyn_python_test_environment}")
+
   # The note on grounded lines that a current moves over a seabed without friction: once per run
   # on every standalone route, and never without a grounded node, a current, or a march, or
   # with friction.

doc/theory.rst

@@ -1834,7 +1834,12 @@ mesh-sequenced route described below first and this route as its fallback):
    take over, because the seabed contact and waterline laws are only piecewise smooth (and the
    curvature of a structured bathymetry grid is not differentiated). Stages that jump
    further than a small fraction of the line are retried at a smaller step; a first stage that
-   does not converge is restarted at full ``EI``.
+   does not converge is restarted at full ``EI``. Such a descent from the seed (and the one past
+   a limit point of the continuation) can run for more than a thousand iterations while the
+   energy falls and the residual rises and falls, and that count scatters with the last digits
+   of the input, such as the line's azimuth. Past its allowance it therefore runs on while the
+   two-norm of the scaled residual keeps halving within an allowance, to one and a half
+   allowances (a longer descent can settle on a folded branch that the other routes avoid).
 3. **Polish and diagnosis** on the deck mesh at the static tolerance. A residual floor above it
    is accepted at :math:`10^{-6}` (or, for a line on the seabed, :math:`10^{-5}`, nodes
    alternating across the 1 µm contact blend) and reported. An element folded beyond

doc/troubleshooting.rst

@@ -568,6 +568,18 @@ but a few geometries are genuinely hard:
   sparse long section is not over-coarsened merely because the deck also contains many short
   elements. Flat-bed installed lines with an already supported anchor retain their grounded-tail
   branch homotopy and are still accepted only after exact-mesh residual and kink checks.
+- **A long, light finite-EI line that solves at one orientation and not at another.** Turning a
+  deck about the vertical (its end points, current and seabed with it) leaves the physics
+  unchanged, but not the arithmetic: the static solve of a long line from its seed is a long
+  descent whose path the last digits of the input can change. Azimuths, declinations and the
+  global frame are defined in :doc:`conventions`; check that the deck's are the intended ones
+  before anything else. A descent from the seed runs on while its residual keeps falling and
+  stops when it has stalled (:doc:`theory`, Static equilibrium), so such a line normally solves
+  at every orientation; one that still fails, or that reaches a different shape at a nearby
+  orientation, is near a branch boundary of its own geometry. Inspect the static profile
+  (``.static.out``: positions, tension and curvature along the line), compare it with the deck at
+  a refined mesh, and judge whether the shape is the one the line would take. CableDyn does not
+  choose between physically admissible branches for you.
 - **A very slack cable on a frictionless flat bed.** If its rest length exceeds the vertical-plus-
   horizontal L-shaped path between the endpoints, a smooth single-touchdown equilibrium may not
   exist: the surplus must fold or form another contact region. CableDyn rejects a resulting

src/CableDyn_HermiteCableStatic.f90

@@ -309,7 +309,7 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_
                                           waterline_z, dry_buoyancy, merit_line_search, pseudo_transient, &
                                           energy_minimization, stable, water_weight, current, residual_out, &
                                           tangent_out, friction, residual_unmet, torsion, energy_out, &
-                                          end_share_hold, f_nodal_twist)
+                                          end_share_hold, f_nodal_twist, stall_limit)
     !! Solve the finite-EI static equilibrium of a Hermite bending-cable line.
     !!
     !! l0(ne), EA(ne), EI(ne), w(ne) : per-element rest length, axial/bending stiffness,
@@ -451,6 +451,15 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_
     !!                                  attachments in a current, evaluated at a twisted state)
     !!                                  then leaves the untwisted state in equilibrium under
     !!                                  f_nodal alone.
+    !! stall_limit (optional, >= 1)    : max_iter becomes the cap of a long descent, which
+    !!                                  stops (not converged) once stall_limit updates have
+    !!                                  not brought the two-norm of the scaled residual below
+    !!                                  half its value at the last such fall (the first: the
+    !!                                  start), so after stall_limit updates without any. A
+    !!                                  descent that is still converging is then not cut at a
+    !!                                  fixed count; the norm, unlike the largest component
+    !!                                  that decides convergence, does not change when the
+    !!                                  line is turned.
     REAL(wp), INTENT(IN) :: l0(:), EA(:), EI(:), w(:), seed(:)
     INTEGER, INTENT(IN) :: fixed_dofs(:)
     REAL(wp), INTENT(IN) :: seabed_z, kn, tol, damping
@@ -489,8 +498,13 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_
     REAL(wp), INTENT(IN), OPTIONAL :: end_share_hold(2)
     TYPE(CD_HermiteTorsionType), INTENT(INOUT), OPTIONAL :: torsion
     REAL(wp), INTENT(IN), OPTIONAL :: f_nodal_twist(:)
+    INTEGER, INTENT(IN), OPTIONAL :: stall_limit
 
     INTEGER :: ne, nn, ndof, c, it, es, i, b, nci, n_buoy, omp_threads, axial_order, bending_order
+    ! stall_limit: the residual merit at the start of an update, its value at the last halving
+    ! of the norm, and the number of updates that had reached that state.
+    REAL(wp) :: merit_update, merit_fall
+    INTEGER :: it_fall
     REAL(wp) :: eifac, rnorm, fscale, bfac, rmerit, ptc_tau, ptc_rprev, pi_total, pi_abs, lm_mu, chord
     REAL(wp), ALLOCATABLE :: q(:), trib(:), w_cur(:), q_base(:)
     REAL(wp) :: endconn_k(2), endconn_d0(3, 2), direction_norm
@@ -615,6 +629,11 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_
         .NOT. CD_Is_Finite(damping) .OR. damping <= CD_ZERO .OR. damping > CD_ONE) THEN
       CALL fail('need n_cont>=1, max_iter>=1, finite tol>0, finite 0<damping<=1'); RETURN
     END IF
+    IF (PRESENT(stall_limit)) THEN
+      IF (stall_limit < 1) THEN
+        CALL fail('need stall_limit >= 1'); RETURN
+      END IF
+    END IF
     IF (PRESENT(residual_history)) THEN
       IF (SIZE(residual_history) < max_iter) THEN
         CALL fail('residual_history must have at least max_iter entries'); RETURN
@@ -1034,10 +1053,20 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_
         drag_now = use_current .AND. b == n_buoy .AND. c == nci
         IF (PRESENT(history_count)) history_count = 0
         IF (PRESENT(residual_history)) residual_history = CD_ZERO
+        it_fall = 0
         DO it = 1, max_iter
           ! newton_step measures the residual at the CURRENT q and updates only when
           ! not yet converged, so it never over-steps a converged state.
           CALL newton_step(rnorm, es, ErrMsg, .TRUE.)
+          IF (PRESENT(stall_limit) .AND. rnorm >= tol .AND. es == CD_HCSTAT_OK) THEN
+            ! The norm has halved when its square (the merit) falls below a quarter.
+            IF (it == 1) THEN
+              merit_fall = merit_update
+            ELSE IF (merit_update < 0.25_wp*merit_fall) THEN
+              merit_fall = merit_update
+              it_fall = it - 1
+            END IF
+          END IF
           ! A stall is sticky by design: every later stage then runs a single step, so the
           ! stage loops still end at the full load and full EI, where the final evaluation
           ! gives the verdict (a stalled solve ends NOCONVERGE).
@@ -1053,6 +1082,9 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_
           IF (PRESENT(residual_history)) residual_history(it) = rnorm
           IF (PRESENT(history_count)) history_count = it
           IF (rnorm < tol) EXIT
+          IF (PRESENT(stall_limit)) THEN
+            IF (it - it_fall >= stall_limit) EXIT
+          END IF
         END DO
       END DO
       q_sec2 = q_sec1
@@ -1450,6 +1482,7 @@ RECURSIVE SUBROUTINE newton_step(rnorm, es, em, do_update)
       IF (.NOT. CD_Is_Finite(rnorm)) THEN
         es = CD_HCSTAT_NOCONVERGE; em = 'non-finite residual'; RETURN
       END IF
+      IF (do_update) merit_update = rmerit
       ! Pure residual evaluation, or already at tolerance: do not step.
       IF (.NOT. do_update .OR. rnorm < tol) RETURN
 
@@ -3532,6 +3565,11 @@ SUBROUTINE CD_HermiteCable_Static_Solve_Continuation(l0, EA, EI, w, seed, fixed_
     !! full-stiffness descent from the last converged state is tried (a limit point where
     !! the equilibrium snaps to another stable configuration); the continuation then stops
     !! with NOCONVERGE, as it does once iteration_budget Newton iterations are spent.
+    !! A full-stiffness descent (that one, or the restart from the seed of a first stage
+    !! that does not converge) has RESTART_ITER_FACTOR stage allowances; its energy pass runs
+    !! on to one and a half of them, within the budget, while its residual keeps falling, and
+    !! stops after that allowance without the norm of the residual halving (stall_limit of
+    !! CD_HermiteCable_Static_Solve).
     !! Intermediate stages use a looser tolerance (MAX(tol, 1e-6)); the final stage uses tol.
     !! The budget makes the failure path bounded: a caller's fallback route receives control
     !! after at most that many Newton iterations.
@@ -3565,7 +3603,8 @@ SUBROUTINE CD_HermiteCable_Static_Solve_Continuation(l0, EA, EI, w, seed, fixed_
     REAL(wp), PARAMETER :: STAGE_MOVE_TOL = 0.02_wp
     ! Every stage (energy solve and pseudo-transient retry) gets this multiple of max_iter.
     INTEGER, PARAMETER :: PTC_ITER_FACTOR = 4
-    ! A full-stiffness restart from the seed descends from far away; its stage allowance.
+    ! A full-stiffness restart from the seed descends from far away: its stage allowance, of
+    ! the iterations without progress in the energy pass.
     INTEGER, PARAMETER :: RESTART_ITER_FACTOR = 10
     INTEGER :: ne, nn, it, es, n_accepted
     REAL(wp) :: lambda, lambda_try, dlog, stage_tol, res
@@ -4324,9 +4363,19 @@ SUBROUTINE stage_solve(ei_x, seed_x, q_x, curv_x, res_x, it_x, es_x, em_x, ptc,
       CHARACTER(*), INTENT(OUT) :: em_x
       LOGICAL, INTENT(IN) :: ptc
       REAL(wp), INTENT(IN), OPTIONAL :: f_frozen(:)
-      INTEGER :: iter_x
+      INTEGER :: iter_x, stall_x
       iter_x = PTC_ITER_FACTOR*max_iter
       IF (restart_stage) iter_x = MAX(iter_x, MIN(iteration_budget - iters_out, RESTART_ITER_FACTOR*iter_x))
+      ! The energy pass of a restart descends from far away (on a long line, through more than
+      ! a thousand iterations in which the residual rises and falls while the energy falls):
+      ! past the restart allowance it runs on while the residual keeps falling, and stops after
+      ! an allowance without progress. A fixed count sat within the scatter of the length of
+      ! such a descent, so the last digits of the input (as from turning the line about the
+      ! vertical) decided whether the statics solved. The extension ends at one and a half
+      ! allowances: a descent still going there can settle on a folded branch that a
+      ! different start avoids (the PTC retry and the caller's other routes follow).
+      stall_x = iter_x
+      IF (restart_stage .AND. .NOT. ptc) iter_x = MAX(iter_x, MIN(iteration_budget - iters_out, (3*iter_x)/2))
       IF (PRESENT(endconn_stiffness)) THEN
         CALL CD_HermiteCable_Static_Solve(l0, EA, ei_x, w, seed_x, fixed_dofs, seabed_z, kn, 1, iter_x, &
                                           stage_tol, CD_ONE, q_x, curv_x, res_x, it_x, es_x, em_x, &
@@ -4340,7 +4389,7 @@ SUBROUTINE stage_solve(ei_x, seed_x, q_x, curv_x, res_x, it_x, es_x, em_x, ptc,
                                           endconn_stiffness=k_stage, endconn_direction=endconn_direction, &
                                           endconn_mode=endconn_mode, waterline_z=waterline_z, &
                                           dry_buoyancy=dry_buoyancy, water_weight=water_weight, current=cur_stage, &
-                                          friction=fric_stage)
+                                          friction=fric_stage, stall_limit=stall_x)
       ELSE
         CALL CD_HermiteCable_Static_Solve(l0, EA, ei_x, w, seed_x, fixed_dofs, seabed_z, kn, 1, iter_x, &
                                           stage_tol, CD_ONE, q_x, curv_x, res_x, it_x, es_x, em_x, &
@@ -4352,7 +4401,8 @@ SUBROUTINE stage_solve(ei_x, seed_x, q_x, curv_x, res_x, it_x, es_x, em_x, ptc,
                                           axial_quadrature_order=axial_quadrature_order, &
                                           bending_quadrature_order=bending_quadrature_order, &
                                           waterline_z=waterline_z, dry_buoyancy=dry_buoyancy, &
-                                          water_weight=water_weight, current=cur_stage, friction=fric_stage)
+                                          water_weight=water_weight, current=cur_stage, friction=fric_stage, &
+                                          stall_limit=stall_x)
       END IF
     END SUBROUTINE stage_solve
 

tests/test_hermite_static_azimuth.py

@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: Apache-2.0
+# Copyright (c) 2026 Jae Hoon Seo, SMI Lab, Inha University
+"""The finite-EI statics of a line do not depend on its azimuth.
+
+A 3.8 km suspended three-section lazy wave on 540 elements, clamped (Rigid) at its fairlead
+along a direction 11 deg out of the vertical plane of its chord, is turned about the vertical
+axis. Its first continuation stage fails on the coarse mesh at the clamp and the statics restart
+at full stiffness from the catenary seed, an energy descent whose length scatters with the last
+digits of the input (1171 to 1661 iterations over 28 azimuths). It was cut at 1600, so the
+statics failed at 328.29 and 329 deg and solved elsewhere; it now runs while its residual keeps
+falling. Every azimuth (seven across 328 to 329.5 deg and four elsewhere) must solve, with the
+same fairlead and anchor tensions to 1e-6.
+
+usage: test_hermite_static_azimuth.py --exe <cabledyn> --output <dir>
+"""
+
+from __future__ import annotations
+
+import argparse
+import math
+import os
+import subprocess
+import sys
+from pathlib import Path
+
+DECK = """--------------------- CableDyn Input File ------------------------------------
+Suspended lazy wave clamped at its fairlead out of its chord plane, turned about the vertical
+--------------------- LINE TYPES ---------------------------------------
+TypeName Diam Mass EA BA EI Cd_n Cd_t Ca_n Ca_t
+(-) (m) (kg/m) (N) (N-s) (N-m2) (-) (-) (-) (-)
+t1 0.13901651567091225 25.009060114295032 644913120.4893594 -1.0 203482.2288474299 1.2 0.01 1.0 0.0
+t2 0.6428183465278386 200.38465492542826 1054908164.887091 -0.8 38575.24712022166 1.2 0.01 1.0 0.0
+t3 0.07107464816032726 25.258028055571554 454644505.00028795 -1.0 2138.24079651583 1.2 0.008 1.0 0.0
+--------------------- POINTS -------------------------------------------
+ID Type X Y Z
+(-) (-) (m) (m) (m)
+1 Coupled {xa!r} {ya!r} -63.39980772163334
+2 Fixed {xb!r} {yb!r} -754.4713708188342
+--------------------- LINES --------------------------------------------
+ID NodeA NodeB Outputs
+(-) (-) (-) (-)
+1 1 2 -
+--------------------- SECTIONS -----------------------------------------
+LineID LineType Length NumSegs
+(-) (-) (m) (-)
+1 t1 1351.43480777781 200
+1 t2 792.1830512136142 140
+1 t3 1678.538077522514 200
+--------------------- END CONNECTIONS -------------------------------
+LineID End Stiffness EzX EzY EzZ
+(-) (-) (N-m/rad) (-) (-) (-)
+1 A Rigid {ex!r} {ey!r} -0.4254613873973588
+--------------------- OPTIONS ------------------------------------------
+0.01 dtM
+0.0 TMax
+--------------------- OUTPUTS ------------------------------------------
+FairTen1
+AnchTen1
+--------------------- need this line -----------------------------------
+"""
+
+A = (14.254910394242287, 60.0359685337248)
+B = (-2330.8134158224498, 2012.8230287080746)
+EZ = (-0.5700516816811296, 0.702868186821989)
+ANGLES = [0.0, 90.0, 200.0, 328.0, 328.29064504314925, 328.5, 328.75, 329.0, 329.25, 329.5, 330.0]
+
+
+def turn(p, th):
+    c, s = math.cos(math.radians(th)), math.sin(math.radians(th))
+    return c * p[0] - s * p[1], s * p[0] + c * p[1]
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--exe", required=True)
+    ap.add_argument("--output", required=True)
+    a = ap.parse_args()
+    out = Path(a.output)
+    out.mkdir(parents=True, exist_ok=True)
+    env = dict(os.environ)
+    env.setdefault("OPENBLAS_NUM_THREADS", "1")
+    failures = []
+    ref = None
+    for th in ANGLES:
+        (xa, ya), (xb, yb), (ex, ey) = turn(A, th), turn(B, th), turn(EZ, th)
+        name = "azimuth_%g" % th
+        deck = out / (name + ".dat")
+        deck.write_text(
+            DECK.format(xa=xa, ya=ya, xb=xb, yb=yb, ex=ex, ey=ey), encoding="utf-8", newline="\n"
+        )
+        root = out / name
+        p = subprocess.run(
+            [a.exe, str(deck), str(root)], capture_output=True, text=True, env=env, cwd=str(out)
+        )
+        if p.returncode != 0:
+            failures.append(
+                "%g deg: exit %d: %s" % (th, p.returncode, (p.stdout + p.stderr)[-300:])
+            )
+            continue
+        rows = Path(str(root) + ".out").read_text(encoding="utf-8").splitlines()
+        vals = [float(v) for v in rows[2].split()[1:3]]
+        if ref is None:
+            ref = vals
+        for v, r, ch in zip(vals, ref, ("FairTen1", "AnchTen1")):
+            if abs(v - r) > 1.0e-6 * abs(r):
+                failures.append("%g deg: %s %.8e differs from %.8e" % (th, ch, v, r))
+    print(
+        "lazy wave turned about the vertical: %d azimuths, %d failures"
+        % (len(ANGLES), len(failures))
+    )
+    for f in failures:
+        print("FAIL:", f)
+    return 1 if failures else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())