diff --git a/.github/workflows/oswec-variable-yaw-250s-step-audit.yml b/.github/workflows/oswec-variable-yaw-250s-step-audit.yml index 49af329f..1c405c48 100644 --- a/.github/workflows/oswec-variable-yaw-250s-step-audit.yml +++ b/.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 diff --git a/tests/matlab/oswec_variable_yaw_250s_step_audit.m b/tests/matlab/oswec_variable_yaw_250s_step_audit.m index d26e93f4..c1984dca 100644 --- a/tests/matlab/oswec_variable_yaw_250s_step_audit.m +++ b/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), ... diff --git a/tests/test_variable_yaw_quarter_step_audit.py b/tests/test_variable_yaw_quarter_step_audit.py index 9ee736b7..257d41be 100644 --- a/tests/test_variable_yaw_quarter_step_audit.py +++ b/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,7 +85,7 @@ 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) @@ -56,43 +93,26 @@ def test_irregular_120s_quarter_step_pair(): 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))