Repository navigation
fix: n_modes-independent dense solves, weight restoring on deck-built floaters, WindIO discovery path scope - #142
Conversation
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.
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.
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.
|
You have reached your Codex usage limits for code reviews. You can see your limits in the Codex usage dashboard. |
|
@codex review |
Codex Review SummaryThis comment shows the latest Codex review activity on this pull request.
ℹ️ About Codex in GitHubYour team has set up Codex to review pull requests in this repo. Reviews are triggered when you
Codex reacts with 👀 while any review is running, comments if it has suggestions, and reacts with 👍 once all reviews finish with no findings. |
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 5266d4fdd8
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| 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)) |
There was a problem hiding this comment.
The project declares numpy>=1.26, but np.trapezoid was introduced in NumPy 2.0. On a supported NumPy 1.26 installation, every deck-backed or WindIO floating constructor that calls _gravitational_restoring now raises AttributeError before returning a model. Use the 1.26-compatible np.trapz or raise the minimum dependency to NumPy 2.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Fixed in 0a2f234. This was wider than the new weight term, since checks.py and the ElastoDyn blade reader on master also called np.trapezoid. Every call site now goes through one version-independent helper in pybmodes._numeric, and tests/test_numeric.py scans the package so the NumPy 2 spelling can't creep back in.
| for pm in bmi.point_masses: | ||
| moment += float(pm.mass) * (base_z + float(pm.height) - ref_z) |
There was a problem hiding this comment.
For a floating model constructed from HydroDyn or WindIO and then modified through the public Tower.add_point_mass() API, this loop never sees the new mass because _gravitational_restoring is evaluated only during construction. The FEM mass matrix includes the added lump, but hydro_K retains the old roll/pitch weight term, producing inconsistent and potentially incorrect rigid-body frequencies; recompute the weight contribution when masses change or immediately before solving.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Fixed in 0a2f234. add_point_mass now adds the lump's weight to the roll and pitch restoring on deck-built floaters, measured the same way as the constructor's term. The OC3 integration test checks it against a full recompute.
…asses 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.
|
@codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 0a2f234d52
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| mu, vecs = eigh( | ||
| gm_sym, gk_sym + shift * gm_sym, | ||
| subset_by_index=(ngd - k, ngd - 1), |
There was a problem hiding this comment.
When an indefinite pencil has a negative eigenvalue much larger in magnitude than later soft modes, the required shift makes their transformed mu values indistinguishable at floating-point precision, so this subset can silently select the wrong eigenpairs. For example, K = diag(-1e16, 1, 2, 3, 4), M = I, and n_modes=2 returns [-1e16, 4], while requesting three modes returns the correct [-1e16, 1, 2]. The residual retry cannot detect this because the selected fourth mode is itself an exact eigenpair with zero residual, leaving the public result dependent on n_modes and omitting lower modes; the shifted path therefore needs ordering verification or a fallback when the shift overwhelms the requested soft spectrum.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Fixed in 882d4fb. When a shifted pencil leaves any mode in the window below the μ resolution, the window now grows until its edge sits in a resolvable gap, and each swamped cluster is re-solved by Rayleigh-Ritz on its own subspace. Your example returns [-1e16, 1, 2, 3, 4][:n] for every n. The shift estimate also no longer depends on the requested window, so a rotated 40-DOF version with two eigenvalues near -1e14 gives bit-identical soft modes for 3 and 9 requested modes. Both cases are in tests/fem/test_dense_inverted_solver.py.
…ution 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.
|
@codex review |
|
Codex Review: Didn't find any major issues. 🎉 Reviewed commit: ℹ️ About Codex in GitHubYour team has set up Codex to review pull requests in this repo. Reviews are triggered when you
If Codex has suggestions, it will comment; otherwise it will react with 👍. Codex can also answer questions or update the PR. Try commenting "@codex address that feedback". |
1 participant
Three defects found while integrating pyBModes 1.18.0 into a downstream engine. There is one commit for each.
1. Dense symmetric solves depended on
n_modesRoot cause. Below the 500-DOF sparse threshold, the dense path solved
eigh(K, M, subset_by_index=...). That call reduces the pencil through a Cholesky factor ofM. A backward-stable symmetric eigensolver gets each eigenvalue of the reduced matrix only to an absolute accuracy of abouteps · λ_max, so the relative error on the lowest mode is abouteps · λ_max / λ_min. The ElastoDyn adapters floor rotary inertia at 1e-6, which puts that ratio near 1e15 (cond(M) ≈ 2e11 on the blade). The low modes therefore carried errors of several percent, and those errors moved with the requested window. The problem is the spectral range, not a bad DOF or missing scaling: diagonal equilibration leaves the eigenvalues unchanged, and lowering the threshold would only have hidden the defect behindeigsh.Fix. Solve the inverted pencil
M x = μ (K + s M) x, withμ = 1/(λ + s). This is the shift-invert transformation the sparse path already uses (Ericsson & Ruhe 1980). It puts the wanted modes at the top of the spectrum, where absolute accuracy becomes relative accuracy.sis 0 whenKis positive definite.sstarts from an estimate of the lowest eigenvalues and grows untilK + s Mfactorises. Rigid-body modes and the negative eigenvalues of a column loaded past buckling are still returned.test_ill_conditioned_mass.pynow runs against that route explicitly.NREL 5MW (bundled deck), first modes in Hz:
Backward error drops from ~5e-2 to ~1e-10.
test_campbell::test_first_flap_lifts_endpoint_to_endpointsat at 2.95 % against its 3 % gate on one BLAS build; it now measures 7.6 %.2. Floating models built from decks lost roll and pitch
Root cause.
Tower.from_elastodyn_with_mooringsethydro_Kto the WAMIT.hst. That matrix holds onlyρ g I_wp + ρ g V z_Band leaves out the body weight, which OpenFAST applies through ElastoDyn gravity. For the ballast-stabilised OC3 spar this givesC44 = C55 = −5.0e9 N·m/rad. The general eigensolver then filtered out the two negative eigenvalues and backfilled the window, so 4 of the 6 rigid-body modes came back with no warning.Fix. The weight restoring
−g Σ mᵢ (zᵢ − z_ref)of the platform, tower (includingsec_massscaling and point masses) and RNA is now added tohydro_K.draft,cm_pform,ref_msl).from_windio_floatingadds the same term in both tiers, so the deck-backed tier stays equivalent tofrom_elastodyn_with_mooring.Tower.runnow warns whenhydro_K + mooring_Khas a negative direction, and names the DOFs. The sign test is on the inertia-scaled matrix (Sylvester's law of inertia).OC3 Hywind from the r-test decks (HydroDyn + MoorDyn + ElastoDyn), in Hz:
Yaw is the catenary-only value. The published 0.1210 Hz includes the report's additional 9.834e7 N·m/rad crowfoot spring, which the r-test MoorDyn deck does not carry. With that spring added, pyBModes gives 0.127 Hz; the remaining ~5 % comes from the adapter's RNA yaw-inertia proxy.
Known and now visible. The bundled
08_nrel5mw_oc4semisample and the OC4 r-test deck still have negative roll/pitch restoring. OC4 keeps its ballast in HydroDyn filled members rather than inPtfmMass, and the platform model does not include those members. Both now raise the new warning instead of silently labelling tower-bending modes as roll/pitch. Handling filled-member ballast is left for a follow-up.3. WindIO discovery with "OpenFAST" in the parent path
discover_windio_inputsexcluded any yaml whose absolute path contained "OpenFAST". As a result, a checkout living under e.g.~/OpenFAST-GUI/found no ontology, andtest_review_hardening::test_discover_skips_bad_sidecar...failed there. The check now looks only at the path below the search root.Validation
tests/fem/test_dense_inverted_solver.py:tests/test_floating_weight_restoring.py:VALIDATION.mdrows and citations updated;scripts/audit_validation_claims.pyis clean.ruffandmypyclean.