Gate unchanged RM3 PTO Extension inputs with physical Python runs

.github/workflows/rm3-pto-extension-published-physical.yml

@@ -0,0 +1,69 @@
+name: RM3 PTO Extension unchanged-input physical comparison
+
+on:
+  pull_request:
+    paths:
+      - '.github/workflows/rm3-pto-extension-published-physical.yml'
+      - 'tests/test_rm3_pto_extension_published_physical.py'
+      - 'tests/fixtures/matlab/rm3_pto_extension_published_physical/**'
+      - 'tests/matlab/reference_model_baseline.m'
+      - 'wecsim/dynamics/caseDynamics.py'
+      - 'wecsim/dynamics/rm3Regular.py'
+      - 'wecsim/hydrodynamics/radiation.py'
+  workflow_dispatch:
+
+jobs:
+  published-source-and-physical:
+    runs-on: ubuntu-latest
+    timeout-minutes: 60
+    steps:
+      - uses: actions/checkout@v7
+      - uses: actions/checkout@v7
+        with:
+          repository: WEC-Sim/WEC-Sim
+          ref: 0753b2e47f2457c078751dcfe5d251d1767b80ab
+          path: matlab-ref
+          filter: blob:none
+          sparse-checkout: |
+            source
+            examples/RM3
+      - uses: actions/checkout@v7
+        with:
+          repository: WEC-Sim/WEC-Sim_Applications
+          ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e
+          path: applications
+          filter: blob:none
+          lfs: true
+          sparse-checkout: |
+            _Common_Input_Files/RM3
+            RM3_PTO_Extension
+      - uses: matlab-actions/setup-matlab@v3
+        with:
+          release: R2025b
+          products: Simulink Simscape Simscape_Multibody
+      - uses: matlab-actions/run-command@v3
+        with:
+          command: addpath('tests/matlab'); reference_model_baseline('RM3_PTO_Extension')
+      - uses: actions/setup-python@v7
+        with:
+          python-version: '3.12'
+          cache: pip
+          cache-dependency-path: pyproject.toml
+      - run: python -m pip install -e ".[test]"
+      - name: Verify fresh unchanged source and independent physical gate
+        run: >-
+          python -m pytest -q -s
+          tests/test_rm3_pto_extension_published_physical.py
+          tests/test_rm3_pto_extension_parity.py
+        env:
+          WEC_SIM_APPLICATIONS_DIR: applications
+          WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output
+      - uses: actions/upload-artifact@v4
+        if: always()
+        with:
+          name: rm3-pto-extension-published-physical
+          path: |
+            matlab-reference-model-output/RM3_PTO_Extension_*_body*.csv
+            matlab-reference-model-output/RM3_PTO_Extension_*_pto1.csv
+            applications/_Common_Input_Files/RM3/hydroData/rm3.h5
+          if-no-files-found: warn

tests/fixtures/matlab/rm3_pto_extension_published_physical/README.md

