diff --git a/CHANGELOG.md b/CHANGELOG.md index 38950b3..3da75d0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,22 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ## [Unreleased] -(nothing yet) +### Added + +- **Material / outfitting overrides on the WindIO tower paths (#133).** + `Tower.from_windio` and `Tower.from_windio_with_monopile` accept optional + `E`, `rho`, `nu` and `outfitting_factor` keywords. Each defaults to + `None` (use the ontology value); pass a number to override it, so a + material or outfitting sensitivity sweep runs straight off a WindIO file + without editing the yaml. On the combined monopile+tower path the + override applies to both segments. Existing callers are unaffected. +- **Per-mode generalised mass and stiffness on `ModalResult` (#134).** + A solve now attaches `generalized_mass` (kg) and `generalized_stiffness` + (N/m) arrays, one entry per mode, normalised to unit tip lateral + displacement so a mode's modal mass and stiffness can be read off and + compared against another tool's modal report. `stiffness = (2π·f)² · + mass` by construction; a rigid-body platform mode (tip barely moves) + yields `NaN`. Not serialised (derivable from the mode shapes). ## [1.16.1] — 2026-07-06 diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 4017961..700780d 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -474,6 +474,10 @@ def read_windio_monopile_tower( thickness_interp: str = "linear", n_nodes: int | None = None, water_depth: float | None = None, + E: float | None = None, + rho: float | None = None, + nu: float | None = None, + outfitting_factor: float | None = None, ) -> WindIOMonopileTower: """Reduce the ``monopile`` and ``tower`` components and splice them into one fixed-bottom cantilever (issue #92). @@ -509,6 +513,10 @@ def read_windio_monopile_tower( present; ``None`` with no ontology value keeps the monopile base as the clamp (correct only when the axis already starts at the mudline). + E, rho, nu, outfitting_factor : optional material / outfitting overrides + (issue #133), each ``None`` by default (use the ontology value). + When given, the override is applied to **both** the monopile and the + tower segment before the tube reduction. Raises ------ @@ -532,6 +540,19 @@ def read_windio_monopile_tower( yaml_path, component=component_tower, thickness_interp=thickness_interp, ) + # Optional material / outfitting overrides (issue #133), applied to both + # freshly-read segments so a single value drives a whole-structure + # sensitivity sweep. mp / tw are local dataclasses, so mutate in place. + for _seg in (mp, tw): + if E is not None: + _seg.E = E + if rho is not None: + _seg.rho = rho + if nu is not None: + _seg.nu = nu + if outfitting_factor is not None: + _seg.outfitting_factor = outfitting_factor + # Rigid fixed-base monopiles are clamped at the mudline, not at the # embedded pile tip. When the monopile reference_axis extends below the # seabed (e.g. IEA-15: axis -75 -> +15 with the mudline at -30), the diff --git a/src/pybmodes/models/_pipeline.py b/src/pybmodes/models/_pipeline.py index 96018ba..ddaf71f 100644 --- a/src/pybmodes/models/_pipeline.py +++ b/src/pybmodes/models/_pipeline.py @@ -31,6 +31,55 @@ from pybmodes.models.result import ModalResult +def _generalised_mass_stiffness( + eigvecs: np.ndarray, + gm: np.ndarray, + active: np.ndarray, + freqs_hz: np.ndarray, + ref_mr: float, +) -> tuple[np.ndarray, np.ndarray]: + """Physical modal (generalised) mass in kg and stiffness in N/m per mode, + normalised to unit tip lateral displacement (issue #134). + + The non-dimensional generalised mass ``phi^T gm phi`` re-dimensionalises + by ``ref_mr = RM * radius``; dividing by the squared tip lateral + displacement DOF applies the unit-tip normalisation (the arbitrary + eigenvector scale cancels). Generalised stiffness follows from + ``omega^2 * mass``. The tip node is global DOFs 0-5 with lateral + displacements at ``v_disp`` (1) and ``w_disp`` (3); a mode whose tip + barely moves laterally (a rigid-body platform mode) yields ``NaN``. + """ + active = np.asarray(active) + + def _row(global_dof: int) -> int | None: + hits = np.nonzero(active == global_dof)[0] + return int(hits[0]) if hits.size else None + + iv, iw = _row(1), _row(3) + n = eigvecs.shape[1] + gen_mass = np.full(n, np.nan) + gen_stiff = np.full(n, np.nan) + for i in range(n): + phi = eigvecs[:, i] + gm_form = float(phi @ gm @ phi) + v = float(eigvecs[iv, i]) if iv is not None else 0.0 + w = float(eigvecs[iw, i]) if iw is not None else 0.0 + tip = float(np.hypot(v, w)) + # Scale-aware guard: a torsion / heave / yaw-dominated mode has only + # numerical lateral tip motion (tip ~ eps of the mode's peak), for + # which the unit-tip normalisation is undefined. An exact-zero check + # would let that tiny value through and divide by it, so gate on the + # tip motion being a real fraction of the mode's largest DOF. + col_scale = float(np.max(np.abs(phi))) + if (col_scale <= 0.0 or tip <= 1.0e-8 * col_scale + or not np.isfinite(gm_form) or gm_form <= 0.0): + continue + m = ref_mr * gm_form / (tip * tip) + gen_mass[i] = m + gen_stiff[i] = (2.0 * np.pi * float(freqs_hz[i])) ** 2 * m + return gen_mass, gen_stiff + + @dataclass class _SectionPropsView: """Minimal struct of pre-nondimensionalised section-property @@ -268,7 +317,12 @@ def run_fem( and np.asarray(bmi.support.distr_m).size > 0: ignored.append("distributed added mass (distr_m)") + gen_mass, gen_stiff = _generalised_mass_stiffness( + eigvecs, gm, active, freqs_hz, nd.ref_mr + ) + return ModalResult( frequencies=freqs_hz, shapes=shapes, mode_labels=mode_labels, ignored_physics=tuple(ignored), diagnostics=diagnostics, + generalized_mass=gen_mass, generalized_stiffness=gen_stiff, ) diff --git a/src/pybmodes/models/result.py b/src/pybmodes/models/result.py index 16b211c..d44077e 100644 --- a/src/pybmodes/models/result.py +++ b/src/pybmodes/models/result.py @@ -79,6 +79,15 @@ class ModalResult: mass-matrix conditioning). Attached by :func:`run_fem`; **not** serialised (it describes the live solve run, not the persistent result data) and excluded from equality. Added 1.14.0. + generalized_mass, generalized_stiffness : optional ``(n_modes,)`` arrays + of the physical modal (generalised) mass in kg and stiffness in N/m + for each mode, normalised to **unit tip lateral displacement** so + they are directly comparable to another tool's modal output (e.g. + an OrcaFlex modal report). ``stiffness = (2π·f)² · mass`` by + construction. Most meaningful for bending modes; for a rigid-body + platform mode (tip displacement ≈ 0) the entry is ``NaN``. Attached + by :func:`run_fem` on a live solve; not serialised (derivable from + the shapes). Added 1.17.0 (issue #134). """ # NOTE: field order is the positional constructor ABI for this @@ -99,6 +108,12 @@ class ModalResult: # Transient solve telemetry — excluded from equality and from both # serialisers (it is about the run, not the persisted result). diagnostics: SolverDiagnostics | None = field(default=None, compare=False) + # Per-mode physical generalised (modal) mass and stiffness, normalised + # to unit tip lateral displacement (issue #134). Populated on a live + # solve; excluded from equality and not serialised (derivable from the + # shapes). ``generalized_mass`` in kg, ``generalized_stiffness`` in N/m. + generalized_mass: np.ndarray | None = field(default=None, compare=False) + generalized_stiffness: np.ndarray | None = field(default=None, compare=False) # ------------------------------------------------------------------ # Shared integrity check diff --git a/src/pybmodes/models/tower.py b/src/pybmodes/models/tower.py index c417923..b8357c2 100644 --- a/src/pybmodes/models/tower.py +++ b/src/pybmodes/models/tower.py @@ -326,6 +326,10 @@ def from_windio( n_nodes: int | None = None, lumped_rna_cal: bool = False, rna_angle_units: str = "auto", + E: float | None = None, + rho: float | None = None, + nu: float | None = None, + outfitting_factor: float | None = None, ) -> Tower: """Build a tower (or monopile) model directly from a **WindIO** ontology ``.yaml`` (issue #35). @@ -378,6 +382,13 @@ def from_windio( the WindIO rad/deg ambiguity by magnitude; pass ``"rad"`` or ``"deg"`` to take the file at its word. Ignored unless ``lumped_rna_cal`` is set. + E, rho, nu, outfitting_factor : optional material / outfitting + overrides (issue #133). Each defaults to ``None``, meaning + "use the value from the ontology"; pass a number to override + it, e.g. for a sensitivity sweep straight off a WindIO file + without editing the yaml. ``E`` (Pa), ``rho`` (kg/m^3), ``nu`` + (-) set the isotropic tower-wall material; ``outfitting_factor`` + scales the non-structural mass (see :meth:`from_geometry`). Notes ----- @@ -427,8 +438,12 @@ def from_windio( g.outer_diameter, g.wall_thickness, flexible_length=g.flexible_length, - E=g.E, rho=g.rho, nu=g.nu, - outfitting_factor=g.outfitting_factor, + E=g.E if E is None else E, + rho=g.rho if rho is None else rho, + nu=g.nu if nu is None else nu, + outfitting_factor=( + g.outfitting_factor if outfitting_factor is None else outfitting_factor + ), hub_conn=hub_conn, tip_mass=tip_mass, n_nodes=n_nodes, @@ -447,6 +462,10 @@ def from_windio_with_monopile( water_depth: float | None = None, lumped_rna_cal: bool = False, rna_angle_units: str = "auto", + E: float | None = None, + rho: float | None = None, + nu: float | None = None, + outfitting_factor: float | None = None, ) -> Tower: """Build a combined **monopile + tower** fixed-bottom cantilever from a WindIO ontology ``.yaml`` (issue #92). @@ -494,6 +513,12 @@ def from_windio_with_monopile( rna_angle_units : how the auto-RNA reads ``cone_angle`` / ``uptilt`` when ``lumped_rna_cal=True`` (``"auto"`` / ``"rad"`` / ``"deg"``); see :meth:`from_windio`. Ignored unless ``lumped_rna_cal`` is set. + E, rho, nu, outfitting_factor : optional material / outfitting + overrides (issue #133), each defaulting to ``None`` (use the + ontology value). When given, the override is applied to **both** + the monopile and the tower segment, so a single value drives a + whole-structure sensitivity sweep off a WindIO file. For + per-segment materials, edit the ontology instead. Notes ----- @@ -528,6 +553,7 @@ def from_windio_with_monopile( thickness_interp=thickness_interp, n_nodes=n_nodes, water_depth=water_depth, + E=E, rho=rho, nu=nu, outfitting_factor=outfitting_factor, ) tip = _coerce_tip_mass(tip_mass) bmi = _build_bmi_skeleton( diff --git a/tests/test_material_override_and_modal_mass.py b/tests/test_material_override_and_modal_mass.py new file mode 100644 index 0000000..05758f2 --- /dev/null +++ b/tests/test_material_override_and_modal_mass.py @@ -0,0 +1,209 @@ +"""Material / outfitting overrides on the WindIO paths (issue #133) and +per-mode generalised mass / stiffness on the result (issue #134). + +Self-contained: builds tiny hand-written WindIO ontologies in ``tmp_path`` +and analytical tower geometry inline, no external data. +""" +from __future__ import annotations + +import pathlib +import textwrap + +import numpy as np +import pytest + +from pybmodes.io.bmi import TipMassProps +from pybmodes.models import Tower + +_MIN_TOWER = textwrap.dedent("""\ + components: + tower: + outer_shape: + outer_diameter: + grid: [0.0, 0.5, 1.0] + values: [8.0, 7.0, 6.0] + structure: + outfitting_factor: 1.1 + layers: + - name: tower_wall + material: steel + thickness: + grid: [0.0, 1.0] + values: [0.05, 0.02] + reference_axis: + z: + grid: [0.0, 1.0] + values: [20.0, 120.0] + materials: + - name: steel + E: 2.0e11 + rho: 7800.0 + nu: 0.3 + G: 7.7e10 + """) + +_MIN_MONOPILE_TOWER = textwrap.dedent("""\ + components: + monopile: + outer_shape: + outer_diameter: + grid: [0.0, 1.0] + values: [9.0, 9.0] + structure: + outfitting_factor: 1.0 + layers: + - name: monopile_wall + material: steel + thickness: + grid: [0.0, 1.0] + values: [0.08, 0.08] + reference_axis: + z: + grid: [0.0, 1.0] + values: [-30.0, 10.0] + tower: + outer_shape: + outer_diameter: + grid: [0.0, 0.5, 1.0] + values: [9.0, 7.5, 6.0] + structure: + outfitting_factor: 1.0 + layers: + - name: tower_wall + material: steel + thickness: + grid: [0.0, 1.0] + values: [0.05, 0.02] + reference_axis: + z: + grid: [0.0, 1.0] + values: [10.0, 110.0] + materials: + - name: steel + E: 2.0e11 + rho: 7850.0 + nu: 0.3 + """) + + +def _tower_yaml(tmp_path: pathlib.Path) -> pathlib.Path: + pytest.importorskip("yaml") + p = tmp_path / "tower.yaml" + p.write_text(_MIN_TOWER, encoding="utf-8") + return p + + +def _f1(model: Tower) -> float: + return float(model.run(n_modes=2, check_model=False).frequencies[0]) + + +# --- issue #133: material / outfitting overrides on from_windio ------------ + +def test_from_windio_none_matches_ontology(tmp_path: pathlib.Path) -> None: + """Passing ``None`` (the default) reproduces the ontology material exactly.""" + p = _tower_yaml(tmp_path) + base = _f1(Tower.from_windio(p)) + same = _f1(Tower.from_windio(p, E=None, rho=None, nu=None, outfitting_factor=None)) + assert same == pytest.approx(base) + + +def test_from_windio_E_override_scales_frequency(tmp_path: pathlib.Path) -> None: + """f scales as sqrt(E) at fixed mass, so overriding E from 2.0e11 to 1.5e11 + lowers the frequency by sqrt(1.5/2.0) (issue #133).""" + p = _tower_yaml(tmp_path) + base = _f1(Tower.from_windio(p)) # ontology E = 2.0e11 + soft = _f1(Tower.from_windio(p, E=1.5e11)) + assert soft / base == pytest.approx(np.sqrt(1.5e11 / 2.0e11), rel=1e-6) + + +def test_from_windio_outfitting_override_scales_frequency( + tmp_path: pathlib.Path, +) -> None: + """Outfitting scales mass, so f ~ 1/sqrt(outfitting). Overriding the + ontology's 1.1 to 1.0 raises the frequency by sqrt(1.1) (issue #133).""" + p = _tower_yaml(tmp_path) + base = _f1(Tower.from_windio(p)) # ontology outfitting = 1.1 + bare = _f1(Tower.from_windio(p, outfitting_factor=1.0)) + assert bare / base == pytest.approx(np.sqrt(1.1), rel=1e-6) + + +def test_from_windio_rho_override_scales_frequency(tmp_path: pathlib.Path) -> None: + """f ~ 1/sqrt(rho) at fixed stiffness (no tip mass).""" + p = _tower_yaml(tmp_path) + base = _f1(Tower.from_windio(p)) # ontology rho = 7800 + heavy = _f1(Tower.from_windio(p, rho=2 * 7800.0)) + assert heavy / base == pytest.approx(1.0 / np.sqrt(2.0), rel=1e-6) + + +def test_from_windio_with_monopile_override(tmp_path: pathlib.Path) -> None: + """The override reaches both segments of the combined monopile+tower path + and softens the whole structure (issue #133).""" + pytest.importorskip("yaml") + p = tmp_path / "mp.yaml" + p.write_text(_MIN_MONOPILE_TOWER, encoding="utf-8") + base = _f1(Tower.from_windio_with_monopile(p, water_depth=30.0)) + soft = _f1(Tower.from_windio_with_monopile(p, water_depth=30.0, E=1.0e11)) + assert soft < base # halving E on both segments lowers the frequency + + +# --- issue #134: per-mode generalised mass / stiffness --------------------- + +def _clamped_tower_with_rna() -> Tower: + z = np.linspace(0.0, 1.0, 11) + od = np.full_like(z, 6.0) + wall = np.full_like(z, 0.03) + rna = TipMassProps( + mass=4.0e5, cm_offset=0.0, cm_axial=3.0, + ixx=2.0e6, iyy=2.0e6, izz=2.0e6, ixy=0.0, izx=0.0, iyz=0.0, + ) + return Tower.from_geometry( + z, od, wall, flexible_length=100.0, E=2.0e11, rho=7850.0, + hub_conn=1, tip_mass=rna, n_nodes=40, + ) + + +def test_generalized_mass_stiffness_populated() -> None: + """The result carries a physical generalised mass (kg) and stiffness + (N/m) per mode (issue #134).""" + res = _clamped_tower_with_rna().run(n_modes=4, check_model=False) + assert res.generalized_mass is not None + assert res.generalized_stiffness is not None + assert res.generalized_mass.shape == res.frequencies.shape + assert res.generalized_stiffness.shape == res.frequencies.shape + # first bending mode: finite, positive, and larger than the tip mass + # (tip translation plus a share of the tower and the CM lever). + assert np.isfinite(res.generalized_mass[0]) + assert res.generalized_mass[0] > 4.0e5 + + +def test_generalized_mass_nan_for_near_zero_tip() -> None: + """A mode whose tip has only numerical lateral motion (torsion / heave / + yaw dominated) yields NaN rather than a huge value from dividing by a tiny + tip displacement (Codex review on #134).""" + from pybmodes.models._pipeline import _generalised_mass_stiffness + + # active DOFs: tip v_disp (global 1), tip w_disp (global 3), one other DOF + active = np.array([1, 3, 99]) + gm = np.eye(3) + # mode 0: genuine lateral tip (w_disp = 1). mode 1: tip is 1e-12 of a large + # non-lateral (twist-like) DOF -> undefined normalisation. + eigvecs = np.array([ + [0.0, 0.0], + [1.0, 1.0e-12], + [0.5, 1.0], + ]) + freqs = np.array([0.5, 2.0]) + m, k = _generalised_mass_stiffness(eigvecs, gm, active, freqs, 10.0) + assert np.isfinite(m[0]) and m[0] > 0.0 + assert np.isnan(m[1]) and np.isnan(k[1]) + + +def test_generalized_mass_stiffness_recovers_frequency() -> None: + """sqrt(K/M)/(2*pi) reproduces each mode's frequency (issue #134).""" + res = _clamped_tower_with_rna().run(n_modes=4, check_model=False) + m = res.generalized_mass + k = res.generalized_stiffness + ok = np.isfinite(m) & np.isfinite(k) & (m > 0) + assert ok.any() + f_from_km = np.sqrt(k[ok] / m[ok]) / (2.0 * np.pi) + assert np.allclose(f_from_km, res.frequencies[ok], rtol=1e-6)