From dd0b21d52eaa2fd645c055f3235b30df5be9a1fe Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Tue, 29 Sep 2026 19:40:26 +0900 Subject: [PATCH 1/5] fix: solve the inverted pencil on the dense symmetric path The dense path reduced K x = lambda M x through a Cholesky factor of M, and a backward-stable symmetric eigensolver returns each eigenvalue of the reduced matrix only to an absolute accuracy of eps * lambda_max. The ElastoDyn adapters floor rotary inertia at 1e-6, which puts lambda_max / lambda_min near 1e15, so the lowest modes carried errors of several percent that moved with the subset_by_index window: the NREL 5MW blade's first flap at 12.1 rpm read 0.7154 / 0.7096 / 0.7661 / 0.7282 Hz for n_modes 4 / 6 / 10 / 20. Solve M x = mu (K + s M) x instead, the shift-invert transformation the sparse path already applies. The wanted modes sit at the top of that spectrum, where the same absolute accuracy is relative. The shift is zero for a definite K and otherwise grown from an estimate of the lowest eigenvalues until K + s M factorises, so rigid-body modes and negative eigenvalues are still returned. Requests reaching into the upper half of the spectrum take those modes from the mass-reduced form, which is the accurate one there. The blade now gives 0.72865 Hz at every window, matching the sparse path; backward error drops from ~5e-2 to ~1e-10. The mass-reduced form stays as the fallback, with the residual retry behind it; the retry tests run against that route explicitly. --- CHANGELOG.md | 40 +++- VALIDATION.md | 5 +- src/pybmodes/fem/solver.py | 208 ++++++++++++++++++- src/pybmodes/options.py | 33 +-- tests/fem/test_dense_inverted_solver.py | 250 +++++++++++++++++++++++ tests/fem/test_gravity_and_point_mass.py | 13 +- tests/fem/test_ill_conditioned_mass.py | 54 +++-- tests/test_campbell.py | 16 +- 8 files changed, 563 insertions(+), 56 deletions(-) create mode 100644 tests/fem/test_dense_inverted_solver.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 23ff0d0..10d3c5c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,45 @@ 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 %. + + 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 %. ## [1.18.0] — 2026-08-13 diff --git a/VALIDATION.md b/VALIDATION.md index ff2a0c9..25c0afa 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,7 @@ 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). - 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). diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 8e6a582..656c9a3 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,195 @@ 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): + rough = eigh( + gk_sym, gm_sym, eigvals_only=True, + subset_by_index=(0, min(k, _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 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) + + +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/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/tests/fem/test_dense_inverted_solver.py b/tests/fem/test_dense_inverted_solver.py new file mode 100644 index 0000000..ff5da28 --- /dev/null +++ b/tests/fem/test_dense_inverted_solver.py @@ -0,0 +1,250 @@ +"""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, + ) + + +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] From 110b073f4b6f63e8ca5a8e46a72c90621c1cd21e Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Tue, 29 Sep 2026 19:40:38 +0900 Subject: [PATCH 2/5] fix: add the weight restoring to deck-built floating models Tower.from_elastodyn_with_mooring took roll / pitch restoring from the WAMIT .hst alone, which holds only rho g I_wp + rho g V z_B. OpenFAST adds the body weight through ElastoDyn gravity; the model did not, so the ballast-stabilised OC3 Hywind spar read C44 = C55 = -5.0e9 N m/rad and the general eigensolver silently dropped both negative-eigenvalue modes, returning four of the six rigid-body modes. Add -g sum m_i (z_i - z_ref) for the platform, tower and RNA to hydro_K (+6.17e9 N m/rad on OC3), measured from the hydro/mooring reference in the MSL-datum conventions, and only when HydroDyn is given, since weight without buoyancy is not a floating body. from_windio_floating adds the same term in both tiers so the deck-backed tier stays equivalent. OC3 Hywind from the r-test decks now returns all six rigid-body modes: surge / sway / heave / roll / pitch / yaw = 0.00806 / 0.00806 / 0.03241 / 0.03383 / 0.03383 / 0.04121 Hz against Jonkman (2010) NREL/TP-500-47535 0.0080 / 0.0080 / 0.0324 / 0.0342 / 0.0343 Hz; yaw is catenary-only, the published 0.1210 Hz including an extra crowfoot spring the deck lacks. Tower.run now warns when hydro_K + mooring_K has a negative direction instead of letting the solver drop those modes unannounced. --- CHANGELOG.md | 31 +++ VALIDATION.md | 4 + docs/conventions.rst | 11 ++ src/pybmodes/io/bmi.py | 15 +- src/pybmodes/io/windio_floating.py | 6 +- src/pybmodes/models/_pipeline.py | 72 +++++++ src/pybmodes/models/_platform.py | 85 +++++++- src/pybmodes/models/tower.py | 29 +++ tests/test_floating_weight_restoring.py | 253 ++++++++++++++++++++++++ tests/test_windio_floating.py | 28 +++ 10 files changed, 528 insertions(+), 6 deletions(-) create mode 100644 tests/test_floating_weight_restoring.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 10d3c5c..8092db3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,6 +47,37 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). `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. ## [1.18.0] — 2026-08-13 diff --git a/VALIDATION.md b/VALIDATION.md index 25c0afa..2d36ce7 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -117,6 +117,7 @@ by title + publisher: - 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). @@ -226,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/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_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..325f453 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,83 @@ 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.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 += float(np.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..e672b08 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, ) @@ -902,6 +903,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,6 +1069,18 @@ 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, + ) + obj = cls.__new__(cls) obj._bmi = bmi obj._sp = sp @@ -1434,6 +1456,13 @@ 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 diff --git a/tests/test_floating_weight_restoring.py b/tests/test_floating_weight_restoring.py new file mode 100644 index 0000000..44db8cc --- /dev/null +++ b/tests/test_floating_weight_restoring.py @@ -0,0 +1,253 @@ +"""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 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_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_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- From 5266d4fdd863c0a95d8de5c7b8e9b933f3b1bffd Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Tue, 29 Sep 2026 19:40:38 +0900 Subject: [PATCH 3/5] fix: scope the WindIO OpenFAST-tree exclusion to the search root discover_windio_inputs skipped any yaml whose absolute path contained "OpenFAST", so a checkout living under a directory such as OpenFAST-GUI found no ontology at all. Examine only the part of the path below the directory being searched; an OpenFAST deck directory inside the tree is still skipped. --- CHANGELOG.md | 5 +++++ src/pybmodes/workflows/windio.py | 18 +++++++++++++++++- tests/test_review_hardening.py | 24 ++++++++++++++++++++++++ 3 files changed, 46 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8092db3..fd82a09 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -78,6 +78,11 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). 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. +- **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/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/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: From 0a2f234d5298727165913a6c33b21e2fa0b71c21 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 30 Sep 2026 16:18:05 +0900 Subject: [PATCH 4/5] fix: keep NumPy 1.26 working and the weight term in step with point masses np.trapezoid only exists from NumPy 2.0, but the package declares numpy>=1.26, so the pre-solve checks, the ElastoDyn blade reader and the floating weight term all raised AttributeError there. They now share one version-independent helper, and a test keeps the NumPy 2 spelling out of the package. Tower.add_point_mass on a deck-built floater now adds the lump's weight to the roll and pitch restoring, so it matches a model rebuilt with the lump in place. Two integration tests pinned an exactly degenerate first tower pair on the IEA-3.4 and UPSCALE 25MW decks. Their RNA roll and pitch inertias differ, which splits the pair by 0.09 % and 1.0 %; the dense, sparse and inverted-pencil solves agree on it. The old degeneracy was solver error, so the tests now pin the split. --- CHANGELOG.md | 20 ++++++++ .../sample_inputs/reference_turbines/build.py | 3 +- src/pybmodes/_numeric.py | 35 ++++++++++++++ src/pybmodes/checks.py | 5 +- src/pybmodes/io/_elastodyn/adapter.py | 3 +- src/pybmodes/io/_elastodyn/types.py | 4 +- src/pybmodes/io/windio.py | 14 +----- src/pybmodes/models/_platform.py | 3 +- src/pybmodes/models/tower.py | 24 +++++++++- tests/test_classifier.py | 45 +++++++++--------- tests/test_elastodyn_reader.py | 13 +++-- tests/test_floating_weight_restoring.py | 36 ++++++++++++++ tests/test_numeric.py | 47 +++++++++++++++++++ 13 files changed, 205 insertions(+), 47 deletions(-) create mode 100644 src/pybmodes/_numeric.py create mode 100644 tests/test_numeric.py diff --git a/CHANGELOG.md b/CHANGELOG.md index fd82a09..e9d4f99 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -78,6 +78,26 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). 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 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/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/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/models/_platform.py b/src/pybmodes/models/_platform.py index 325f453..877caca 100644 --- a/src/pybmodes/models/_platform.py +++ b/src/pybmodes/models/_platform.py @@ -194,6 +194,7 @@ def _gravitational_restoring( """ import numpy as np + from pybmodes._numeric import trapezoid from pybmodes.io.bmi import PlatformSupport ps = bmi.support @@ -211,7 +212,7 @@ def _gravitational_restoring( 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 += float(np.trapezoid(mass_den * (z - ref_z), 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 diff --git a/src/pybmodes/models/tower.py b/src/pybmodes/models/tower.py index e672b08..63d41b7 100644 --- a/src/pybmodes/models/tower.py +++ b/src/pybmodes/models/tower.py @@ -212,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: @@ -1078,12 +1084,14 @@ def from_elastodyn_with_mooring( # 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, + 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 @@ -1467,6 +1475,7 @@ def _tower_grid_od_wt( obj = cls.__new__(cls) obj._bmi = bmi obj._sp = sp + obj._weight_g = float(g) obj.coeff_validation = None return obj @@ -1761,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)) @@ -1777,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/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 index 44db8cc..68d6153 100644 --- a/tests/test_floating_weight_restoring.py +++ b/tests/test_floating_weight_restoring.py @@ -148,6 +148,31 @@ def _oc3_sample(): 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.""" @@ -244,6 +269,17 @@ def test_independent_of_the_requested_mode_count(self): 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.""" 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) + ) From 882d4fb2b7a690b6a12998c153da3d0447c4a290 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 30 Sep 2026 19:24:04 +0900 Subject: [PATCH 5/5] =?UTF-8?q?fix:=20recover=20soft=20modes=20a=20large?= =?UTF-8?q?=20indefinite=20shift=20pushes=20below=20=CE=BC=20resolution?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With one strongly negative eigenvalue the definite shift is so large that every soft mode maps to the same μ, so their order was arbitrary and the result depended on n_modes: K = diag(-1e16, 1, 2, 3, 4) gave [-1e16, 4] for two modes. The window now grows until its edge falls in a resolvable gap, and each swamped cluster is re-solved by Rayleigh-Ritz on its own subspace. The shift estimate no longer depends on the requested window either. --- CHANGELOG.md | 8 ++ src/pybmodes/fem/solver.py | 100 +++++++++++++++++++++++- tests/fem/test_dense_inverted_solver.py | 28 +++++++ 3 files changed, 135 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e9d4f99..83b7465 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,14 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). 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 diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 656c9a3..0815a7c 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -939,9 +939,11 @@ def _definite(shift: float) -> bool: 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(k, _SHIFT_SCALE_MODES) - 1), + 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: @@ -960,6 +962,10 @@ def _definite(shift: float) -> bool: ) 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( @@ -987,6 +993,98 @@ def _inverted_window( 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, diff --git a/tests/fem/test_dense_inverted_solver.py b/tests/fem/test_dense_inverted_solver.py index ff5da28..255858c 100644 --- a/tests/fem/test_dense_inverted_solver.py +++ b/tests/fem/test_dense_inverted_solver.py @@ -225,6 +225,34 @@ def test_negative_eigenvalues_are_kept(self): 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):