@@ -0,0 +1,53 @@
+# Unchanged RM3 PTO Extension physical comparison
+
+This fixture archives the **unchanged published** 30 s, 0.1 s `ode4`
+`RM3_PTO_Extension/float` and `/spar` source runs. The six arrays in
+`source.npz` are the complete 301-row, 25-column MATLAB body and PTO CSVs,
+loaded as `float64` without interpolation or rounding. The exact generated
+`rm3.h5` used by those runs is also archived. Both inputs have no wave, zero
+PTO spring and damping, and a 5 m initial PTO extension.
+
+- Source: [successful R2025b MATLAB reference run 38282150325](https://github.com/cmudrc/wec-sim-python/actions/runs/38282150325), repository commit `7c403d3e2e07b9940e4a2a280d375000b8985fea`.
+- Pinned WEC-Sim core: `0753b2e47f2457c078751dcfe5d251d1767b80ab`.
+- Pinned Applications: `d53d4d4c9eda2581f04204f5d394a6ef84bb099e`.
+- `rm3.h5` SHA-256: `dac35aa42321e1ade04e8935b38680d7b70e7b247edac9b5b800d01142fcc259`.
+- `source.npz` SHA-256: `535eaabbf0ec68949ffffff215548003696a0ad99029214a396b37d6fa67e00c`.
+
+The nine time, active heave/speed, and PTO columns match the repository's
+earlier `rm3_pto_extension_step.npz` published 0.1 s arms exactly. The source
+HDF5's active center-of-gravity, volume, restoring, added-mass, and radiation
+IRF datasets equal those in `examples/data/rm3_end_stops.h5` by value. The
+portable test checks both identities and can compare all six arrays to a
+fresh MATLAB reference output when `WEC_SIM_MATLAB_MODEL_OUTPUT_DIR` is set.
+The [dedicated source workflow](../../../../.github/workflows/rm3-pto-extension-published-physical.yml)
+regenerates both unchanged inputs from the pinned checkouts and runs that
+fresh comparison together with the independent Python gate.
+
+`tests/test_rm3_pto_extension_published_physical.py` advances a public
+`run_case` using physical implicit added mass, the exact source HDF5, no
+wave, the original 0.1 s grid and 30 s duration, and the published initial
+coordinates. It then checks all six positions and velocities of both bodies,
+PTO stroke and speed, and zero force and power against the complete MATLAB
+records. The limits are **inherited unchanged** from
+`tests/test_rm3_pto_extension_parity.py`: active float heave/stroke 90 mm
+and speed 210 mm/s; active spar 18 mm and 8 mm/s; the original limits on
+all other channels also apply. The physical run receives no MATLAB states,
+accelerations, forces, or solver steps. An opt-in source-`ode4` run is used
+only as a negative control: it nearly reproduces the published numerical
+history, whereas the physical trajectory is measurably distinct.
+
+On local arm64, the physical float's largest heave/stroke gap is 85.849 mm
+and speed gap 202.600 mm/s, leaving margins of 4.151 mm and 7.400 mm/s.
+The spar gaps are 14.877 mm and 6.522 mm/s, leaving margins of 3.123 mm
+and 1.478 mm/s. These are formal passes under the inherited gates with
+limited margin. The original limits preceded this source-step investigation
+but were not prospective thresholds chosen blind to all earlier errors.
+
+The published source remains numerically sensitive. Changing only its
+Simulink `FixedStep` from 0.1 to 0.0125 s removes over 99% of the float's
+and over 99.99% of the spar's squared heave/speed gap to a refined source
+run; see `../rm3_pto_extension_fixedstep_only/README.md`. This test therefore
+does not claim that the coarse published trajectory is source-converged or
+that Python reproduces it pointwise. The separate corrected fine-step
+baseline and the original published comparison retain distinct fixtures.
+Case status and badge interpretation require separate review.

tests/test_rm3_pto_extension_published_physical.py

@@ -0,0 +1,220 @@
+"""Portable physical gates for the unchanged RM3 PTO Extension inputs."""
+
+import hashlib
+import os
+from pathlib import Path
+
+import h5py
+import numpy as np
+import pytest
+
+from wecsim.dynamics.caseDynamics import run_case
+
+
+ROOT = Path(__file__).resolve().parents[1]
+FIXTURE = ROOT / "tests/fixtures/matlab/rm3_pto_extension_published_physical"
+SOURCE = FIXTURE / "source.npz"
+HYDRO = FIXTURE / "rm3.h5"
+STEP = ROOT / "tests/fixtures/matlab/rm3_pto_extension_step.npz"
+FRESH_SOURCE = os.environ.get("WEC_SIM_MATLAB_MODEL_OUTPUT_DIR")
+FRESH_APPLICATIONS = os.environ.get("WEC_SIM_APPLICATIONS_DIR")
+
+# Retain the limits declared in test_rm3_pto_extension_parity.py before the
+# source-step investigation. Only the qualifying solver changes here: this
+# test uses the public physical implicit-mass path, not source-style ode4.
+CASES = (
+    ("float", {"float_heave": 5}, 0, 0.090, 0.210, 0.050),
+    ("spar", {"spar_heave": -5}, 1, 0.018, 0.008, 0.008),
+)
+
+
+def _bound(actual, expected, limit, label):
+    assert actual.shape == expected.shape, label
+    assert np.isfinite(actual).all() and np.isfinite(expected).all(), label
+    error = float(np.max(np.abs(actual - expected)))
+    print(f"{label}: {error:.9g} (limit {limit:.9g})")
+    assert error < limit, label
+    return error
+
+
+def _case(coordinate, *, numerical=False):
+    simulation = {
+        "dt": 0.1, "end_time": 30, "radiation_memory": 60,
+        "added_mass_scheme": "simulink_delay" if numerical else "implicit",
+    }
+    if numerical:
+        simulation["solver"] = "wecsim_ode4"
+    return {
+        "simulation": simulation,
+        "wave": {"type": "none"},
+        "bodies": [
+            {"hydro_file": str(HYDRO), "hydro_body": 1,
+             "mass": "equilibrium", "pitch_inertia": 21_306_090.66},
+            {"hydro_file": str(HYDRO), "hydro_body": 2,
+             "mass": "equilibrium", "pitch_inertia": 94_407_091.24},
+        ],
+        "constraint": {"kind": "floating_joint", "location": [0, 0, 0],
+                       "initial_coordinate": coordinate},
+        "pto": {"kind": "relative_heave", "damping": 0},
+    }
+
+
+def test_published_physical_fixture_matches_pinned_source():
+    assert hashlib.sha256(HYDRO.read_bytes()).hexdigest() == (
+        "dac35aa42321e1ade04e8935b38680d7b70e7b247edac9b5b800d01142fcc259"
+    )
+    assert hashlib.sha256(SOURCE.read_bytes()).hexdigest() == (
+        "535eaabbf0ec68949ffffff215548003696a0ad99029214a396b37d6fa67e00c"
+    )
+    active = (
+        "properties/cg", "properties/disp_vol",
+        "hydro_coeffs/linear_restoring_stiffness",
+        "hydro_coeffs/added_mass/inf_freq",
+        "hydro_coeffs/radiation_damping/impulse_response_fun/K",
+        "hydro_coeffs/radiation_damping/impulse_response_fun/t",
+    )
+    if FRESH_APPLICATIONS:
+        apps = Path(FRESH_APPLICATIONS)
+        for name in ("float", "spar"):
+            input_file = apps / "RM3_PTO_Extension" / name / "wecSimInputFile.m"
+            assert input_file.is_file(), input_file
+            contents = input_file.read_text()
+            for setting in (
+                "simu.endTime = 30;", "simu.solver = 'ode4';",
+                "simu.dt = 0.1;", "waves = waveClass('noWaveCIC');",
+                "pto(1).stiffness = 0;", "pto(1).damping = 0;",
+                "pto(1).extension.PositionTargetValue = 5;",
+            ):
+                assert setting in contents, (name, setting)
+    with h5py.File(HYDRO) as source, h5py.File(
+        ROOT / "examples/data/rm3_end_stops.h5"
+    ) as bundled:
+        for body in ("body1", "body2"):
+            for suffix in active:
+                key = f"{body}/{suffix}"
+                np.testing.assert_array_equal(source[key], bundled[key])
+        if FRESH_APPLICATIONS:
+            fresh_hydro = (Path(FRESH_APPLICATIONS) /
+                           "_Common_Input_Files/RM3/hydroData/rm3.h5")
+            with h5py.File(fresh_hydro) as fresh:
+                for body in ("body1", "body2"):
+                    for suffix in active:
+                        key = f"{body}/{suffix}"
+                        if suffix.endswith("impulse_response_fun/K"):
+                            np.testing.assert_allclose(
+                                source[key], fresh[key], rtol=0, atol=1e-10,
+                            )
+                        else:
+                            np.testing.assert_array_equal(source[key], fresh[key])
+
+    with np.load(SOURCE, allow_pickle=False) as saved, \
+            np.load(STEP, allow_pickle=False) as prior:
+        assert set(saved.files) == {
+            f"{name}_{part}" for name in ("float", "spar")
+            for part in ("body1", "body2", "pto1")
+        }
+        for name in ("float", "spar"):
+            body1, body2, pto = (
+                saved[f"{name}_{part}"] for part in ("body1", "body2", "pto1")
+            )
+            for array in (body1, body2, pto):
+                assert array.shape == (301, 25) and np.isfinite(array).all()
+                np.testing.assert_allclose(
+                    array[:, 0], np.arange(301) * .1, rtol=0, atol=1e-10,
+                )
+            assert np.max(np.abs(pto[:, 15])) < 1e-7
+            assert np.max(np.abs(pto[:, 21])) < 1e-7
+            compact = np.column_stack((
+                body1[:, 0], body1[:, 3], body1[:, 9],
+                body2[:, 3], body2[:, 9],
+                pto[:, 3], pto[:, 9], pto[:, 15], pto[:, 21],
+            ))
+            np.testing.assert_array_equal(compact, prior[f"{name}_0p1"])
+            if FRESH_SOURCE:
+                for part, archived in (("body1", body1),
+                                       ("body2", body2), ("pto1", pto)):
+                    fresh = np.loadtxt(
+                        Path(FRESH_SOURCE) /
+                        f"RM3_PTO_Extension_{name}_{part}.csv",
+                        delimiter=",", ndmin=2,
+                    )
+                    np.testing.assert_allclose(fresh, archived, rtol=1e-10,
+                                               atol=1e-6)
+
+
+@pytest.mark.parametrize(
+    "name,coordinate,active_body,heave_limit,speed_limit,control_gap",
+    CASES,
+)
+def test_unchanged_published_input_with_independent_physical_motion(
+    name, coordinate, active_body, heave_limit, speed_limit, control_gap,
+):
+    # MATLAB states, accelerations, forces, and solver steps are read only
+    # after the independent physical run has completed.
+    response = run_case(_case(coordinate))
+    assert response.wave_elevation is None
+    assert response.body_position.shape == response.body_velocity.shape == (
+        301, 2, 6,
+    )
+    with np.load(SOURCE, allow_pickle=False) as saved:
+        source = {part: saved[f"{name}_{part}"].copy()
+                  for part in ("body1", "body2", "pto1")}
+
+    for index in range(2):
+        expected = source[f"body{index + 1}"]
+        _bound(response.time, expected[:, 0], 1e-10,
+               f"{name} body{index + 1} time")
+        for axis in range(6):
+            if axis == 2:
+                limits = ((heave_limit, speed_limit) if index == active_body
+                          else (1e-8, 1e-8))
+            elif axis == 0:
+                limits = (3e-5, 1e-5)
+            elif axis == 4:
+                limits = (1e-6, 1e-6)
+            else:
+                limits = (1e-9, 1e-9)
+            _bound(response.body_position[:, index, axis],
+                   expected[:, 1 + axis], limits[0],
+                   f"{name} body{index + 1} axis {axis} position")
+            _bound(response.body_velocity[:, index, axis],
+                   expected[:, 7 + axis], limits[1],
+                   f"{name} body{index + 1} axis {axis} speed")
+
+    with h5py.File(HYDRO) as hydro:
+        center_gap = (
+            np.asarray(hydro["body1/properties/cg"]).ravel()[2]
+            - np.asarray(hydro["body2/properties/cg"]).ravel()[2]
+        )
+    pitch = response.body_position[:, 0, 4]
+    pitch_speed = response.body_velocity[:, 0, 4]
+    cosine = np.cos(pitch)
+    stroke = ((response.body_position[:, 0, 2]
+               - response.body_position[:, 1, 2]) / cosine - center_gap)
+    speed = ((response.body_velocity[:, 0, 2]
+              - response.body_velocity[:, 1, 2]
+              + (center_gap + stroke) * np.sin(pitch) * pitch_speed) / cosine)
+    pto = source["pto1"]
+    _bound(response.time, pto[:, 0], 1e-10, f"{name} PTO time")
+    assert abs(pto[0, 3] - 5) < 1e-10
+    assert abs(stroke[0] - 5) < 1e-8
+    assert abs(speed[0]) < 1e-8
+    _bound(stroke, pto[:, 3], heave_limit, f"{name} PTO stroke")
+    _bound(speed, pto[:, 9], speed_limit, f"{name} PTO speed")
+    _bound(response.pto_force, pto[:, 15], 1e-7, f"{name} PTO zero force")
+    _bound(response.pto_force * speed, pto[:, 21], 1e-7,
+           f"{name} PTO zero power")
+
+    # Strict negative control: the opt-in source-numerical ode4 path nearly
+    # reproduces the archived trace, while the qualifying implicit trajectory
+    # is distinct. A switch of the qualifying run to that route fails here.
+    numerical = run_case(_case(coordinate, numerical=True))
+    reference_heave = source[f"body{active_body + 1}"][:, 3]
+    _bound(numerical.body_position[:, active_body, 2], reference_heave,
+           1e-8, f"{name} source-numerical control")
+    separation = float(np.max(np.abs(
+        response.body_position[:, active_body, 2]
+        - numerical.body_position[:, active_body, 2]
+    )))
+    print(f"{name} physical/source-numerical separation: {separation:.9g}")
+    assert separation > control_gap