Keep original OWC physical motion failure explicit

tests/test_owc_orifice_corrected_full.py

@@ -9,6 +9,7 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
 
 from wecsim import OrificePTO, PMWave, WEC
 
@@ -21,6 +22,15 @@
 ORIGINAL = FIXTURES / "owc_orifice_corrected_piston/source.npz"
 PUBLISHED = FIXTURES / "owc_orifice_130s/source.npz"
 HYDRO = FIXTURES / "owc_orifice.h5"
+# Full-duration physical limits declared before the corrected-source runs;
+# retain them for both the corrected gate and the original-input diagnostic.
+PHYSICAL_MOTION_LIMITS = np.array([.020, .010, .001, .002])
+PHYSICAL_WORK_LIMIT = .01
+MOTION_NAMES = ("surge", "heave", "pitch", "flexible displacement")
+
+
+class PublishedPhysicalMotionGap(AssertionError):
+    """The physical Python path exceeds a declared original-input limit."""
 
 
 def _peak(actual, expected, limit, name):
@@ -30,7 +40,31 @@ def _peak(actual, expected, limit, name):
     assert error <= limit, f"{name}: {error:.8g} > {limit:.8g}"
 
 
-def test_corrected_owc_full_motion_and_air_system(tmp_path):
+@pytest.fixture(scope="module")
+def physical_owc_130s(tmp_path_factory):
+    """Advance one independent physical 130 s trajectory for both comparisons."""
+    with np.load(PUBLISHED, allow_pickle=False) as saved:
+        phases = saved["phases"].copy()
+        parameters = saved["orifice_parameters"].copy()
+    phase_file = tmp_path_factory.mktemp("owc_phases") / "phases.csv"
+    np.savetxt(phase_file, phases, delimiter=",")
+    wec = WEC("OWC physical piston-path comparison")
+    floating = wec.body("OWC", HYDRO, inertia=(99.28, 11.04, 99.2))
+    wec.floating_gbm(
+        floating, orifice=OrificePTO(*parameters),
+        heave_linear_damping=100, mode_linear_damping=100,
+        heave_drag_cd=1.2, heave_drag_area=8,
+        pitch_drag_cd=1.2, pitch_drag_area=8,
+    )
+    result = wec.run(
+        PMWave(height=1, period=4, phase_file=phase_file),
+        dt=.005, end_time=130, ramp_time=10, radiation_memory=15,
+    )
+    assert result.case["constraint"]["orifice_force_path"] == "coupled"
+    return result
+
+
+def test_corrected_owc_full_motion_and_air_system(physical_owc_130s):
     with np.load(SOURCE, allow_pickle=False) as saved:
         source = {name: saved[name] for name in saved.files}
     with np.load(ORIGINAL, allow_pickle=False) as saved:
@@ -63,20 +97,7 @@ def test_corrected_owc_full_motion_and_air_system(tmp_path):
     _peak(source_air.absorbed_power / 1000, air[:, 5], 1e-7,
           "source power law")
 
-    phase_file = tmp_path / "phases.csv"
-    np.savetxt(phase_file, phases, delimiter=",")
-    wec = WEC("OWC corrected piston-path gate")
-    floating = wec.body("OWC", HYDRO, inertia=(99.28, 11.04, 99.2))
-    wec.floating_gbm(
-        floating, orifice=orifice,
-        heave_linear_damping=100, mode_linear_damping=100,
-        heave_drag_cd=1.2, heave_drag_area=8,
-        pitch_drag_cd=1.2, pitch_drag_area=8,
-    )
-    result = wec.run(
-        PMWave(height=1, period=4, phase_file=phase_file),
-        dt=.005, end_time=130, ramp_time=10, radiation_memory=15,
-    )
+    result = physical_owc_130s
     _peak(result.time, time, 1e-10, "Python time")
     _peak(result.wave_elevation, source["wave"][:, 1], 1e-10,
           "Python incident wave")
@@ -92,17 +113,17 @@ def test_corrected_owc_full_motion_and_air_system(tmp_path):
     ))
     source_reproduction = np.max(np.abs(
         source_motion - existing["fine_motion"]), axis=0)
