diff --git a/CHANGELOG.md b/CHANGELOG.md index d54a840..7cc4c05 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,15 @@ All notable changes to CableDyn are recorded in this file. The format follows `EI = 0` + finite-EI route) the finite-EI cable steps with their residual, tangent and solve counts and times. Every output file is the same with and without it. See the manual's command-line page (Run-time profile). +- A note when a finite-EI line starts its dynamics from a static state that is unstable out of + its vertical plane. A line solved in that plane (free in torsion, no current or end direction + across it) has its out-of-plane DOFs held in the statics and released in the dynamics. On + elements longer than half the bending length the touchdown can make that state unstable out of + the plane, which the planar solve does not see. One stability test with the dynamic + constraints now runs at the hand-off (0.35 ms on an 87-element lazy wave, 6 ms on 2048 + elements; on a line with `ATTACHMENTS`, on the discrete equilibrium after the attachment + polish) and the note names the line, its mesh and the factor that meets the bending length, + and recommends a finer mesh or `adaptive_mesh True`. Results do not change. - Torsion of finite-EI lines in coupled OpenFAST (`CompMooring = 5`) and FAST.Farm runs. The End A frame of a line restrained in torsion turns with the platform orientation, the torque moment is returned to the platform, `PtfmInit` twists the line in the static initialisation, @@ -60,6 +69,19 @@ All notable changes to CableDyn are recorded in this file. The format follows `MD_Shared` farm (`validation/scripts/fastfarm_shared_ab.py`): tension means agree to 0.02 %, extremes to 0.25 %, and platform motions to millimetres. The migration guide now says that on a line with no `Fixed` end `FairTen` and `AnchTen` name the opposite ends from MoorDyn's. +- The power cable of `examples/iea15mw_umaine_mixed_cabledyn.dat` is meshed to the + half-bending-length rule: 248 elements (2 m, 2.5 m and 2.5 m in its three sections) instead + of 87 (4 m, 4 m and 10 m). On the 10 m grounded elements the touchdown was under-resolved and + the static state, solved in the vertical plane of the cable, was unstable out of that plane. + The static fairlead tension of the cable changes from 11 898 N to 11 922 N (11 921 N on a + 1 m mesh), its anchor tension from 1218 N to 1230 N, and the touchdown curvature from + 0.058 1/m to 0.079 1/m (0.080 1/m on the 1 m mesh). The curvature output moves with the + mesh: `Curv4N20`/`BendMom4N20` (76 m from the hang-off) become `Curv4N38`/`BendMom4N38` + (74 m, the peak of the sag bend). In the one-hour full-turbine run of the validation record + the hang-off tension statistics are unchanged (11.5 ± 1.1 kN), the touchdown curvature peaks + at 0.109 1/m instead of 0.075 1/m, and the anchor tension stays at its static 1.23 kN, where + the coarse mesh showed stick-slip peaks up to 4.1 kN. The validation record and the tutorials + are updated. ### Fixed @@ -79,9 +101,9 @@ All notable changes to CableDyn are recorded in this file. The format follows out of its plane. The seabed contact acts at the nodes, which carry no moment, so at an anchor resting on the floor the end moment of the end element's weight turned the free end tangent into the floor: the element dipped below the seabed between its nodes (27 mm for the - power cable of `examples/iea15mw_umaine_mixed_cabledyn.dat`) and carried axial forces from - −47 kN to +127 kN about a line tension of 1.2 kN. Once the hang-off swayed, the element - buckled sideways; with seabed friction the anchor tension then locked at about 2.5 times its + power cable of `examples/iea15mw_umaine_mixed_cabledyn.dat` on its former 10 m grounded + elements) and carried axial forces from −47 kN to +127 kN about a line tension of 1.2 kN. + Once the hang-off swayed, the element buckled sideways; with seabed friction the anchor tension then locked at about 2.5 times its static value for the rest of the run (1.22 kN to 3.1 kN in a coupled OpenFAST run while the platform drifted toward the anchor), and without friction the end nodes wandered by up to 0.18 m. The seabed now carries the end element's weight with its end moment at a held line end @@ -102,6 +124,20 @@ All notable changes to CableDyn are recorded in this file. The format follows grounded run turned sideways on the seabed: when the refinement of the catenary continuation stops on a compression in a current, its failure now names the current limit, so the driver retries out of the plane. +- A finite-EI line solved in 3D (restrained in torsion, or with an end connection direction or + a current out of the vertical plane of its chord) on a deck mesh too coarse at touchdown or at a section + change no longer stops with "the converged state is an unstable equilibrium". The stable + equilibrium had read as a saddle: elements longer than the bending length cannot follow + the bend, the axial force at their Gauss points oscillates into compression, and that + compression lowered the out-of-plane stiffness until the stability test failed. Such a + state is now refined, from the mesh the route returned, to the mesh whose elements are at + most half the bending length: the catenary-continuation route is repeated there, as for a + failed polish, and the state of the mesh-sequenced route is prolongated there and polished. + This also covers a mesh that an adaptive or safety refinement left too coarse; a mesh that + already meets the rule is not refined again. A note names both meshes. An equilibrium that + is still unstable on the refined mesh stops as before. Decks that solved before give the + same results. `tests/data/lazy_wave_touchdown_torsion_coarse.dat` is the regression case, + solved with each route first. ## [0.1.1] - 2026-10-09 diff --git a/CMakeLists.txt b/CMakeLists.txt index b3d6db2..7d28801 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2014,6 +2014,65 @@ add_test(NAME driver_clamped_cable_pinned_rod -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_clamped_cable_pinned_rod -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) set_tests_properties(driver_clamped_cable_pinned_rod PROPERTIES LABELS "fortran;driver;integration") +# A lazy-wave cable restrained in torsion (3D static solve) whose deck mesh is too coarse at +# touchdown: the deck-mesh equilibrium reads as unstable, the refined one is stable and used. +add_test(NAME driver_lazy_wave_torsion_coarse_refined + COMMAND ${CMAKE_COMMAND} -DDRIVER=$ -DEXPECT_RC=0 + "-DARGS=${CMAKE_CURRENT_SOURCE_DIR}/tests/data/lazy_wave_touchdown_torsion_coarse.dat|coarse" + "-DONCE=failed the stability test with elements longer than half the bending length" + "-DMUST_MATCH=stable on the refined mesh of 348 elements;torsion solved" + -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_lazy_wave_torsion_coarse_refined + -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) +set_tests_properties(driver_lazy_wave_torsion_coarse_refined PROPERTIES + LABELS "fortran;driver;integration;static;torsion") +# The same deck with the mesh-sequenced route first: its state is prolongated to the refined mesh. +add_test(NAME driver_lazy_wave_torsion_coarse_sequenced + COMMAND ${CMAKE_COMMAND} -DDRIVER=$ -DEXPECT_RC=0 + "-DARGS=${CMAKE_CURRENT_SOURCE_DIR}/tests/data/lazy_wave_touchdown_torsion_coarse_sequenced.dat|coarse" + "-DONCE=failed the stability test with elements longer than half the bending length" + "-DMUST_MATCH=stable on the refined mesh of 348 elements;torsion solved" + -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_lazy_wave_torsion_coarse_sequenced + -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) +set_tests_properties(driver_lazy_wave_torsion_coarse_sequenced PROPERTIES + LABELS "fortran;driver;integration;static;torsion") +# The same cable free in torsion (planar static): the coarse planar state is unstable out of its +# plane and the hand-off names it once; on a mesh that resolves the bending length it does not. +add_test(NAME driver_lazy_wave_planar_coarse_note + COMMAND ${CMAKE_COMMAND} -DDRIVER=$ -DEXPECT_RC=0 + "-DARGS=${CMAKE_CURRENT_SOURCE_DIR}/tests/data/lazy_wave_touchdown_planar_coarse.dat|coarse" + "-DONCE=starts the dynamics from a static state solved in its vertical plane that is unstable" + "-DMUST_MATCH=mesh of 87 elements, the coarsest 10.00 m long;4x meets that length" + -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_lazy_wave_planar_coarse_note + -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) +set_tests_properties(driver_lazy_wave_planar_coarse_note PROPERTIES + LABELS "fortran;driver;integration;static") +add_test(NAME driver_lazy_wave_planar_refined_quiet + COMMAND ${CMAKE_COMMAND} -DDRIVER=$ -DEXPECT_RC=0 + "-DARGS=${CMAKE_CURRENT_SOURCE_DIR}/tests/data/lazy_wave_touchdown_planar_refined.dat|refined" + "-DMUST_NOT=unstable out of" + -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_lazy_wave_planar_refined_quiet + -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) +set_tests_properties(driver_lazy_wave_planar_refined_quiet PROPERTIES + LABELS "fortran;driver;integration;static") +# The same pair with the buoyant stretch built from discrete ATTACHMENTS: the hand-off check runs +# on the discrete equilibrium the attachment polish gives. +add_test(NAME driver_lazy_wave_planar_coarse_att_note + COMMAND ${CMAKE_COMMAND} -DDRIVER=$ -DEXPECT_RC=0 + "-DARGS=${CMAKE_CURRENT_SOURCE_DIR}/tests/data/lazy_wave_touchdown_planar_coarse_att.dat|coarse_att" + "-DONCE=starts the dynamics from a static state solved in its vertical plane that is unstable" + "-DMUST_MATCH=mesh of 87 elements, the coarsest 10.00 m long;4x meets that length;25 attachment.s. lumped" + -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_lazy_wave_planar_coarse_att_note + -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) +set_tests_properties(driver_lazy_wave_planar_coarse_att_note PROPERTIES + LABELS "fortran;driver;integration;static") +add_test(NAME driver_lazy_wave_planar_refined_att_quiet + COMMAND ${CMAKE_COMMAND} -DDRIVER=$ -DEXPECT_RC=0 + "-DARGS=${CMAKE_CURRENT_SOURCE_DIR}/tests/data/lazy_wave_touchdown_planar_refined_att.dat|refined_att" + "-DMUST_NOT=unstable out of" + -DWORKDIR=${CMAKE_CURRENT_BINARY_DIR}/driver_lazy_wave_planar_refined_att_quiet + -P ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_driver_output.cmake) +set_tests_properties(driver_lazy_wave_planar_refined_att_quiet PROPERTIES + LABELS "fortran;driver;integration;static") if(WIN32 AND CABLEDYN_RUNTIME_DLLS) add_test( diff --git a/VALIDATION.md b/VALIDATION.md index 2fcc76c..6349554 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -317,10 +317,10 @@ scripts, results, and source provenance are in ~3% std. Tension agreement with MoorDyn is scored in the settled case above and in the DLC matrix below. - *Mixed deck*: the shipped example `iea15mw_umaine_mixed_cabledyn.dat` (three - chain moorings and an 87-element finite-EI power cable that touches down on the + chain moorings and a 248-element finite-EI power cable that touches down on the 200 m seabed) runs to completion under full coupled physics with finite, bounded - channels. Coupled touchdown is covered by the 40 Hz `openfast_touchdown` test; no - full-DLC touchdown result is claimed. + channels. Coupled touchdown is covered by the 40 Hz `openfast_touchdown` test and by + the one-hour full-turbine run in the coupled finite-EI touchdown row below. - *Wall-clock*: see *Coupled throughput* under [L5](#wall-clock-time). - *`PtfmInit`*: MoorDyn-F's convention (rigid transform of Coupled/Vessel points, undisplaced reference `Position` + `TranslationDisp`); tested by @@ -734,7 +734,7 @@ stop with a named error. | Active tension (`CtrlChan`) | Pass verified | CONTROL parsing, last-segment length/rate modulation, and exact unit tests; coupled A/B in which a ServoDyn cable-control DLL drives the same pay-out/haul-in `DeltaL` schedule into `=5` (CableDyn) and `=3` (stock MoorDyn) in one binary: hold tensions agree to <0.03% and the tension responses to <0.6% (controlled fairleads) / <0.7% (anchors) ([`validation/scripts/activetension/`](https://github.com/SMI-Lab-Inha/CableDyn/tree/main/validation/scripts/activetension), [`validation/scripts/openfast_activetension_ab.py`](https://github.com/SMI-Lab-Inha/CableDyn/blob/main/validation/scripts/openfast_activetension_ab.py)) | | Coupled `Rigid6` bodies + `Point3` buoys | Pass verified | 6-DOF `Rigid6` motion with hydrostatic (`C33`/`C44`/`C55`) restoring, checkpoint state (bit-identical `openfast -restart` of a 13.9 MN-swinging tethered body), and SeaState drag / Froude–Krylov / added mass. `Point3` buoys take buoyancy and drag/added-mass/FK; a submerged volume-only buoy takes FK. Tests: `case_coupled_rigid6`/`_mirror`/`_fluid`, `case_coupled_point3_buoy`/`_fluid` | | Coupled rigid rods | Pass verified | host-driven `RodA/B` endpoints, correction rewind, bit-identical 24-scalar-per-rod checkpoint state, and two Gauss fluid nodes per rod with SeaState drag/Froude–Krylov/added mass; `case_coupled_rod`, `_mirror`, and `_fluid` pass, and the OpenFAST module compiles with rod decks enabled | -| Coupled finite-EI touchdown/contact | Pass verified | flat `WtrDpth` and equivalent constant structured bathymetry initialise and step with matching loads; friction and checkpoint restore are tested. The IEA-15 mixed example has 38 grounded nodes and completes a 20 s aggregate trajectory at `dtM=0.025 s` with 2 m surge, 1 m heave, and 0.5 m/s host current, without fallback or stall recovery, at most 14 Newton iterations. A full-turbine run of the same model (one hour after a 200 s transient; TurbSim normal turbulence, 10.6 m/s, 17 % TI; JONSWAP Hs 8 m, Tp 12 s, aligned with the wind and the cable; AeroDyn with the rotor held at 7.55 rpm, no controller, rigid blades and tower; SeaState grid covering the lines) completes without recoveries or non-finite values. The cable's hang-off tension is 11.5 ± 1.1 kN (static 11.9 kN, maximum 17.9 kN), the same 38 nodes stay grounded, and the peak curvature, 0.160 1/m at the hog bend, is mesh-converged (0.156 1/m on a refined mesh); the 10 m elements of the laid section under-resolve the touchdown curvature by about 30 %. The chain-only twin of the run matches stock MoorDyn over 1200 s to 0.3 % in platform motion and fairlead tension statistics. The anchor tension of the laid cable is 1.39 kN (median, static 1.22 kN) with wave-frequency stick-slip peaks up to 4.1 kN lasting under 1 s; before the anchor-end correction of this release it locked at about 2.75 kN for the whole run | +| Coupled finite-EI touchdown/contact | Pass verified | flat `WtrDpth` and equivalent constant structured bathymetry initialise and step with matching loads; friction and checkpoint restore are tested. The IEA-15 mixed example (its power cable on 248 elements, each at most half the bending length) has 149 grounded nodes and completes a 20 s aggregate trajectory at `dtM=0.025 s` with 2 m surge, 1 m heave, and 0.5 m/s host current, without fallback or stall recovery, at most 12 Newton iterations. A full-turbine run of the same model (one hour after a 200 s transient; TurbSim normal turbulence, 10.6 m/s, 17 % TI; JONSWAP Hs 8 m, Tp 12 s, aligned with the wind and the cable; AeroDyn with the rotor held at 7.55 rpm, no controller, rigid blades and tower; SeaState grid covering the lines) completes without recoveries or non-finite values. The cable's hang-off tension is 11.5 ± 1.1 kN (static 11.9 kN, maximum 17.9 kN). The 149 nodes grounded at the static state stay grounded, joined by one more once the platform has drifted (the touchdown moves 2.5 m toward the platform). The peak curvature is 0.156 1/m at the sag bend, 72 m from the hang-off, and 0.109 1/m at the touchdown. On the former 87-element mesh (10 m elements along the laid section) the peak was 0.160 1/m and the touchdown curvature 0.075 1/m, under-resolved by about 30 %. The chain-only twin of the run matches stock MoorDyn over 1200 s to 0.3 % in platform motion and fairlead tension statistics. The anchor tension of the laid cable stays at its static 1.23 kN (maximum 1.24 kN): the seabed friction of the 370 m laid section holds the tension changes at the touchdown. On the former mesh it read 1.39 kN (median) with wave-frequency stick-slip peaks up to 4.1 kN lasting under 1 s, which the resolved mesh does not show, and before the anchor-end correction of this release it locked at about 2.75 kN for the whole run. The run takes 27.1 min of wall time (18.8 min on the former mesh) | | `WaterKin` from file | Pass verified | See [WaterKin and SeaState modes](#waterkin-and-seastate-modes) | | Syrope polyester (`SYROPE:` syntax; [Falkenberg et al., 2017](https://cabledyn.readthedocs.io/en/latest/references.html#ref-falkenberg2017)) | Pass verified | MoorDyn-F/C `BA_s|BA_d` rate and force equations; analytic position/velocity tangents with finite-difference checks; projected, bracketed backward-Euler slow-state update with a strong-unloading boundary test; OWC running maximum; `SYROPE IC` `Tmax0/Tmean0` history with range/order checks; both states carried through every lifecycle stage and the coupled checkpoint. `BA_s` is kept separate from line axial damping (`BA`) and tested against double counting. Composite, hydrodynamic, finite-EI, current, and wave combinations stop with a named error | | VIV (`compVIV`) | Not supported | outside CableDyn's scope | diff --git a/doc/modeling.rst b/doc/modeling.rst index e3249cd..ecaaf2d 100644 --- a/doc/modeling.rst +++ b/doc/modeling.rst @@ -86,6 +86,22 @@ refine regions with high curvature or rapidly varying load: hang-off and bend-st buoyancy-module transitions, sag bends, touchdown, and short weighted attachments. Avoid using a material boundary solely as a substitute for mesh refinement. +**Starting mesh for a finite-EI line.** At touchdown, at a sag or hog bend and on each side of a +section change, keep the element length at most half the bending length +:math:`(EI/|w|)^{1/3}` of the section, with :math:`w` its submerged weight per metre. For a +power cable with :math:`EI = 2\times10^4` N m² and :math:`w = 158` N/m that is 2.5 m. A longer +element cannot follow the bend: it stays stiff in bending, and its stiff axial term absorbs the +mismatch, so the pointwise axial force oscillates about the element mean, by tens of kN on a +lazy wave of 1–2 kN tension at 10 m elements. The ``.t.out`` segment tensions, evaluated at the +element midpoints, show these spikes and negative values near touchdown and at section changes; +they shrink as the mesh is refined, and the element-mean force (what enters equilibrium, see +``tensile_safety`` in :doc:`theory`) is far smoother. The compressed Gauss points also lower the +out-of-plane stiffness: a line solved in 3D (restrained in torsion, or with an end connection +direction or a current out of its vertical plane) can read as an unstable equilibrium on such a +mesh. The static solve then refines that state to the rule above, whichever route reached it and +whether or not an adaptive refinement preceded, and reports it in a note (see +:doc:`troubleshooting`); meshing to the rule avoids the extra solve. + Perform a mesh study on engineering observables, not node identity. At minimum compare endpoint tension and the peak/position of curvature or touchdown response. Refine until the change is below the project tolerance. Cross-code comparisons should interpolate by normalised arc length as diff --git a/doc/openfast.rst b/doc/openfast.rst index 1e70136..2d7b687 100644 --- a/doc/openfast.rst +++ b/doc/openfast.rst @@ -201,7 +201,7 @@ Mixed moorings and power cables One deck may combine ``EI = 0`` chains and finite-EI (``EI > 0``) power cables; each line uses its own element. In ``iea15mw_umaine_mixed_cabledyn.dat`` line 4 runs from a platform hang-off to a fixed termination on the 200 m seabed; its bare/buoyant/bare sections form the lazy wave and -a grounded tail (38 of 88 nodes on the seabed at equilibrium). ``WtrDpth``, ``kBot``, ``cBot``, +a grounded tail (149 of 249 nodes on the seabed at equilibrium). ``WtrDpth``, ``kBot``, ``cBot``, and ``frictionMu`` define contact; a structured ``bathymetryFile`` may replace the flat floor. CableDyn finds touchdown during static initialisation. Inspect ``.CD.static.out`` (curvature, grounded length, section transitions, minimum bend radius) before the time history. diff --git a/doc/troubleshooting.rst b/doc/troubleshooting.rst index 4646882..fbfe63b 100644 --- a/doc/troubleshooting.rst +++ b/doc/troubleshooting.rst @@ -463,6 +463,17 @@ but a few geometries are genuinely hard: of the twist cannot be trusted beyond it): reduce ``dtM`` or slow the imposed roll. ``tangent turns by more than 120 degrees`` names an element or node pair where the twist of the centreline is no longer reliable: refine the mesh there. +- **A finite-EI equilibrium reported unstable.** ``found no physical equilibrium ... the + converged state is an unstable equilibrium (tangent stiffness not positive definite)`` means + the state that balanced its residual is a saddle of the potential energy: a compressed strut, a + buckled arch or an inverted fold. A coarse mesh can also produce it. On elements longer than + half the bending length :math:`(EI/|w|)^{1/3}` the axial force oscillates into compression at + Gauss points near touchdown or a section change, which lowers the out-of-plane stiffness of a + line solved in 3D. Each static route therefore refines such a state to a mesh that resolves + the bending length before it fails: the catenary-continuation route is repeated there, and + the mesh-sequenced state is prolongated and polished there. If the message names that refined + mesh as well (``still unstable on the refined mesh``), the instability is physical: check the + layout, the line length and the buoyancy. - **A cold start at full load.** The static solve continues the load in stages; if you have hand-tuned a case into a bad basin, remove the tuning and let the default continuation run from the catenary/arch seed. @@ -527,6 +538,29 @@ Warnings and notes that do not stop a run Newton system loses the precision to converge. An 80 m lazy wave (6.5 mm limit) solves at 8 mm elements and not at 4 to 6 mm, where the solve stops within seconds and names the mesh. Coarsen ``NumSegs`` to at least the length the note gives. + * - ``Note: line ...: the equilibrium on the deck mesh (... elements) failed the stability test + with elements longer than half the bending length; stable on the refined mesh of ... + elements`` (or ``the equilibrium on the mesh of ... elements (refined from the deck mesh of + ...) failed ...``) + - on the deck mesh, or on a refinement of it that was still too coarse, the static state + read as unstable, an artefact of elements longer than half the bending length + :math:`(EI/|w|)^{1/3}` (see :doc:`modeling`, step 5). The run uses + the refined mesh, statics and dynamics alike, and its outputs have that many segments. + Mesh the touchdown and section changes to that rule to keep the deck mesh. + * - ``Note: line ... starts the dynamics from a static state solved in its vertical plane that + is unstable out of that plane on its mesh of ... elements, the coarsest ... m long where + half the bending length (EI/|w|)^(1/3) is ... m; ...`` + - a finite-EI line free in torsion is solved in its vertical plane, and the dynamics release + its out-of-plane DOFs. On elements longer than half the bending length the axial force at + touchdown or a section change oscillates into compression and makes that planar state + unstable out of the plane, so a lateral disturbance (a cross current, a sway motion) can + move the line from it in a way the mesh, not the physics, decides. On a line with + ``ATTACHMENTS`` the test runs on the discrete equilibrium, after the attachment polish, + and the bending length uses the smeared weight of the buoyant stretch. The results are + those of the deck as written. Refine the mesh by the factor the note gives, or set + ``adaptive_mesh True``; a line restrained in torsion is solved in 3D and refined + automatically (see :doc:`modeling`, step 5). When ``adaptive_mesh True`` is already set, + its own refinement criteria stopped short of the bending length: refine the deck mesh. * - ``Note: line ... mesh too coarse for its bends`` - the element-mean axial force of a finite-EI line, averaged with its neighbours, is compressive, and the pointwise resultant oscillates well beyond it. Refine ``NumSegs`` in diff --git a/doc/tutorial_openfast.rst b/doc/tutorial_openfast.rst index 50fb7b8..a04334e 100644 --- a/doc/tutorial_openfast.rst +++ b/doc/tutorial_openfast.rst @@ -315,15 +315,18 @@ seabed. Copy it into the case folder and point a copy of the ``.fst`` at it: Created CableDyn model: 4 line object(s), 8 point(s), 6 section(s) [EI=0: 3, finite-EI: 1]. ... - Line 4 fairlead effective tension: 11898 N - force [Fx, Fy, Fz]: [1218.2, 0, -11835] N, inclination=84.123 deg - line tangent: inclination=83.927 deg, declination=173.93 deg, azimuth=0 deg + Line 4 fairlead effective tension: 11922 N + force [Fx, Fy, Fz]: [1230, 0, -11858] N, inclination=84.078 deg + line tangent: inclination=83.943 deg, declination=173.94 deg, azimuth=0 deg ... - Curv4N20 = 9.15231E-02 (1/m) - BendMom4N20 = 1821.3 (N.m) + Curv4N38 = 9.06020E-02 (1/m) + BendMom4N38 = 1803 (N.m) CableDyn found the cable's touchdown and lazy-wave shape by itself: in -``Mixed_CableDyn.CD.static.out``, 38 of line 4's 88 nodes lie on the seabed. The chain moorings +``Mixed_CableDyn.CD.static.out``, 149 of line 4's 249 nodes lie on the seabed. Node 38, 74 m +from the hang-off, is the peak of the sag bend. The cable's elements are at most half its bending +length :math:`(EI/|w|)^{1/3}` (the deck's comments give the numbers), which resolves the bends +and the touchdown. The chain moorings use the cable element and the power cable the cubic-Hermite bending element, in one deck and one coupled solve. Plot line 4's ``Curvature`` against ``ArcLength`` before looking at the time history, and request ``Curv4N``/``BendMom4N`` channels at the hang-off, sag bend, arch, and diff --git a/examples/iea15mw_umaine_mixed_cabledyn.dat b/examples/iea15mw_umaine_mixed_cabledyn.dat index 157284e..2bd92fe 100644 --- a/examples/iea15mw_umaine_mixed_cabledyn.dat +++ b/examples/iea15mw_umaine_mixed_cabledyn.dat @@ -34,9 +34,14 @@ LineID LineType Length NumSegs 1 main 850.0 50 2 main 850.0 50 3 main 850.0 50 -4 bare 80.0 20 -4 buoy 100.0 25 -4 bare 420.0 42 +4 bare 80.0 40 +4 buoy 100.0 40 +4 bare 420.0 168 +-- Line 4 mesh: every element at most half the bending length (EI/|w|)^(1/3) of its section, so +-- the touchdown, the sag and hog bends and the section changes are resolved: bare cable +-- (w = 157.8 N/m) 5.02 m -> 2.51 m, 2.0 m and 2.5 m elements; buoyancy section (|w| = 80.1 N/m) +-- 6.28 m -> 3.14 m, 2.5 m elements. A coarser mesh (10 m along the grounded run) reads the +-- planar equilibrium as unstable out of its plane at touchdown (see doc/modeling.rst, step 5). --------------------------- OPTIONS ------------------------------------------- 9.80665 g - Gravitational acceleration (m/s^2); overridden by OpenFAST [default: 9.80665] 1025.0 rhoW - Water density (kg/m^3); overridden by OpenFAST [default: 1025] @@ -71,7 +76,7 @@ none waves - Wave model, height, period, and direction (m, s, d "AnchTen4" "FairIncl4" "AnchIncl4" -"Curv4N20" -"BendMom4N20" +"Curv4N38" +"BendMom4N38" END --------------------------- need this line ------------------------------------ diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index 6872c54..f51d51c 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -1543,13 +1543,18 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, ! A planar equilibrium in a current that is a saddle out of its plane (see route_pass). REAL(wp), ALLOCATABLE :: l0_saddle(:), ea_saddle(:), ei_saddle(:), w_saddle(:), q_saddle(:), curv_saddle(:) INTEGER, ALLOCATABLE :: fixed_saddle(:) - LOGICAL :: saddle_tol_note, saddle_kink + LOGICAL :: saddle_tol_note, saddle_kink, saddle_stab + INTEGER :: saddle_stab_ne REAL(wp) :: saddle_res CHARACTER(320) :: saddle_notes LOGICAL :: primary_done, primary_inadmissible, tolerance_note, stable, kink_refined, sequence_done LOGICAL :: primary_walked, walk_rejected ! Why the bow behind the fairlead was not reached (bow_route), blank when not tried. - CHARACTER(1024) :: bow_msg + CHARACTER(1024) :: bow_msg, refine_msg + ! Set when the deck-mesh equilibrium failed the stability test and the refined one passed. + LOGICAL :: stability_refined + ! Elements of the mesh whose equilibrium failed the stability test before that refinement. + INTEGER :: stability_from_ne LOGICAL :: tensile_error CHARACTER(400) :: note_text CHARACTER(512) :: limit_text @@ -1993,11 +1998,15 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, resolution_note = '' tolerance_note = .FALSE. kink_refined = .FALSE. + stability_refined = .FALSE. + stability_from_ne = 0 sym_broken = .FALSE. planar_unstable = .FALSE. walk_rejected = .FALSE. saddle_tol_note = .FALSE. saddle_kink = .FALSE. + saddle_stab = .FALSE. + saddle_stab_ne = 0 saddle_res = CD_ZERO saddle_notes = '' ! In a current the traced path can end at a fold or a pitchfork of the planar branch @@ -2014,6 +2023,9 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, END IF END IF route_pass: DO pass = 1, 3 + ! set by the catenary-continuation route only: a state the sequenced route produced + ! was not walked + primary_walked = .FALSE. IF (primary_done) EXIT route_pass ! A planar equilibrium that is unstable out of its plane is not repeated by the ! mesh-sequenced route (also planar): the symmetry-breaking pass follows. @@ -2345,6 +2357,13 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, IF (.NOT. stable) THEN sequence_done = .FALSE. sequence_msg = 'the converged state is an unstable equilibrium (tangent stiffness not positive definite)' + ! As for the catenary-continuation route below, on the mesh that resolves the + ! bending length; the sequenced state is prolongated there, not re-solved. + IF (bending_refine_scale(l0_am, ei_am, w_am) >= 2) THEN + CALL sequenced_refined_polish(sequence_done, refine_msg) + IF (.NOT. sequence_done .AND. LEN_TRIM(refine_msg) > 0) & + sequence_msg = TRIM(sequence_msg)//'; '//TRIM(short_message(refine_msg, 240)) + END IF END IF END IF ELSE @@ -2405,6 +2424,16 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, IF (.NOT. stable) THEN primary_done = .FALSE. primary_msg = 'the converged state is an unstable equilibrium (tangent stiffness not positive definite)' + ! A mesh longer than half the bending length can turn a stable equilibrium into a + ! discrete saddle (see stability_refine in hermite_static_primary): the route is + ! repeated once on the mesh refined from the returned one to resolve that length, + ! and its state audited the same way. A returned mesh that already resolves it, + ! adaptively refined or not, is not refined again: its instability stands. + IF (.NOT. primary_walked .AND. bending_refine_scale(l0_am, ei_am, w_am) >= 2) THEN + CALL stability_refined_route(primary_done, refine_msg) + IF (.NOT. primary_done .AND. LEN_TRIM(refine_msg) > 0) & + primary_msg = TRIM(primary_msg)//'; '//TRIM(short_message(refine_msg, 240)) + END IF END IF END IF ! The equilibrium the anchor walk reached (the line too long to hang as a catenary) @@ -2452,6 +2481,14 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, ! toward its end can hold in the plane only as a saddle). Otherwise the out-of-plane ! (symmetry-breaking) pass finds the stable equilibrium. IF ((sequence_done .OR. primary_done) .AND. pass <= 2) THEN + ! On a mesh longer than half the bending length the out-of-plane negative mode can be + ! the mesh's, not a pitchfork of the line: the planar state is first carried to the + ! mesh that resolves that length and polished there (sequenced_refined_polish), and + ! kept if it is not a saddle there; otherwise the coarse state is restored and + ! classified as before. + IF (.NOT. primary_walked .AND. bending_refine_scale(l0_am, ei_am, w_am) >= 2) THEN + IF (planar_current_saddle()) CALL saddle_refined_retry() + END IF IF (planar_current_saddle()) THEN primary_done = .FALSE. sequence_done = .FALSE. @@ -2463,6 +2500,8 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, saddle_notes = compression_note saddle_tol_note = tolerance_note saddle_kink = kink_refined + saddle_stab = stability_refined + saddle_stab_ne = stability_from_ne saddle_res = res CALL MOVE_ALLOC(l0_am, l0_saddle) CALL MOVE_ALLOC(ea_am, ea_saddle) @@ -2474,6 +2513,7 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, compression_note = '' tolerance_note = .FALSE. kink_refined = .FALSE. + stability_refined = .FALSE. IF (ALLOCATED(sequence_iters)) DEALLOCATE (sequence_iters) IF (ALLOCATED(q_sequence)) DEALLOCATE (q_sequence) IF (ALLOCATED(q_coupled_tol)) DEALLOCATE (q_coupled_tol) @@ -2503,6 +2543,8 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, compression_note = saddle_notes tolerance_note = saddle_tol_note kink_refined = saddle_kink + stability_refined = saddle_stab + stability_from_ne = saddle_stab_ne res = saddle_res primary_done = .TRUE. saddle_kept = .TRUE. @@ -2558,6 +2600,19 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, CALL report_note('line '//TRIM(int_to_str(ln%id))//': '//TRIM(compression_note)) END IF END IF + IF (stability_refined .AND. (sequence_done .OR. primary_done)) THEN + IF (stability_from_ne == SIZE(l0)) THEN + WRITE (note_text, '(A,I0,A,I0,A,I0,A)') 'line ', ln%id, ': the equilibrium on the deck mesh (', SIZE(l0), & + ' elements) failed the stability test with elements longer than half the bending length; '// & + 'stable on the refined mesh of ', SIZE(l0_am), ' elements' + ELSE + WRITE (note_text, '(A,I0,A,I0,A,I0,A,I0,A)') 'line ', ln%id, ': the equilibrium on the mesh of ', & + stability_from_ne, ' elements (refined from the deck mesh of ', SIZE(l0), ') failed the stability '// & + 'test with elements longer than half the bending length; stable on the refined mesh of ', SIZE(l0_am), & + ' elements' + END IF + CALL report_note(note_text) + END IF IF (kink_refined .AND. primary_done) THEN WRITE (note_text, '(A,I0,A,I0,A,I0,A)') 'line ', ln%id, ': the deck mesh (', SIZE(l0), & ' elements) could not resolve the curvature of the equilibrium; refined to ', SIZE(l0_am), ' elements' @@ -2616,6 +2671,8 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, ', xi=', diag%axial_strain_min_xi, ', strain=', diag%axial_strain_min CALL fail_solve(ErrStat, ErrMsg, TRIM(em)); RETURN END IF + ! A line with ATTACHMENTS is checked on its discrete equilibrium, after the attachment polish. + IF (.NOT. has_att) CALL handoff_lateral_note(l0, l0_am, ea_am, ei_am, w_am, w_am, diam, dry_b, fixed_am, .FALSE.) ! Build the coarse Morison arrays, then map every per-element quantity onto the ! refined reference-arc mesh. This covers both global refinement and the local ! axial-oscillation refinement used by tensile safety. Every child inherits the @@ -2668,6 +2725,10 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & endconn_mode, opts, static_tol, ErrStat, ErrMsg, bathymetry, cur_obj) IF (ErrStat /= CD_DECKDRV_OK) RETURN + ! On the final mesh, with the bare weights the polished state carries; the bending-length + ! rule of the note uses the smeared section weights w. + CALL handoff_lateral_note(l0, l0, ea, ei, w_bare, w, diam, opts%rhoW*opts%g*QUARTER_PI*diam**2, fixed_dofs, & + .TRUE.) w = w_bare WRITE (note_text, '(A,I0,A,I0,A,ES10.3,A)') 'line ', ln%id, ': ', SIZE(att_node), & ' attachment(s) lumped at the nearest nodes (largest arc offset ', att_offset, ' m)' @@ -2932,6 +2993,214 @@ SUBROUTINE bow_route(done, why) END IF END SUBROUTINE bow_route + SUBROUTINE stability_refined_route(done, why) + !! The catenary-continuation route repeated on the mesh that resolves the bending length, + !! refined from the mesh the route returned (the deck mesh, or a refinement of it that is + !! still too coarse), after that state failed the stability test (hermite_static_primary + !! with stability_refine). done is set with the audited, stable refined state in the *_am + !! outputs; otherwise why says what stopped it. + LOGICAL, INTENT(OUT) :: done + CHARACTER(*), INTENT(OUT) :: why + REAL(wp), ALLOCATABLE :: l0_s(:), ea_s(:), ei_s(:), w_s(:), diam_s(:), dry_s(:) + INTEGER, ALLOCATABLE :: fixed_s(:) + INTEGER :: es_r + CHARACTER(1024) :: em_r + LOGICAL :: inadm_r, walked_r, tol_r, kink_r + done = .FALSE. + why = '' + stability_from_ne = SIZE(l0_am) + ! The returned mesh is the input of the repetition, whose outputs overwrite the *_am arrays. + CALL MOVE_ALLOC(l0_am, l0_s) + CALL MOVE_ALLOC(ea_am, ea_s) + CALL MOVE_ALLOC(ei_am, ei_s) + CALL MOVE_ALLOC(w_am, w_s) + CALL MOVE_ALLOC(fixed_am, fixed_s) + diam_s = diam_on(l0, l0_s, diam) + dry_s = dry_b + IF (SIZE(l0_s) /= SIZE(l0)) CALL remap_elem_by_reference_arc(dry_s, l0, l0_s) + IF (PRESENT(bathymetry) .AND. .NOT. bathymetry_flat) THEN + CALL hermite_static_primary(anchor, fairlead, l0_s, ea_s, ei_s, w_s, diam_s, dry_s, fixed_s, & + contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & + endconn_mode, opts, static_tol, l0_am, ea_am, ei_am, w_am, q_static, curv, & + fixed_am, applied_scale, diag, res, iters, es_r, em_r, inadm_r, tol_r, kink_r, & + bathymetry=bathymetry, cur_obj=cur_obj, walked=walked_r, stability_refine=.TRUE.) + ELSE + CALL hermite_static_primary(anchor, fairlead, l0_s, ea_s, ei_s, w_s, diam_s, dry_s, fixed_s, & + contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & + endconn_mode, opts, static_tol, l0_am, ea_am, ei_am, w_am, q_static, curv, & + fixed_am, applied_scale, diag, res, iters, es_r, em_r, inadm_r, tol_r, kink_r, & + cur_obj=cur_obj, walked=walked_r, stability_refine=.TRUE.) + END IF + IF (es_r /= CD_HCSTAT_OK) THEN + why = em_r + RETURN + END IF + ! A walked refined state is not substituted for an unwalked one. + IF (walked_r) THEN + why = 'the refined route reached its state by walking the anchor' + RETURN + END IF + CALL audit_refined_state(done, why) + IF (done) tolerance_note = tol_r + END SUBROUTINE stability_refined_route + + SUBROUTINE sequenced_refined_polish(done, why) + !! The mesh-sequenced route's counterpart of stability_refined_route: its state, which + !! failed the stability test on a mesh longer than half the bending length, is carried to + !! the mesh that resolves that length by cubic-Hermite prolongation and polished there, as + !! the route's own refinement does, so the branch the sequencing selected is kept and the + !! fine mesh is never cold-started. done is set with the audited, stable refined state in + !! the *_am outputs; otherwise why says what stopped it. + LOGICAL, INTENT(OUT) :: done + CHARACTER(*), INTENT(OUT) :: why + REAL(wp), ALLOCATABLE :: l0_r(:), q_r(:), ea_r(:), ei_r(:), w_r(:), diam_r(:), dry_r(:) + INTEGER, ALLOCATABLE :: fixed_r(:) + INTEGER :: es_r, it_r, scale, k, parent, k_tol + REAL(wp) :: tol_r(3) + REAL(wp), PARAMETER :: CONTACT_FLOOR_TOL = 1.0e-5_wp + CHARACTER(1024) :: em_r + done = .FALSE. + why = '' + stability_from_ne = SIZE(l0_am) + scale = bending_refine_scale(l0_am, ei_am, w_am) + IF (scale < 2) RETURN + CALL CD_HermiteCable_Refine_Mesh(l0_am, q_static, fixed_am, scale, l0_r, q_r, fixed_r, es_r, em_r, & + inherit_dofs=[2, 5]) + IF (es_r /= CD_HCSTAT_OK) THEN + why = em_r + RETURN + END IF + ALLOCATE (ea_r(SIZE(l0_r)), ei_r(SIZE(l0_r)), w_r(SIZE(l0_r))) + DO k = 1, SIZE(l0_r) + parent = (k - 1)/scale + 1 + ea_r(k) = ea_am(parent) + ei_r(k) = ei_am(parent) + w_r(k) = w_am(parent) + END DO + diam_r = diam_on(l0, l0_r, diam) + dry_r = dry_b + CALL remap_elem_by_reference_arc(dry_r, l0, l0_r) + ! The tolerance ladder of the routes (the tight target, the coupled static tolerance, the + ! contact floor): the state refined here may have been accepted only at one of the + ! relaxed tolerances, which the refined polish must be allowed to reach as well. + tol_r = [static_tol, MAX(static_tol, COUPLED_STATIC_TOL), & + MERGE(MAX(static_tol, CONTACT_FLOOR_TOL), CD_ZERO, contact_enabled)] + DO k_tol = 1, 3 + IF (k_tol > 1) THEN + IF (.NOT. (tol_r(k_tol) > MAXVAL(tol_r(1:k_tol - 1)))) CYCLE + IF (INDEX(em_r, 'Newton did not converge') == 0) EXIT + END IF + IF (ALLOCATED(l0_am)) DEALLOCATE (l0_am) + IF (ALLOCATED(ea_am)) DEALLOCATE (ea_am) + IF (ALLOCATED(ei_am)) DEALLOCATE (ei_am) + IF (ALLOCATED(w_am)) DEALLOCATE (w_am) + IF (ALLOCATED(q_static)) DEALLOCATE (q_static) + IF (ALLOCATED(curv)) DEALLOCATE (curv) + IF (ALLOCATED(fixed_am)) DEALLOCATE (fixed_am) + IF (PRESENT(bathymetry) .AND. .NOT. bathymetry_flat) THEN + CALL CD_HermiteCable_Static_Solve_AutoMesh(l0_r, ea_r, ei_r, w_r, q_r, fixed_r, contact_z, CD_ZERO, & + 8, 80, tol_r(k_tol), 0.7_wp, CD_HC_HKAPPA_TARGET, CD_ZERO, & + CD_HC_REFINE_CAP, [2, 5], l0_am, ea_am, ei_am, w_am, & + q_static, curv, fixed_am, applied_scale, diag, res, it_r, & + es_r, em_r, contact_diameter=diam_r, contact_kbot=opts%kbot, & + bathymetry=bathymetry, contact_frame_cs=frame_cs, & + initial_polish=.TRUE., require_resolved=opts%adaptive_mesh, & + axial_quadrature_order=opts%axial_quadrature_order, & + bending_quadrature_order=opts%bending_quadrature_order, & + require_tensile=tensile_error, & + compression_mode=CD_HC_COMPRESSION_STRICT, & + refine_non_axial=opts%adaptive_mesh, & + endconn_stiffness=endconn_k, endconn_direction=endconn_d0_local, & + endconn_mode=endconn_mode, waterline_z=CD_ZERO, & + dry_buoyancy=dry_r, water_weight=opts%rhoW*opts%g, current=cur_obj) + ELSE IF (contact_enabled) THEN + CALL CD_HermiteCable_Static_Solve_AutoMesh(l0_r, ea_r, ei_r, w_r, q_r, fixed_r, contact_z, CD_ZERO, & + 8, 80, tol_r(k_tol), 0.7_wp, CD_HC_HKAPPA_TARGET, CD_ZERO, & + CD_HC_REFINE_CAP, [2, 5], l0_am, ea_am, ei_am, w_am, & + q_static, curv, fixed_am, applied_scale, diag, res, it_r, & + es_r, em_r, contact_diameter=diam_r, contact_kbot=opts%kbot, & + contact_frame_cs=frame_cs, & + initial_polish=.TRUE., require_resolved=opts%adaptive_mesh, & + axial_quadrature_order=opts%axial_quadrature_order, & + bending_quadrature_order=opts%bending_quadrature_order, & + require_tensile=tensile_error, & + compression_mode=CD_HC_COMPRESSION_STRICT, & + refine_non_axial=opts%adaptive_mesh, & + endconn_stiffness=endconn_k, endconn_direction=endconn_d0_local, & + endconn_mode=endconn_mode, waterline_z=CD_ZERO, & + dry_buoyancy=dry_r, water_weight=opts%rhoW*opts%g, current=cur_obj) + ELSE + CALL CD_HermiteCable_Static_Solve_AutoMesh(l0_r, ea_r, ei_r, w_r, q_r, fixed_r, contact_z, CD_ZERO, & + 8, 80, tol_r(k_tol), 0.7_wp, CD_HC_HKAPPA_TARGET, CD_ZERO, & + CD_HC_REFINE_CAP, [2, 5], l0_am, ea_am, ei_am, w_am, & + q_static, curv, fixed_am, applied_scale, diag, res, it_r, & + es_r, em_r, & + initial_polish=.TRUE., require_resolved=opts%adaptive_mesh, & + axial_quadrature_order=opts%axial_quadrature_order, & + bending_quadrature_order=opts%bending_quadrature_order, & + require_tensile=tensile_error, & + compression_mode=CD_HC_COMPRESSION_STRICT, & + refine_non_axial=opts%adaptive_mesh, & + endconn_stiffness=endconn_k, endconn_direction=endconn_d0_local, & + endconn_mode=endconn_mode, waterline_z=CD_ZERO, & + dry_buoyancy=dry_r, water_weight=opts%rhoW*opts%g, current=cur_obj) + END IF + iters = iters + it_r + IF (es_r == CD_HCSTAT_OK) THEN + ! reported whenever the accepted tolerance is above the target, also when the outer + ! tolerance had already been relaxed before this polish + IF (tol_r(k_tol) > static_tol_target) tolerance_note = .TRUE. + EXIT + END IF + END DO + IF (es_r /= CD_HCSTAT_OK) THEN + WRITE (why, '(A,I0,A)') 'on the ', scale, 'x refined mesh: ' + why = TRIM(why)//' '//TRIM(em_r) + RETURN + END IF + CALL audit_refined_state(done, why) + END SUBROUTINE sequenced_refined_polish + + SUBROUTINE audit_refined_state(done, why) + !! The audit of a state refined after a failed stability test (stability_refined_route, + !! sequenced_refined_polish): the physical-branch audit and the stability test of the + !! routes, on the refined mesh. Sets stability_refined when both pass. + LOGICAL, INTENT(OUT) :: done + CHARACTER(*), INTENT(OUT) :: why + INTEGER :: es_r + CHARACTER(1024) :: em_r + LOGICAL :: stable_r + CALL audit_hermite_static(l0_am, q_static, ea_am, ei_am, w_am, diam_on(l0, l0_am, diam), endconn_mode, & + done, why, compression_note) + IF (.NOT. done) THEN + why = 'on the refined mesh: '//TRIM(why) + compression_note = '' + RETURN + END IF + IF (PRESENT(bathymetry) .AND. .NOT. bathymetry_flat) THEN + CALL hermite_static_stable(l0, l0_am, ea_am, ei_am, w_am, diam, dry_b, q_static, fixed_am, & + contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & + endconn_mode, opts, stable_r, es_r, em_r, bathymetry=bathymetry, cur_obj=cur_obj) + ELSE + CALL hermite_static_stable(l0, l0_am, ea_am, ei_am, w_am, diam, dry_b, q_static, fixed_am, & + contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & + endconn_mode, opts, stable_r, es_r, em_r, cur_obj=cur_obj) + END IF + IF (es_r /= CD_DECKDRV_OK) THEN + done = .FALSE. + why = em_r + ELSE IF (.NOT. stable_r) THEN + done = .FALSE. + WRITE (why, '(A,I0,A)') 'still unstable on the refined mesh of ', SIZE(l0_am), ' elements' + END IF + IF (done) THEN + stability_refined = .TRUE. + why = '' + ELSE + compression_note = '' + END IF + END SUBROUTINE audit_refined_state + SUBROUTINE still_water_compression(text) !! On a failed solve in a current: whether the same line without the current is !! already in axial compression (an element-mean force below -1e-3 of the peak @@ -3003,6 +3272,54 @@ LOGICAL FUNCTION symmetric_current_retry() RESULT(retry) retry = .TRUE. END FUNCTION symmetric_current_retry + SUBROUTINE saddle_refined_retry() + !! A planar state in an in-plane current that is a saddle out of its plane on a mesh + !! longer than half the bending length: polished on the refined mesh that resolves that + !! length (sequenced_refined_polish, which audits it with the planar constraints) and + !! kept when it is no saddle there; otherwise the coarse state and its notes are restored. + REAL(wp), ALLOCATABLE :: l0_c(:), ea_c(:), ei_c(:), w_c(:), q_c(:), curv_c(:) + INTEGER, ALLOCATABLE :: fixed_c(:) + TYPE(CD_HermiteResolutionType) :: diag_c + REAL(wp) :: res_c + INTEGER :: scale_c, ne_c + LOGICAL :: stab_c, done_r + CHARACTER(320) :: note_c + CHARACTER(1024) :: why_r + LOGICAL :: tol_c + l0_c = l0_am + ea_c = ea_am + ei_c = ei_am + w_c = w_am + q_c = q_static + curv_c = curv + fixed_c = fixed_am + diag_c = diag + res_c = res + scale_c = applied_scale + ne_c = stability_from_ne + stab_c = stability_refined + note_c = compression_note + tol_c = tolerance_note + CALL sequenced_refined_polish(done_r, why_r) + IF (done_r) THEN + IF (.NOT. planar_current_saddle()) RETURN + END IF + CALL MOVE_ALLOC(l0_c, l0_am) + CALL MOVE_ALLOC(ea_c, ea_am) + CALL MOVE_ALLOC(ei_c, ei_am) + CALL MOVE_ALLOC(w_c, w_am) + CALL MOVE_ALLOC(q_c, q_static) + CALL MOVE_ALLOC(curv_c, curv) + CALL MOVE_ALLOC(fixed_c, fixed_am) + diag = diag_c + res = res_c + applied_scale = scale_c + stability_from_ne = ne_c + stability_refined = stab_c + compression_note = note_c + tolerance_note = tol_c + END SUBROUTINE saddle_refined_retry + SUBROUTINE symmetric_current_polish(done) !! Pass 3 solved the line under the perturbed current: restore the exact current and !! polish from that state on its mesh, then audit it again. As on the @@ -3095,7 +3412,14 @@ SUBROUTINE symmetric_current_polish(done) CALL hermite_static_stable(l0, l0_am, ea_am, ei_am, w_am, diam, dry_b, q_static, fixed_am, & contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & endconn_mode, opts, done, es_s, em_s, cur_obj=cur_obj) - IF (es_s /= CD_DECKDRV_OK) done = .FALSE. + IF (es_s /= CD_DECKDRV_OK) THEN + done = .FALSE. + RETURN + END IF + ! A polished 3D state that fails the stability test on a mesh longer than half the + ! bending length is carried to the mesh that resolves that length and polished there, + ! as the mesh-sequenced route's state is, before the planar saddle is fallen back on. + IF (.NOT. done) CALL sequenced_refined_polish(done, em_s) END SUBROUTINE symmetric_current_polish LOGICAL FUNCTION in_plane_current() RESULT(ok) @@ -3130,6 +3454,117 @@ LOGICAL FUNCTION planar_current_saddle() RESULT(saddle) saddle = es_p == CD_DECKDRV_OK .AND. .NOT. stable_3d END FUNCTION planar_current_saddle + SUBROUTINE handoff_lateral_note(l0_ref, l0_n, ea_n, ei_n, w_n, w_rule, diam_ref, dry_ref, fixed_n, after_polish) + !! A diagnostic of the hand-off, never a change of state: a static state solved in the + !! vertical plane of the line (out-of-plane DOFs held) starts dynamics in which only the + !! end positions are held. On elements longer than half the bending length the axial + !! force at Gauss points near touchdown or a section change oscillates into compression, + !! which can make that state unstable out of its plane (a line solved in 3D is refined + !! instead, see stability_refined_route). One stability test with the dynamic + !! constraints, a single tangent assembly and banded factorisation, names it in a note. + !! A test that cannot be evaluated prints nothing. Not run where another check covers + !! it: an in-plane current (planar_current_saddle) and a coarse pass of that current. + !! after_polish: the state is the discrete equilibrium of an attachment polish, which + !! no other stability test has seen (planar_current_saddle and the 3D audit run on the + !! smeared state), so neither an in-plane current nor a line solved in 3D skips it; a + !! line solved in 3D is tested with the same end constraints. The lumped attachment drag, + !! a fixed-point load of the polish, is not part of the tangent. The arrays describe the + !! state q_static is checked on: the final mesh l0_n with its + !! element properties and weights w_n (the bare weights after an attachment polish, + !! whose lumped nodal loads are dead loads and leave the tangent unchanged), diam_ref + !! and dry_ref on the mesh l0_ref, and the static constraints fixed_n. The bending-length + !! rule of the note uses the section weights w_rule (smeared over the attachments). + REAL(wp), INTENT(IN) :: l0_ref(:), l0_n(:), ea_n(:), ei_n(:), w_n(:), w_rule(:), diam_ref(:), dry_ref(:) + INTEGER, INTENT(IN) :: fixed_n(:) + LOGICAL, INTENT(IN) :: after_polish + ! Largest factor the note names (a bending length below 1e-6 of the element length is + ! named as at least this factor). + INTEGER, PARAMETER :: SCALE_NOTE_MAX = 2**20 + INTEGER, ALLOCATABLE :: fixed_3d(:) + INTEGER :: nd, es_h, k, scale_h + LOGICAL :: stable_h + REAL(wp) :: h_max, rule, need_h, ratio_h + CHARACTER(160) :: fix_h + CHARACTER(1024) :: em_h + CHARACTER(600) :: text_h + CHARACTER(160) :: state_h, tail_h + CHARACTER(9) :: at_least_h + LOGICAL :: solved_3d + IF (saddle_kept) RETURN + IF (.NOT. after_polish .AND. (require_3d_static .OR. in_plane_current())) RETURN + IF (PRESENT(coarse_pass)) THEN + IF (coarse_pass) RETURN + END IF + IF (.NOT. ALLOCATED(q_static)) RETURN + ! Only the end translations held: the state was solved in 3D (require_3d_static, or a + ! symmetry-broken retry of a planar saddle), which the static audit covered before an + ! attachment polish but not after it. + solved_3d = require_3d_static .OR. SIZE(fixed_n) <= 6 + IF (solved_3d .AND. .NOT. after_polish) RETURN + IF (solved_3d) THEN + state_h = 'its discrete attachment equilibrium, which is unstable' + tail_h = 'a disturbance can move the line away from it' + ELSE + state_h = 'a static state solved in its vertical plane that is unstable out of that plane' + tail_h = 'a lateral disturbance can move the line out of the plane' + END IF + nd = SIZE(q_static) + fixed_3d = [1, 2, 3, nd - 5, nd - 4, nd - 3] + IF (PRESENT(bathymetry) .AND. .NOT. bathymetry_flat) THEN + CALL hermite_static_stable(l0_ref, l0_n, ea_n, ei_n, w_n, diam_ref, dry_ref, q_static, fixed_3d, & + contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & + endconn_mode, opts, stable_h, es_h, em_h, bathymetry=bathymetry, cur_obj=cur_obj) + ELSE + CALL hermite_static_stable(l0_ref, l0_n, ea_n, ei_n, w_n, diam_ref, dry_ref, q_static, fixed_3d, & + contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0_local, & + endconn_mode, opts, stable_h, es_h, em_h, cur_obj=cur_obj) + END IF + IF (es_h /= CD_DECKDRV_OK .OR. stable_h) RETURN + ! The advice names the uncapped power of two that meets the rule (the automatic + ! refinement is capped, the deck mesh is not). + need_h = bending_refine_need(l0_n, ei_n, w_rule) + scale_h = 1 + DO WHILE (REAL(scale_h, wp) < need_h .AND. scale_h < SCALE_NOTE_MAX) + scale_h = 2*scale_h + END DO + at_least_h = '' + IF (REAL(scale_h, wp) < need_h) at_least_h = 'at least' + IF (need_h > CD_ONE) THEN + ! The coarsest element against its own half bending length (the element of the + ! largest ratio, which sets need_h). + h_max = CD_ZERO + rule = HUGE(CD_ONE) + ratio_h = CD_ZERO + DO k = 1, SIZE(l0_n) + IF (ei_n(k) > CD_ZERO .AND. ABS(w_rule(k)) > CD_ZERO) THEN + IF (l0_n(k)/(0.5_wp*(ei_n(k)/ABS(w_rule(k)))**(CD_ONE/3.0_wp)) > ratio_h) THEN + rule = 0.5_wp*(ei_n(k)/ABS(w_rule(k)))**(CD_ONE/3.0_wp) + h_max = l0_n(k) + ratio_h = h_max/rule + END IF + END IF + END DO + ! AutoMesh refines by its own criteria, which need not meet the bending length: a run + ! that already sets adaptive_mesh is told to refine the deck mesh. + IF (opts%adaptive_mesh) THEN + fix_h = 'x meets that length; adaptive_mesh True did not)' + ELSE + fix_h = 'x meets that length) or set adaptive_mesh True (a line restrained in torsion is '// & + 'refined automatically)' + END IF + WRITE (text_h, '(A,I0,A,I0,A,F0.2,A,F0.2,A,I0,A)', IOSTAT=k) 'line ', ln%id, & + ' starts the dynamics from '//TRIM(state_h)//' on its mesh of ', SIZE(l0_n), & + ' elements, the coarsest ', h_max, ' m long where half the bending length (EI/|w|)^(1/3) is ', rule, & + ' m; refine the deck mesh ('//TRIM(at_least_h)//REPEAT(' ', MIN(1, LEN_TRIM(at_least_h))), scale_h, & + TRIM(fix_h) + ELSE + WRITE (text_h, '(A,I0,A,I0,A)', IOSTAT=k) 'line ', ln%id, & + ' starts the dynamics from '//TRIM(state_h)//' on its mesh of ', SIZE(l0_n), & + ' elements, which resolves the bending length; '//TRIM(tail_h) + END IF + CALL report_note(text_h) + END SUBROUTINE handoff_lateral_note + LOGICAL FUNCTION coarse_current_eligible() RESULT(ok) !! coarse_current_first applies to an in-plane current (in_plane_current) without !! seabed friction or attachments, on a deck build that is not itself the coarse one, @@ -3621,12 +4056,43 @@ SUBROUTINE attachment_static_drag(q, node, att, opts, frame_cs, f, ErrStat, ErrM END DO END SUBROUTINE attachment_static_drag + PURE INTEGER FUNCTION bending_refine_scale(l0, ei, w) RESULT(scale) + !! The power-of-two refinement factor (at most CD_HC_REFINE_CAP, and at most + !! REFINED_MAX_ELEMENTS elements after it) that brings every element of a finite-EI mesh to + !! at most half its bending length (EI/|w|)^(1/3); 1 when the mesh already meets it or the + !! caps allow no doubling. The resolution rule of the refined repetition of the static route + !! and of the stability retry. A capped factor can fall short of the rule: whether the mesh + !! meets it is bending_refine_need(...) <= 1. + REAL(wp), INTENT(IN) :: l0(:), ei(:), w(:) + ! Element cap of a refined repetition of the static route. + INTEGER, PARAMETER :: REFINED_MAX_ELEMENTS = 4096 + REAL(wp) :: need + need = bending_refine_need(l0, ei, w) + scale = 1 + DO WHILE (REAL(scale, wp) < need .AND. scale < CD_HC_REFINE_CAP .AND. SIZE(l0)*2*scale <= REFINED_MAX_ELEMENTS) + scale = 2*scale + END DO + END FUNCTION bending_refine_scale + + PURE REAL(wp) FUNCTION bending_refine_need(l0, ei, w) RESULT(need) + !! The uncapped refinement a finite-EI mesh needs for every element to be at most half its + !! bending length (EI/|w|)^(1/3): the largest ratio of element length to that length, at + !! least 1 (1 when the mesh meets the rule). + REAL(wp), INTENT(IN) :: l0(:), ei(:), w(:) + INTEGER :: k + need = CD_ONE + DO k = 1, SIZE(l0) + IF (ei(k) > CD_ZERO .AND. ABS(w(k)) > CD_ZERO) & + need = MAX(need, l0(k)/(0.5_wp*(ei(k)/ABS(w(k)))**(CD_ONE/3.0_wp))) + END DO + END FUNCTION bending_refine_need + RECURSIVE SUBROUTINE hermite_static_primary(anchor, fairlead, l0, ea, ei, w, diam, dry_b, fixed_dofs, & contact_enabled, contact_z, frame_cs, endconn_k, endconn_d0, & endconn_mode, opts, static_tol, l0_am, ea_am, ei_am, w_am, q_static, & curv, fixed_am, applied_scale, diag, res, iters, ErrStat, ErrMsg, & inadmissible, tolerance_note, kink_refined, bathymetry, refined_pass, & - cur_obj, walked, fold_walk) + cur_obj, walked, fold_walk, stability_refine) !! Primary finite-EI static route of a deck cable (the OrcaFlex/MoorPy practice of !! "catenary statics, then full statics"): !! @@ -3663,6 +4129,15 @@ RECURSIVE SUBROUTINE hermite_static_primary(anchor, fairlead, l0, ea, ei, w, dia !! refined-mesh polish of its (unresolved) solution: the whole route is then repeated !! once on a uniformly refined mesh whose elements are at most half the local bending !! length (refined_pass marks that repetition), which also sets kink_refined. + !! stability_refine (optional): run that refined repetition directly, refined from l0. The + !! caller passes it, with l0 the mesh the route returned (the deck mesh or a refinement of + !! it), when the equilibrium on that mesh failed the stability test: the pointwise axial + !! resultant of a stiff-EA element longer than half the bending length oscillates about + !! its element mean, and its compressed Gauss points lower the geometric stiffness of the + !! transverse DOFs, so an under-resolved touchdown or section change can read as a saddle. + !! A physical instability persists on the refined mesh; this artefact does not. ErrStat + !! is CD_HCSTAT_NOCONVERGE when l0 already resolves the bending length or the + !! refined route fails (ErrMsg says which). REAL(wp), INTENT(IN) :: anchor(3), fairlead(3), l0(:), ea(:), ei(:), w(:), diam(:), dry_b(:) INTEGER, INTENT(IN) :: fixed_dofs(:) LOGICAL, INTENT(IN) :: contact_enabled @@ -3685,6 +4160,8 @@ RECURSIVE SUBROUTINE hermite_static_primary(anchor, fairlead, l0, ea, ei, w, dia LOGICAL, INTENT(OUT), OPTIONAL :: walked ! Trace the bow behind the fairlead (see fold_walk above). LOGICAL, INTENT(IN), OPTIONAL :: fold_walk + ! Go straight to the refined repetition (see stability_refine above). + LOGICAL, INTENT(IN), OPTIONAL :: stability_refine ! Static tolerance of the coupled route; the floor accepted when the tighter ! dynamic-consistent tolerance is below the attainable residual. @@ -3698,8 +4175,6 @@ RECURSIVE SUBROUTINE hermite_static_primary(anchor, fairlead, l0, ea, ei, w, dia ! Continuation floor on a mesh finer than sqrt(EI/EA) (see budget_cont). INTEGER, PARAMETER :: FINE_MIN_ITERATIONS = 20 INTEGER, PARAMETER :: NFIT = 33 - ! Element cap of the refined repetition of the route. - INTEGER, PARAMETER :: REFINED_MAX_ELEMENTS = 4096 INTEGER :: nnode, ne, i, es, cont_iters, polish_iters, budget_cont, cmode, k_walk ! Anchor retreats tried for a line too long for its span (see seed_walk). INTEGER, PARAMETER :: WALK_MAX_TRIES = 6 @@ -3728,6 +4203,17 @@ RECURSIVE SUBROUTINE hermite_static_primary(anchor, fairlead, l0, ea, ei, w, dia res = CD_ZERO ne = SIZE(l0) nnode = ne + 1 + IF (PRESENT(stability_refine)) THEN + IF (stability_refine) THEN + ! refined_restart sets ErrStat only on success and appends why a refined mesh failed. + ErrStat = CD_HCSTAT_NOCONVERGE + ErrMsg = 'unstable on the returned mesh' + CALL refined_restart() + IF (ErrStat /= CD_HCSTAT_OK .AND. TRIM(ErrMsg) == 'unstable on the returned mesh') & + ErrMsg = 'unstable on the returned mesh, which already resolves the bending length' + RETURN + END IF + END IF origin = CD_ZERO IF (.NOT. PRESENT(bathymetry)) origin = anchor(1:2) @@ -4422,21 +4908,11 @@ RECURSIVE SUBROUTINE refined_restart() !! refined equilibrium is returned in the *_am outputs like an adaptive refinement. REAL(wp), ALLOCATABLE :: l0_r(:), q_r(:), q_ref_r(:), ea_r(:), ei_r(:), w_r(:), diam_r(:), dry_r(:) INTEGER, ALLOCATABLE :: fixed_r(:) - REAL(wp) :: need INTEGER :: scale, scale_need, k, es_r, parent CHARACTER(512) :: em_r LOGICAL :: inadm_r, walked_sub IF (PRESENT(refined_pass) .OR. fold) RETURN - need = CD_ONE - DO k = 1, ne - IF (ei(k) > CD_ZERO .AND. ABS(w(k)) > CD_ZERO) & - need = MAX(need, l0(k)/(0.5_wp*(ei(k)/ABS(w(k)))**(CD_ONE/3.0_wp))) - END DO - scale_need = 1 - DO WHILE (REAL(scale_need, wp) < need .AND. scale_need < CD_HC_REFINE_CAP .AND. & - ne*2*scale_need <= REFINED_MAX_ELEMENTS) - scale_need = 2*scale_need - END DO + scale_need = bending_refine_scale(l0, ei, w) IF (scale_need < 2) RETURN ! The needed factor first; when its continuation exhausts its element-scaled iteration ! budget, the coarser powers of two down to 2x (a halved mesh, still resolving the diff --git a/tests/data/lazy_wave_touchdown_planar_coarse.dat b/tests/data/lazy_wave_touchdown_planar_coarse.dat new file mode 100644 index 0000000..79068b9 --- /dev/null +++ b/tests/data/lazy_wave_touchdown_planar_coarse.dat @@ -0,0 +1,44 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable free in torsion, deck mesh too coarse at touchdown +-- The cable of lazy_wave_touchdown_torsion_coarse.dat free in torsion, so its static solve runs in +-- the vertical plane of the line. On the coarse deck mesh (10 m elements along the grounded run, +-- half the bending length is 2.5 m) that planar state is unstable out of its plane, and the dynamics +-- would start from it with the line free to leave the plane: the driver names it in a note, once, +-- and solves the deck unchanged. +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +buoy 0.29 59.53 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 20 +4 buoy 100.0 25 +4 bare 420.0 42 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Free 0 1 0 0 +4 B Rigid 1 0 0 Free 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.025 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +--------------------- END ---------------------------------------------------- diff --git a/tests/data/lazy_wave_touchdown_planar_coarse_att.dat b/tests/data/lazy_wave_touchdown_planar_coarse_att.dat new file mode 100644 index 0000000..1330ddc --- /dev/null +++ b/tests/data/lazy_wave_touchdown_planar_coarse_att.dat @@ -0,0 +1,48 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable with discrete buoyancy modules, deck mesh too coarse at touchdown +-- The cable of lazy_wave_touchdown_planar_coarse.dat with its 100 m buoyant stretch built from the +-- bare cable plus 25 discrete buoyancy modules at a 4 m pitch, one on each node from 80 m to 176 m +-- (each carries the extra mass and displaced volume of 4 m of the buoy section). Its static solve +-- runs in the vertical plane of the line and ends with the attachment polish. On the coarse deck +-- mesh (10 m elements along the grounded run) the discrete planar state is unstable out of its +-- plane: the driver names it in a note, once, after the polish, and solves the deck unchanged. +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 20 +4 bare 100.0 25 +4 bare 420.0 42 +--------------------- ATTACHMENTS ------------------------------------------- +LineID ArcLength Mass Volume CdA Ca CdAx +(-) (m) (kg) (m^3) (m^2) (-) (m^2) +4 80.0:4.0:176.0 91.32 0.18378 1.392 1.0 0.0 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Free 0 1 0 0 +4 B Rigid 1 0 0 Free 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.025 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +--------------------- END ---------------------------------------------------- diff --git a/tests/data/lazy_wave_touchdown_planar_refined.dat b/tests/data/lazy_wave_touchdown_planar_refined.dat new file mode 100644 index 0000000..c7007f0 --- /dev/null +++ b/tests/data/lazy_wave_touchdown_planar_refined.dat @@ -0,0 +1,41 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable free in torsion, mesh resolving the bending length +-- lazy_wave_touchdown_planar_coarse.dat with every element at most 2.5 m, half the bending +-- length of the bare cable: the planar state is stable out of its plane and no note is printed. +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +buoy 0.29 59.53 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 40 +4 buoy 100.0 50 +4 bare 420.0 168 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Free 0 1 0 0 +4 B Rigid 1 0 0 Free 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.025 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +--------------------- END ---------------------------------------------------- diff --git a/tests/data/lazy_wave_touchdown_planar_refined_att.dat b/tests/data/lazy_wave_touchdown_planar_refined_att.dat new file mode 100644 index 0000000..062bb09 --- /dev/null +++ b/tests/data/lazy_wave_touchdown_planar_refined_att.dat @@ -0,0 +1,45 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable with discrete buoyancy modules, mesh resolving the bending length +-- The cable of lazy_wave_touchdown_planar_coarse_att.dat (25 discrete buoyancy modules at a 4 m +-- pitch) on 2 m, 2 m and 2.5 m elements, each at most half the bending length: its discrete planar +-- state is stable out of its plane and the hand-off prints no note. +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 40 +4 bare 100.0 50 +4 bare 420.0 168 +--------------------- ATTACHMENTS ------------------------------------------- +LineID ArcLength Mass Volume CdA Ca CdAx +(-) (m) (kg) (m^3) (m^2) (-) (m^2) +4 80.0:4.0:176.0 91.32 0.18378 1.392 1.0 0.0 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Free 0 1 0 0 +4 B Rigid 1 0 0 Free 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.025 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +--------------------- END ---------------------------------------------------- diff --git a/tests/data/lazy_wave_touchdown_torsion_coarse.dat b/tests/data/lazy_wave_touchdown_torsion_coarse.dat new file mode 100644 index 0000000..c20db89 --- /dev/null +++ b/tests/data/lazy_wave_touchdown_torsion_coarse.dat @@ -0,0 +1,47 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable restrained in torsion, deck mesh too coarse at touchdown +-- A finite-EI lazy-wave cable clamped at both ends and restrained in torsion (Rigid, zero +-- pretwist), so its static solve runs in 3D. The deck mesh is 4 m up to the end of the +-- buoyancy section and 10 m along the grounded run, longer than half the bending length +-- (EI/w)^(1/3) = 5.0 m of the bare cable. On that mesh the pointwise axial resultant of the +-- under-resolved touchdown elements oscillates into compression at Gauss points, and the +-- tangent of the stable equilibrium shows a spurious out-of-plane negative eigenvalue there. +-- The static solve repeats the route on the refined mesh, which is stable, and reports it once. +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +buoy 0.29 59.53 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 20 +4 buoy 100.0 25 +4 bare 420.0 42 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Rigid 0 1 0 0 +4 B Rigid 1 0 0 Rigid 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.025 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +Torq4N1 +--------------------- END ---------------------------------------------------- diff --git a/tests/data/lazy_wave_touchdown_torsion_coarse_sequenced.dat b/tests/data/lazy_wave_touchdown_torsion_coarse_sequenced.dat new file mode 100644 index 0000000..abba947 --- /dev/null +++ b/tests/data/lazy_wave_touchdown_torsion_coarse_sequenced.dat @@ -0,0 +1,44 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable restrained in torsion, deck mesh too coarse at touchdown +-- lazy_wave_touchdown_torsion_coarse.dat with the mesh-sequenced static route first +-- (OPTIONS cable_statics sequenced): its deck-mesh state fails the stability test the same way, +-- and is prolongated to the mesh that resolves the bending length and polished there. +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +buoy 0.29 59.53 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 20 +4 buoy 100.0 25 +4 bare 420.0 42 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Rigid 0 1 0 0 +4 B Rigid 1 0 0 Rigid 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.025 dtM +0.0 TMax +sequenced cable_statics +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +Torq4N1 +--------------------- END ---------------------------------------------------- diff --git a/tests/data/lozon_coupled_harmonics.dat b/tests/data/lozon_coupled_harmonics.dat index 627d49d..7d5738e 100644 --- a/tests/data/lozon_coupled_harmonics.dat +++ b/tests/data/lozon_coupled_harmonics.dat @@ -1,5 +1,5 @@ # L4-1c hang-off harmonic table (cos convention: A cos(2 pi t/T + phase_deg)) -# source: CableDyn coupled OpenFAST run of the IEA-15-240-RWT UMaine mixed example (Main_mixed_600.out); +# source: CableDyn coupled OpenFAST run of the IEA-15-240-RWT UMaine mixed example, 87-element cable mesh (Main_mixed_600.out); # offset (30.0,-20.0); window 232.4s; period<=40s; capture 85.4% 16 4.131042 0.131826 # period_s surge_amp surge_phase_deg heave_amp heave_phase_deg diff --git a/tests/test_cable_statics_matrix.py b/tests/test_cable_statics_matrix.py index 7ebbb1b..769be72 100644 --- a/tests/test_cable_statics_matrix.py +++ b/tests/test_cable_statics_matrix.py @@ -59,9 +59,13 @@ from pathlib import Path # The per-deck bound is on the solver's CPU time, which does not grow with machine load; the -# wall-clock limit only stops a hung process and is scaled far above any loaded-machine run. +# wall-clock limit only stops a hung process and is scaled far above any loaded-machine run: +# at least DECK_WALL_LIMIT and WALL_PER_CPU times the configuration's CPU bound (--cpu-limit), +# so a Debug build, whose CPU bound is ten times larger, is not stopped as hung while it is +# still inside that bound. DECK_CPU_LIMIT = 60.0 DECK_WALL_LIMIT = 600.0 +WALL_PER_CPU = 4.0 def _windows_cpu_seconds(handle) -> float: @@ -190,7 +194,7 @@ def self_intersects(xyz: list) -> bool: return False -def run_case(exe: str, root: Path, case: dict) -> dict: +def run_case(exe: str, root: Path, case: dict, wall_limit: float = DECK_WALL_LIMIT) -> dict: cdir = root / case["id"] cdir.mkdir(parents=True, exist_ok=True) if case.get("bathy"): @@ -201,7 +205,7 @@ def run_case(exe: str, root: Path, case: dict) -> dict: env.setdefault("OMP_NUM_THREADS", "1") t0 = time.time() rc, log, cpu = run_measured([exe, str(cdir / "deck.dat"), str(cdir / "out")], str(cdir), env, - DECK_WALL_LIMIT) + wall_limit) return {"id": case["id"], "rc": rc, "log": log, "time": time.time() - t0, "cpu": cpu, "dir": cdir} @@ -314,16 +318,17 @@ def main() -> int: cases = json.loads(Path(a.cases).read_text(encoding="utf-8")) root = Path(a.output) root.mkdir(parents=True, exist_ok=True) + wall_limit = max(DECK_WALL_LIMIT, WALL_PER_CPU * a.cpu_limit) t0 = time.time() with ThreadPoolExecutor(a.workers) as ex: - results = list(ex.map(lambda c: run_case(a.exe, root, c), cases)) + results = list(ex.map(lambda c: run_case(a.exe, root, c, wall_limit), cases)) byid = {r["id"]: r for r in results} failures = 0 for case in cases: res = byid[case["id"]] errs = check_case(case, res) if res["rc"] == "TIMEOUT": - errs.append("no exit within the %.0f s wall-clock hang guard" % DECK_WALL_LIMIT) + errs.append("no exit within the %.0f s wall-clock hang guard" % wall_limit) # a deck's own bound (a fast failure) scales with the configuration's bound elif not res["cpu"] <= case.get("cpu_limit", DECK_CPU_LIMIT) * a.cpu_limit / DECK_CPU_LIMIT: errs.append("CPU time %.1f s beyond the %.0f s per-deck bound" % ( diff --git a/tests/test_deck_hermite_cable.f90 b/tests/test_deck_hermite_cable.f90 index 8d6e6a7..6ab5778 100644 --- a/tests/test_deck_hermite_cable.f90 +++ b/tests/test_deck_hermite_cable.f90 @@ -553,14 +553,14 @@ SUBROUTINE case_coupled_seabed_contact() END SUBROUTINE case_coupled_seabed_contact SUBROUTINE case_grounded_anchor_end() - !! A lazy-wave power cable laid on a frictional flat bed up to its pinned anchor on the - !! floor (Line 4 of examples/iea15mw_umaine_mixed_cabledyn.dat: 200 m depth, kBot 3e6, - !! frictionMu 0.35, 10 m grounded elements). The seabed carries the weight of the anchor + !! A lazy-wave power cable laid on a frictional flat bed up to its pinned anchor on the floor + !! (Line 4 of examples/iea15mw_umaine_mixed_cabledyn.dat on its former mesh: 200 m depth, kBot + !! 3e6, frictionMu 0.35, 10 m grounded elements). The seabed carries the weight of the anchor !! element, so its free end tangent stays along the floor; turned into the floor by the - !! unbalanced end moment of the weight, the element dipped about 27 mm below the bed, - !! carried axial forces from -47 kN to +127 kN about the 1.2 kN line tension, and - !! buckled sideways once the hang-off swayed, after which the friction locked the anchor - !! tension at about 2.5 times its static value while the line slackened. + !! unbalanced end moment of the weight, the element dipped about 27 mm below the bed, carried + !! axial forces from -47 kN to +127 kN about the 1.2 kN line tension, and buckled sideways once + !! the hang-off swayed, after which the friction locked the anchor tension at about 2.5 times + !! its static value while the line slackened. TYPE(CD_HFMF_ModuleType), ALLOCATABLE :: cable INTEGER :: cn, es, base, k CHARACTER(256) :: em diff --git a/tests/test_openfast_aggregate.f90 b/tests/test_openfast_aggregate.f90 index b02aaba..fe58f15 100644 --- a/tests/test_openfast_aggregate.f90 +++ b/tests/test_openfast_aggregate.f90 @@ -4253,11 +4253,11 @@ SUBROUTINE write_lozon_semitaut_deck(path) END SUBROUTINE write_lozon_semitaut_deck SUBROUTINE case_touchdown_cable_end_tension() - !! The IEA-15MW touchdown power cable of the mixed example (line 2 beside one mooring - !! chain): bare / buoyant / bare sections, stiff EA = 4.69e8 N on 4-10 m elements, - !! ending on a 3e6 Pa/m seabed. Statics carry no horizontal load, so the anchor end sees - !! the grounded run's tension H (~1.2 kN). The node-sampled stretch EA*(|dr/ds| - 1) at - !! the held anchor read 78 kN; the element end force equals H. + !! The IEA-15MW touchdown power cable of the mixed example on its former mesh (line 2 beside one + !! mooring chain): bare / buoyant / bare sections, stiff EA = 4.69e8 N on 4-10 m elements, + !! ending on a 3e6 Pa/m seabed. Statics carry no horizontal load, so the anchor end sees the + !! grounded run's tension H (~1.2 kN). The node-sampled stretch EA*(|dr/ds| - 1) at the held + !! anchor read 78 kN; the element end force equals H. TYPE(CD_AGG_ModuleType) :: agg INTEGER :: u, ios, es, k REAL(wp) :: v(3) diff --git a/validation/OPENFAST_DLC_VERIFICATION.md b/validation/OPENFAST_DLC_VERIFICATION.md index e04cd1e..8dab5f4 100644 --- a/validation/OPENFAST_DLC_VERIFICATION.md +++ b/validation/OPENFAST_DLC_VERIFICATION.md @@ -169,14 +169,17 @@ Tensions and platform dynamic response agree to within 0.6 %, apart from the Ptf ### Finite-EI power cable, severe sea (Hs=11) The cases above use the all-chain mooring. This case uses the mixed deck -`examples/iea15mw_umaine_mixed_cabledyn.dat`: three chains plus an 87-element (20 + 25 + 42) -finite-EI cubic-Hermite lazy-wave cable that runs from the platform (point 7, `Coupled`) to a -touchdown on the 200 m seabed and ends at an anchor (point 8, `Fixed`), in the full turbine under -Hs=11 waves. The run completed (exit 0, all 2400 steps, 0 NaN/Inf), and its output contains the -finite-EI curvature and bending-moment channels `Curv4N20` and `BendMom4N20`, which confirms that -the cable ran on the finite-EI element. At the static initial state the standalone driver gives -`Curv4N20` = 0.0915 1/m and `BendMom4N20` = 1821 N·m on this deck. This case has no MoorDyn-F -twin and establishes completion only; cable accuracy is established against OrcaFlex in cases L3-4, L3-5, L3-6 and +`examples/iea15mw_umaine_mixed_cabledyn.dat`: three chains plus a finite-EI cubic-Hermite +lazy-wave cable that runs from the platform (point 7, `Coupled`) to a touchdown on the 200 m +seabed and ends at an anchor (point 8, `Fixed`), in the full turbine under Hs=11 waves. The run +used the deck's earlier 87-element (20 + 25 + 42) mesh of the cable; the example now meshes it +with 248 elements to resolve its bending length, and the one-hour full-turbine run of the +coupled finite-EI touchdown row in [`VALIDATION.md`](../VALIDATION.md) uses that mesh. The run +completed (exit 0, all 2400 steps, 0 NaN/Inf), and its output contains the finite-EI curvature +and bending-moment channels of the sag bend, which confirms that the cable ran on the finite-EI +element. At the static initial state the standalone driver gives 0.0915 1/m and 1821 N·m on +that mesh (`Curv4N20`, `BendMom4N20`), and 0.0906 1/m and 1803 N·m on the current one +(`Curv4N38`, `BendMom4N38`). This case has no MoorDyn-F twin and establishes completion only; cable accuracy is established against OrcaFlex in cases L3-4, L3-5, L3-6 and L4-1c of [`VALIDATION.md`](../VALIDATION.md). ### Coupled rod between two fairleads, still water