diff --git a/CHANGELOG.md b/CHANGELOG.md index 41a39e8e..7a298876 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 687b3541..ee3a8387 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,21 @@ 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. 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. 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 ! 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 0dd1b2a5..d3d4aae8 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 @@ -42,6 +43,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 +316,77 @@ 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 + LOGICAL :: divided + 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') + 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() + !! 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