-    assert np.all(source_reproduction <= [.020, .010, .001, .002]), (
+    assert np.all(source_reproduction <= PHYSICAL_MOTION_LIMITS), (
         source_reproduction
     )
     assert abs(float(source["summary"][0])
                - float(existing["fine_work_internal"])) / float(
                    existing["fine_work_internal"]) <= .005
-    _peak(pose[:, 0], body[:, 1], .020, "surge")
-    _peak(pose[:, 2], body[:, 3], .010, "heave")
+    _peak(pose[:, 0], body[:, 1], PHYSICAL_MOTION_LIMITS[0], "surge")
+    _peak(pose[:, 2], body[:, 3], PHYSICAL_MOTION_LIMITS[1], "heave")
     _peak(result.coordinates["pitch_unwrapped"].position,
-          source_pitch, .001, "pitch")
-    _peak(flexible.position[:, 0], flex[:, 1], .002,
+          source_pitch, PHYSICAL_MOTION_LIMITS[2], "pitch")
+    _peak(flexible.position[:, 0], flex[:, 1], PHYSICAL_MOTION_LIMITS[3],
           "flexible displacement")
     _peak(velocity[:, 0], body[:, 7], .020, "surge speed")
     _peak(velocity[:, 2], body[:, 9], .020, "heave speed")
@@ -125,4 +146,55 @@ def test_corrected_owc_full_motion_and_air_system(tmp_path):
     assert not np.any(python_air.compressibility_flag)
     source_work = float(source["summary"][0])
     python_work = float(np.trapezoid(pto.absorbed_power, time))
-    assert abs(python_work - source_work) / source_work <= .01
+    assert abs(python_work - source_work) / source_work <= PHYSICAL_WORK_LIMIT
+
+
+@pytest.mark.xfail(
+    strict=True, raises=PublishedPhysicalMotionGap,
+    reason="the published flexible state omits the piston reaction",
+)
+def test_original_published_owc_physical_motion(physical_owc_130s):
+    """Keep the original published-input physical failure explicit in CI.
+
+    The passing corrected-source gate above uses this same Python trajectory.
+    A separate source-only same-step control in
+    ``test_owc_orifice_corrected_piston.py`` isolates the changed piston-force
+    connection. This comparison uses the untouched published source output.
+    """
+    with np.load(PUBLISHED, allow_pickle=False) as saved:
+        source = {name: saved[name] for name in saved.files}
+    result = physical_owc_130s
+    np.testing.assert_allclose(result.time, source["time"], rtol=0, atol=1e-10)
+    np.testing.assert_allclose(result.wave_elevation,
+                               source["wave_elevation"], rtol=0, atol=1e-10)
+    assert source["phases"].shape == (500,)
+    source_motion = np.column_stack((
+        source["body_position"], source["flex_position"],
+    ))
+    python_motion = np.column_stack((
+        result.bodies["OWC"].position[:, 0],
+        result.bodies["OWC"].position[:, 2],
+        result.coordinates["pitch_unwrapped"].position,
+        result.flexible_modes["OWC"].position[:, 0],
+    ))
+    gap = np.max(np.abs(python_motion - source_motion), axis=0)
+    source_work = float(np.trapezoid(source["orifice_power"], source["time"]))
+    python_work = float(np.trapezoid(
+        result.ptos["orifice"].absorbed_power, result.time,
+    ))
+    assert source_work > 0 and python_work > 0
+    work_gap = abs(python_work - source_work) / source_work
+    print("published OWC physical motion gap:", gap,
+          "work gap:", work_gap)
+    violations = [
+        f"{name}: {error:.8g} > {limit:.8g}"
+        for name, error, limit in zip(
+            MOTION_NAMES, gap, PHYSICAL_MOTION_LIMITS,
+        ) if error > limit
+    ]
+    if work_gap > PHYSICAL_WORK_LIMIT:
+        violations.append(
+            f"PTO work: {work_gap:.8g} > {PHYSICAL_WORK_LIMIT:.8g}"
+        )
+    if violations:
+        raise PublishedPhysicalMotionGap("; ".join(violations))