diff --git a/.github/workflows/rm3-pto-extension-published-physical.yml b/.github/workflows/rm3-pto-extension-published-physical.yml new file mode 100644 index 00000000..9df98b08 --- /dev/null +++ b/.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 diff --git a/tests/fixtures/matlab/rm3_pto_extension_published_physical/README.md b/tests/fixtures/matlab/rm3_pto_extension_published_physical/README.md new file mode 100644 index 00000000..57218196 --- /dev/null +++ b/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. diff --git a/tests/fixtures/matlab/rm3_pto_extension_published_physical/rm3.h5 b/tests/fixtures/matlab/rm3_pto_extension_published_physical/rm3.h5 new file mode 100644 index 00000000..675ea7a5 Binary files /dev/null and b/tests/fixtures/matlab/rm3_pto_extension_published_physical/rm3.h5 differ diff --git a/tests/fixtures/matlab/rm3_pto_extension_published_physical/source.npz b/tests/fixtures/matlab/rm3_pto_extension_published_physical/source.npz new file mode 100644 index 00000000..3c97083e Binary files /dev/null and b/tests/fixtures/matlab/rm3_pto_extension_published_physical/source.npz differ diff --git a/tests/test_rm3_pto_extension_published_physical.py b/tests/test_rm3_pto_extension_published_physical.py new file mode 100644 index 00000000..c59ca3e9 --- /dev/null +++ b/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