diff --git a/CMakeLists.txt b/CMakeLists.txt index 54d8a9d..e7cf48b 100644 --- a/CMakeLists.txt +++ b/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 $ + --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. diff --git a/doc/theory.rst b/doc/theory.rst index 6d84938..ee8f8b6 100644 --- a/doc/theory.rst +++ b/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 diff --git a/doc/troubleshooting.rst b/doc/troubleshooting.rst index 203f9e7..fd141c6 100644 --- a/doc/troubleshooting.rst +++ b/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 diff --git a/src/CableDyn_HermiteCableStatic.f90 b/src/CableDyn_HermiteCableStatic.f90 index f7b4c8b..6451254 100644 --- a/src/CableDyn_HermiteCableStatic.f90 +++ b/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= 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 diff --git a/tests/test_hermite_static_azimuth.py b/tests/test_hermite_static_azimuth.py new file mode 100644 index 0000000..b777e30 --- /dev/null +++ b/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 --output +""" + +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())