Gate bounded 140 s OSWEC selected-heading quarter-step prefix

.github/workflows/oswec-variable-yaw-250s-step-audit.yml

@@ -16,7 +16,7 @@ permissions:
 jobs:
   source:
     runs-on: ubuntu-latest
-    timeout-minutes: 240
+    timeout-minutes: 330
     steps:
       - uses: actions/checkout@v7
       - uses: actions/checkout@v7

tests/matlab/oswec_variable_yaw_250s_step_audit.m

@@ -1,5 +1,5 @@
 function oswec_variable_yaw_250s_step_audit()
-% Refine the first 120 s of the pinned 250 s two-step artifact without
+% Refine the first 140 s of the pinned 250 s two-step artifact without
 % rerunning its costly simulations. Do not infer full-duration parity.
 repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath'))));
 addpath(genpath(fullfile(repoRoot, 'matlab-ref', 'source')));
@@ -30,7 +30,7 @@ function oswec_variable_yaw_250s_step_audit()
 sourceDir = fullfile(repoRoot, 'applications', ...
     'Variable_Hydro', 'Passive_Yaw');
 caseDir = fullfile(repoRoot, 'applications', 'Variable_Hydro', ...
-    'paired_Passive_Yaw_irregular_120s');
+    'paired_Passive_Yaw_irregular_140s');
 [copied, copyMessage] = copyfile(sourceDir, caseDir);
 assert(copied, copyMessage);
 bankScript = fullfile(caseDir, 'hydroData', 'bemio.m');
@@ -49,22 +49,22 @@ function oswec_variable_yaw_250s_step_audit()
     contains(original, 'waves.phaseSeed = 1;'), ...
     'The pinned irregular variable-yaw input changed');
 quarterInput = strrep(original, 'simu.endTime = 600;', ...
-    'simu.endTime = 120;');
+    'simu.endTime = 140;');
 quarterInput = strrep(quarterInput, 'simu.dt = 0.01;', ...
     'simu.dt = 0.0025;');
 write_input(inputFile, quarterInput);
 cd(caseDir);
 waveFlag = 'irregular';
 wecSim;
-check_run(simu, waves, bemDirections, output, 0.0025, 48001);
+check_run(simu, waves, bemDirections, output, 0.0025, 56001);
 quarter = trajectory(output, bemDirections);
 quarterWave = waves.waveAmpTime;
 quarterComponents = [waves.omega(:), waves.amplitude(:), ...
     waves.dOmega(:), waves.phase(:)];
