From e0d09d1ad94255b0f8945e434bc016d374049848 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 10:12:16 -0400 Subject: [PATCH 1/4] Add nearest-heading yaw bank for variable hydro application --- .../workflows/reference-model-baselines.yml | 12 +- tests/matlab/reference_model_baseline.m | 55 ++++++++ tests/test_variable_yaw_heading_bank.py | 117 ++++++++++++++++++ wecsim/api.py | 19 ++- wecsim/caseDynamics.py | 16 ++- wecsim/passiveYaw.py | 39 +++++- 6 files changed, 249 insertions(+), 9 deletions(-) create mode 100644 tests/test_variable_yaw_heading_bank.py diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index c50b067..8b30dc4 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -84,7 +84,7 @@ on: model: description: Run all models or a selected reference application type: choice - options: [all, RM3_MCR_ARRAY, SPHERE_ELEVATION_IMPORT, MORISON_FIXED, SPHERE_MOVING_MORISON, SPHERE_MOVING_MORISON_WAVE, RM3_DD_PTO, SPHERE_MPC, GBM_BARGE, SPHERE_REACTIVE_DDPTO, SPHERE_VARIABLE_MASS, ELLIPSOID_NLH_REG, ELLIPSOID_NLH_CIC, ELLIPSOID_NLH_ODE45, OSWEC_PASSIVE_YAW_IRR_CONT, OSWEC_PASSIVE_YAW_IRR, OSWEC_PASSIVE_YAW, OSWEC_FULL_DIR, OSWEC_MULTI_WAVE, RM3_B2B, RM3_MOORING_MATRIX, SPHERE_MEAN_DRIFT, RM3_MCR_SEASTATE, RM3_END_STOPS, RM3_END_STOPS_STEP, RM3_END_STOPS_STEP_FINE, RM3_END_STOPS_STEP_FINER, RM3_END_STOPS_FULL_FINE, RM3_END_STOPS_FULL_FINER, RM3_END_STOPS_FULL_PAIR] + options: [all, RM3_MCR_ARRAY, SPHERE_ELEVATION_IMPORT, MORISON_FIXED, SPHERE_MOVING_MORISON, SPHERE_MOVING_MORISON_WAVE, RM3_DD_PTO, SPHERE_MPC, GBM_BARGE, SPHERE_REACTIVE_DDPTO, SPHERE_VARIABLE_MASS, ELLIPSOID_NLH_REG, ELLIPSOID_NLH_CIC, ELLIPSOID_NLH_ODE45, OSWEC_PASSIVE_YAW_IRR_CONT, OSWEC_PASSIVE_YAW_IRR, OSWEC_PASSIVE_YAW, OSWEC_VARIABLE_YAW_2DEG, OSWEC_FULL_DIR, OSWEC_MULTI_WAVE, RM3_B2B, RM3_MOORING_MATRIX, SPHERE_MEAN_DRIFT, RM3_MCR_SEASTATE, RM3_END_STOPS, RM3_END_STOPS_STEP, RM3_END_STOPS_STEP_FINE, RM3_END_STOPS_STEP_FINER, RM3_END_STOPS_FULL_FINE, RM3_END_STOPS_FULL_FINER, RM3_END_STOPS_FULL_PAIR] default: all jobs: @@ -94,7 +94,7 @@ jobs: strategy: fail-fast: false matrix: - model: ${{ fromJSON(inputs.model == 'RM3_MCR_ARRAY' && '["RM3_MCR_ARRAY"]' || inputs.model == 'SPHERE_ELEVATION_IMPORT' && '["SPHERE_ELEVATION_IMPORT"]' || inputs.model == 'MORISON_FIXED' && '["MORISON_FIXED"]' || inputs.model == 'SPHERE_MOVING_MORISON' && '["SPHERE_MOVING_MORISON"]' || inputs.model == 'SPHERE_MOVING_MORISON_WAVE' && '["SPHERE_MOVING_MORISON_WAVE"]' || inputs.model == 'RM3_DD_PTO' && '["RM3_DD_PTO"]' || inputs.model == 'SPHERE_MPC' && '["SPHERE_MPC"]' || inputs.model == 'GBM_BARGE' && '["GBM_BARGE"]' || inputs.model == 'SPHERE_REACTIVE_DDPTO' && '["SPHERE_REACTIVE_DDPTO"]' || inputs.model == 'SPHERE_VARIABLE_MASS' && '["SPHERE_VARIABLE_MASS"]' || inputs.model == 'ELLIPSOID_NLH_REG' && '["ELLIPSOID_NLH_REG"]' || inputs.model == 'ELLIPSOID_NLH_CIC' && '["ELLIPSOID_NLH_CIC"]' || inputs.model == 'ELLIPSOID_NLH_ODE45' && '["ELLIPSOID_NLH_ODE45"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' && '["OSWEC_PASSIVE_YAW_IRR_CONT"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR' && '["OSWEC_PASSIVE_YAW_IRR"]' || inputs.model == 'OSWEC_PASSIVE_YAW' && '["OSWEC_PASSIVE_YAW"]' || inputs.model == 'OSWEC_FULL_DIR' && '["OSWEC_FULL_DIR"]' || inputs.model == 'OSWEC_MULTI_WAVE' && '["OSWEC_MULTI_WAVE"]' || inputs.model == 'RM3_B2B' && '["RM3_B2B"]' || inputs.model == 'RM3_MOORING_MATRIX' && '["RM3_MOORING_MATRIX"]' || inputs.model == 'SPHERE_MEAN_DRIFT' && '["SPHERE_MEAN_DRIFT"]' || inputs.model == 'RM3_MCR_SEASTATE' && '["RM3_MCR_SEASTATE"]' || inputs.model == 'RM3_END_STOPS' && '["RM3_END_STOPS"]' || inputs.model == 'RM3_END_STOPS_STEP' && '["RM3_END_STOPS_STEP"]' || inputs.model == 'RM3_END_STOPS_STEP_FINE' && '["RM3_END_STOPS_STEP_FINE"]' || inputs.model == 'RM3_END_STOPS_STEP_FINER' && '["RM3_END_STOPS_STEP_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_FINE' && '["RM3_END_STOPS_FULL_FINE"]' || inputs.model == 'RM3_END_STOPS_FULL_FINER' && '["RM3_END_STOPS_FULL_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_PAIR' && '["RM3_END_STOPS_FULL_FINE", "RM3_END_STOPS_FULL_FINER"]' || '["SPHERE_ELEVATION_IMPORT", "MORISON_FIXED", "SPHERE_MOVING_MORISON", "SPHERE_MOVING_MORISON_WAVE", "RM3_DD_PTO", "RM3", "OSWEC", "OSWEC_Nonhydro", "GBM_BARGE", "SPHERE_VARIABLE_MASS", "SPHERE_REACTIVE_DDPTO", "OSWEC_FULL_DIR", "OSWEC_MULTI_WAVE", "OSWEC_PASSIVE_YAW", "Sphere", "SPHERE_MEAN_DRIFT", "RM3_B2B", "RM3_END_STOPS", "RM3_PTO_Extension", "RM3_Radiation_Options", "Sphere_Passive", "Sphere_PTO_Config", "Sphere_Reactive_PI", "Sphere_Declutching", "Sphere_Latching", "RM3_MCR", "RM3_MCR_ARRAY", "RM3_MCR_EXCEL", "RM3_MCR_MAT", "RM3_MCR_SEASTATE", "ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", "ELLIPSOID_NLH_ODE45"]') }} + model: ${{ fromJSON(inputs.model == 'RM3_MCR_ARRAY' && '["RM3_MCR_ARRAY"]' || inputs.model == 'SPHERE_ELEVATION_IMPORT' && '["SPHERE_ELEVATION_IMPORT"]' || inputs.model == 'MORISON_FIXED' && '["MORISON_FIXED"]' || inputs.model == 'SPHERE_MOVING_MORISON' && '["SPHERE_MOVING_MORISON"]' || inputs.model == 'SPHERE_MOVING_MORISON_WAVE' && '["SPHERE_MOVING_MORISON_WAVE"]' || inputs.model == 'RM3_DD_PTO' && '["RM3_DD_PTO"]' || inputs.model == 'SPHERE_MPC' && '["SPHERE_MPC"]' || inputs.model == 'GBM_BARGE' && '["GBM_BARGE"]' || inputs.model == 'SPHERE_REACTIVE_DDPTO' && '["SPHERE_REACTIVE_DDPTO"]' || inputs.model == 'SPHERE_VARIABLE_MASS' && '["SPHERE_VARIABLE_MASS"]' || inputs.model == 'ELLIPSOID_NLH_REG' && '["ELLIPSOID_NLH_REG"]' || inputs.model == 'ELLIPSOID_NLH_CIC' && '["ELLIPSOID_NLH_CIC"]' || inputs.model == 'ELLIPSOID_NLH_ODE45' && '["ELLIPSOID_NLH_ODE45"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' && '["OSWEC_PASSIVE_YAW_IRR_CONT"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR' && '["OSWEC_PASSIVE_YAW_IRR"]' || inputs.model == 'OSWEC_PASSIVE_YAW' && '["OSWEC_PASSIVE_YAW"]' || inputs.model == 'OSWEC_VARIABLE_YAW_2DEG' && '["OSWEC_VARIABLE_YAW_2DEG"]' || inputs.model == 'OSWEC_FULL_DIR' && '["OSWEC_FULL_DIR"]' || inputs.model == 'OSWEC_MULTI_WAVE' && '["OSWEC_MULTI_WAVE"]' || inputs.model == 'RM3_B2B' && '["RM3_B2B"]' || inputs.model == 'RM3_MOORING_MATRIX' && '["RM3_MOORING_MATRIX"]' || inputs.model == 'SPHERE_MEAN_DRIFT' && '["SPHERE_MEAN_DRIFT"]' || inputs.model == 'RM3_MCR_SEASTATE' && '["RM3_MCR_SEASTATE"]' || inputs.model == 'RM3_END_STOPS' && '["RM3_END_STOPS"]' || inputs.model == 'RM3_END_STOPS_STEP' && '["RM3_END_STOPS_STEP"]' || inputs.model == 'RM3_END_STOPS_STEP_FINE' && '["RM3_END_STOPS_STEP_FINE"]' || inputs.model == 'RM3_END_STOPS_STEP_FINER' && '["RM3_END_STOPS_STEP_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_FINE' && '["RM3_END_STOPS_FULL_FINE"]' || inputs.model == 'RM3_END_STOPS_FULL_FINER' && '["RM3_END_STOPS_FULL_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_PAIR' && '["RM3_END_STOPS_FULL_FINE", "RM3_END_STOPS_FULL_FINER"]' || '["SPHERE_ELEVATION_IMPORT", "MORISON_FIXED", "SPHERE_MOVING_MORISON", "SPHERE_MOVING_MORISON_WAVE", "RM3_DD_PTO", "RM3", "OSWEC", "OSWEC_Nonhydro", "GBM_BARGE", "SPHERE_VARIABLE_MASS", "SPHERE_REACTIVE_DDPTO", "OSWEC_FULL_DIR", "OSWEC_MULTI_WAVE", "OSWEC_PASSIVE_YAW", "Sphere", "SPHERE_MEAN_DRIFT", "RM3_B2B", "RM3_END_STOPS", "RM3_PTO_Extension", "RM3_Radiation_Options", "Sphere_Passive", "Sphere_PTO_Config", "Sphere_Reactive_PI", "Sphere_Declutching", "Sphere_Latching", "RM3_MCR", "RM3_MCR_ARRAY", "RM3_MCR_EXCEL", "RM3_MCR_MAT", "RM3_MCR_SEASTATE", "ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", "ELLIPSOID_NLH_ODE45"]') }} name: ${{ matrix.model }} steps: - uses: actions/checkout@v7 @@ -108,7 +108,7 @@ jobs: source examples/${{ matrix.model }} - uses: actions/checkout@v7 - if: matrix.model == 'SPHERE_ELEVATION_IMPORT' || matrix.model == 'MORISON_FIXED' || matrix.model == 'SPHERE_MOVING_MORISON' || matrix.model == 'SPHERE_MOVING_MORISON_WAVE' || matrix.model == 'RM3_DD_PTO' || matrix.model == 'SPHERE_MPC' || matrix.model == 'GBM_BARGE' || matrix.model == 'SPHERE_REACTIVE_DDPTO' || matrix.model == 'SPHERE_VARIABLE_MASS' || matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45' || matrix.model == 'Sphere' || matrix.model == 'SPHERE_MEAN_DRIFT' || matrix.model == 'Sphere_Passive' || matrix.model == 'Sphere_PTO_Config' || matrix.model == 'Sphere_Reactive_PI' || matrix.model == 'Sphere_Declutching' || matrix.model == 'Sphere_Latching' || matrix.model == 'RM3_B2B' || matrix.model == 'RM3_MOORING_MATRIX' || matrix.model == 'RM3_END_STOPS' || matrix.model == 'RM3_END_STOPS_STEP' || matrix.model == 'RM3_END_STOPS_STEP_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINER' || matrix.model == 'RM3_END_STOPS_STEP_FINER' || matrix.model == 'RM3_PTO_Extension' || matrix.model == 'RM3_Radiation_Options' || matrix.model == 'RM3_MCR' || matrix.model == 'RM3_MCR_ARRAY' || matrix.model == 'RM3_MCR_EXCEL' || matrix.model == 'RM3_MCR_MAT' || matrix.model == 'RM3_MCR_SEASTATE' || matrix.model == 'OSWEC_Nonhydro' || matrix.model == 'OSWEC_MULTI_WAVE' || matrix.model == 'OSWEC_FULL_DIR' || matrix.model == 'OSWEC_PASSIVE_YAW' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' + if: matrix.model == 'SPHERE_ELEVATION_IMPORT' || matrix.model == 'MORISON_FIXED' || matrix.model == 'SPHERE_MOVING_MORISON' || matrix.model == 'SPHERE_MOVING_MORISON_WAVE' || matrix.model == 'RM3_DD_PTO' || matrix.model == 'SPHERE_MPC' || matrix.model == 'GBM_BARGE' || matrix.model == 'SPHERE_REACTIVE_DDPTO' || matrix.model == 'SPHERE_VARIABLE_MASS' || matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45' || matrix.model == 'Sphere' || matrix.model == 'SPHERE_MEAN_DRIFT' || matrix.model == 'Sphere_Passive' || matrix.model == 'Sphere_PTO_Config' || matrix.model == 'Sphere_Reactive_PI' || matrix.model == 'Sphere_Declutching' || matrix.model == 'Sphere_Latching' || matrix.model == 'RM3_B2B' || matrix.model == 'RM3_MOORING_MATRIX' || matrix.model == 'RM3_END_STOPS' || matrix.model == 'RM3_END_STOPS_STEP' || matrix.model == 'RM3_END_STOPS_STEP_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINER' || matrix.model == 'RM3_END_STOPS_STEP_FINER' || matrix.model == 'RM3_PTO_Extension' || matrix.model == 'RM3_Radiation_Options' || matrix.model == 'RM3_MCR' || matrix.model == 'RM3_MCR_ARRAY' || matrix.model == 'RM3_MCR_EXCEL' || matrix.model == 'RM3_MCR_MAT' || matrix.model == 'RM3_MCR_SEASTATE' || matrix.model == 'OSWEC_Nonhydro' || matrix.model == 'OSWEC_MULTI_WAVE' || matrix.model == 'OSWEC_FULL_DIR' || matrix.model == 'OSWEC_PASSIVE_YAW' || matrix.model == 'OSWEC_VARIABLE_YAW_2DEG' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' with: repository: WEC-Sim/WEC-Sim_Applications ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e @@ -128,6 +128,7 @@ jobs: Passive_Yaw/PassiveYawOFF Passive_Yaw/PassiveYawON Passive_Yaw/PassiveYawRegression + Variable_Hydro/Passive_Yaw Body-to-Body_Interactions/B2B_Case1 Body-to-Body_Interactions/B2B_Case2 Body-to-Body_Interactions/B2B_Case3 @@ -402,6 +403,11 @@ jobs: env: WEC_SIM_SPHERE_H5: applications/_Common_Input_Files/Sphere/hydroData/sphere.h5 WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output + - run: python -m pytest -q tests/test_variable_yaw_heading_bank.py + if: matrix.model == 'OSWEC_VARIABLE_YAW_2DEG' + env: + WEC_SIM_APPLICATIONS_DIR: applications + WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - run: python -m pytest -q tests/test_sphere_pto_config_parity.py if: matrix.model == 'Sphere_PTO_Config' env: diff --git a/tests/matlab/reference_model_baseline.m b/tests/matlab/reference_model_baseline.m index b100e57..7dffd65 100644 --- a/tests/matlab/reference_model_baseline.m +++ b/tests/matlab/reference_model_baseline.m @@ -315,6 +315,33 @@ function reference_model_baseline(model) end cases = ["PassiveYawOFF", "PassiveYawON"]; caseDirs = fullfile(repoRoot, 'applications', 'Passive_Yaw', cases); + case "OSWEC_VARIABLE_YAW_2DEG" + commonHydro = fullfile(repoRoot, 'applications', ... + '_Common_Input_Files', 'OSWEC', 'hydroData'); + cd(commonHydro); + if ~isfile('oswec.h5') + bemio; + end + sourceDir = fullfile(repoRoot, 'applications', ... + 'Variable_Hydro', 'Passive_Yaw'); + caseDir = fullfile(repoRoot, 'applications', ... + 'Variable_Hydro', 'paired_Passive_Yaw_2deg'); + [copied, copyMessage] = copyfile(sourceDir, caseDir); + assert(copied, copyMessage); + bankScript = fullfile(caseDir, 'hydroData', 'bemio.m'); + contents = fileread(bankScript); + oldGrid = 'newDirs = -40:0.05:40;'; + assert(contains(contents, oldGrid), ... + 'The pinned variable-yaw bank generation changed'); + contents = strrep(contents, oldGrid, 'newDirs = -40:2:40;'); + fid = fopen(bankScript, 'w'); + assert(fid ~= -1, 'Could not write the two-degree bank script'); + fprintf(fid, '%s', contents); + fclose(fid); + cd(fullfile(caseDir, 'hydroData')); + bemio; + cases = "regular_2deg"; + caseDirs = string(caseDir); case {"OSWEC_PASSIVE_YAW_IRR", "OSWEC_PASSIVE_YAW_IRR_CONT"} hydroDir = fullfile(repoRoot, 'applications', '_Common_Input_Files', ... 'OSWEC', 'hydroData'); @@ -484,6 +511,13 @@ function reference_model_baseline(model) for iCase = 1:numel(cases) cd(caseDirs(iCase)); + if string(model) == "OSWEC_VARIABLE_YAW_2DEG" + waveFlag = 'regular'; + bemDirections = -40:2:40; + assert(numel(bemDirections) == 41 && ... + isfile('hydroData/oswec_10.h5'), ... + 'The two-degree hydrodynamic bank is incomplete'); + end if string(model) == "SPHERE_MPC" instrument_sphere_mpc(); end @@ -910,6 +944,27 @@ function reference_model_baseline(model) writematrix(waves.waveAmpTime, fullfile(outDir, sprintf( ... 'OSWEC_PASSIVE_YAW_%s_wave.csv', cases(iCase)))); end + if string(model) == "OSWEC_VARIABLE_YAW_2DEG" + assert(simu.dt == 0.01 && simu.endTime == 600 && ... + simu.rampTime == 100 && simu.cicEndTime == 40 && ... + strcmp(waves.type, 'regular') && waves.height == 2.5 && ... + waves.period == 8 && waves.direction == 10 && ... + body(1).mass == 12700 && body(2).mass == 999 && ... + body(1).variableHydro.option == 1 && ... + pto(1).damping == 120000 && ... + isequal(pto(1).location, [0 0 -8.9]) && ... + isequal(constraint(1).location, [0 0 -10]), ... + 'The pinned variable-yaw settings changed'); + index = output.bodies(1).hydroForceIndex(:); + assert(all(index >= 1 & index <= numel(bemDirections)) && ... + all(index == round(index)), ... + 'The variable-yaw selected indices are invalid'); + writematrix([output.bodies(1).time(:), index, ... + bemDirections(index(:))], fullfile(outDir, ... + 'OSWEC_VARIABLE_YAW_2DEG_selected_heading.csv')); + writematrix(waves.waveAmpTime, fullfile(outDir, ... + 'OSWEC_VARIABLE_YAW_2DEG_wave.csv')); + end if any(string(model) == ["OSWEC_PASSIVE_YAW_IRR", ... "OSWEC_PASSIVE_YAW_IRR_CONT"]) expectedThreshold = double(string(model) == "OSWEC_PASSIVE_YAW_IRR"); diff --git a/tests/test_variable_yaw_heading_bank.py b/tests/test_variable_yaw_heading_bank.py new file mode 100644 index 0000000..d3d0ed5 --- /dev/null +++ b/tests/test_variable_yaw_heading_bank.py @@ -0,0 +1,117 @@ +"""Nearest-heading selection used by the published variable-hydro yaw case.""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +from wecsim import RegularWave, WEC, WorldPoint +from wecsim.passiveYaw import NearestHeadingExcitation, PassiveYawExcitation + + +def _model(): + headings = np.arange(0.0, 360.0, 10.0) + real = np.zeros((6, len(headings))) + real[5] = headings + return PassiveYawExcitation( + headings, real, np.zeros_like(real), + incident_direction=10, omega=1, amplitude=1, ramp_time=0, + ) + + +def test_nearest_bank_uses_source_relative_angle_and_first_tie(): + bank = NearestHeadingExcitation(_model(), np.arange(-40, 41, 2)) + assert bank.heading(0) == 10 + assert bank.heading(np.deg2rad(9)) == 0 + assert bank.heading(np.deg2rad(-31)) == 40 + assert bank.heading(np.deg2rad(51)) == -40 + assert bank.heading(np.deg2rad(9)) == 0 # a one-degree tie chooses lower + + yaw = np.deg2rad(9.2) + actual = bank.force(0, yaw) + expected = _model().force(0, yaw, coefficient_heading=0) + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-13) + assert actual[5] == 0 + + +def test_nearest_bank_rejects_invalid_grid(): + model = _model() + for headings in ([0], [0, 0], [2, 1], [-181, 0], [0, np.nan]): + with pytest.raises(ValueError, match="heading bank"): + NearestHeadingExcitation(model, headings) + + +@pytest.mark.skipif( + not (os.environ.get("WEC_SIM_APPLICATIONS_DIR") + and os.environ.get("WEC_SIM_MATLAB_MODEL_OUTPUT_DIR")), + reason="fresh variable-yaw MATLAB output not provided", +) +def test_two_degree_bank_against_pinned_matlab(): + applications = Path(os.environ["WEC_SIM_APPLICATIONS_DIR"]) + source = Path(os.environ["WEC_SIM_MATLAB_MODEL_OUTPUT_DIR"]) + prefix = "OSWEC_VARIABLE_YAW_2DEG" + flap = np.loadtxt(source / f"{prefix}_regular_2deg_body1.csv", delimiter=",") + base = np.loadtxt(source / f"{prefix}_regular_2deg_body2.csv", delimiter=",") + pto = np.loadtxt(source / f"{prefix}_regular_2deg_pto1.csv", delimiter=",") + selected = np.loadtxt(source / f"{prefix}_selected_heading.csv", delimiter=",") + wave = np.loadtxt(source / f"{prefix}_wave.csv", delimiter=",") + hydro = applications / "_Common_Input_Files/OSWEC/hydroData/oswec.h5" + + wec = WEC("OSWEC two-degree variable yaw") + moving = wec.body( + "flap", hydro, mass=12700, inertia=(1.85e6,) * 3, + passive_yaw=True, yaw_heading_bank=np.arange(-40, 41, 2), + ) + wec.body("base", hydro, mass=999, inertia=(999,) * 3) + yaw = wec.coordinate("yaw", moving.move("yaw", pivot=WorldPoint(0, 0, -8.9))) + wec.rotational_pto("hinge", yaw, damping=120000) + result = wec.run(RegularWave(2.5, 8, 10), dt=0.01, + end_time=600, ramp_time=100) + + assert flap.shape == base.shape == pto.shape == (60001, 25) + assert selected.shape == (60001, 3) + assert wave.shape == (60001, 2) + np.testing.assert_allclose(result.time, flap[:, 0], rtol=0, atol=1e-10) + np.testing.assert_allclose(result.wave_elevation, wave[:, 1], rtol=0, atol=1e-12) + np.testing.assert_allclose(result.bodies["base"].position, + base[:, 1:7], rtol=0, atol=1e-10) + np.testing.assert_allclose(result.bodies["base"].velocity, + base[:, 7:13], rtol=0, atol=1e-10) + + model = PassiveYawExcitation.from_hydro_data( + _hydro_data(hydro), omega=2 * np.pi / 8, + incident_direction=10, amplitude=1.25, ramp_time=100, + rho=1000, g=9.81, + ) + bank = NearestHeadingExcitation(model, np.arange(-40, 41, 2)) + source_heading = np.array([bank.heading(angle) for angle in flap[:, 6]]) + heading_error = np.max(np.abs(source_heading - selected[:, 2])) + print(f"source heading selection max error: {heading_error:.6g} deg") + assert heading_error <= 2.0 # source output may record the prior accepted step + + source_path_force = np.stack([ + bank.force(t, angle) for t, angle in zip(result.time, flap[:, 6]) + ]) + force_error = np.max(np.abs(source_path_force[:, 5] - flap[:, 24])) + yaw_error = np.max(np.abs(result.bodies["flap"].position[:, 5] - flap[:, 6])) + speed_error = np.max(np.abs(result.bodies["flap"].velocity[:, 5] - flap[:, 12])) + torque_error = np.max(np.abs(result.ptos["hinge"].force - pto[:, 17])) + print(f"source-path yaw excitation max error: {force_error:.6g} N m") + print(f"trajectory max errors: {yaw_error:.6g} rad, " + f"{speed_error:.6g} rad/s, {torque_error:.6g} N m") + assert np.isfinite([force_error, yaw_error, speed_error, torque_error]).all() + assert force_error < 1000 + assert yaw_error < 0.02 + assert speed_error < 0.005 + assert torque_error < 600 + + +def _hydro_data(path): + from wecsim.bodyClass import BodyClass + + body = BodyClass(str(path)) + body.bodyNumber = 1 + body.bodyTotal = 2 + body.readH5file() + return body.hydroData diff --git a/wecsim/api.py b/wecsim/api.py index c445d63..5396410 100644 --- a/wecsim/api.py +++ b/wecsim/api.py @@ -86,6 +86,7 @@ class Body: drag_area: float = 0.0 variable_hydro: VariableHydro | None = None passive_yaw_threshold: float = 0.0 + yaw_heading_bank: tuple[float, ...] | None = None fixed: bool = False center_gravity: tuple[float, float, float] | None = None volume: float = 0.0 @@ -376,6 +377,7 @@ def body(self, name: str, hydro_file: str | Path, *, mean_drift: str = "none", passive_yaw: bool = False, passive_yaw_threshold: float = 0.0, + yaw_heading_bank: Sequence[float] | None = None, geometry_file: str | Path | None = None, nonlinear_hydro: str | None = None, drag_coefficient: float = 0.0, @@ -388,10 +390,20 @@ def body(self, name: str, hydro_file: str | Path, *, or passive_yaw_threshold < 0 or (passive_yaw_threshold > 0 and passive_yaw is not True)): raise ValueError("passive_yaw_threshold needs a nonnegative degree value and passive_yaw=True") + if yaw_heading_bank is not None: + headings = np.asarray(yaw_heading_bank, dtype=float) + if (not passive_yaw or passive_yaw_threshold + or headings.ndim != 1 or len(headings) < 2 + or not np.isfinite(headings).all() + or not np.all(np.diff(headings) > 0) + or headings[0] < -180 or headings[-1] > 180): + raise ValueError("yaw_heading_bank needs passive_yaw, no threshold, and ordered directions in [-180, 180]") + yaw_heading_bank = tuple(float(value) for value in headings) body = Body(name, hydro_file, mass, tuple(inertia), hydro_body, mean_drift, passive_yaw, geometry_file, nonlinear_hydro, drag_coefficient, drag_area, - passive_yaw_threshold=passive_yaw_threshold) + passive_yaw_threshold=passive_yaw_threshold, + yaw_heading_bank=yaw_heading_bank) self.bodies.append(body) return body @@ -633,6 +645,8 @@ def to_case( body_case["passive_yaw"] = True if body.passive_yaw_threshold: body_case["passive_yaw_threshold"] = body.passive_yaw_threshold + if body.yaw_heading_bank is not None: + body_case["yaw_heading_bank"] = list(body.yaw_heading_bank) if body.geometry_file is not None: body_case["geometry_file"] = str(body.geometry_file) if body.nonlinear_hydro is not None: @@ -736,7 +750,8 @@ def _floating_joint_case(self, wave, simulation, or body.mean_drift != "none" or body.passive_yaw or body.geometry_file is not None or body.nonlinear_hydro is not None or body.variable_hydro is not None or body.drag_coefficient - or body.drag_area or body.passive_yaw_threshold): + or body.drag_area or body.passive_yaw_threshold + or body.yaw_heading_bank is not None): raise ValueError("floating_joint needs equilibrium-mass hydrodynamic bodies without extra force models") inertia = np.asarray(body.inertia, dtype=float) if inertia.shape != (3,) or not np.isfinite(inertia).all() or inertia[1] <= 0: diff --git a/wecsim/caseDynamics.py b/wecsim/caseDynamics.py index e8d1372..b1cae58 100644 --- a/wecsim/caseDynamics.py +++ b/wecsim/caseDynamics.py @@ -38,7 +38,8 @@ ) from .nonlinearHydro import HeaveMeshHydro from .passiveYaw import ( - HeldPassiveYawExcitation, PassiveYawExcitation, SampledPassiveYawExcitation, + HeldPassiveYawExcitation, NearestHeadingExcitation, + PassiveYawExcitation, SampledPassiveYawExcitation, ) from .ptoConnections import build_linear_ptos from .rm3Regular import solve_rm3_regular @@ -133,7 +134,8 @@ def _hydro_file(body, base_dir): _section(body, "body", {"hydro_file"}, {"hydro_file", "hydro_body", "mass", "pitch_inertia", "inertia", "coordinate_map", "name", "mean_drift", "fixed", - "passive_yaw", "passive_yaw_threshold", "geometry_file", "nonlinear_hydro", + "passive_yaw", "passive_yaw_threshold", "yaw_heading_bank", + "geometry_file", "nonlinear_hydro", "drag_coefficient", "drag_area", "variable_hydro", "morison_elements"}) raw = body["hydro_file"] if not isinstance(raw, str) or not raw: @@ -1348,6 +1350,12 @@ def moving_terms(at_time, coordinate, speed): or wave["type"] != "pm") for index, threshold in enumerate(yaw_thresholds)): raise ValueError("positive passive_yaw_threshold needs PM passive yaw") + yaw_banks = [spec.get("yaw_heading_bank") for spec in bodies] + if any(bank is not None and (not bodies[index].get("passive_yaw", False) + or yaw_thresholds[index] + or wave["type"] != "regular") + for index, bank in enumerate(yaw_banks)): + raise ValueError("yaw_heading_bank needs regular-wave passive yaw without a threshold") if passive_indices: yaw_map = np.zeros((6, 1)) yaw_map[5, 0] = 1 @@ -1596,6 +1604,10 @@ def excitation(at_time): incident_direction=direction, amplitude=height / 2, ramp_time=ramp_time, rho=rho, g=g, ) + if yaw_banks[index - 1] is not None: + passive_model = NearestHeadingExcitation( + passive_model, yaw_banks[index - 1], + ) if isinstance(passive_model, HeldPassiveYawExcitation): state_excitation = passive_model diff --git a/wecsim/passiveYaw.py b/wecsim/passiveYaw.py index 1832b6a..6526182 100644 --- a/wecsim/passiveYaw.py +++ b/wecsim/passiveYaw.py @@ -60,9 +60,13 @@ def from_hydro_data(cls, hydro_data, *, omega, incident_direction, return cls(headings, real, imaginary, incident_direction, omega, amplitude, ramp_time) - def force(self, time: float, yaw: float) -> np.ndarray: + def force(self, time: float, yaw: float, *, + coefficient_heading: float | None = None) -> np.ndarray: """Return the six-component excitation in the world frame.""" - relative_heading = (self.incident_direction - np.degrees(yaw)) % 360 + relative_heading = ( + self.incident_direction - np.degrees(yaw) + if coefficient_heading is None else coefficient_heading + ) % 360 real = np.array([ np.interp(relative_heading, self.directions, row, period=360) for row in self.real @@ -85,6 +89,37 @@ def force(self, time: float, yaw: float) -> np.ndarray: return world +@dataclass(frozen=True) +class NearestHeadingExcitation: + """Select the nearest BEM heading as in the variable-hydro yaw example. + + The published direction-bank files share mass, restoring, and radiation + coefficients; their excitation coefficients differ by heading. This + selector uses the full-direction HDF5 input for those coefficients. + """ + + model: PassiveYawExcitation + headings: np.ndarray + + def __post_init__(self): + headings = np.asarray(self.headings, dtype=float) + if (headings.ndim != 1 or headings.size < 2 + or not np.isfinite(headings).all() + or not np.all(np.diff(headings) > 0) + or headings[0] < -180 or headings[-1] > 180): + raise ValueError("heading bank needs ordered directions in [-180, 180]") + object.__setattr__(self, "headings", headings) + + def heading(self, yaw: float) -> float: + relative = self.model.incident_direction - np.degrees(yaw) + return float(self.headings[np.argmin(np.abs(self.headings - relative))]) + + def force(self, time: float, yaw: float) -> np.ndarray: + return self.model.force( + time, yaw, coefficient_heading=self.heading(yaw), + ) + + @dataclass(frozen=True) class SampledPassiveYawExcitation: """Broadband force at each BEM heading, interpolated at the body's yaw. From 5cac416137c2f8be1a2920c0a27e2d9e2b892356 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 10:17:50 -0400 Subject: [PATCH 2/4] Write selected yaw headings as a column --- tests/matlab/reference_model_baseline.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/matlab/reference_model_baseline.m b/tests/matlab/reference_model_baseline.m index 7dffd65..d07216b 100644 --- a/tests/matlab/reference_model_baseline.m +++ b/tests/matlab/reference_model_baseline.m @@ -959,8 +959,9 @@ function reference_model_baseline(model) assert(all(index >= 1 & index <= numel(bemDirections)) && ... all(index == round(index)), ... 'The variable-yaw selected indices are invalid'); + selectedDirection = bemDirections(index(:)); writematrix([output.bodies(1).time(:), index, ... - bemDirections(index(:))], fullfile(outDir, ... + selectedDirection(:)], fullfile(outDir, ... 'OSWEC_VARIABLE_YAW_2DEG_selected_heading.csv')); writematrix(waves.waveAmpTime, fullfile(outDir, ... 'OSWEC_VARIABLE_YAW_2DEG_wave.csv')); From 943ac3cf9ee2244f87b7b0989c93bb13967bed82 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 10:29:23 -0400 Subject: [PATCH 3/4] Gate variable yaw force and pre-event motion, document trajectory gap --- .github/workflows/reference-model-baselines.yml | 2 +- PARITY.md | 1 + README.md | 7 +++++++ tests/test_variable_yaw_heading_bank.py | 16 ++++++++++++---- 4 files changed, 21 insertions(+), 5 deletions(-) diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index 8b30dc4..e4d45a3 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -403,7 +403,7 @@ jobs: env: WEC_SIM_SPHERE_H5: applications/_Common_Input_Files/Sphere/hydroData/sphere.h5 WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - - run: python -m pytest -q tests/test_variable_yaw_heading_bank.py + - run: python -m pytest -q -s tests/test_variable_yaw_heading_bank.py if: matrix.model == 'OSWEC_VARIABLE_YAW_2DEG' env: WEC_SIM_APPLICATIONS_DIR: applications diff --git a/PARITY.md b/PARITY.md index 7ef188b..c677cd0 100644 --- a/PARITY.md +++ b/PARITY.md @@ -57,6 +57,7 @@ the production Python code. The live wave comparison passed on 6 October | OSWEC `Multiple_Wave_Spectra` with two independent PM seas | Pinned MATLAB Applications `Multiple_Wave_Spectra` input and BEMIO-generated OSWEC HDF5 | The public `run_case` interface sums separate 2 m, 0° and 1 m, 90° seas with saved independent phase realizations over 100 s. Combined elevation agrees within `1e-11` m and both bodies' six excitation-force components within `1e-6` N or N m. The fixed hydrodynamic base remains stationary. With the ordinary implicit added-mass default, maximum flap differences are 106 mm surge, 10.8 mm heave, and 0.0213 rad pitch; PTO torque differs by at most 7.25 kN m. The explicit source-delay comparison reduces these to 10.3 mm, 1.20 mm, 0.00207 rad, and 818 N m. Both dynamics paths have separate paired gates. The source enables passive-yaw preprocessing, but its hinge has zero yaw motion; the source's spline frequency interpolation is selected explicitly for this case. Moving-yaw dynamics and fixed-base reactions remain unsupported. | | OSWEC `Full_Directional_Waves` imported spectrum | Pinned MATLAB Applications `Full_Directional_Waves` input and BEMIO-generated OSWEC HDF5 | The imported 47-frequency, 180-heading spectrum and realized phases reproduce all 8,001 wave-elevation samples within `9e-14` m. The source force block omits heading-bin width although its wave-elevation block includes it: logged excitation is about `1/sqrt(2° in radians) = 5.352` times the physically integrated force. Python keeps heading integration as its default. Explicit `matlab_omitted` quadrature plus `spline_frequency` interpolation matches all six logged excitation components for both bodies within `4e-7` N or N m. With that source forcing and ordinary implicit added mass, the 400 s flap differs by at most 12.8 mm surge, 3.52 mm heave, 0.00258 rad pitch, and 56.6 N m PTO torque. The separate opt-in Simulink-delay comparison reduces these to 1.53 mm, 0.293 mm, 0.000306 rad, and 5.54 N m. Source-compatible agreement does not establish that the omitted-width force is physically correct. The fixed base remains stationary; its ground reactions are not computed. | | OSWEC passive yaw off and on | Pinned MATLAB Applications `PassiveYawOFF` and `PassiveYawON` inputs, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectories | The public `WEC` interface defines a yawing flap, fixed hydrodynamic base, and 120,000 N m s/rad rotational PTO. Both published 600 s cases use 2.5 m, 8 s waves at 10°. With passive yaw off, Python and MATLAB flap yaw agree within `4e-14` rad. With passive yaw on, excitation follows the body-relative wave heading and turns the flap toward the 10° incident heading; yaw and angular-speed differences stay below 0.00149 rad and 0.000193 rad/s. PTO torque and absorbed power differ by at most 23.2 N m and 0.128 W, accounting for MATLAB's negative absorbed-power sign. On the saved MATLAB yaw trajectory, all six reconstructed excitation components differ by at most 100 N or N m except yaw moment, which differs by 467 N m. The source holds heading coefficients until yaw changes by 0.01°, while Python interpolates continuously. The fixed base is stationary; full three-dimensional joint mechanics remain unsupported. | +| OSWEC variable-hydrodynamics passive yaw, 2° bank | Pinned MATLAB Applications `Variable_Hydro/Passive_Yaw` with its published regular-wave, 600 s, 0.01 s case and the first published `-40:2:40` direction bank; fresh R2025b [paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37791308393) | The Python `WEC.body(..., passive_yaw=True, yaw_heading_bank=range(-40, 41, 2))` selects the nearest body-relative BEM heading using the source rule. The selected heading matches MATLAB at all 60,001 saved MATLAB states; yaw excitation evaluated on that path differs by at most 867 N m. Before the first heading-event split at 64.57 s, yaw, yaw speed, and PTO torque differ by at most `3.89e-5` rad, `3.32e-5` rad/s, and 3.98 N m. At 64.57 s the paths differ by only `3.7e-5` rad but select different headings for one sample. Discrete events amplify that offset; over the full run the maximum yaw, speed, and PTO torque differences are 0.164 rad, 0.0221 rad/s, and 2.65 kN m. Full-trajectory parity is **unestablished**. This mode covers the published pure-yaw/fixed-base regular-wave bank, not general variable hydrodynamics. | | OSWEC irregular passive yaw | Pinned MATLAB Applications `PassiveYawRegression` input, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectory | The public `PMWave` replays the published 500-frequency phase realization over 250 s with 40 s radiation memory. Frequency bins, spectral amplitudes, and widths agree within `6e-15`; elevation agrees within `1e-13` m. The public `passive_yaw_threshold=1` selects the published coefficient hold and nearby BEM-heading snap. On MATLAB's saved yaw path, that opt-in force law matches all six logged excitation components within `1e-7` N or N m. The source's yaw radiation convolution reconstructs from its saved speed and pinned HDF5 within `5.4e-9` N m; added mass, hydrodynamic force balance, and rigid inertia balance close within `1.4e-8` N m. In a native Python run a heading update occurs one 0.01 s sample early at 53.45 s: MATLAB's relative heading is 10.99980° and Python's is 11.00004°, on opposite sides of the 11° update boundary. The yaw difference there is `4.3e-6` rad, and replaying MATLAB's logged force produces nearly the same pre-event drift (`4.29e-6` rad). Later event differences accumulate to 0.203 rad maximum yaw and 0.0277 rad/s speed differences over 250 s. With continuous heading interpolation, Python differs from the published trajectory by up to 0.452 rad yaw and 0.101 rad/s yaw speed. Replaying MATLAB's logged six-component force through the Python dynamics reduces these maximum differences to 0.00180 rad and 0.0000935 rad/s. The source force laws and balances are paired, but published-case trajectory parity with native excitation remains **unestablished** because the threshold amplifies small integration differences. | | OSWEC irregular passive yaw, continuous-heading control | Same pinned MATLAB input and phase realization, with only `body(1:2).yaw.threshold` changed from 1° to 0° in a temporary application copy | MATLAB and Python then agree over all 25,001 time samples: maximum differences are `1.1e-7` N or N m across the six excitation components on the MATLAB yaw path, `3.16e-5` rad flap/PTO angle, `7.30e-6` rad/s yaw/PTO speed, `0.875` N m PTO torque, and `0.0123` W source-signed power. Wave elevation differs by less than `1e-13` m and the fixed base remains stationary. This control validates the continuous-heading Python dynamics; it does not erase the published 1° threshold difference. | | OSWEC fixed nonhydrodynamic base | Pinned MATLAB Applications `Nonhydro_Body` case and its BEMIO-generated OSWEC HDF5 | The regular-wave solver reports the stationary base and nonlinear flap motion about the PTO hinge. Against a fresh 400 s MATLAB R2025b run (4,001 samples), maximum flap position differences are 20.8 mm surge, 7.9 mm heave, and 0.00444 rad pitch; velocity differences are 15.6 mm/s surge, 6.8 mm/s heave, and 0.00328 rad/s pitch. All six excitation-force components differ by less than 1 N, the base position and velocity agree exactly, and zero PTO torque differs only by MATLAB numerical noise below `3.3e-7` N m. The base's ground-constraint reaction forces are not calculated. | diff --git a/README.md b/README.md index 27d9158..6af8778 100644 --- a/README.md +++ b/README.md @@ -207,6 +207,13 @@ keeps continuous interpolation. The sampled setting is available for PM waves with one incident direction; small trajectory differences can change its update sample and accumulate over long runs. +For the published `Variable_Hydro/Passive_Yaw` regular-wave case, pass +`yaw_heading_bank=range(-40, 41, 2)` alongside `passive_yaw=True` on the +flap. This chooses the nearest 2° BEM heading from the wave direction relative +to yaw. The source-path force and heading selection are paired to MATLAB; +the long trajectory remains sensitive to one-sample heading changes (see +[`PARITY.md`](PARITY.md)). + The published `Morison_Element/morisonElement` application has a fixed monopile and tower without HDF5 hydrodynamic bodies. Its nonlinear drag and fluid-inertia force can be configured with a body-local point: diff --git a/tests/test_variable_yaw_heading_bank.py b/tests/test_variable_yaw_heading_bank.py index d3d0ed5..672e7a2 100644 --- a/tests/test_variable_yaw_heading_bank.py +++ b/tests/test_variable_yaw_heading_bank.py @@ -88,7 +88,7 @@ def test_two_degree_bank_against_pinned_matlab(): source_heading = np.array([bank.heading(angle) for angle in flap[:, 6]]) heading_error = np.max(np.abs(source_heading - selected[:, 2])) print(f"source heading selection max error: {heading_error:.6g} deg") - assert heading_error <= 2.0 # source output may record the prior accepted step + assert heading_error == 0 source_path_force = np.stack([ bank.force(t, angle) for t, angle in zip(result.time, flap[:, 6]) @@ -102,9 +102,17 @@ def test_two_degree_bank_against_pinned_matlab(): f"{speed_error:.6g} rad/s, {torque_error:.6g} N m") assert np.isfinite([force_error, yaw_error, speed_error, torque_error]).all() assert force_error < 1000 - assert yaw_error < 0.02 - assert speed_error < 0.005 - assert torque_error < 600 + + # Before the first heading event separates the two paths, the motion + # remains paired. The source's discontinuous selector magnifies the + # subsequent one-sample event offset; full-run trajectory parity is open. + before_split = result.time < 64.57 + assert np.max(np.abs(result.bodies["flap"].position[before_split, 5] + - flap[before_split, 6])) < 5e-5 + assert np.max(np.abs(result.bodies["flap"].velocity[before_split, 5] + - flap[before_split, 12])) < 5e-5 + assert np.max(np.abs(result.ptos["hinge"].force[before_split] + - pto[before_split, 17])) < 5 def _hydro_data(path): From eb45677c28c94ebfd498ca0c3a9d3c25cdc714ed Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 10:37:44 -0400 Subject: [PATCH 4/4] Point variable yaw parity note to passing paired run --- PARITY.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PARITY.md b/PARITY.md index c677cd0..f718b21 100644 --- a/PARITY.md +++ b/PARITY.md @@ -57,7 +57,7 @@ the production Python code. The live wave comparison passed on 6 October | OSWEC `Multiple_Wave_Spectra` with two independent PM seas | Pinned MATLAB Applications `Multiple_Wave_Spectra` input and BEMIO-generated OSWEC HDF5 | The public `run_case` interface sums separate 2 m, 0° and 1 m, 90° seas with saved independent phase realizations over 100 s. Combined elevation agrees within `1e-11` m and both bodies' six excitation-force components within `1e-6` N or N m. The fixed hydrodynamic base remains stationary. With the ordinary implicit added-mass default, maximum flap differences are 106 mm surge, 10.8 mm heave, and 0.0213 rad pitch; PTO torque differs by at most 7.25 kN m. The explicit source-delay comparison reduces these to 10.3 mm, 1.20 mm, 0.00207 rad, and 818 N m. Both dynamics paths have separate paired gates. The source enables passive-yaw preprocessing, but its hinge has zero yaw motion; the source's spline frequency interpolation is selected explicitly for this case. Moving-yaw dynamics and fixed-base reactions remain unsupported. | | OSWEC `Full_Directional_Waves` imported spectrum | Pinned MATLAB Applications `Full_Directional_Waves` input and BEMIO-generated OSWEC HDF5 | The imported 47-frequency, 180-heading spectrum and realized phases reproduce all 8,001 wave-elevation samples within `9e-14` m. The source force block omits heading-bin width although its wave-elevation block includes it: logged excitation is about `1/sqrt(2° in radians) = 5.352` times the physically integrated force. Python keeps heading integration as its default. Explicit `matlab_omitted` quadrature plus `spline_frequency` interpolation matches all six logged excitation components for both bodies within `4e-7` N or N m. With that source forcing and ordinary implicit added mass, the 400 s flap differs by at most 12.8 mm surge, 3.52 mm heave, 0.00258 rad pitch, and 56.6 N m PTO torque. The separate opt-in Simulink-delay comparison reduces these to 1.53 mm, 0.293 mm, 0.000306 rad, and 5.54 N m. Source-compatible agreement does not establish that the omitted-width force is physically correct. The fixed base remains stationary; its ground reactions are not computed. | | OSWEC passive yaw off and on | Pinned MATLAB Applications `PassiveYawOFF` and `PassiveYawON` inputs, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectories | The public `WEC` interface defines a yawing flap, fixed hydrodynamic base, and 120,000 N m s/rad rotational PTO. Both published 600 s cases use 2.5 m, 8 s waves at 10°. With passive yaw off, Python and MATLAB flap yaw agree within `4e-14` rad. With passive yaw on, excitation follows the body-relative wave heading and turns the flap toward the 10° incident heading; yaw and angular-speed differences stay below 0.00149 rad and 0.000193 rad/s. PTO torque and absorbed power differ by at most 23.2 N m and 0.128 W, accounting for MATLAB's negative absorbed-power sign. On the saved MATLAB yaw trajectory, all six reconstructed excitation components differ by at most 100 N or N m except yaw moment, which differs by 467 N m. The source holds heading coefficients until yaw changes by 0.01°, while Python interpolates continuously. The fixed base is stationary; full three-dimensional joint mechanics remain unsupported. | -| OSWEC variable-hydrodynamics passive yaw, 2° bank | Pinned MATLAB Applications `Variable_Hydro/Passive_Yaw` with its published regular-wave, 600 s, 0.01 s case and the first published `-40:2:40` direction bank; fresh R2025b [paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37791308393) | The Python `WEC.body(..., passive_yaw=True, yaw_heading_bank=range(-40, 41, 2))` selects the nearest body-relative BEM heading using the source rule. The selected heading matches MATLAB at all 60,001 saved MATLAB states; yaw excitation evaluated on that path differs by at most 867 N m. Before the first heading-event split at 64.57 s, yaw, yaw speed, and PTO torque differ by at most `3.89e-5` rad, `3.32e-5` rad/s, and 3.98 N m. At 64.57 s the paths differ by only `3.7e-5` rad but select different headings for one sample. Discrete events amplify that offset; over the full run the maximum yaw, speed, and PTO torque differences are 0.164 rad, 0.0221 rad/s, and 2.65 kN m. Full-trajectory parity is **unestablished**. This mode covers the published pure-yaw/fixed-base regular-wave bank, not general variable hydrodynamics. | +| OSWEC variable-hydrodynamics passive yaw, 2° bank | Pinned MATLAB Applications `Variable_Hydro/Passive_Yaw` with its published regular-wave, 600 s, 0.01 s case and the first published `-40:2:40` direction bank; fresh R2025b [paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37792909630) | The Python `WEC.body(..., passive_yaw=True, yaw_heading_bank=range(-40, 41, 2))` selects the nearest body-relative BEM heading using the source rule. The selected heading matches MATLAB at all 60,001 saved MATLAB states; yaw excitation evaluated on that path differs by at most 867 N m. Before the first heading-event split at 64.57 s, yaw, yaw speed, and PTO torque differ by at most `3.89e-5` rad, `3.32e-5` rad/s, and 3.98 N m. At 64.57 s the paths differ by only `3.7e-5` rad but select different headings for one sample. Discrete events amplify that offset; over the full run the maximum yaw, speed, and PTO torque differences are 0.164 rad, 0.0221 rad/s, and 2.65 kN m. Full-trajectory parity is **unestablished**. This mode covers the published pure-yaw/fixed-base regular-wave bank, not general variable hydrodynamics. | | OSWEC irregular passive yaw | Pinned MATLAB Applications `PassiveYawRegression` input, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectory | The public `PMWave` replays the published 500-frequency phase realization over 250 s with 40 s radiation memory. Frequency bins, spectral amplitudes, and widths agree within `6e-15`; elevation agrees within `1e-13` m. The public `passive_yaw_threshold=1` selects the published coefficient hold and nearby BEM-heading snap. On MATLAB's saved yaw path, that opt-in force law matches all six logged excitation components within `1e-7` N or N m. The source's yaw radiation convolution reconstructs from its saved speed and pinned HDF5 within `5.4e-9` N m; added mass, hydrodynamic force balance, and rigid inertia balance close within `1.4e-8` N m. In a native Python run a heading update occurs one 0.01 s sample early at 53.45 s: MATLAB's relative heading is 10.99980° and Python's is 11.00004°, on opposite sides of the 11° update boundary. The yaw difference there is `4.3e-6` rad, and replaying MATLAB's logged force produces nearly the same pre-event drift (`4.29e-6` rad). Later event differences accumulate to 0.203 rad maximum yaw and 0.0277 rad/s speed differences over 250 s. With continuous heading interpolation, Python differs from the published trajectory by up to 0.452 rad yaw and 0.101 rad/s yaw speed. Replaying MATLAB's logged six-component force through the Python dynamics reduces these maximum differences to 0.00180 rad and 0.0000935 rad/s. The source force laws and balances are paired, but published-case trajectory parity with native excitation remains **unestablished** because the threshold amplifies small integration differences. | | OSWEC irregular passive yaw, continuous-heading control | Same pinned MATLAB input and phase realization, with only `body(1:2).yaw.threshold` changed from 1° to 0° in a temporary application copy | MATLAB and Python then agree over all 25,001 time samples: maximum differences are `1.1e-7` N or N m across the six excitation components on the MATLAB yaw path, `3.16e-5` rad flap/PTO angle, `7.30e-6` rad/s yaw/PTO speed, `0.875` N m PTO torque, and `0.0123` W source-signed power. Wave elevation differs by less than `1e-13` m and the fixed base remains stationary. This control validates the continuous-heading Python dynamics; it does not erase the published 1° threshold difference. | | OSWEC fixed nonhydrodynamic base | Pinned MATLAB Applications `Nonhydro_Body` case and its BEMIO-generated OSWEC HDF5 | The regular-wave solver reports the stationary base and nonlinear flap motion about the PTO hinge. Against a fresh 400 s MATLAB R2025b run (4,001 samples), maximum flap position differences are 20.8 mm surge, 7.9 mm heave, and 0.00444 rad pitch; velocity differences are 15.6 mm/s surge, 6.8 mm/s heave, and 0.00328 rad/s pitch. All six excitation-force components differ by less than 1 N, the base position and velocity agree exactly, and zero PTO torque differs only by MATLAB numerical noise below `3.3e-7` N m. The base's ground-constraint reaction forces are not calculated. |