diff --git a/CHANGELOG.md b/CHANGELOG.md index 23ff0d0..83b7465 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,109 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ## [Unreleased] -(nothing yet) +### Fixed + +- **Dense symmetric solves depended on how many modes were requested.** + The dense path reduced `K x = λ M x` through a Cholesky factor of the + mass matrix, and a backward-stable symmetric eigensolver gets each + eigenvalue of the reduced matrix only to an *absolute* accuracy of + about `eps · λ_max`. The ElastoDyn adapters floor rotary inertia at + 1e-6, which puts `λ_max / λ_min` near 1e15, so the lowest modes carried + errors of several percent that moved with the `subset_by_index` + window. On the bundled NREL 5MW blade at 12.1 rpm the first flap mode + read **0.7154 / 0.7096 / 0.7661 / 0.7282 Hz** for `n_modes` = 4 / 6 / + 10 / 20, and first edge 1.1221 / 1.1184 / 1.0884 / 1.1303 Hz; the tower + fore-aft / side-side pair moved and swapped. Every model below the + 500-DOF sparse threshold took this path — every ElastoDyn blade and + tower. + + The dense path now solves the inverted pencil `M x = μ (K + s M) x`, + `μ = 1 / (λ + s)`, which puts the wanted modes at the top of the + spectrum where the same absolute accuracy is a relative one — the + shift-invert transformation the sparse path already used. The shift is + zero when `K` is positive definite and is otherwise grown from an + estimate of the lowest eigenvalues until `K + s M` factorises, so + rigid-body modes and the negative eigenvalues of a column loaded past + buckling are still returned. A request reaching into the upper half of + the spectrum takes those modes from the mass-reduced form, which is the + accurate one there. Results are now independent of `n_modes` and agree + with the sparse path to roundoff: blade flap1 / edge1 / flap2 = + **0.72865 / 1.09754 / 2.00842 Hz** at every window, tower 1st pair + 0.33483 / 0.33561 Hz, backward error ~1e-10 (was ~5e-2). Frequencies of + every dense-path model change by up to the old error — on the NREL 5MW + blade up to ~5 % on the lowest flap / edge modes, on the tower ~0.3 %. + + When one strongly negative eigenvalue forces a shift much larger than + the soft end of the spectrum, those soft modes all map to the same `μ` + to working precision. The window is then grown until its edge falls in + a resolvable gap, and each such cluster is re-solved by Rayleigh-Ritz + on its own subspace, so the soft modes come back in order and + independent of `n_modes` (`K = diag(−1e16, 1, 2, 3, 4)` returns + `[−1e16, 1]` for two modes, not `[−1e16, 4]`). + + The mass-reduced form is kept as the fallback for a pencil no shift + makes definite, and the residual retry of 1.18.0 still stands behind + it; the near-singular-mass case it was built for now solves directly + on the dense path without a retry or a warning. + `tests/test_campbell.py::test_first_flap_lifts_endpoint_to_endpoint`, + which sat at 2.95 % against its 3 % gate on some BLAS builds, now + measures 7.6 %. +- **Floating models built from OpenFAST decks lost their roll and pitch + modes.** `Tower.from_elastodyn_with_mooring` took the restoring from + the WAMIT `.hst`, which holds only the hydrostatic part + `ρ g I_wp + ρ g V z_B`. OpenFAST adds the body weight through + ElastoDyn's gravity; the model did not, so a ballast-stabilised spar + read as unstable — OC3 Hywind `C44 = C55 = −5.0e9 N·m/rad` — and the + general eigensolver silently dropped the two negative-eigenvalue modes, + returning four of the six rigid-body modes. The constructor now adds the + weight restoring `−g Σ mᵢ (zᵢ − z_ref)` of the platform, the tower and + the RNA to `hydro_K` (+6.17e9 N·m/rad on OC3: +6.58e9 of ballast, + −0.41e9 of tower and RNA). It is added only when a HydroDyn deck is + given, since weight without the buoyancy that balances it does not + describe a floating body. `Tower.from_windio_floating` adds the same + term in both of its tiers, so the deck-backed tier stays equivalent to + `from_elastodyn_with_mooring`. OC3 Hywind from the r-test decks now + gives surge / sway / heave / roll / pitch / yaw = 0.00806 / 0.00806 / + 0.03241 / **0.03383 / 0.03383** / 0.04121 Hz against Jonkman (2010) + 0.0080 / 0.0080 / 0.0324 / 0.0342 / 0.0343 Hz (roll and pitch were + absent before). Yaw is the catenary-only value; the published 0.1210 Hz + includes an additional 9.834e7 N·m/rad crowfoot spring the r-test + MoorDyn deck does not carry. The coupled first tower pair moves from + 0.4887 / 0.4903 to 0.4922 / 0.4938 Hz. +- **A negative-stiffness rigid-body mode is no longer dropped silently.** + `Tower.run` on a free-free floating model now warns (`UserWarning`) + when `hydro_K + mooring_K` has a negative direction, naming the + dominant DOFs: the platform is statically unstable there and the + general solver leaves those modes out of the result. The bundled + `08_nrel5mw_oc4semi` sample triggers it — the OC4 DeepCwind deck keeps + its ballast in HydroDyn filled members rather than `PtfmMass`, which + the platform model does not include, so its roll / pitch restoring is + negative and the modes labelled roll / pitch there are tower bending. +- **`Tower.add_point_mass` on a deck-built floater now updates the + weight restoring.** The weight term is formed when the model is + constructed, so a lump added afterwards reached the FEM mass matrix but + not the roll / pitch restoring. The lump's `−m g (z − z_ref)` is now + added to `hydro_K` as well, which matches a model rebuilt with the lump + in place. +- **NumPy 1.26 support.** Several code paths called `np.trapezoid`, + which only exists from NumPy 2.0, so on the supported `numpy>=1.26` + every `Tower.run()` (through the pre-solve checks), every ElastoDyn + blade read and every floating constructor failed with + `AttributeError`. They now share one version-independent trapezoid + helper, and a test keeps the NumPy-2-only spelling out of the package. +- **The first tower pair of the IEA-3.4 and IFE UPSCALE 25MW decks is no + longer degenerate.** Both towers are symmetric in stiffness, but the + lumped RNA is not (pitch inertia 1.9× and 3.0× roll inertia, plus an + `ixz` product), which splits the first FA / SS pair by 0.09 % and 1.0 %. + The dense solver above, the sparse path and an `eigh(M, K)` solve all + agree on the split. The old dense path returned the pair degenerate + only because its error on the lowest modes was larger than the split, + and two integration tests had pinned that; they now pin the split. +- **WindIO discovery ignored every yaml under a path containing + "OpenFAST".** The exclusion meant for the `OpenFAST/` deck directory of + an RWT layout matched the absolute path, so a checkout living under + e.g. `~/OpenFAST-GUI/` found no ontology at all. Only the part of the + path below the search root is examined now. ## [1.18.0] — 2026-08-13 diff --git a/VALIDATION.md b/VALIDATION.md index ff2a0c9..2d36ce7 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -102,7 +102,9 @@ metrics: | Discrete mid-span point mass (issue #35) | cantilever carrying one lump at station $a$: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1) | 1st frequency; mesh-position independence | < 0.5 %; coarse-vs-fine < 0.2 % | (within tol) | [`tests/fem/test_gravity_and_point_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_gravity_and_point_mass.py) | no | | Distributed Winkler soil bed vs the lumped mudline condensation (issue #118) | Psaroudakis et al. (2021) / Yu & Amdahl (2023) Eq. 25 is the exact static condensation of a constant-EI pile on a bed of rate $k = D_P E_{SO}$ | coupled 1st frequency, distributed bed vs lumped springs | < 1 % | 0.4 % (the embedded pile inertia the condensed form drops) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | Distributed Winkler bed converges on the rigid mudline clamp (issue #118) | limit $E_{SO} \to \infty$; residual compliance scales as the elastic length $(4EI/k)^{1/4}$ | 1st frequency at $E_{SO} \times 10^8$ vs the clamped model | < 1 % | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | -| Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes | < 0.5 % | (within tol; the unguarded symmetric solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | +| Dense symmetric solve is independent of the requested mode count | NREL 5MW ElastoDyn blade (12.1 rpm) and tower, bundled deck; shift-invert spectral transformation (Ericsson & Ruhe 1980) | lowest modes for 4 vs 20 requested; dense vs sparse shift-invert | `rtol = 1e-6` | identical to roundoff (blade flap1 0.72865 Hz at every window; was 0.715 / 0.710 / 0.766 / 0.728 Hz for 4 / 6 / 10 / 20) | [`tests/fem/test_dense_inverted_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_dense_inverted_solver.py) | no | +| Near-singular mass solved directly on the dense path | cantilever carrying one lump: $f = \sqrt{3EI/(m L^3)}/2\pi$ (Blevins 1979, Table 8-1), 4e5:1 lump-to-beam mass ratio | 1st frequency, no residual retry, backward error | < 0.5 %; residual < 1e-6 | (within tol) | [`tests/fem/test_dense_inverted_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_dense_inverted_solver.py) | no | +| Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes, through the residual retry behind the mass-reduced fallback | < 0.5 % | (within tol; the unguarded mass-reduced solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | | Rigid-body modes are not lost to the solver guard | construction (rank-deficient $K$, well-conditioned $M$; requested subsets of 1, 3, 6 and 10 modes) | zero modes retained, requested mode count returned | exact | (within tol) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | | Deck-reader `n_nodes` refinement (issue #58) | Euler-Bernoulli closed form, modes 1-3, plus self-convergence $n{=}100$ vs $200$ | bending frequencies | < 1 %; self-convergence < 0.2 % | (within tol) | [`tests/test_refine_mesh.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_refine_mesh.py) | no | @@ -114,6 +116,8 @@ by title + publisher: - Blevins (1979). *Formulas for Natural Frequency and Mode Shape*. Krieger Publishing. (Textbook; no DOI.) - Karnovsky & Lebed (2001). *Formulas for Structural Dynamics*. McGraw-Hill. (Textbook; no DOI.) - Timoshenko & Gere (1961). *Theory of Elastic Stability*, 2nd ed. McGraw-Hill. (Textbook; no DOI.) +- Ericsson & Ruhe (1980). *The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems*. *Mathematics of Computation* 35(152), 1251–1268. DOI [10.1090/S0025-5718-1980-0583502-2](https://doi.org/10.1090/S0025-5718-1980-0583502-2). +- Faltinsen (1990). *Sea Loads on Ships and Offshore Structures*. Cambridge University Press. (Textbook; no DOI.) - Greenhill (1881). *Determination of the greatest height consistent with stability that a vertical pole or mast can be made*. *Proceedings of the Cambridge Philosophical Society* 4, 65–73. (No DOI; tabulated in Timoshenko & Gere 1961, Art. 2.16.) - Shadlou & Bhattacharya (2016). *Dynamic stiffness of monopiles supporting offshore wind turbine generators*. *Soil Dynamics and Earthquake Engineering* 88, 15–32. DOI [10.1016/j.soildyn.2016.04.002](https://doi.org/10.1016/j.soildyn.2016.04.002). - Yu & Amdahl (2023). *A Rayleigh-Ritz solution for high order natural frequencies and eigenmodes of monopile supported offshore wind turbines considering tapered towers and soil-pile interactions*. *Marine Structures* 92, 103482. DOI [10.1016/j.marstruc.2023.103482](https://doi.org/10.1016/j.marstruc.2023.103482). @@ -223,6 +227,9 @@ is behavioural / contract-style. | `MooringSystem.stiffness_matrix` — OC3 yaw stiffness (catenary-only) | OC3 r-test MoorDyn `.dat` | `K[5,5] ≈ 1.156e7 N·m/rad` (catenary lines only; OC3's delta-line crowfoot adds the bulk of the published 9.83e7) | `rel = 5%` | (within tol) | [`tests/test_mooring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_mooring.py) | yes | | `MooringSystem.stiffness_matrix` — 3-fold-symmetry signatures | OC3 r-test MoorDyn `.dat` | K[0,0] = K[1,1], K[3,3] = K[4,4], K[0,4]·K[1,3] < 0 | `rel = 1e-4` | (within tol) | [`tests/test_mooring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_mooring.py) | yes | | `Tower.from_elastodyn_with_mooring` — end-to-end OC3 coupled solve | OC3 ElastoDyn + r-test MoorDyn `.dat` | hub\_conn=2, PlatformSupport populated, 1st tower-bending in 0.4–0.6 Hz | `rel = 5%` | within 1.2 % of Jonkman 2010's 0.482 Hz | [`tests/test_mooring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_mooring.py) | yes | +| Floating weight restoring $-g \sum m_i (z_i - z_{ref})$ in roll / pitch (platform + tower + RNA) | closed form; classical $C_{44} = \rho g V (z_B - z_G) + \rho g I_{wp}$ (Faltinsen 1990, ch. 2) | 6×6 weight matrix, reference-point shift, mass scaling, point masses | `rel = 1e-9` | exact | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | no | +| Negative platform restoring is reported, not silently dropped | construction (OC3 Hywind sample with the weight term removed) | `UserWarning` naming roll / pitch | always | (within tol) | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | no | +| **`Tower.from_elastodyn_with_mooring` — OC3 Hywind rigid-body frequencies** (WAMIT `.hst` + weight of platform, tower and RNA) | Jonkman (2010) NREL/TP-500-47535, FAST full-system natural frequencies: surge/sway 0.0080, heave 0.0324, roll 0.0342, pitch 0.0343 Hz | six labelled rigid-body modes | surge/sway/roll/pitch 3 %, heave 2 % | surge/sway 0.8 %, heave 0.0 %, roll 1.1 %, pitch 1.4 % (roll/pitch were missing before the weight term) | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | yes | | `MudlineFoundation.from_soil_properties` — Yu & Amdahl (2023) Table 9 DTU 10 MW flexible pile, E\_SO = 30 MPa | Yu and Amdahl (2023) Marine Structures 92, 103482, Table 9 row 2 | K\_hh, K\_hr, K\_rr via Shadlou and Bhattacharya (2016) formulas | `rel = 10%` | within 3.0 % on all three springs | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | `MudlineFoundation.from_soil_properties` — Yu & Amdahl (2023) Table 9 DTU 10 MW rigid pile, E\_SO = 30 MPa | Yu and Amdahl (2023) Table 9 row 1 | K\_hh, K\_hr, K\_rr via Shadlou and Bhattacharya (2016) rigid-pile formulas | `rel = 5%` | within 2.9 % on all three springs | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | `MudlineFoundation.as_mooring_K` — cross-coupling sign matches Jonkman 2010 OC3 convention | construction | K\[0,4\] < 0, K\[1,3\] > 0, K\[0,4\]·K\[1,3\] < 0 | always | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | diff --git a/docs/conventions.rst b/docs/conventions.rst index d61fdc2..af9345b 100644 --- a/docs/conventions.rst +++ b/docs/conventions.rst @@ -289,6 +289,17 @@ Common pitfalls * **Leaving ``mooring_K`` in the platform/hull frame.** Like hydro, it must be referenced to the tower axis. +* **Passing a WAMIT ``.hst`` as the whole restoring.** The ``.hst`` is + hydrostatic only (``ρ g I_wp + ρ g V z_B``); the body weight's + ``−m g z_G`` in roll / pitch is applied by OpenFAST through ElastoDyn + gravity. ``hydro_K + mooring_K`` must carry both. The deck + constructors add the weight of platform, tower and RNA — each mass + at its height above the ``ref_msl`` reference — for you; a + hand-built ``PlatformSupport`` has to include it itself. Without it + a ballast-stabilised spar is negative in roll / pitch (OC3 Hywind: + ``−5.0e9 N·m/rad``), and ``Tower.run`` warns that those rigid-body + modes are unstable rather than returning them. + * **Confusing OpenFAST's two frames.** ``PtfmCMzt`` (CM) and ``PtfmRefzt`` (reference) are both in the tower-base *t*-frame; pyBmodes' ``ref_msl`` is ``PtfmRefzt``. Don't add one to the other. diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/build.py b/src/pybmodes/_examples/sample_inputs/reference_turbines/build.py index 7133353..d233d10 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/build.py +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/build.py @@ -89,6 +89,7 @@ if SRC_DIR.is_dir() and str(SRC_DIR) not in sys.path: sys.path.insert(0, str(SRC_DIR)) +from pybmodes._numeric import trapezoid # noqa: E402 from pybmodes.io.bmi import PlatformSupport # noqa: E402 from pybmodes.io.elastodyn_reader import ( # noqa: E402 read_elastodyn_blade, @@ -172,7 +173,7 @@ def _blade_total_mass(blade) -> float: span = np.asarray(blade.bl_fract, dtype=float) rho = np.asarray(blade.b_mass_den, dtype=float) # bl_fract is normalized [0, 1]; multiply by blade length elsewhere. - return float(np.trapezoid(rho, span)) + return trapezoid(rho, span) def _emit_section_properties_table( diff --git a/src/pybmodes/_numeric.py b/src/pybmodes/_numeric.py new file mode 100644 index 0000000..6762dc5 --- /dev/null +++ b/src/pybmodes/_numeric.py @@ -0,0 +1,35 @@ +# Copyright 2024-2026 Jae Hoon Seo +# Marine Structural Mechanics and Integrity Lab (SMI Lab), Inha University +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Small numerical helpers shared across the package.""" + +from __future__ import annotations + +import numpy as np +from numpy.typing import ArrayLike + +__all__ = ["trapezoid"] + + +def trapezoid(y: ArrayLike, x: ArrayLike) -> float: + """Trapezoidal integral of ``y`` over ``x``. + + Implemented directly rather than through ``np.trapezoid``, which only + exists from NumPy 2.0, or ``np.trapz``, which NumPy 2 deprecates, so + the package works across the whole advertised ``numpy>=1.26`` range. + """ + y = np.asarray(y, dtype=float) + x = np.asarray(x, dtype=float) + return float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(x))) diff --git a/src/pybmodes/checks.py b/src/pybmodes/checks.py index 376a2c9..a588f6a 100644 --- a/src/pybmodes/checks.py +++ b/src/pybmodes/checks.py @@ -102,6 +102,7 @@ import numpy as np +from pybmodes._numeric import trapezoid from pybmodes.options import DEFAULT_CHECK_OPTIONS as _CHECK_OPTIONS if TYPE_CHECKING: @@ -269,7 +270,7 @@ def _check_section_properties_finite( ) -> None: """Flag any non-finite (NaN / ±Inf) entry in the numeric section- property fields. ERROR-severity because every downstream - consumer (``np.trapezoid``, ``np.linalg.eigh``, the FE assembly) + consumer (the trapezoid mass integral, ``np.linalg.eigh``, the FE assembly) silently produces NaN-filled outputs on NaN inputs. The per-field checks below this one (``_check_span_monotonic`` / @@ -395,7 +396,7 @@ def _check_rna_vs_tower_mass( mass_den = np.asarray(sp.mass_den, dtype=float) if span_phys.size < 2 or not np.all(np.diff(span_phys) > 0): return # other checks will flag the bad span; skip cleanly here - tower_mass = float(np.trapezoid(mass_den, span_phys)) + tower_mass = trapezoid(mass_den, span_phys) if tower_mass <= 0.0: return # mass-density check will catch this rna = float(bmi.tip_mass.mass) diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 8e6a582..0815a7c 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -28,8 +28,11 @@ dense LAPACK solve for the few-lowest-modes case on a 500+ DOF tower mesh. 2. **Dense symmetric** — ``scipy.linalg.eigh`` on the symmetrised - matrices. Path used for small / mid-size symmetric problems and - when the sparse path fails to converge (logged as a warning). + matrices, applied to the *inverted* pencil ``M x = μ (K + s M) x`` + so that the lowest modes are resolved to relative rather than + absolute accuracy (see :func:`_solve_dense_inverted`). Path used for + small / mid-size symmetric problems and when the sparse path fails + to converge (logged as a warning). 3. **Dense general** — ``scipy.linalg.eig`` for genuinely asymmetric systems (offshore decks where the rigid-arm transformation makes the platform-support block non-symmetric). Matches BModes JJ. Also @@ -39,12 +42,18 @@ The residual retry ------------------ -``scipy.linalg.eigh`` reduces ``K x = λ M x`` through a Cholesky factor -of the **mass** matrix, and that reduction degrades once ``M`` is nearly -singular — a very light beam carrying a very heavy lump. The failure is -silent: LAPACK returns confidently wrong low modes rather than raising. -On the case in ``tests/fem/test_ill_conditioned_mass.py`` it reported -0.103 Hz against a true 0.0436 Hz. +``scipy.linalg.eigh(K, M)`` reduces ``K x = λ M x`` through a Cholesky +factor of the **mass** matrix, and that reduction degrades once ``M`` is +nearly singular — a very light beam carrying a very heavy lump. The +failure is silent: LAPACK returns confidently wrong low modes rather +than raising. On the case in ``tests/fem/test_ill_conditioned_mass.py`` +it reported 0.103 Hz against a true 0.0436 Hz. + +The dense symmetric path no longer reduces that way first: it solves +the inverted pencil, which factorises ``K`` and gets that case right +directly (``tests/fem/test_dense_inverted_solver.py``). The mass-reduced +form remains as the fallback for a pencil no shift makes definite, and +the retry below still stands behind it. :func:`solve_modes` therefore measures the backward error of a dense symmetric solve and, above @@ -821,18 +830,293 @@ def _solve_sparse_shift_invert( def _solve_dense_symmetric( gk: np.ndarray, gm: np.ndarray, n_modes: int | None, ) -> tuple[np.ndarray, np.ndarray]: - """Dense LAPACK eigh on the symmetrised matrices. ``n_modes=None`` - requests the full spectrum; otherwise a subset slice is taken.""" + """Dense symmetric solve on the symmetrised matrices. ``n_modes=None`` + requests the full spectrum; otherwise the lowest ``n_modes``. + + Routed through :func:`_solve_dense_inverted`, which factorises the + stiffness side the way the sparse shift-invert path does, so the two + paths answer to the same precision. The mass-reduced standard form is + kept only as the fallback for a pencil no shift can make definite + (:func:`_solve_dense_mass_reduced`); the residual retry in + :func:`solve_modes` still stands behind that route. + """ gk_sym = 0.5 * (gk + gk.T) gm_sym = 0.5 * (gm + gm.T) + try: + return _solve_dense_inverted(gk_sym, gm_sym, n_modes) + except np.linalg.LinAlgError as exc: + _log.info( + "solve_modes: no definite shift found (%r); using the " + "mass-reduced dense eigh", exc, + ) + return _solve_dense_mass_reduced(gk_sym, gm_sym, n_modes) + + +def _solve_dense_mass_reduced( + gk_sym: np.ndarray, gm_sym: np.ndarray, n_modes: int | None, +) -> tuple[np.ndarray, np.ndarray]: + """LAPACK ``eigh(K, M)``: reduction through a Cholesky factor of ``M``. + + Every eigenvalue it returns carries an absolute error of order + ``eps * lambda_max``. That is harmless when the spectrum is compact + and ruinous when it is not — see :func:`_solve_dense_inverted`. + """ if n_modes is not None: - subset = (0, min(n_modes, gk.shape[0]) - 1) + subset = (0, min(n_modes, gk_sym.shape[0]) - 1) eigvals, eigvecs = eigh(gk_sym, gm_sym, subset_by_index=subset) else: eigvals, eigvecs = eigh(gk_sym, gm_sym) return np.asarray(eigvals), np.asarray(eigvecs) +# How many of the lowest eigenvalues set the scale of a shift. Enough to +# see past the six rigid-body modes of a free-free beam to the first +# elastic ones, few enough that the scale stays that of the low end. +_SHIFT_SCALE_MODES = 12 +# Attempts at a definite shift, growing tenfold each time, before the +# mass-reduced route is used instead. +_SHIFT_ATTEMPTS = 20 + + +def _solve_dense_inverted( + gk_sym: np.ndarray, gm_sym: np.ndarray, n_modes: int | None, +) -> tuple[np.ndarray, np.ndarray]: + """Lowest modes of ``K x = λ M x`` from the inverted pencil. + + The standard reduction ``eigh(K, M)`` solves ``L^-1 K L^-T y = λ y`` + with ``M = L L^T``, and a backward-stable symmetric eigensolver gets + each ``λ`` of that matrix to an **absolute** accuracy of about + ``eps * λ_max``. The error on the lowest mode is therefore + ``eps * λ_max / λ_min`` relative — the spectral range of the pencil, + not ``cond(M)`` alone. A beam with a small rotary-inertia floor (the + ElastoDyn adapters floor it at 1e-6) has ``λ_max / λ_min`` near 1e15, + so the lowest modes come back with an error of order one, and one + that moves with the ``subset_by_index`` window: the NREL 5MW blade's + first flap mode at 12.1 rpm read 0.715 / 0.710 / 0.766 / 0.728 Hz for + 4 / 6 / 10 / 20 requested modes. A very heavy lump on a very light + beam is the same failure reached through a near-singular ``M``. + + Solving ``M x = μ (K + s M) x`` instead, with ``μ = 1 / (λ + s)``, + puts the wanted modes at the *top* of the spectrum, where the same + absolute accuracy ``eps * μ_max`` is relative accuracy — the dense + counterpart of the shift-invert the sparse path already uses. The + returned modes then no longer depend on how many were requested. + + The shift ``s`` is zero when ``K`` is positive definite. Otherwise — + rigid-body modes, or the negative eigenvalues a column loaded past + its buckling weight has — it is set from a first estimate of the + lowest eigenvalues and grown until ``K + s M`` factorises. The error + on mode ``i`` is then about ``eps (λ_i + s)^2 / (λ_1 + s)``, and with + ``s`` at least twice the magnitude of the lowest eigenvalues that + stays relative to the modes in the window. A ``K`` that factorises + despite rigid-body modes is left unshifted: their eigenvalues are + roundoff of order ``eps λ_max``, which leaves the relative error of a + low elastic mode near ``λ_i / λ_max``. + + The inverted form is relatively accurate at the low end only; the + mass-reduced form at the high end. A request reaching into the upper + half of the spectrum takes each mode from whichever form is more + accurate there, so a full-spectrum request stays accurate throughout. + + Raises :class:`numpy.linalg.LinAlgError` when no shift makes the + pencil definite (a singular ``M`` along a direction ``K`` does not + stiffen), for the caller to fall back on. + """ + from scipy.linalg import cho_factor + + ngd = gk_sym.shape[0] + k = ngd if n_modes is None else min(n_modes, ngd) + if k <= 0: + return np.empty(0), np.empty((ngd, 0)) + + def _definite(shift: float) -> bool: + try: + cho_factor(gk_sym + shift * gm_sym, lower=True, + check_finite=False) + except np.linalg.LinAlgError: + return False + return True + + shift = 0.0 + if not _definite(shift): + # The same estimate whatever ``k`` is, so the shift, and with it + # every returned mode, does not depend on the requested window. + rough = eigh( + gk_sym, gm_sym, eigvals_only=True, + subset_by_index=(0, min(ngd, _SHIFT_SCALE_MODES) - 1), + ) + scale = float(np.max(np.abs(rough))) + if not np.isfinite(scale) or scale <= 0.0: + scale = float(np.finfo(float).eps) * ( + float(np.abs(np.diag(gk_sym)).max()) + / max(float(np.abs(np.diag(gm_sym)).max()), 1.0e-300) + ) + shift = 2.0 * scale + for _ in range(_SHIFT_ATTEMPTS): + if _definite(shift): + break + shift *= 10.0 + else: + raise np.linalg.LinAlgError( + "no shift makes K + s M positive definite" + ) + + eigvals, eigvecs = _inverted_window(gk_sym, gm_sym, k, shift) + if shift > 0.0: + eigvals, eigvecs = _resolve_swamped_modes( + gk_sym, gm_sym, k, shift, eigvals, eigvecs, + ) + + if 2 * k > ngd: + eigvals, eigvecs = _splice_upper_spectrum( + gk_sym, gm_sym, k, shift, eigvals, eigvecs, + ) + return eigvals, eigvecs + + +def _inverted_window( + gk_sym: np.ndarray, gm_sym: np.ndarray, k: int, shift: float, +) -> tuple[np.ndarray, np.ndarray]: + """The ``k`` largest ``μ`` of ``M x = μ (K + s M) x``, as ascending + ``λ = 1/μ - s``. ``μ <= 0`` (a massless direction) maps to ``+inf`` + and so sorts to the top, where the splice replaces it.""" + ngd = gk_sym.shape[0] + mu, vecs = eigh( + gm_sym, gk_sym + shift * gm_sym, + subset_by_index=(ngd - k, ngd - 1), + ) + mu = np.asarray(mu)[::-1] + vecs = np.asarray(vecs)[:, ::-1] + with np.errstate(divide="ignore"): + lam = np.where(mu > 0.0, 1.0 / np.where(mu > 0.0, mu, 1.0) - shift, + np.inf) + return lam, np.ascontiguousarray(vecs) + + +# A backward-stable solve of the inverted pencil fixes each μ to about +# ``ngd * eps * μ_max``; this factor is the safety margin on that bound. +_MU_RESOLUTION_ULPS = 100.0 +# A mode is swamped by the shift when the uncertainty that μ resolution +# leaves on ``λ = 1/μ - s`` exceeds this fraction of ``|λ|``. +_SWAMPED_RTOL = 1.0e-8 + + +def _resolve_swamped_modes( + gk_sym: np.ndarray, + gm_sym: np.ndarray, + k: int, + shift: float, + eigvals: np.ndarray, + eigvecs: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Recover modes a large shift has pushed below the μ resolution. + + With ``μ = 1/(λ + s)`` the uncertainty on ``λ`` is about + ``δμ (λ + s)^2``. When ``s`` is much larger than the soft end of the + spectrum (one strongly negative eigenvalue forces a large shift), + every soft mode maps to the same ``μ`` to working precision: their + order is arbitrary, a window edge falling among them picks any of + them, and ``1/μ - s`` returns noise. ``K = diag(-1e16, 1, 2, 3, 4)`` + with ``M = I`` gave ``[-1e16, 4]`` for two modes. + + The subspace such a cluster spans is still well determined, because + it is separated from the rest of the μ spectrum. So the window is + grown until its edge falls in a resolvable gap, and each cluster that + holds a swamped mode is re-solved by Rayleigh-Ritz on ``(K, M)`` + restricted to that subspace. A Rayleigh quotient is only quadratically + sensitive to the cluster's mixing with the well-separated modes, so + this recovers the soft modes to the accuracy ``(K, M)`` itself allows. + The lowest ``k`` of the result are returned. + """ + ngd = gk_sym.shape[0] + eps = float(np.finfo(float).eps) + + def swamped(lam: np.ndarray, tol_mu: float) -> np.ndarray: + finite = np.isfinite(lam) + out = np.zeros(lam.shape, dtype=bool) + dlam = tol_mu * (lam[finite] + shift) ** 2 + out[finite] = dlam > _SWAMPED_RTOL * np.abs(lam[finite]) + return out + + def mu_of(lam: np.ndarray) -> np.ndarray: + with np.errstate(divide="ignore"): + return np.where(np.isfinite(lam), 1.0 / (lam + shift), 0.0) + + mu_max = float(mu_of(eigvals[:1])[0]) + if not np.isfinite(mu_max) or mu_max <= 0.0: + return eigvals, eigvecs + tol_mu = _MU_RESOLUTION_ULPS * ngd * eps * mu_max + if not swamped(eigvals, tol_mu).any(): + return eigvals, eigvecs + + # Grow the window until the μ just beyond it is resolvably below the + # last μ inside it, so no cluster is cut in two. + m = k + while m < ngd: + lam_w, vec_w = _inverted_window(gk_sym, gm_sym, m + 1, shift) + mu_w = mu_of(lam_w) + if mu_w[m - 1] - mu_w[m] > tol_mu: + lam_w, vec_w = lam_w[:m], vec_w[:, :m] + break + m = min(ngd, 2 * m) + else: + lam_w, vec_w = _inverted_window(gk_sym, gm_sym, ngd, shift) + mu_w = mu_of(lam_w) + flagged = swamped(lam_w, tol_mu) + + lam_out = np.array(lam_w, dtype=float) + vec_out = np.array(vec_w, dtype=float) + start = 0 + for stop in range(1, len(mu_w) + 1): + if stop < len(mu_w) and mu_w[stop - 1] - mu_w[stop] <= tol_mu: + continue + cols = np.arange(start, stop) + cols = cols[np.isfinite(lam_w[cols])] + if cols.size and flagged[cols].any(): + v = vec_w[:, cols] + kr = v.T @ gk_sym @ v + mr = v.T @ gm_sym @ v + vals, y = eigh(0.5 * (kr + kr.T), 0.5 * (mr + mr.T)) + lam_out[cols] = vals + vec_out[:, cols] = v @ y + start = stop + + order = np.argsort(lam_out, kind="stable")[:k] + return lam_out[order], np.ascontiguousarray(vec_out[:, order]) + + +def _splice_upper_spectrum( + gk_sym: np.ndarray, + gm_sym: np.ndarray, + k: int, + shift: float, + eigvals: np.ndarray, + eigvecs: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Take each mode from whichever form resolves it more accurately. + + The inverted form's error on mode ``i`` is about + ``eps (λ_i + s)^2 / (λ_1 + s)``; the mass-reduced form's is about + ``eps λ_max``. Both are monotone in ``i``, so there is one crossover, + and below it the two forms describe the same modes in the same order. + """ + full_vals, full_vecs = eigh(gk_sym, gm_sym) + full_vals = np.asarray(full_vals) + lam_max = float(np.max(np.abs(full_vals))) + base = float(eigvals[0]) + shift + with np.errstate(over="ignore", invalid="ignore"): + inverted_err = (eigvals + shift) ** 2 / base + use_mass_reduced = ~(inverted_err < lam_max) + if not use_mass_reduced.any(): + return eigvals, eigvecs + cut = int(np.argmax(use_mass_reduced)) + vals = np.concatenate([eigvals[:cut], full_vals[cut:k]]) + vecs = np.concatenate( + [eigvecs[:, :cut], np.asarray(full_vecs)[:, cut:k]], axis=1, + ) + return vals, vecs + + def _solve_dense_general( gk: np.ndarray, gm: np.ndarray, n_modes: int | None, ) -> tuple[np.ndarray, np.ndarray]: diff --git a/src/pybmodes/io/_elastodyn/adapter.py b/src/pybmodes/io/_elastodyn/adapter.py index a088c38..9eb9314 100644 --- a/src/pybmodes/io/_elastodyn/adapter.py +++ b/src/pybmodes/io/_elastodyn/adapter.py @@ -44,6 +44,7 @@ import numpy as np +from pybmodes._numeric import trapezoid from pybmodes.io._elastodyn.types import ( ElastoDynBlade, ElastoDynMain, @@ -374,7 +375,7 @@ def _tower_top_assembly_mass( bl_len = main.tip_rad - main.hub_rad bmd = blade.b_mass_den * blade.adj_bl_ms s = blade.bl_fract * bl_len - m_bl_each = float(np.trapezoid(bmd, s)) + m_bl_each = trapezoid(bmd, s) else: m_bl_each = 0.0 m_blades = main.num_bl * m_bl_each diff --git a/src/pybmodes/io/_elastodyn/types.py b/src/pybmodes/io/_elastodyn/types.py index a15035d..d18f85b 100644 --- a/src/pybmodes/io/_elastodyn/types.py +++ b/src/pybmodes/io/_elastodyn/types.py @@ -31,6 +31,8 @@ import numpy as np +from pybmodes._numeric import trapezoid + @dataclass class ElastoDynMain: @@ -112,7 +114,7 @@ def compute_rot_mass(self, blade: ElastoDynBlade) -> float: bl_len = self.tip_rad - self.hub_rad s = blade.bl_fract * bl_len bl_mass_per_blade = float( - blade.adj_bl_ms * np.trapezoid(blade.b_mass_den, s) + blade.adj_bl_ms * trapezoid(blade.b_mass_den, s) ) return self.hub_mass + self.num_bl * bl_mass_per_blade diff --git a/src/pybmodes/io/bmi.py b/src/pybmodes/io/bmi.py index 8747b3c..abcd119 100644 --- a/src/pybmodes/io/bmi.py +++ b/src/pybmodes/io/bmi.py @@ -206,9 +206,18 @@ class PlatformSupport: Infinite-frequency hydrodynamic added mass ``A_inf``, 6×6, about the ``ref_msl`` reference, OpenFAST DOF order. hydro_K : np.ndarray - Hydrostatic restoring ``C_hst``, 6×6, about the ``ref_msl`` - reference. May be (legitimately) negative-definite in - roll/pitch for an unstable spar stabilised by mooring. + Hydrostatic restoring, 6×6, about the ``ref_msl`` reference. + ``hydro_K + mooring_K`` must hold **all** the restoring the + rigid-body modes see, including the weight term ``−m g z_G`` in + roll / pitch: a WAMIT ``.hst`` excludes it (OpenFAST applies it + through ElastoDyn gravity), so the deck constructors + (:meth:`~pybmodes.models.Tower.from_elastodyn_with_mooring`, + :meth:`~pybmodes.models.Tower.from_windio_floating`) add the + weight of platform, tower and RNA here. A BModes ``.bmi`` deck + such as OC3 Hywind folds it into ``mooring_K`` instead. A + ballast-stabilised spar's ``.hst`` alone is negative in + roll / pitch; a negative total is reported by ``Tower.run`` as + a ``UserWarning``. mooring_K : np.ndarray Linearised mooring stiffness, 6×6, about the ``ref_msl`` reference, OpenFAST DOF order. diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 0cd67a9..9b3b574 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -50,22 +50,12 @@ import numpy as np +from pybmodes._numeric import trapezoid from pybmodes.io.sec_props import SectionProperties if TYPE_CHECKING: from pybmodes.io.bmi import TipMassProps -def _trapezoid(y: np.ndarray, x: np.ndarray) -> float: - """Trapezoidal integral of ``y`` over ``x``. - - Implemented directly rather than via ``np.trapezoid`` (NumPy >= 2.0) / - ``np.trapz`` (the deprecated <2.0 spelling) so the reader works across - the advertised ``numpy>=1.26`` floor without a version shim. - """ - y = np.asarray(y, dtype=float) - x = np.asarray(x, dtype=float) - return float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(x))) - def _second_moment(w: np.ndarray, c: np.ndarray, x: np.ndarray) -> float: """Exact ``integral w(x)*c(x)**2 dx`` for piecewise-linear ``w`` and ``c``. @@ -1021,7 +1011,7 @@ def _blade_span_mass_inertia( xyz = np.vstack(coords).T seg = np.linalg.norm(np.diff(xyz, axis=0), axis=1) s = np.concatenate([[0.0], np.cumsum(seg)]) - mass = _trapezoid(mpl, s) + mass = trapezoid(mpl, s) # Resolve each section onto the rotor axes (r = in-plane radial from the # shaft, t = in-plane tangential, a = along the shaft). The reference axis # gives z spanwise, x prebend (flapwise, out of the rotor plane) and y diff --git a/src/pybmodes/io/windio_floating.py b/src/pybmodes/io/windio_floating.py index 0b6c143..99cc43a 100644 --- a/src/pybmodes/io/windio_floating.py +++ b/src/pybmodes/io/windio_floating.py @@ -268,8 +268,10 @@ def hydrostatic_restoring( ) -> np.ndarray: """6×6 hydrostatic restoring (DOF order surge, sway, heave, roll, pitch, yaw) from member geometry — the WAMIT/`.hst` buoyancy + - waterplane convention (no gravitational term; that enters via the - body mass elsewhere). For a freely-floating semi the heave / + waterplane convention (no gravitational term: the ``−m g z_G`` + weight restoring is added by :meth:`Tower.from_windio_floating` + through :func:`pybmodes.models._platform._gravitational_restoring`, + which also counts the tower and RNA). For a freely-floating semi the heave / roll / pitch entries are geometry-exact, so this matches a potential-flow `.hst` closely (integration anchor).""" S = Sx = Sy = Sxx = Syy = Sxy = 0.0 diff --git a/src/pybmodes/models/_pipeline.py b/src/pybmodes/models/_pipeline.py index 0059586..58cbc8f 100644 --- a/src/pybmodes/models/_pipeline.py +++ b/src/pybmodes/models/_pipeline.py @@ -17,6 +17,7 @@ from __future__ import annotations +import warnings from dataclasses import dataclass import numpy as np @@ -191,6 +192,75 @@ def _gravity_axial_force( return axf_g, grav_w +_RIGID_DOF_NAMES = ("surge", "sway", "heave", "roll", "pitch", "yaw") +# Relative size below which a negative restoring eigenvalue is noise +# rather than an unstable direction. Mooring stiffness is a finite +# difference of the catenary solution, so a direction with no restoring +# at all (the yaw of a single column with centreline fairleads) comes back +# at a few parts per million of the stiffest one, either sign. A real +# instability is of the order of the restoring itself. +_NEGATIVE_RESTORING_RTOL = 1.0e-4 + + +def _warn_on_negative_platform_restoring(plat: PlatformSupport) -> None: + """Warn when the platform restoring has a negative-stiffness direction. + + On a free-free floating model the rigid-body modes are resisted by + ``hydro_K + mooring_K`` alone — the flexible tower adds nothing to a + rigid motion — so a negative direction in that matrix is a + rigid-body mode with a negative eigenvalue. The general dense solver + filters such eigenvalues out, and before this check the mode simply + went missing from the result: OC3 Hywind built from its WAMIT + ``.hst`` without the body weight came back with four of its six + rigid-body modes and no indication why. + + The sign test is on the matrix scaled by the platform's own inertia, + ``D^-1/2 K D^-1/2``, which keeps the count of negative eigenvalues + (Sylvester's law of inertia) while putting translations and rotations + on a common scale for the roundoff tolerance. + """ + k = np.asarray(plat.hydro_K, dtype=float) + np.asarray( + plat.mooring_K, dtype=float, + ) + if k.shape != (6, 6) or not np.all(np.isfinite(k)): + return + k = 0.5 * (k + k.T) + d = np.diag( + np.asarray(plat.i_matrix, dtype=float) + + np.asarray(plat.hydro_M, dtype=float) + ).copy() + # A DOF with no inertia of its own (a yaw inertia left at zero) takes + # the largest of its group, so its roundoff is not magnified. + for group in (slice(0, 3), slice(3, 6)): + part = d[group] + valid = np.isfinite(part) & (part > 0.0) + part[~valid] = float(part[valid].max()) if valid.any() else 1.0 + scale = 1.0 / np.sqrt(d) + lam, vec = np.linalg.eigh(k * np.outer(scale, scale)) + worst = float(np.max(np.abs(lam), initial=0.0)) + negative = lam < -_NEGATIVE_RESTORING_RTOL * worst + if not negative.any(): + return + names = sorted({ + _RIGID_DOF_NAMES[int(np.argmax(np.abs(vec[:, i])))] + for i in np.flatnonzero(negative) + }, key=_RIGID_DOF_NAMES.index) + warnings.warn( + f"the platform restoring (hydro_K + mooring_K) has " + f"{int(negative.sum())} negative-stiffness direction(s), dominated " + f"by {', '.join(names)}: the platform is statically unstable " + f"there, and the corresponding rigid-body mode(s) have negative " + f"eigenvalues, which the general (non-symmetric) solver path " + f"leaves out of the returned modes. A WAMIT " + f".hst excludes the body weight (the -m g z_G term), and ballast " + f"carried outside PtfmMass — e.g. HydroDyn filled members — is not " + f"in the model either; check that the restoring includes every " + f"weight and buoyancy contribution.", + UserWarning, + stacklevel=4, # past run_fem and Tower.run, to the caller + ) + + def run_fem( bmi: BMIFile, n_modes: int = 20, @@ -280,6 +350,8 @@ def run_fem( # hand-built BMIFile instances follow the same path. platform_nd = None if isinstance(bmi.support, PlatformSupport): + if bmi.hub_conn == 2: + _warn_on_negative_platform_restoring(bmi.support) platform_nd = nondim_platform(bmi.support, nd) # Embedded tension wires within the platform block plat_wires = bmi.support.wires diff --git a/src/pybmodes/models/_platform.py b/src/pybmodes/models/_platform.py index ee8179a..877caca 100644 --- a/src/pybmodes/models/_platform.py +++ b/src/pybmodes/models/_platform.py @@ -13,7 +13,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -"""Platform-scalar parsing + inertia-matrix assembly helpers. +"""Platform-scalar parsing, inertia-matrix and weight-restoring helpers. Used by :meth:`pybmodes.models.Tower.from_elastodyn_with_mooring` and :meth:`pybmodes.models.Tower.from_windio_floating` (deck-fallback @@ -36,6 +36,9 @@ if TYPE_CHECKING: import numpy as np + from pybmodes.io.bmi import BMIFile + from pybmodes.io.sec_props import SectionProperties + def _scan_platform_fields(dat_path: pathlib.Path) -> dict[str, float]: """Scan an ElastoDyn ``.dat`` for the platform scalars used to @@ -137,3 +140,84 @@ def _platform_inertia_matrix(ptfm: dict[str, float]) -> np.ndarray: i_mat[4, 4] = ptfm["PtfmPIner"] # pitch inertia about CM (DOF 4) i_mat[5, 5] = ptfm["PtfmYIner"] # yaw inertia about CM return i_mat + + +#: Standard gravity (m/s²), the OpenFAST ``Gravity`` default. +STANDARD_GRAVITY = 9.80665 + + +def _gravitational_restoring( + bmi: BMIFile, + sp: SectionProperties, + *, + g: float = STANDARD_GRAVITY, +) -> np.ndarray: + """Roll / pitch restoring from the weight of the whole floating system. + + A WAMIT ``.hst`` holds only the hydrostatic part of the restoring, + ``ρ g I_wp + ρ g V z_B`` — the weight of the body is deliberately left + out, because OpenFAST adds it through ElastoDyn's gravity loading of + the platform, tower and rotor-nacelle assembly. A model assembled from + those decks has to add it back, or a ballast-stabilised spar reads as + hydrostatically unstable: on OC3 Hywind the ``.hst`` alone gives + ``C44 = C55 = −5.0e9 N·m/rad``. + + For a rigid rotation ``θ`` about a horizontal axis through the + hydro/mooring reference point, a mass ``m`` whose centre sits at + height ``z`` above that point drops by ``z θ² / 2``, so its weight + contributes ``−m g z`` to both ``C44`` and ``C55`` — the ``−m g z_G`` + half of the classical restoring ``ρ g V (z_B − z_G) + ρ g I_wp`` + (Faltinsen 1990, *Sea Loads on Ships and Offshore Structures*, ch. 2), + and the half that Jonkman (2007) NREL/TP-500-41958 leaves to the + structural model rather than the hydrodynamic one. The sum runs over + + - the platform: ``mass_pform`` at ``z = −cm_pform``; + - the flexible tower: ``mass_den`` (times the ``sec_mass`` scaling) + integrated from the base at ``z = −draft`` to the top at + ``z = radius``, plus any discrete point masses; + - the tower-top lump: ``tip_mass`` at ``z = radius + cm_axial``, + + all measured from MSL (see :doc:`/conventions`) and taken relative to + the reference point at ``z = −ref_msl``. Mass above the reference + destabilises, mass below stabilises. + + Only the rigid-body term is added. The tower's own geometric + (P-delta) softening under its weight is a separate effect that a + free-base model does not include (``Tower.run(gravity=...)`` refuses + it there). Horizontal CM offsets are not given a weight coupling + either: the conventions assume a platform trimmed with its centre of + gravity over its centre of buoyancy, where those terms cancel against + their buoyancy counterparts. + + Returns a 6×6 matrix in OpenFAST DOF order with only ``[3, 3]`` and + ``[4, 4]`` set. + """ + import numpy as np + + from pybmodes._numeric import trapezoid + from pybmodes.io.bmi import PlatformSupport + + ps = bmi.support + if not isinstance(ps, PlatformSupport): + raise TypeError( + "weight restoring needs a PlatformSupport; got " + f"{type(ps).__name__}" + ) + base_z = -float(ps.draft) + top_z = float(bmi.radius) + ref_z = -float(ps.ref_msl) + + span = np.asarray(sp.span_loc, dtype=float) + mass_den = np.asarray(sp.mass_den, dtype=float) * float(bmi.scaling.sec_mass) + z = base_z + span * (top_z - base_z) + # First moment of each part about the reference point, in kg·m. + moment = float(ps.mass_pform) * (-float(ps.cm_pform) - ref_z) + moment += trapezoid(mass_den * (z - ref_z), z) + for pm in bmi.point_masses: + moment += float(pm.mass) * (base_z + float(pm.height) - ref_z) + tip = bmi.tip_mass + moment += float(tip.mass) * (top_z + float(tip.cm_axial) - ref_z) + + c = np.zeros((6, 6)) + c[3, 3] = c[4, 4] = -g * moment + return c diff --git a/src/pybmodes/models/tower.py b/src/pybmodes/models/tower.py index d758eb3..63d41b7 100644 --- a/src/pybmodes/models/tower.py +++ b/src/pybmodes/models/tower.py @@ -40,6 +40,7 @@ from pybmodes.io.sec_props import SectionProperties from pybmodes.models._pipeline import run_fem from pybmodes.models._platform import ( + _gravitational_restoring, _platform_inertia_matrix, _scan_platform_fields, ) @@ -211,6 +212,12 @@ class Tower: # deck-derived model, where those numbers don't exist. _construction: ConstructionInputs | None = None + # Gravity used for the weight restoring folded into ``hydro_K`` by + # the deck-built floating constructors, so that ``add_point_mass`` + # can add the new lump's share. ``None`` when no weight term was + # added (every other constructor, or a floater without HydroDyn). + _weight_g: float | None = None + def __init__( self, bmi_path: str | pathlib.Path, *, n_nodes: int | None = None, ) -> None: @@ -902,6 +909,15 @@ def from_elastodyn_with_mooring( to the WAMIT ``.1`` and ``.hst`` files). Optional — if ``hydrodyn_dat_path`` is omitted, both default to zero, so the resulting model couples only mooring + platform inertia. + - **Gravitational restoring** of the platform, tower and RNA, + ``−m g z`` in roll and pitch about the reference point, added + to ``hydro_K`` alongside ``C_hst``. The ``.hst`` excludes the + body weight, which OpenFAST applies through ElastoDyn's + gravity; without it a ballast-stabilised spar has negative + roll / pitch stiffness. Added only when a HydroDyn deck is + given, since weight without the buoyancy that balances it does + not describe a floating body. See + :func:`pybmodes.models._platform._gravitational_restoring`. - **Platform inertia** from the ``PtfmMass`` / ``PtfmRIner`` / ``PtfmPIner`` / ``PtfmYIner`` / ``PtfmCM*`` / ``PtfmRefzt`` scalars in the ElastoDyn main file. The 6 × 6 ``i_matrix`` is @@ -1059,9 +1075,23 @@ def from_elastodyn_with_mooring( bmi.tow_support = 1 bmi.support = platform_support + # The WAMIT ``.hst`` excludes the body weight, which OpenFAST + # supplies through ElastoDyn's gravity loading of the platform, + # tower and RNA. Without it a ballast-stabilised spar reads as + # unstable in roll and pitch (OC3 Hywind: C44 = C55 = −5.0e9 + # N·m/rad from the .hst alone), and those two rigid-body modes + # were lost. Only paired with hydrostatics: weight without the + # buoyancy that balances it is not a floating body. + if hydrodyn_dat_path is not None: + platform_support.hydro_K = C_hst + _gravitational_restoring( + bmi, sp, g=STANDARD_GRAVITY, + ) + obj = cls.__new__(cls) obj._bmi = bmi obj._sp = sp + if hydrodyn_dat_path is not None: + obj._weight_g = STANDARD_GRAVITY if n_nodes is not None: obj.refine_mesh(n_nodes) return obj @@ -1434,10 +1464,18 @@ def _tower_grid_od_wt( bmi.hub_conn = 2 bmi.tow_support = 1 bmi.support = platform_support + # Both hydrostatic tiers — the WAMIT ``.hst`` and the member + # waterplane integration — are buoyancy-only; add the weight of + # platform, tower and RNA exactly as from_elastodyn_with_mooring + # does, so the deck-backed tier stays equivalent to it. + platform_support.hydro_K = C_hst + _gravitational_restoring( + bmi, sp, g=g, + ) obj = cls.__new__(cls) obj._bmi = bmi obj._sp = sp + obj._weight_g = float(g) obj.coeff_validation = None return obj @@ -1732,6 +1770,8 @@ def add_point_mass( would also change the centrifugal tension distribution, which is a separate modelling track. """ + import numpy as np + from pybmodes.io.bmi import PlatformSupport, PointMass pm = PointMass(height=float(height), mass=float(mass)) @@ -1748,6 +1788,17 @@ def add_point_mass( f"the very top." ) self._bmi.point_masses = (*self._bmi.point_masses, pm) + # A floater built from decks carries the system weight in hydro_K; + # the new lump's weight belongs there too, measured the same way + # as in _gravitational_restoring (beam base at z = -draft, relative + # to the reference point at z = -ref_msl). + ps = self._bmi.support + if self._weight_g is not None and isinstance(ps, PlatformSupport): + z_rel = -float(ps.draft) + pm.height + float(ps.ref_msl) + k = np.array(ps.hydro_K, dtype=float, copy=True) + k[3, 3] -= self._weight_g * pm.mass * z_rel + k[4, 4] -= self._weight_g * pm.mass * z_rel + ps.hydro_K = k return self def run( diff --git a/src/pybmodes/options.py b/src/pybmodes/options.py index 09ddfff..f059ed5 100644 --- a/src/pybmodes/options.py +++ b/src/pybmodes/options.py @@ -73,18 +73,21 @@ class SolverOptions: residual_retry_threshold : float, default 0.1 Largest per-mode relative residual ``||K x - λ M x|| / ||K x||`` a symmetric solve may return before - the general dense path is tried as well. ``scipy.linalg.eigh`` - reduces the generalised problem through a Cholesky factor of the - **mass** matrix, which loses accuracy once that matrix is nearly - singular — a very light beam carrying a very heavy lump — and - returns confidently wrong low modes rather than failing. + the general dense path is tried as well. The mass-reduced + ``scipy.linalg.eigh(K, M)`` reduces the generalised problem + through a Cholesky factor of the **mass** matrix, which loses + accuracy once that matrix is nearly singular — a very light beam + carrying a very heavy lump — and returns confidently wrong low + modes rather than failing. The dense symmetric path now solves + the inverted pencil first, which does not have this failure, so + the retry only guards the mass-reduced fallback. The band this sits in is narrower than it looks. Ordinary solves - land at or below ~1e-3; a real deck whose adapter leaves ``M`` - genuinely ill-conditioned (the bundled NREL 5MW land tower, cond - ~4e10) reaches ~2e-2 and is *not* meant to trigger; the degraded - regime starts around 0.7. The default splits the last two gaps - with roughly 5x either side. + land at or below ~1e-3; on the mass-reduced route a real deck + whose adapter leaves ``M`` genuinely ill-conditioned (the bundled + NREL 5MW land tower, cond ~4e10) reaches ~1e-2 and is *not* + meant to trigger; the degraded regime starts around 0.7. The + default splits the last two gaps with roughly 5x either side. residual_retry_resolved : float, default 1e-3 Backward error the candidate must reach on a mode before that mode can justify a swap. Guards the case the ratio alone cannot: @@ -124,11 +127,11 @@ class SolverOptions: default sits in the middle with roughly 60x margin on the noise side and 100x on the rescue side. - A marginal win is refused for a second reason as well: on the - bundled NREL 5MW land deck the general path is 1.4x better while - *breaking* a physically real degenerate fore-aft / side-side pair - the symmetric solver resolves exactly, which the downstream FA / - SS classifier depends on. + A marginal win is refused for a second reason as well: accepting + replaces the whole spectrum, and a general solve that is only + slightly better on one mode is no evidence that it is better on + the others — on the mass-reduced route the bundled NREL 5MW land + deck showed a 1.4x win that moved its fore-aft / side-side pair. """ sparse_ndof_threshold: int = 500 diff --git a/src/pybmodes/workflows/windio.py b/src/pybmodes/workflows/windio.py index 01b0fda..12d5940 100644 --- a/src/pybmodes/workflows/windio.py +++ b/src/pybmodes/workflows/windio.py @@ -200,6 +200,22 @@ def _doc_is_floating(doc: dict | None) -> bool: return isinstance(comps, dict) and "floating_platform" in comps +def _inside_openfast_tree(p: pathlib.Path, root: pathlib.Path) -> bool: + """Whether ``p`` sits in an OpenFAST deck tree *below* ``root``. + + The RWT layout keeps OpenFAST inputs (whose yaml sidecars are not + ontologies) under an ``OpenFAST`` directory next to the ontology, so + those are skipped. Only the part of the path below the search root is + examined: a checkout that merely lives somewhere under a directory + whose name contains "OpenFAST" (``~/OpenFAST-GUI/...``) must not have + every yaml in it excluded, which is what matching on the absolute + path did. + """ + return any( + "openfast" in part.lower() for part in p.relative_to(root).parts + ) + + def discover_windio_inputs( path: str | pathlib.Path, ) -> WindioDiscovery: @@ -229,7 +245,7 @@ def discover_windio_inputs( elif path.is_dir(): cands = sorted( p for p in path.rglob("*.yaml") - if "OpenFAST" not in str(p) and "openfast" not in str(p) + if not _inside_openfast_tree(p, path) and _load_windio_doc(p) is not None ) if not cands: diff --git a/tests/fem/test_dense_inverted_solver.py b/tests/fem/test_dense_inverted_solver.py new file mode 100644 index 0000000..255858c --- /dev/null +++ b/tests/fem/test_dense_inverted_solver.py @@ -0,0 +1,278 @@ +"""The dense symmetric path solves the inverted pencil, so its lowest +modes do not depend on how many were requested. + +``scipy.linalg.eigh(K, M)`` reduces the pencil through a Cholesky factor +of ``M`` and returns each eigenvalue to an *absolute* accuracy of about +``eps * lambda_max`` (the Weyl bound for a backward-stable symmetric +eigensolver; Golub and Van Loan 2013, *Matrix Computations*, 4th ed., +ch. 8). The relative error on the lowest mode is therefore about +``eps * lambda_max / lambda_min``. The ElastoDyn adapters floor rotary +inertia at 1e-6, which pushes that ratio to ~1e15 on the NREL 5MW blade, +so the reduction returned first-flap frequencies of 0.715 / 0.710 / +0.766 / 0.728 Hz for 4 / 6 / 10 / 20 requested modes. + +Solving ``M x = mu (K + s M) x`` with ``mu = 1 / (lambda + s)`` puts the +wanted modes at the top of the spectrum, where the same absolute error +is a relative one — the standard shift-invert spectral transformation +(Ericsson and Ruhe 1980, *Math. Comp.* 35(152), 1251-1268; Bathe 2014, +*Finite Element Procedures*, 2nd ed., ch. 11), and the transformation the +sparse ``eigsh(sigma=0)`` path already applies. + +References for the physics: a cantilever whose beam mass is negligible +next to a tip lump is a spring-mass oscillator on the static tip +stiffness ``3 EI / L^3`` (Blevins 1979, *Formulas for Natural Frequency +and Mode Shape*, Table 8-1). +""" + +from __future__ import annotations + +import warnings + +import numpy as np +import pytest + +from pybmodes.cli import _resolve_examples_root +from pybmodes.fem import solver +from pybmodes.fem.assembly import assemble +from pybmodes.fem.nondim import RM, ROMG, make_params, nondim_tip_mass +from pybmodes.fem.solver import eigvals_to_hz, solve_modes +from pybmodes.io.bmi import TipMassProps + +NREL5MW_LAND = ( + _resolve_examples_root() / "reference_decks" / "nrel5mw_land" + / "NRELOffshrBsline5MW_Onshore_ElastoDyn.dat" +) +_needs_deck = pytest.mark.skipif( + not NREL5MW_LAND.is_file(), + reason=f"bundled NREL 5MW deck not present at {NREL5MW_LAND}", +) + +L = 100.0 +EI = 1.0e10 +M_TIP = 4.0e5 + + +def _cantilever_with_tip_lump( + nselt: int, mass_den: float, +) -> tuple[np.ndarray, np.ndarray]: + """Same construction as ``tests/fem/test_ill_conditioned_mass.py``.""" + nd = make_params(radius=L, hub_rad=0.0, rot_rpm=0.0) + eiy = EI / nd.ref4 + eli = 1.0 / nselt + el = np.full(nselt, eli) + xb = np.array([1.0 - (i + 1) * eli for i in range(nselt)]) + tip = nondim_tip_mass( + TipMassProps(mass=M_TIP, cm_offset=0.0, cm_axial=0.0, ixx=0.0, + iyy=0.0, izz=0.0, ixy=0.0, izx=0.0, iyz=0.0), + nd, beam_type=2, id_form=1, hub_conn=1, + ) + gk, gm, _ = assemble( + nselt=nselt, el=el, xb=xb, cfe=np.zeros(nselt), + eiy=np.full(nselt, eiy), eiz=np.full(nselt, eiy), + gj=np.full(nselt, 1.0e3 * eiy), eac=np.full(nselt, 100.0), + rmas=np.full(nselt, mass_den / RM), + skm1=np.full(nselt, 1.0e-5), skm2=np.full(nselt, 1.0e-5), + eg=np.zeros(nselt), ea=np.zeros(nselt), omega2=0.0, + sec_loc=np.array([0.0, 1.0]), str_tw=np.zeros(2), hub_conn=1, + tip_mass=tip, + ) + return gk, gm + + +def _frequencies(model_factory, n_modes: int) -> np.ndarray: + with warnings.catch_warnings(): + # The bundled blade carries two expected stiffness-jump warnings. + warnings.simplefilter("ignore", UserWarning) + return np.asarray( + model_factory().run(n_modes=n_modes).frequencies, dtype=float, + ) + + +def _blade(): + from pybmodes.models import RotatingBlade + + return RotatingBlade.from_elastodyn(NREL5MW_LAND) + + +def _tower(): + from pybmodes.models import Tower + + return Tower.from_elastodyn(NREL5MW_LAND) + + +@_needs_deck +class TestNModesIndependence: + """The regression: the lowest modes must not move with ``n_modes``.""" + + @pytest.mark.parametrize("factory", [_blade, _tower], ids=["blade", "tower"]) + def test_lowest_modes_agree_for_4_and_20_requested(self, factory): + f4 = _frequencies(factory, 4) + f20 = _frequencies(factory, 20) + np.testing.assert_allclose(f20[:4], f4, rtol=1.0e-6) + + def test_blade_first_flap_across_every_window(self): + """The four windows that spread 8 % before the fix.""" + flap = [_frequencies(_blade, n)[0] for n in (4, 6, 10, 20)] + assert max(flap) - min(flap) < 1.0e-6 * min(flap) + + @pytest.mark.parametrize("factory", [_blade, _tower], ids=["blade", "tower"]) + def test_dense_agrees_with_sparse(self, factory, monkeypatch): + """Both paths are shift-invert at heart, so they must agree.""" + dense = _frequencies(factory, 10) + monkeypatch.setattr(solver, "_SPARSE_NDOF_THRESHOLD", 0) + sparse = _frequencies(factory, 10) + np.testing.assert_allclose(dense, sparse, rtol=1.0e-6) + + def test_the_dense_path_actually_ran(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + diag = _blade().run(n_modes=10).diagnostics + assert diag is not None + assert diag.path == "dense_symmetric" + assert diag.residual_fallback is False + # The mass-reduced route left this deck at a backward error of + # ~5e-2; the inverted one solves it to roundoff. + assert diag.max_residual < 1.0e-6 + + def test_the_mass_reduced_route_really_did_depend_on_the_window(self): + """Pins the failure the fix removes, on the same matrices.""" + from scipy.linalg import eigh + + import pybmodes.models._pipeline as pipeline + + captured = {} + original = pipeline.solve_modes + + def _capture(gk, gm, n_modes=None, **kwargs): + captured["pair"] = (gk.copy(), gm.copy()) + return original(gk, gm, n_modes, **kwargs) + + mp = pytest.MonkeyPatch() + mp.setattr(pipeline, "solve_modes", _capture) + try: + _frequencies(_blade, 4) + finally: + mp.undo() + gk, gm = captured["pair"] + gk = 0.5 * (gk + gk.T) + gm = 0.5 * (gm + gm.T) + low = [ + eigh(gk, gm, subset_by_index=(0, n - 1), eigvals_only=True)[0] + for n in (4, 10) + ] + assert abs(low[0] - low[1]) / low[1] > 1.0e-2 + + +class TestNearSingularMassSolvesDirectly: + """The heavy-lump case the residual retry was built for no longer + needs it: factorising K is unaffected by a near-singular M.""" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_analytic_frequency_without_a_retry(self, nselt): + gk, gm = _cantilever_with_tip_lump(nselt, 1.0e-2) + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + f = float(eigvals_to_hz(eigvals, ROMG)[0]) + f_ref = float(np.sqrt(3.0 * EI / (M_TIP * L**3)) / (2.0 * np.pi)) + assert f == pytest.approx(f_ref, rel=5.0e-3) + assert diag.path == "dense_symmetric" + assert diag.residual_fallback is False + assert diag.max_residual < 1.0e-6 + + +class TestShiftSelection: + """``K`` is not always definite; the shift has to cope.""" + + @staticmethod + def _free_free(seed: int, n_zero: int) -> tuple[np.ndarray, np.ndarray]: + n = 14 + rng = np.random.default_rng(seed) + a = rng.normal(size=(n, n)) + m = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n - n_zero)) + k = b @ b.T + return 0.5 * (k + k.T), 0.5 * (m + m.T) + + @pytest.mark.parametrize("n_zero", [1, 3, 6]) + def test_rigid_body_modes_and_elastic_modes(self, n_zero): + from scipy.linalg import eigh + + gk, gm = self._free_free(11, n_zero) + ref = eigh(gk, gm, eigvals_only=True) + vals, vecs = solver._solve_dense_symmetric(gk, gm, 8) + scale = float(np.max(np.abs(ref))) + assert np.all(np.abs(vals[:n_zero]) < 1.0e-9 * scale) + np.testing.assert_allclose(vals[n_zero:], ref[n_zero:8], rtol=1.0e-9) + residual = gk @ vecs - (gm @ vecs) * vals + assert np.max(np.abs(residual)) < 1.0e-8 * np.max(np.abs(gk)) + + def test_negative_eigenvalues_are_kept(self): + """An indefinite K — a column loaded past its buckling weight — + must come back with its negative modes, as ``eigh`` would.""" + from scipy.linalg import eigh + + rng = np.random.default_rng(3) + q, _r = np.linalg.qr(rng.normal(size=(10, 10))) + gk = q @ np.diag([-2.0, -0.5, 0.3, 1.0, 2.0, 3.0, 5.0, 8.0, 13.0, + 21.0]) @ q.T + gk = 0.5 * (gk + gk.T) + gm = np.eye(10) + vals, _v = solver._solve_dense_symmetric(gk, gm, 5) + np.testing.assert_allclose( + vals, eigh(gk, gm, eigvals_only=True)[:5], rtol=1.0e-10, + ) + + @pytest.mark.parametrize("n_modes", [1, 2, 3, 4, 5]) + def test_large_negative_mode_does_not_swamp_the_soft_ones(self, n_modes): + """One strongly negative eigenvalue forces a shift so large that + every soft mode maps to the same ``μ = 1/(λ + s)``. Their order was + then arbitrary and ``n_modes = 2`` returned ``[-1e16, 4]``.""" + gk = np.diag([-1.0e16, 1.0, 2.0, 3.0, 4.0]) + vals, vecs = solver._solve_dense_symmetric(gk, np.eye(5), n_modes) + np.testing.assert_allclose( + vals, [-1.0e16, 1.0, 2.0, 3.0, 4.0][:n_modes], rtol=1.0e-12, + ) + residual = gk @ vecs - vecs * vals + assert np.max(np.abs(residual)) < 1.0e-6 + + def test_swamped_modes_do_not_depend_on_the_window(self): + """A rotated version, where the soft cluster is no longer aligned + with the coordinates. The soft modes can only be as accurate as + ``eps * 1e14`` allows, but they must not move with ``n_modes``.""" + rng = np.random.default_rng(0) + q, _r = np.linalg.qr(rng.normal(size=(40, 40))) + d = np.concatenate([[-1.0e14, -3.0e13], np.arange(1.0, 39.0)]) + gk = q @ np.diag(d) @ q.T + gk = 0.5 * (gk + gk.T) + gm = np.eye(40) + v3 = solver._solve_dense_symmetric(gk, gm, 3)[0] + v9 = solver._solve_dense_symmetric(gk, gm, 9)[0] + np.testing.assert_array_equal(v9[:3], v3) + np.testing.assert_allclose(v9, d[:9], rtol=1.0e-10, atol=1.0e-2) + + +class TestFullSpectrum: + def test_full_request_matches_a_subset_at_the_bottom(self): + gk, gm = _cantilever_with_tip_lump(27, 1.0e-2) + low, _v = solve_modes(gk, gm, n_modes=6) + full, _w = solve_modes(gk, gm, n_modes=None) + assert full.size == gk.shape[0] + # Relative accuracy of mode i is about eps * lambda_i / lambda_1; + # this pencil spans 3e7 across the six, so 1e-6 is roundoff. + np.testing.assert_allclose(full[:6], low, rtol=1.0e-6) + + def test_full_request_keeps_the_top_of_the_spectrum(self): + """The inverted form is weakest at the top, where the splice takes + the mass-reduced form instead; the largest eigenvalue must match + that form, which is accurate there.""" + from scipy.linalg import eigh + + gk, gm = _cantilever_with_tip_lump(27, 1.0e-2) + full, _w = solve_modes(gk, gm, n_modes=None) + top = eigh(0.5 * (gk + gk.T), 0.5 * (gm + gm.T), eigvals_only=True) + assert full[-1] == pytest.approx(top[-1], rel=1.0e-9) + assert np.all(np.isfinite(full)) + assert np.all(np.diff(full) >= 0.0) diff --git a/tests/fem/test_gravity_and_point_mass.py b/tests/fem/test_gravity_and_point_mass.py index 4f4f5fe..a57186e 100644 --- a/tests/fem/test_gravity_and_point_mass.py +++ b/tests/fem/test_gravity_and_point_mass.py @@ -27,6 +27,8 @@ from __future__ import annotations +import warnings + import numpy as np import pytest @@ -277,13 +279,14 @@ def test_lump_matches_static_cantilever_stiffness(self, tmp_path): a = 60.0 m_lump = 5.0e5 # A beam mass this far below the lump makes the generalised mass - # matrix nearly singular, which the symmetric eigensolvers handle - # badly. The solver detects that from the backward error and - # redoes the solve on the general path, so a coarse mesh is fine - # here — see tests/fem/test_ill_conditioned_mass.py. + # matrix nearly singular. The dense solver factorises the + # stiffness side instead of the mass, so this solves directly — + # no residual retry, no warning — and a coarse mesh is fine here; + # see tests/fem/test_dense_inverted_solver.py. tower = _synthetic_tower(tmp_path, mass_den=1.0e-2, n_elements=25) tower.add_point_mass(a, m_lump) - with pytest.warns(RuntimeWarning, match="backward error"): + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) f = tower.run(6, check_model=False).frequencies[0] f_ref = np.sqrt(3.0 * EI_PHYS / (m_lump * a**3)) / (2.0 * np.pi) assert f == pytest.approx(f_ref, rel=5.0e-3) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py index 7bd0b62..fba21be 100644 --- a/tests/fem/test_ill_conditioned_mass.py +++ b/tests/fem/test_ill_conditioned_mass.py @@ -1,13 +1,21 @@ -"""The symmetric eigensolvers degrade silently on a near-singular mass -matrix, and the solver has to notice. - -``scipy.linalg.eigh`` reduces ``K x = lambda M x`` through a Cholesky -factor of the **mass** matrix. When that matrix is nearly singular — a -very light beam carrying a very heavy lump — the reduction loses -accuracy, and LAPACK returns confidently wrong low modes rather than -raising. On the case pinned below the dense symmetric path reported -0.103 Hz against a true 0.0436 Hz, a factor of 2.4, with no error and no -warning. +"""The mass-reduced symmetric eigensolve degrades silently on a +near-singular mass matrix, and the solver has to notice. + +``scipy.linalg.eigh(K, M)`` reduces ``K x = lambda M x`` through a +Cholesky factor of the **mass** matrix. When that matrix is nearly +singular — a very light beam carrying a very heavy lump — the reduction +loses accuracy, and LAPACK returns confidently wrong low modes rather +than raising. On the case pinned below that route reported 0.103 Hz +against a true 0.0436 Hz, a factor of 2.4, with no error and no warning. + +The dense symmetric path no longer takes that route first: it solves the +inverted pencil, which factorises ``K`` and gets this case right on its +own (``tests/fem/test_dense_inverted_solver.py``). The mass-reduced form +survives as the fallback for a pencil no shift makes definite, and the +residual retry pinned here still stands behind it. Every test in this +module therefore runs with the inverted route switched off (the +``_mass_reduced_route`` fixture below), so each one keeps exercising the +guard it was written for. The sparse path is exempt and must stay exempt: ``eigsh(sigma=0, mode='normal')`` factorises ``K`` instead, so a near-singular ``M`` does @@ -19,11 +27,12 @@ be wrong. Three things it does *not* establish, each learned the hard way and each pinned below: -- **A large error is not evidence of a breakdown.** The bundled NREL 5MW - land deck sits at ~2e-2 because its adapter leaves ``M`` at cond ~4e10; - swapping there churns a validated frequency by 0.84 % and splits a - degenerate fore-aft / side-side pair the symmetric solver resolves - exactly, which the FA / SS classifier depends on. +- **A large error is not by itself a case for the general path.** On + the mass-reduced route the bundled NREL 5MW land deck sits at ~1e-2, + because its adapter leaves ``M`` at cond ~4e10 and its pencil spans + some fifteen decades; the general path shares that range and is no + reliable cure, so a swap there is churn rather than a rescue. The + inverted route removes the error at the source (backward error ~1e-10). - **A small error is not evidence of a rescue.** A rigid-body mode's residual divides one roundoff quantity by another, so its value is arbitrary: 12.4, 0.79 and 0.076 have all been measured on healthy @@ -64,6 +73,18 @@ REALISTIC = 1.0e3 +@pytest.fixture(autouse=True) +def _mass_reduced_route(monkeypatch): + """Route every dense symmetric solve through the mass-reduced form, + the only route the residual retry now stands behind.""" + import pybmodes.fem.solver as solvermod + + def _no_definite_shift(*_args, **_kwargs): + raise np.linalg.LinAlgError("inverted route disabled for this module") + + monkeypatch.setattr(solvermod, "_solve_dense_inverted", _no_definite_shift) + + def _analytic() -> float: return float(np.sqrt(3.0 * EI / (M_TIP * L**3)) / (2.0 * np.pi)) @@ -178,8 +199,7 @@ class TestMarginalImprovementIsRefused: A synthetic pair where the symmetric solve is above the threshold but the general path cannot do materially better must keep the symmetric - result — the behaviour that protects the bundled land deck's - degenerate fore-aft / side-side pair. + result. """ def test_marginal_gain_keeps_the_symmetric_result(self, monkeypatch): diff --git a/tests/test_campbell.py b/tests/test_campbell.py index 0de0117..595220e 100644 --- a/tests/test_campbell.py +++ b/tests/test_campbell.py @@ -152,12 +152,16 @@ class TestCentrifugalStiffening: def test_first_flap_lifts_endpoint_to_endpoint( self, spec_sweep: CampbellResult ) -> None: - """Compare parked vs rated only — the per-step noise floor on the - ill-conditioned ElastoDyn-blade FEM is ~5 % on the lowest mode, - so a 3-rpm-point monotonicity test is unreliable. The Wright - (1982) Southwell estimate ``ω² = ω₀² + K·Ω²`` with K ≈ 1.12 (NREL - 5MW) predicts a ~4–5 % frequency lift over 0..12.1 rpm; gate at - 3 % to clear the noise floor. + """Compare parked vs rated. The Wright (1982) Southwell estimate + ``ω² = ω₀² + K·Ω²`` with K ≈ 1.12 (NREL 5MW) predicts a lift of + several percent over 0..12.1 rpm; gate at 3 %. + + The gate used to be marginal (2.95 % measured on one BLAS + build): the mass-reduced dense solve returned the ElastoDyn + blade's first flap mode with an error of several percent that + moved with the requested mode count. The inverted dense solve + resolves it to roundoff (0.677 to 0.729 Hz, +7.6 %); see + ``tests/fem/test_dense_inverted_solver.py``. """ slot = _first_label_index(spec_sweep, "flap") f0 = spec_sweep.frequencies[0, slot] diff --git a/tests/test_classifier.py b/tests/test_classifier.py index 149c123..c640de4 100644 --- a/tests/test_classifier.py +++ b/tests/test_classifier.py @@ -237,26 +237,28 @@ def test_iea34_no_degeneracy_warning(): """The IEA-3.4 modes 1-2 should resolve to clean FA/SS shapes with no warning emitted. - The deck's tower is symmetric in stiffness; with a sufficiently large - eigenvalue request scipy/LAPACK returns the degenerate pair at - *identical* frequencies but in an arbitrarily-rotated basis (≈ 52/48 - FA/SS mix). We use ``n_modes=10`` here to match what the case-study - script in ``cases/iea3mw_land/run.py`` exercises and to put the - eigensolver on the full-solve path; smaller subset requests produce - slightly-different numerical answers where the eigensolver lifts the - degeneracy artificially and the modes come back already-separated. + The tower is symmetric in stiffness, but the lumped RNA is not: its + roll and pitch inertias differ by almost a factor of two (ixx 1.41e6, + iyy 2.70e6 kg·m²) and it carries an ixz product. That splits the first + pair by about 0.09 %, which the dense, sparse and inverted-pencil + solves all agree on. Solvers before 1.19 returned the pair exactly + degenerate in a mixed (≈ 52/48) basis only because their error on the + lowest modes exceeded the split; this test used to assert that + artefact. What it gates now is the user-facing contract that either + way the resolver stays silent and the first two modes come back + FA-pure and SS-pure. """ tower = Tower.from_elastodyn(_IEA34_MAIN) modal = tower.run(n_modes=10) - # Sanity: this deck is expected to produce a degenerate first pair - # (gap exactly zero on a full solve — the structural model is symmetric - # in stiffness and the rigid-RNA c.m. offset is small enough not to - # lift the degeneracy at the FEM precision we operate at). - assert _is_degenerate_pair(modal.shapes[0], modal.shapes[1]), ( - f"expected a degenerate first pair from the full solve; got " - f"freqs {modal.shapes[0].freq_hz:.6f}, {modal.shapes[1].freq_hz:.6f}" + f0, f1 = modal.shapes[0].freq_hz, modal.shapes[1].freq_hz + assert 1e-4 < (f1 - f0) / f0 < 5e-3, ( + f"first pair should be split by the RNA inertia asymmetry; got " + f"{f0:.6f}, {f1:.6f} Hz" ) + # The split is a property of the model, not of the requested window. + f4 = tower.run(n_modes=4).frequencies[:2] + np.testing.assert_allclose(f4, [f0, f1], rtol=1e-6) with warnings.catch_warnings(record=True) as captured: warnings.simplefilter("always") @@ -267,14 +269,11 @@ def test_iea34_no_degeneracy_warning(): + "\n ".join(str(w.message) for w in runtime_warnings) ) - # First rotated shape should be FA-pure (was ~ 0.52 pre-rotation). - p_fa, _ = _shape_participation(rotated[0]) - assert p_fa > 0.99, ( - f"After rotation, FA-aligned mode has p_FA = {p_fa:.4f} " - f"(expected > 0.99). Pre-rotation values were ~ 0.522 / 0.478." - ) - # And the second SS-pure. - _, p_ss = _shape_participation(rotated[1]) + # One of the first two modes is FA-pure and the other SS-pure. + parts = [_shape_participation(s) for s in rotated[:2]] + p_fa = max(p[0] for p in parts) + p_ss = max(p[1] for p in parts) + assert p_fa > 0.99, f"FA-aligned mode has p_FA = {p_fa:.4f}" assert p_ss > 0.99, f"SS-aligned mode has p_SS = {p_ss:.4f}" diff --git a/tests/test_elastodyn_reader.py b/tests/test_elastodyn_reader.py index 3b8f858..c6a42b8 100644 --- a/tests/test_elastodyn_reader.py +++ b/tests/test_elastodyn_reader.py @@ -410,12 +410,15 @@ def test_upscale25_tower_duplicate_pair_stations(): assert np.all(f > 0.0), f"non-physical zero modes returned: {f}" assert np.all(np.diff(f) >= -1e-9) - # 25 MW on a 167 m flexible tower — 1st bending pair near 0.27 Hz. - # FA and SS stiffness are identical in the deck so the pair is - # numerically degenerate. + # 25 MW on a 167 m flexible tower, 1st bending pair near 0.27 Hz. FA + # and SS stiffness are identical in the deck, but the lumped RNA's + # pitch inertia is three times its roll inertia (plus an ixz product), + # which splits the pair by about 1 %. The dense and sparse solves + # agree on that split; the old dense path returned it degenerate only + # because its error on the lowest modes was larger than the split. assert 0.20 < f[0] < 0.35, f"1st FA/SS pair out of band: {f[0]:.4f} Hz" - assert abs(f[1] - f[0]) / f[0] < 1e-3, ( - f"1st FA/SS pair should be degenerate: {f[0]:.4f}, {f[1]:.4f} Hz" + assert abs(f[1] - f[0]) / f[0] < 2e-2, ( + f"1st FA/SS pair should stay close: {f[0]:.4f}, {f[1]:.4f} Hz" ) diff --git a/tests/test_floating_weight_restoring.py b/tests/test_floating_weight_restoring.py new file mode 100644 index 0000000..68d6153 --- /dev/null +++ b/tests/test_floating_weight_restoring.py @@ -0,0 +1,289 @@ +"""Weight restoring on floating models built from OpenFAST decks. + +A WAMIT ``.hst`` holds only the hydrostatic restoring ``ρ g I_wp + +ρ g V z_B``; OpenFAST adds the weight of the platform, tower and RNA +through ElastoDyn's gravity. The classical restoring of a floating body +is ``C44 = C55 = ρ g V (z_B − z_G) + ρ g I_wp`` (Faltinsen 1990, *Sea +Loads on Ships and Offshore Structures*, ch. 2), so a model assembled +from the decks has to add ``−m g z_G`` back. Without it the OC3 Hywind +spar, which is stabilised by ballast rather than by its waterplane, read +``C44 = C55 = −5.0e9 N·m/rad`` and its roll and pitch modes silently +disappeared from the result. + +The end-to-end reference is the OC3 Hywind full-system linearisation in +Jonkman (2010), *Definition of the Floating System for Phase IV of OC3*, +NREL/TP-500-47535, full-system natural frequencies from the FAST +linearisation: surge / sway 0.0080 Hz, heave 0.0324 Hz, roll 0.0342 Hz, +pitch 0.0343 Hz and yaw 0.1210 Hz — the last with the additional +9.834e7 N·m/rad yaw spring the report's mooring-system definition uses to +stand in for the delta-line crowfoot. The r-test MoorDyn deck models +plain catenaries without it. +""" + +from __future__ import annotations + +import copy +import pathlib +import warnings +from types import SimpleNamespace + +import numpy as np +import pytest + +from pybmodes.io.bmi import PlatformSupport, TipMassProps +from pybmodes.models._platform import STANDARD_GRAVITY, _gravitational_restoring + +REPO_ROOT = pathlib.Path(__file__).resolve().parents[1] +SAMPLES = ( + REPO_ROOT / "src" / "pybmodes" / "_examples" / "sample_inputs" + / "reference_turbines" +) +OC3_SAMPLE = ( + SAMPLES / "07_nrel5mw_oc3hywind_spar" + / "07_nrel5mw_oc3hywind_spar_tower.bmi" +) +OC3_DECKS = ( + REPO_ROOT / "external" / "OpenFAST_files" / "r-test" / "glue-codes" + / "openfast" / "5MW_OC3Spar_DLL_WTurb_WavesIrr" +) +OC3_ELASTODYN = OC3_DECKS / "NRELOffshrBsline5MW_OC3Hywind_ElastoDyn.dat" +OC3_MOORDYN = OC3_DECKS / "NRELOffshrBsline5MW_OC3Hywind_MoorDyn.dat" +OC3_HYDRODYN = OC3_DECKS / "NRELOffshrBsline5MW_OC3Hywind_HydroDyn.dat" + +# Jonkman (2010) NREL/TP-500-47535, FAST full-system natural frequencies. +JONKMAN_2010_HZ = { + "surge": 0.0080, "sway": 0.0080, "heave": 0.0324, + "roll": 0.0342, "pitch": 0.0343, "yaw": 0.1210, +} +# Jonkman (2010) mooring-system definition: additional yaw spring that +# stands in for the delta-line crowfoot. +OC3_ADDITIONAL_YAW_STIFFNESS = 9.834e7 + +G = STANDARD_GRAVITY + + +def _support(**overrides) -> PlatformSupport: + fields = dict( + draft=-10.0, cm_pform=90.0, mass_pform=7.0e6, + i_matrix=np.diag([7.0e6, 7.0e6, 7.0e6, 4.0e9, 4.0e9, 1.6e8]), + ref_msl=0.0, hydro_M=np.zeros((6, 6)), hydro_K=np.zeros((6, 6)), + mooring_K=np.zeros((6, 6)), + distr_m_z=np.zeros(0), distr_m=np.zeros(0), + distr_k_z=np.zeros(0), distr_k=np.zeros(0), + ) + fields.update(overrides) + return PlatformSupport(**fields) + + +def _model(support: PlatformSupport, *, mass_den: float = 3000.0, + tip: float = 3.5e5, cm_axial: float = 2.0, + sec_mass: float = 1.0, point_masses=()): + """A hand-built BMI / section-property pair: uniform tower from z = 10 + to z = 90 above MSL, with the floating ``radius = tower-top z`` + convention the deck constructors use.""" + bmi = SimpleNamespace( + support=support, radius=90.0, + scaling=SimpleNamespace(sec_mass=sec_mass), + tip_mass=TipMassProps( + mass=tip, cm_offset=0.0, cm_axial=cm_axial, + ixx=0.0, iyy=0.0, izz=0.0, ixy=0.0, izx=0.0, iyz=0.0, + ), + point_masses=tuple(point_masses), + ) + sp = SimpleNamespace( + span_loc=np.linspace(0.0, 1.0, 11), + mass_den=np.full(11, mass_den), + ) + return bmi, sp + + +class TestClosedForm: + """``−g Σ m_i (z_i − z_ref)`` evaluated by hand.""" + + def test_platform_tower_and_rna(self): + bmi, sp = _model(_support()) + c = _gravitational_restoring(bmi, sp) + # Platform CM 90 m below MSL; tower 80 m long, centroid at 50 m; + # RNA 2 m above the 90 m tower top. + moment = 7.0e6 * -90.0 + 3000.0 * 80.0 * 50.0 + 3.5e5 * 92.0 + assert c[3, 3] == pytest.approx(-G * moment, rel=1e-12) + assert c[4, 4] == pytest.approx(-G * moment, rel=1e-12) + mask = np.ones((6, 6), dtype=bool) + mask[3, 3] = mask[4, 4] = False + assert np.all(c[mask] == 0.0) + + def test_ballast_below_stabilises_and_topside_destabilises(self): + deep, sp = _model(_support(), mass_den=0.0, tip=0.0) + assert _gravitational_restoring(deep, sp)[3, 3] > 0.0 + top, sp = _model(_support(mass_pform=0.0)) + assert _gravitational_restoring(top, sp)[3, 3] < 0.0 + + def test_lever_is_taken_about_the_reference_point(self): + """The ``.hst`` is about ``ref_msl``; so must the weight be. Moving + the reference down by ``d`` adds ``g · m_total · d``.""" + bmi0, sp = _model(_support()) + bmi5, _ = _model(_support(ref_msl=5.0)) + m_total = 7.0e6 + 3000.0 * 80.0 + 3.5e5 + delta = (_gravitational_restoring(bmi5, sp)[3, 3] + - _gravitational_restoring(bmi0, sp)[3, 3]) + assert delta == pytest.approx(-G * m_total * 5.0, rel=1e-9) + + def test_mass_scaling_and_point_masses_count(self): + from pybmodes.io.bmi import PointMass + + base, sp = _model(_support()) + scaled, _ = _model(_support(), sec_mass=1.1, + point_masses=(PointMass(height=40.0, mass=2.0e4),)) + delta = (_gravitational_restoring(scaled, sp)[3, 3] + - _gravitational_restoring(base, sp)[3, 3]) + # 10 % more tower (centroid 50 m) and 20 t at 10 + 40 = 50 m. + assert delta == pytest.approx( + -G * (0.1 * 3000.0 * 80.0 * 50.0 + 2.0e4 * 50.0), rel=1e-9, + ) + + +def _oc3_sample(): + from pybmodes.models import Tower + + return Tower(OC3_SAMPLE) + + +class TestPointMassKeepsWeightInStep: + """``add_point_mass`` after construction must reach the weight term too, + or the FEM mass and the roll / pitch restoring disagree.""" + + def test_lump_adds_its_own_weight(self): + tower = _oc3_sample() + tower._weight_g = G # as the deck-built constructors set it + ps = tower._bmi.support + before = np.array(ps.hydro_K, copy=True) + tower.add_point_mass(30.0, 5.0e4) + z_rel = -ps.draft + 30.0 + ps.ref_msl + delta = ps.hydro_K - before + assert delta[3, 3] == pytest.approx(-G * 5.0e4 * z_rel, rel=1e-12) + assert delta[4, 4] == delta[3, 3] + mask = np.ones((6, 6), dtype=bool) + mask[3, 3] = mask[4, 4] = False + assert np.all(delta[mask] == 0.0) + + def test_model_without_a_weight_term_is_untouched(self): + tower = _oc3_sample() + before = np.array(tower._bmi.support.hydro_K, copy=True) + tower.add_point_mass(30.0, 5.0e4) + np.testing.assert_array_equal(tower._bmi.support.hydro_K, before) + + +class TestNegativeRestoringIsNeverSilent: + """A negative-stiffness rigid-body mode must be reported, not dropped.""" + + def test_stable_platform_is_quiet(self): + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + res = _oc3_sample().run(n_modes=9, check_model=False) + assert set(res.mode_labels[:6]) == { + "surge", "sway", "heave", "roll", "pitch", "yaw", + } + + def test_unstable_roll_and_pitch_warn(self): + tower = _oc3_sample() + ps = copy.deepcopy(tower._bmi.support) + # Strip the weight the BModes deck folds into mooring_K, leaving + # hydrostatics + catenary only — the pre-fix deck-built state. + ps.mooring_K[3, 3] -= 7.46633e6 * G * 89.9155 + ps.mooring_K[4, 4] -= 7.46633e6 * G * 89.9155 + tower._bmi.support = ps + with pytest.warns(UserWarning, match=r"negative-stiffness.*roll, pitch"): + res = tower.run(n_modes=9, check_model=False) + assert "roll" not in res.mode_labels + assert "pitch" not in res.mode_labels + + +_needs_oc3 = pytest.mark.skipif( + not (OC3_ELASTODYN.is_file() and OC3_MOORDYN.is_file() + and OC3_HYDRODYN.is_file()), + reason=f"OC3 Hywind r-test decks not present under {OC3_DECKS}", +) + + +def _oc3_from_decks(): + from pybmodes.models import Tower + + return Tower.from_elastodyn_with_mooring( + OC3_ELASTODYN, OC3_MOORDYN, OC3_HYDRODYN, + ) + + +def _by_label(res) -> dict[str, float]: + return { + lbl: float(f) + for lbl, f in zip(res.mode_labels, res.frequencies) + if lbl is not None + } + + +@pytest.mark.integration +@_needs_oc3 +class TestOC3HywindFromDecks: + def test_all_six_rigid_body_modes_are_returned(self): + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + res = _oc3_from_decks().run(n_modes=9, check_model=False) + assert set(res.mode_labels[:6]) == set(JONKMAN_2010_HZ) + + def test_weight_is_added_to_the_hst(self): + from pybmodes.io.wamit_reader import HydroDynReader + + tower = _oc3_from_decks() + c_hst = HydroDynReader(OC3_HYDRODYN).read_platform_matrices().C_hst + c_w = tower._bmi.support.hydro_K - c_hst + # Platform ballast (+6.58e9) outweighs tower + RNA (−4.1e8). + assert c_hst[3, 3] == pytest.approx(-4.999e9, rel=1e-3) + assert c_w[3, 3] == pytest.approx(6.171e9, rel=1e-3) + assert c_w[4, 4] == c_w[3, 3] + + @pytest.mark.parametrize( + "dof, rel", + [("surge", 0.03), ("sway", 0.03), ("heave", 0.02), + ("roll", 0.03), ("pitch", 0.03)], + ) + def test_matches_jonkman_2010(self, dof, rel): + f = _by_label(_oc3_from_decks().run(n_modes=9, check_model=False)) + assert f[dof] == pytest.approx(JONKMAN_2010_HZ[dof], rel=rel) + + def test_yaw_with_the_published_additional_spring(self): + """The deck carries no crowfoot spring, so yaw is the catenary-only + 0.041 Hz. Adding the report's spring brings it to the published + value; the remaining ~5 % is the RNA yaw inertia, which the + ElastoDyn adapter assembles from a slender-body nacelle proxy.""" + tower = _oc3_from_decks() + assert _by_label(tower.run(n_modes=9, check_model=False))["yaw"] == ( + pytest.approx(0.0412, rel=0.02) + ) + tower._bmi.support.mooring_K[5, 5] += OC3_ADDITIONAL_YAW_STIFFNESS + f = _by_label(tower.run(n_modes=9, check_model=False)) + assert f["yaw"] == pytest.approx(JONKMAN_2010_HZ["yaw"], rel=0.08) + + def test_independent_of_the_requested_mode_count(self): + tower = _oc3_from_decks() + f9 = tower.run(n_modes=9, check_model=False).frequencies + f20 = tower.run(n_modes=20, check_model=False).frequencies + np.testing.assert_allclose(f20[:9], f9, rtol=1e-6) + + def test_point_mass_matches_a_full_recompute(self): + from pybmodes.io.wamit_reader import HydroDynReader + + tower = _oc3_from_decks().add_point_mass(40.0, 1.0e5) + c_hst = HydroDynReader(OC3_HYDRODYN).read_platform_matrices().C_hst + np.testing.assert_allclose( + tower._bmi.support.hydro_K - c_hst, + _gravitational_restoring(tower._bmi, tower._sp), + rtol=1e-12, atol=1e-3, + ) + + def test_mooring_only_model_gets_no_weight_term(self): + """Without HydroDyn there is no buoyancy to balance the weight, so + adding it would not describe a floating body.""" + from pybmodes.models import Tower + + tower = Tower.from_elastodyn_with_mooring(OC3_ELASTODYN, OC3_MOORDYN) + assert np.all(tower._bmi.support.hydro_K == 0.0) diff --git a/tests/test_numeric.py b/tests/test_numeric.py new file mode 100644 index 0000000..96f0bc7 --- /dev/null +++ b/tests/test_numeric.py @@ -0,0 +1,47 @@ +"""Shared numerical helpers, and the NumPy-version guard they exist for. + +``pyproject.toml`` declares ``numpy>=1.26``. ``np.trapezoid`` only exists +from NumPy 2.0 and ``np.trapz`` is deprecated there, so any direct call +breaks one end of the supported range. CI installs a current NumPy and +would not notice, hence the source scan below. +""" + +from __future__ import annotations + +import pathlib +import re + +import numpy as np +import pytest + +from pybmodes._numeric import trapezoid + +SRC = pathlib.Path(__file__).resolve().parents[1] / "src" / "pybmodes" + + +def test_linear_integrand_is_exact_on_a_nonuniform_grid(): + x = np.array([0.0, 0.3, 1.1, 2.0, 4.5]) + assert trapezoid(3.0 * x + 2.0, x) == pytest.approx( + 1.5 * 4.5 ** 2 + 2.0 * 4.5, rel=1e-14) + + +def test_single_point_integrates_to_zero(): + assert trapezoid([5.0], [1.0]) == 0.0 + + +def test_returns_a_python_float(): + assert type(trapezoid([1, 2, 3], [0, 1, 2])) is float + + +def test_no_version_specific_numpy_integrators_in_the_package(): + pattern = re.compile(r"\bnp\.(trapezoid|trapz)\s*\(") + offenders = [ + f"{path.relative_to(SRC)}:{n}" + for path in SRC.rglob("*.py") + for n, line in enumerate(path.read_text(encoding="utf-8").splitlines(), 1) + if pattern.search(line) + ] + assert not offenders, ( + "use pybmodes._numeric.trapezoid instead (numpy>=1.26 support): " + + ", ".join(offenders) + ) diff --git a/tests/test_review_hardening.py b/tests/test_review_hardening.py index 88098d6..2be3756 100644 --- a/tests/test_review_hardening.py +++ b/tests/test_review_hardening.py @@ -106,6 +106,30 @@ def test_discover_skips_bad_sidecar_and_finds_valid_ontology( assert disc.yaml.name == "good.yaml" +def test_discover_ignores_openfast_in_the_path_above_the_search_root( + tmp_path: pathlib.Path, +) -> None: + """Regression: the OpenFAST-deck exclusion matched the absolute path, + so a checkout living under a directory such as ``OpenFAST-GUI`` + excluded every yaml in it and discovery found nothing. Only the part + of the path below the search root may be examined — an ``OpenFAST`` + deck directory *inside* the tree is still skipped.""" + pytest.importorskip("yaml") + from pybmodes.workflows.windio import discover_windio_inputs + + root = tmp_path / "OpenFAST-GUI" / "turbine" + (root / "OpenFAST").mkdir(parents=True) + _write(root / "ontology.yaml", "name: ok\ncomponents:\n tower: {}\n") + # An ontology-shaped yaml under the turbine's own OpenFAST tree sorts + # first by name, so it would win if the exclusion stopped working. + _write( + root / "OpenFAST" / "a_sidecar.yaml", + "name: deck\ncomponents:\n tower: {}\n", + ) + disc = discover_windio_inputs(root) + assert disc.yaml == root / "ontology.yaml" + + def test_load_windio_doc_propagates_missing_pyyaml( tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch, ) -> None: diff --git a/tests/test_windio_floating.py b/tests/test_windio_floating.py index 4518bce..a2b445d 100644 --- a/tests/test_windio_floating.py +++ b/tests/test_windio_floating.py @@ -763,6 +763,34 @@ def test_from_windio_floating_yaml_only_modal_smoke(tmp_path) -> None: assert f[6] > 5.0 * f[5] +def test_from_windio_floating_adds_weight_to_the_hydrostatics(tmp_path) -> None: + """The member-waterplane restoring is buoyancy-only, like a WAMIT + ``.hst``; the assembled ``hydro_K`` must add ``−m g z`` for the + platform, tower and RNA exactly as ``from_elastodyn_with_mooring`` + does (Faltinsen 1990, *Sea Loads on Ships and Offshore Structures*, + ch. 2: ``C44 = ρ g V (z_B − z_G) + ρ g I_wp``).""" + pytest.importorskip("yaml") + from pybmodes.io.windio_floating import ( + hydrostatic_restoring, + read_windio_floating, + ) + from pybmodes.models import Tower + from pybmodes.models._platform import _gravitational_restoring + + p = tmp_path / "fowt.yaml" + p.write_text(_FLOAT_TURBINE, encoding="utf-8") + with pytest.warns(UserWarning, match="SCREENING-fidelity"): + tower = Tower.from_windio_floating(p, water_depth=200.0) + c_hst = hydrostatic_restoring(read_windio_floating(p)) + c_w = _gravitational_restoring(tower._bmi, tower._sp) + np.testing.assert_allclose( + tower._bmi.support.hydro_K, c_hst + c_w, rtol=1e-12, atol=1e-6, + ) + # A 20 m-deep column carrying a 128 m tower: the weight term is + # destabilising and material against the waterplane term. + assert c_w[3, 3] < -0.05 * c_hst[3, 3] + + def test_from_windio_floating_screening_honors_rna_tip(tmp_path) -> None: """Regression (issue #83): in the screening tier (no companion ElastoDyn deck, no injected ``platform_support``) the caller-