-fine120 = fine(1:24001, :);
-fine120Wave = fineWave(1:24001, :);
-assert(max(abs(quarter(1:2:end, 1) - fine120(:, 1))) < 1e-10 && ...
-    max(abs(quarterWave(1:2:end, :) - fine120Wave), [], 'all') < 1e-12 && ...
+fine140 = fine(1:28001, :);
+fine140Wave = fineWave(1:28001, :);
+assert(max(abs(quarter(1:2:end, 1) - fine140(:, 1))) < 1e-10 && ...
+    max(abs(quarterWave(1:2:end, :) - fine140Wave), [], 'all') < 1e-12 && ...
     max(abs(quarterComponents - coarseComponents), [], 'all') < 1e-12, ...
     'Quarter-step refinement changed the time grid or incident sea');
 writematrix(quarter, fullfile(outDir, 'quarter_trajectory.csv'));
@@ -73,10 +73,10 @@ function oswec_variable_yaw_250s_step_audit()
 coarseWork = -trapz(coarse(:, 1), coarse(:, 7));
 fineWork = -trapz(fine(:, 1), fine(:, 7));
 quarterWork = -trapz(quarter(:, 1), quarter(:, 7));
-fine120Work = -trapz(fine120(:, 1), fine120(:, 7));
+fine140Work = -trapz(fine140(:, 1), fine140(:, 7));
 mismatch = coarse(:, 4) ~= fine(1:2:end, 4);
 firstMismatch = find(mismatch, 1);
-quarterMismatch = fine120(:, 4) ~= quarter(1:2:end, 4);
+quarterMismatch = fine140(:, 4) ~= quarter(1:2:end, 4);
 firstQuarterMismatch = find(quarterMismatch, 1);
 metrics = struct( ...
     'maxWaveErrorM', max(abs(fineWave(1:2:end, 2) - coarseWave(:, 2))), ...
@@ -89,17 +89,17 @@ function oswec_variable_yaw_250s_step_audit()
     'fineHeadingSwitches', nnz(diff(fine(:, 4))), ...
     'headingMismatchSamples', nnz(mismatch));
 metrics.fineToQuarterMaxYawErrorRad = ...
-    max(abs(fine120(:, 2) - quarter(1:2:end, 2)));
+    max(abs(fine140(:, 2) - quarter(1:2:end, 2)));
 metrics.fineToQuarterMaxYawRateErrorRadPerS = ...
-    max(abs(fine120(:, 3) - quarter(1:2:end, 3)));
+    max(abs(fine140(:, 3) - quarter(1:2:end, 3)));
 metrics.fineToQuarterMaxPtoTorqueErrorNm = ...
-    max(abs(fine120(:, 6) - quarter(1:2:end, 6)));
-metrics.quarterDurationS = 120;
-metrics.fine120PtoWorkJ = fine120Work;
+    max(abs(fine140(:, 6) - quarter(1:2:end, 6)));
+metrics.quarterDurationS = 140;
+metrics.fine140PtoWorkJ = fine140Work;
 metrics.quarterPtoWorkJ = quarterWork;
 metrics.fineToQuarterRelativePtoWorkChange = ...
-    (quarterWork - fine120Work) / fine120Work;
-metrics.fine120HeadingSwitches = nnz(diff(fine120(:, 4)));
+    (quarterWork - fine140Work) / fine140Work;
+metrics.fine140HeadingSwitches = nnz(diff(fine140(:, 4)));
 metrics.quarterHeadingSwitches = nnz(diff(quarter(:, 4)));
 metrics.fineToQuarterHeadingMismatchSamples = nnz(quarterMismatch);
 if isempty(firstMismatch)
@@ -111,13 +111,13 @@ function oswec_variable_yaw_250s_step_audit()
     metrics.firstFineToQuarterHeadingMismatchS = -1;
 else
     metrics.firstFineToQuarterHeadingMismatchS = ...
-        fine120(firstQuarterMismatch, 1);
+        fine140(firstQuarterMismatch, 1);
 end
 fid = fopen(fullfile(outDir, 'metrics.json'), 'w');
 assert(fid ~= -1, 'Could not write source step metrics');
 fprintf(fid, '%s\n', jsonencode(metrics));
 fclose(fid);
-fprintf('Variable-yaw 250 s source with 120 s quarter-step audit: %s\n', ...
+fprintf('Variable-yaw 250 s source with 140 s quarter-step audit: %s\n', ...
     jsonencode(metrics));
 end
 
@@ -130,7 +130,7 @@ function write_input(inputFile, contents)
 
 function check_run(simu, waves, bemDirections, output, expectedDt, expectedSamples)
 assert(strcmp(simu.solver, 'ode4') && simu.dt == expectedDt && ...
-    simu.endTime == 120 && simu.rampTime == 100 && ...
+    simu.endTime == 140 && simu.rampTime == 100 && ...
     strcmp(waves.type, 'irregular') && waves.phaseSeed == 1 && ...
     waves.height == 2.5 && waves.period == 8 && ...
     isequal(bemDirections, -30:0.25:30), ...

tests/test_variable_yaw_quarter_step_audit.py

@@ -1,4 +1,4 @@
-"""Gate the early selected-heading motion against a refined source path."""
+"""Gate the stable 120 s prefix and audit a longer selected-heading path."""
 
 import os
 from pathlib import Path
@@ -17,7 +17,43 @@
 )
 
 
-def test_irregular_120s_quarter_step_pair():
+def _errors(fine, quarter, result):
+    """Compare a quarter-step prefix with the same half-step source window."""
+    fine = fine[:(len(quarter) + 1) // 2]
+    motion = result.bodies["flap"]
+    pto = result.ptos["hinge"]
+    source_work = -np.trapezoid(quarter[:, 6], quarter[:, 0])
+    fine_work = -np.trapezoid(fine[:, 6], fine[:, 0])
+    python_work = np.trapezoid(pto.absorbed_power[:len(quarter)],
+                               result.time[:len(quarter)])
+    assert source_work > 0 and fine_work > 0 and python_work > 0
+    source_step = (
+        np.max(np.abs(fine[:, 1] - quarter[::2, 1])),
+        np.max(np.abs(fine[:, 2] - quarter[::2, 2])),
+        np.max(np.abs(fine[:, 5] - quarter[::2, 5])),
+        abs(source_work - fine_work) / source_work,
+    )
+    python_source = (
+        np.max(np.abs(motion.position[:len(quarter), 5] - quarter[:, 1])),
+        np.max(np.abs(motion.velocity[:len(quarter), 5] - quarter[:, 2])),
+        np.max(np.abs(pto.force[:len(quarter)] - quarter[:, 5])),
+        abs(python_work - source_work) / source_work,
+    )
+    return source_step, python_source
+
+
+def _report(duration, source_step, python_source):
+    print(
+        f"{duration} s quarter-step source self-change: "
+        f"yaw={source_step[0]:.6g} rad, speed={source_step[1]:.6g} rad/s, "
+        f"torque={source_step[2]:.6g} Nm, work={source_step[3]:.3%}; "
+        f"Python/source yaw={python_source[0]:.6g} rad, "
+        f"speed={python_source[1]:.6g} rad/s, "
+        f"torque={python_source[2]:.6g} Nm, work={python_source[3]:.3%}"
+    )
+
+
+def test_irregular_140s_quarter_step_audit():
     source = Path(SOURCE)
     quarter_file = source / "quarter_trajectory.csv"
     assert quarter_file.is_file(), "quarter-step source run is absent"
@@ -26,15 +62,16 @@ def test_irregular_120s_quarter_step_pair():
     fine_wave = np.loadtxt(source / "fine_wave.csv", delimiter=",")
     quarter_wave = np.loadtxt(source / "quarter_wave.csv", delimiter=",")
     assert fine.shape == (50001, 7)
-    assert quarter.shape == (48001, 7)
+    assert quarter.shape == (56001, 7)
     assert fine_wave.shape == (50001, 2)
-    assert quarter_wave.shape == (48001, 2)
+    assert quarter_wave.shape == (56001, 2)
     assert all(np.isfinite(array).all() for array in
                (fine, quarter, fine_wave, quarter_wave))
-    fine120 = fine[:24001]
-    fine120_wave = fine_wave[:24001]
-    np.testing.assert_allclose(quarter[::2, 0], fine120[:, 0], rtol=0, atol=1e-10)
-    np.testing.assert_allclose(quarter_wave[::2], fine120_wave, rtol=0,
+    fine140 = fine[:28001]
+    fine140_wave = fine_wave[:28001]
+    np.testing.assert_allclose(quarter[::2, 0], fine140[:, 0],
+                               rtol=0, atol=1e-10)
+    np.testing.assert_allclose(quarter_wave[::2], fine140_wave, rtol=0,
                                atol=1e-12)
 
     headings = np.arange(-30, 30.25, .25)
@@ -48,51 +85,34 @@ def test_irregular_120s_quarter_step_pair():
 
     hydro = (Path(APPLICATIONS) /
              "_Common_Input_Files/OSWEC/hydroData/oswec.h5")
-    wec = WEC("OSWEC quarter-step selected-heading diagnostic")
+    wec = WEC("OSWEC 140 s selected-heading diagnostic")
     flap = wec.body("flap", hydro, mass=12700, inertia=(1.85e6,) * 3,
                     passive_yaw=True, yaw_heading_bank=headings)
     wec.body("base", hydro, mass=999, inertia=(999,) * 3)
     yaw = wec.coordinate("yaw", flap.move("yaw", pivot=WorldPoint(0, 0, -8.9)))
     wec.rotational_pto("hinge", yaw, damping=120000)
     result = wec.run(
         PMWave(2.5, 8, direction=10, seed=1, phase_generator="matlab"),
-        dt=.0025, end_time=120, ramp_time=100, radiation_memory=40,
+        dt=.0025, end_time=140, ramp_time=100, radiation_memory=40,
     )
     np.testing.assert_allclose(result.time, quarter[:, 0], rtol=0, atol=1e-10)
     np.testing.assert_allclose(result.wave_elevation, quarter_wave[:, 1],
                                rtol=0, atol=1e-12)
 
-    source_work = -np.trapezoid(quarter[:, 6], quarter[:, 0])
-    python_work = np.trapezoid(result.ptos["hinge"].absorbed_power, result.time)
-    assert source_work > 0 and python_work > 0
-    source_step_yaw = np.max(np.abs(fine120[:, 1] - quarter[::2, 1]))
-    source_step_speed = np.max(np.abs(fine120[:, 2] - quarter[::2, 2]))
-    source_step_torque = np.max(np.abs(fine120[:, 5] - quarter[::2, 5]))
-    source_step_work = -np.trapezoid(fine120[:, 6], fine120[:, 0])
-    motion = result.bodies["flap"]
-    python_source_yaw = np.max(np.abs(motion.position[:, 5] - quarter[:, 1]))
-    python_source_speed = np.max(np.abs(motion.velocity[:, 5] - quarter[:, 2]))
-    python_source_torque = np.max(np.abs(
-        result.ptos["hinge"].force - quarter[:, 5]
-    ))
-    source_step_work_error = abs(source_work - source_step_work) / source_work
-    python_work_error = abs(python_work - source_work) / source_work
-
     # These limits exceed the measured 120 s 0.01-to-0.005 s source
     # self-change but still reject a new event-driven trajectory split.
     limits = (.0006, .00015, 20, .001)
-    print(
-        f"120 s quarter-step source self-change: yaw={source_step_yaw:.6g} rad, "
-        f"speed={source_step_speed:.6g} rad/s, "
-        f"torque={source_step_torque:.6g} Nm, work={source_step_work_error:.3%}; "
-        f"Python/source yaw={python_source_yaw:.6g} rad, "
-        f"speed={python_source_speed:.6g} rad/s, "
-        f"torque={python_source_torque:.6g} Nm, work={python_work_error:.3%}"
-    )
-    for measured in (
-        (source_step_yaw, source_step_speed, source_step_torque,
-         source_step_work_error),
-        (python_source_yaw, python_source_speed, python_source_torque,
-         python_work_error),
-    ):
+    source120, python120 = _errors(fine, quarter[:48001], result)
+    _report(120, source120, python120)
+    for measured in (source120, python120):
         assert all(value < limit for value, limit in zip(measured, limits))
+
+    # The pinned 0.01-to-0.005 s source pair stays within these same limits
+    # through 140 s, then exceeds the yaw limit by 145 s. The quarter-step
+    # pair must establish its own stability before extending the gate.
+    source140, python140 = _errors(fine, quarter, result)
+    _report(140, source140, python140)
+    source_stable = all(value < limit for value, limit in zip(source140, limits))
+    print(f"140 s MATLAB step stability under 120 s limits: {source_stable}")
+    assert source_stable, "140 s MATLAB source is step-sensitive"
+    assert all(value < limit for value, limit in zip(python140, limits))