From 79ce13e67684803d3020046cb9b4be3d63212cfe Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 26 Oct 2026 16:36:26 +0900 Subject: [PATCH 1/3] fix(statics): converge weightless taut finite-EI lines The Hermite static solve judges convergence by the residual relative to a force scale: the line's total submerged weight, with a 1 N floor. A weightless or near-neutral taut line (a neutrally buoyant cable, a tendon) carries tension far above its weight, and the round-off of its axial forces EA (|r'| - 1) then lies far above tol x scale: a neutrally buoyant 60 m cable at 61 kN stopped at ||R|| 1e-4 N against a 1e-6 N test and failed every static route, even lying straight along its chord. A taut line, its ends farther apart than its length L, carries at least (chord/L - 1) times its series axial stiffness whatever its weight; the force scale now includes that tension. A slack line (chord <= L) keeps its weight scale, so every such solve is unchanged: all 99 example output files are identical to dev. The weightless cable clamped 30 deg off its chord now converges and follows an independent solution of the planar elastica within 0.14 % (root moment) and 0.07 % (end tension). Signed-off-by: Jae Hoon Seo --- CHANGELOG.md | 7 +++++++ src/CableDyn_HermiteCableStatic.f90 | 12 ++++++++++- tests/test_hermite_cable_static.f90 | 31 +++++++++++++++++++++++++++++ 3 files changed, 49 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 41a39e8..7a29887 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,13 @@ All notable changes to CableDyn are recorded in this file. The format follows range graphs and element file) report the cable's own moment, the line type's EI times the curvature; lines without stiffeners are unchanged. +### Fixed + +- A weightless or near-neutral taut finite-EI line (a neutrally buoyant cable, a tendon) failed + its static solve, even straight along its chord: the convergence test scaled the residual by + the line's weight alone, far below the round-off of its axial forces. A taut line's scale now + includes the tension its stretched chord implies; slack lines are unchanged. + ## [0.3.0] - 2026-10-23 ### Added diff --git a/src/CableDyn_HermiteCableStatic.f90 b/src/CableDyn_HermiteCableStatic.f90 index 687b354..582aeeb 100644 --- a/src/CableDyn_HermiteCableStatic.f90 +++ b/src/CableDyn_HermiteCableStatic.f90 @@ -467,7 +467,7 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_ TYPE(CD_HermiteTorsionType), INTENT(INOUT), OPTIONAL :: torsion INTEGER :: ne, nn, ndof, c, it, es, i, b, nci, n_buoy, omp_threads, axial_order, bending_order - REAL(wp) :: eifac, rnorm, fscale, bfac, rmerit, ptc_tau, ptc_rprev, pi_total, pi_abs, lm_mu + 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 REAL(wp) :: rigid_basis(3, 3, 2), projected_tangent(3) @@ -892,6 +892,16 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_ fscale = MAX(fscale, MAXVAL(ABS(f_nodal(6*(i - 1) + 4:6*(i - 1) + 6)))/trib(i)) END DO END IF + ! The tension a taut line must carry: with its ends farther apart than its length L, a line + ! is stretched by at least chord/L - 1 over its series axial stiffness, so it carries at + ! least that force whatever its weight. A line that carries little weight for its tension + ! -- a weightless or near-neutral taut line, a tendon -- otherwise had a weight scale far + ! below the round-off of its axial forces (EA (|r'| - 1) is formed from coordinate + ! differences), and Newton stagnated at the round-off floor of an equilibrium it had + ! reached: a neutrally buoyant 60 m cable at 61 kN stopped at ||R|| 1e-4 N against tol x + ! 1 N. A slack line (chord <= L) keeps its weight scale. + chord = NORM2(q(6*(nn - 1) + 1:6*(nn - 1) + 3) - q(1:3)) + IF (chord > SUM(l0)) fscale = MAX(fscale, (chord/SUM(l0) - CD_ONE)*SUM(l0)/SUM(l0/EA)) ! Buoyancy load continuation. A net-buoyant section (w < 0, e.g. a lazy-wave buoyancy ! module) makes the equilibrium non-unique: alongside the smooth buoyant arch there are diff --git a/tests/test_hermite_cable_static.f90 b/tests/test_hermite_cable_static.f90 index 0dd1b2a..7185132 100644 --- a/tests/test_hermite_cable_static.f90 +++ b/tests/test_hermite_cable_static.f90 @@ -42,6 +42,7 @@ PROGRAM test_hermite_cable_static CALL check_seabed_mesh_jump() CALL check_scalar_contact_on_slope() CALL check_taut_line_anchored_on_floor() + CALL check_weightless_taut_line() CALL check_laid_line_node_recovery() CALL check_end_support_energy_gradient() CALL check_frictionless_slope_run() @@ -314,6 +315,36 @@ SUBROUTINE check_seabed_mesh_jump() CALL require(dev < 0.1_wp*pen0, 'mesh-jump penetration stays near w/(kBot d)') END SUBROUTINE check_seabed_mesh_jump + SUBROUTINE check_weightless_taut_line() + !! A weightless taut line converges. The relative residual took its force scale from the + !! line's weight alone (1 N floor), far below the round-off of the 61 kN axial force of a + !! stiff line (EA 5e8 N), and Newton stagnated at ||R|| ~ 1e-4 N of an equilibrium it had + !! reached. A 60 m line stretched over a 60.0073 m chord on 212 elements (0.05 m over its + !! first 8 m) must converge straight, with the axial strain of the chord, C/L - 1. + INTEGER, PARAMETER :: NA = 160, NB = 52, NE = NA + NB, NN = NE + 1 + REAL(wp) :: l0(NE), EA(NE), EI(NE), w(NE), seed(6*NN), q(6*NN), curv(NN), res, chord(3), s, stretch + INTEGER :: fixed(6), i, es, it + CHARACTER(512) :: em + l0(1:NA) = 8.0_wp/NA; l0(NA + 1:) = 1.0_wp + EA = 5.0e8_wp; EI = 2.0e4_wp; w = 0.0_wp + chord = [51.97_wp, 0.0_wp, -30.0_wp] + ! a straight seed along the chord, unit tangents + s = 0.0_wp + DO i = 1, NN + IF (i > 1) s = s + l0(i - 1) + seed(6*(i - 1) + 1:6*(i - 1) + 3) = chord*s/60.0_wp + seed(6*(i - 1) + 4:6*(i - 1) + 6) = chord/NORM2(chord) + END DO + fixed = [1, 2, 3, 6*NE + 1, 6*NE + 2, 6*NE + 3] + CALL CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed, -1.0e4_wp, 0.0_wp, 4, 200, 1.0e-6_wp, & + 1.0_wp, q, curv, res, it, es, em) + CALL require(es == CD_HCSTAT_OK, 'weightless taut line converges: '//TRIM(em)) + IF (es /= CD_HCSTAT_OK) RETURN + stretch = NORM2(q(6*(NE/2) + 4:6*(NE/2) + 6)) + CALL require(ABS(stretch - NORM2(chord)/60.0_wp) <= 1.0e-9_wp, 'weightless taut line takes the chord strain') + CALL require(MAXVAL(ABS(curv)) <= 1.0e-8_wp, 'weightless taut line stays straight') + END SUBROUTINE check_weightless_taut_line + SUBROUTINE check_taut_line_anchored_on_floor() !! The floor carries the end element's weight moment only where the element lies on it. !! A taut line anchored on the floor rises into the water from its anchor: its end element From d610ae5d3529b4691cce5569b56bd2a255e7c5e0 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 26 Oct 2026 17:04:53 +0900 Subject: [PATCH 2/3] fix(statics): use the chord tension only with both ends held The taut-chord tension is a lower bound on the line's tension only while the solve holds both end positions. With an end free to move (a cantilever) the seed's chord implies no tension at equilibrium, and an overstretched seed set a force scale that the solve never revised: a cantilever with EA 1e12 N seeded 1 % long took a 1e10 N scale and read its 1000 N weight as converged before bending. The bound now applies only when the end positions are both fixed; such a cantilever reaches the beam deflection w L^4/(8 EI) from either seed. Signed-off-by: Jae Hoon Seo --- src/CableDyn_HermiteCableStatic.f90 | 10 ++++++--- tests/test_hermite_cable_static.f90 | 32 +++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 3 deletions(-) diff --git a/src/CableDyn_HermiteCableStatic.f90 b/src/CableDyn_HermiteCableStatic.f90 index 582aeeb..e3612dc 100644 --- a/src/CableDyn_HermiteCableStatic.f90 +++ b/src/CableDyn_HermiteCableStatic.f90 @@ -899,9 +899,13 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_ ! below the round-off of its axial forces (EA (|r'| - 1) is formed from coordinate ! differences), and Newton stagnated at the round-off floor of an equilibrium it had ! reached: a neutrally buoyant 60 m cable at 61 kN stopped at ||R|| 1e-4 N against tol x - ! 1 N. A slack line (chord <= L) keeps its weight scale. - chord = NORM2(q(6*(nn - 1) + 1:6*(nn - 1) + 3) - q(1:3)) - IF (chord > SUM(l0)) fscale = MAX(fscale, (chord/SUM(l0) - CD_ONE)*SUM(l0)/SUM(l0/EA)) + ! 1 N. A slack line (chord <= L) keeps its weight scale. The bound holds only while the + ! solve holds both end positions, which keeps the chord; with an end free to move (a + ! cantilever) the seed's chord implies no tension at equilibrium. + IF (.NOT. (ANY(freemask(1:3)) .OR. ANY(freemask(6*(nn - 1) + 1:6*(nn - 1) + 3)))) THEN + chord = NORM2(q(6*(nn - 1) + 1:6*(nn - 1) + 3) - q(1:3)) + IF (chord > SUM(l0)) fscale = MAX(fscale, (chord/SUM(l0) - CD_ONE)*SUM(l0)/SUM(l0/EA)) + END IF ! Buoyancy load continuation. A net-buoyant section (w < 0, e.g. a lazy-wave buoyancy ! module) makes the equilibrium non-unique: alongside the smooth buoyant arch there are diff --git a/tests/test_hermite_cable_static.f90 b/tests/test_hermite_cable_static.f90 index 7185132..de9702b 100644 --- a/tests/test_hermite_cable_static.f90 +++ b/tests/test_hermite_cable_static.f90 @@ -343,8 +343,40 @@ SUBROUTINE check_weightless_taut_line() stretch = NORM2(q(6*(NE/2) + 4:6*(NE/2) + 6)) CALL require(ABS(stretch - NORM2(chord)/60.0_wp) <= 1.0e-9_wp, 'weightless taut line takes the chord strain') CALL require(MAXVAL(ABS(curv)) <= 1.0e-8_wp, 'weightless taut line stays straight') + CALL check_free_end_preload() END SUBROUTINE check_weightless_taut_line + SUBROUTINE check_free_end_preload() + !! The taut-chord tension enters the force scale only while both end positions are held. A + !! cantilever (clamped End A, free End B; EA 1e12 N, EI 1e8 N m^2, 100 N/m over 10 m) seeded + !! 1 % overstretched would otherwise take a 1e10 N scale from its seed, and its 1000 N weight + !! would read as converged before the line bent; it must reach the beam deflection + !! w L^4/(8 EI) = 1.25 mm, as from an unstretched seed. + INTEGER, PARAMETER :: NE = 10, NN = NE + 1 + REAL(wp) :: l0(NE), EA(NE), EI(NE), w(NE), seed_a(6*NN), seed_b(6*NN), qa(6*NN), qb(6*NN), curv(NN), res + REAL(wp) :: tip + INTEGER :: i, es, it + CHARACTER(512) :: em + l0 = 1.0_wp; EA = 1.0e12_wp; EI = 1.0e8_wp; w = 100.0_wp + tip = 100.0_wp*10.0_wp**4/(8.0_wp*1.0e8_wp) + seed_a = 0.0_wp; seed_b = 0.0_wp + DO i = 1, NN + seed_a(6*(i - 1) + 1) = REAL(i - 1, wp) + seed_b(6*(i - 1) + 1) = 1.01_wp*REAL(i - 1, wp) + seed_a(6*(i - 1) + 4) = 1.0_wp + seed_b(6*(i - 1) + 4) = 1.0_wp + END DO + CALL CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed_a, [1, 2, 3, 4, 5, 6], -1.0e4_wp, 0.0_wp, 4, 200, & + 1.0e-6_wp, 1.0_wp, qa, curv, res, it, es, em) + CALL require(es == CD_HCSTAT_OK, 'cantilever from an unstretched seed converges: '//TRIM(em)) + CALL CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed_b, [1, 2, 3, 4, 5, 6], -1.0e4_wp, 0.0_wp, 4, 200, & + 1.0e-6_wp, 1.0_wp, qb, curv, res, it, es, em) + CALL require(es == CD_HCSTAT_OK, 'cantilever from an overstretched seed converges: '//TRIM(em)) + CALL require(ABS(-qa(6*NE + 3) - tip) <= 1.0e-3_wp*tip, 'cantilever from an unstretched seed: w L^4/(8 EI)') + CALL require(ABS(-qb(6*NE + 3) - tip) <= 1.0e-3_wp*tip, & + 'cantilever from an overstretched seed: w L^4/(8 EI) (no seed preload in the force scale)') + END SUBROUTINE check_free_end_preload + SUBROUTINE check_taut_line_anchored_on_floor() !! The floor carries the end element's weight moment only where the element lies on it. !! A taut line anchored on the floor rises into the water from its anchor: its end element From 73f492efce106c30b2391a0f3d6190b52373140f Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 26 Oct 2026 17:14:44 +0900 Subject: [PATCH 3/3] fix(statics): skip the chord tension with an EA = 0 element EA = 0 is a valid element input, and the series compliance sum(l0/EA) divided by zero on a held taut seed. An element without axial stiffness leaves the series stiffness, and so the tension bound, at zero: the bound is skipped. The static test checks that no division by zero is raised. Signed-off-by: Jae Hoon Seo --- src/CableDyn_HermiteCableStatic.f90 | 5 +++-- tests/test_hermite_cable_static.f90 | 12 +++++++++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/CableDyn_HermiteCableStatic.f90 b/src/CableDyn_HermiteCableStatic.f90 index e3612dc..ee3a838 100644 --- a/src/CableDyn_HermiteCableStatic.f90 +++ b/src/CableDyn_HermiteCableStatic.f90 @@ -901,8 +901,9 @@ SUBROUTINE CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed_dofs, seabed_ ! reached: a neutrally buoyant 60 m cable at 61 kN stopped at ||R|| 1e-4 N against tol x ! 1 N. A slack line (chord <= L) keeps its weight scale. The bound holds only while the ! solve holds both end positions, which keeps the chord; with an end free to move (a - ! cantilever) the seed's chord implies no tension at equilibrium. - IF (.NOT. (ANY(freemask(1:3)) .OR. ANY(freemask(6*(nn - 1) + 1:6*(nn - 1) + 3)))) THEN + ! cantilever) the seed's chord implies no tension at equilibrium. An element without axial + ! stiffness (EA = 0) leaves the series stiffness, and so the bound, at zero. + IF (.NOT. (ANY(freemask(1:3)) .OR. ANY(freemask(6*(nn - 1) + 1:6*(nn - 1) + 3))) .AND. ALL(EA > CD_ZERO)) THEN chord = NORM2(q(6*(nn - 1) + 1:6*(nn - 1) + 3) - q(1:3)) IF (chord > SUM(l0)) fscale = MAX(fscale, (chord/SUM(l0) - CD_ONE)*SUM(l0)/SUM(l0/EA)) END IF diff --git a/tests/test_hermite_cable_static.f90 b/tests/test_hermite_cable_static.f90 index de9702b..d3d4aae 100644 --- a/tests/test_hermite_cable_static.f90 +++ b/tests/test_hermite_cable_static.f90 @@ -25,7 +25,8 @@ PROGRAM test_hermite_cable_static USE CableDyn_Line, ONLY: CD_Nodal_Seabed_Stiffness, CD_LINE_OK USE CableDyn_Bathymetry, ONLY: CD_BathymetryType, CD_Init_Bathymetry USE, INTRINSIC :: IEEE_ARITHMETIC, ONLY: IEEE_IS_FINITE - USE, INTRINSIC :: IEEE_EXCEPTIONS, ONLY: IEEE_USUAL, IEEE_GET_HALTING_MODE, IEEE_SET_HALTING_MODE, IEEE_SET_FLAG + USE, INTRINSIC :: IEEE_EXCEPTIONS, ONLY: IEEE_USUAL, IEEE_GET_HALTING_MODE, IEEE_SET_HALTING_MODE, IEEE_SET_FLAG, & + IEEE_GET_FLAG, IEEE_DIVIDE_BY_ZERO IMPLICIT NONE LOGICAL :: fp_halt(3) ! saved IEEE halting modes around deliberately non-finite inputs @@ -323,6 +324,7 @@ SUBROUTINE check_weightless_taut_line() !! first 8 m) must converge straight, with the axial strain of the chord, C/L - 1. INTEGER, PARAMETER :: NA = 160, NB = 52, NE = NA + NB, NN = NE + 1 REAL(wp) :: l0(NE), EA(NE), EI(NE), w(NE), seed(6*NN), q(6*NN), curv(NN), res, chord(3), s, stretch + LOGICAL :: divided INTEGER :: fixed(6), i, es, it CHARACTER(512) :: em l0(1:NA) = 8.0_wp/NA; l0(NA + 1:) = 1.0_wp @@ -344,6 +346,14 @@ SUBROUTINE check_weightless_taut_line() CALL require(ABS(stretch - NORM2(chord)/60.0_wp) <= 1.0e-9_wp, 'weightless taut line takes the chord strain') CALL require(MAXVAL(ABS(curv)) <= 1.0e-8_wp, 'weightless taut line stays straight') CALL check_free_end_preload() + ! An element with EA = 0 (a valid input) leaves the series stiffness at zero: the bound is + ! skipped, not computed as l0/EA (no division by zero on a held taut seed). + EA(NA + 1) = 0.0_wp + CALL IEEE_SET_FLAG(IEEE_DIVIDE_BY_ZERO, .FALSE.) + CALL CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed, -1.0e4_wp, 0.0_wp, 1, 1, 1.0e-6_wp, & + 1.0_wp, q, curv, res, it, es, em) + CALL IEEE_GET_FLAG(IEEE_DIVIDE_BY_ZERO, divided) + CALL require(.NOT. divided, 'a held taut line with an EA = 0 element raises no division by zero') END SUBROUTINE check_weightless_taut_line SUBROUTINE check_free_end_preload()