Audit held-yaw step instability and gate continuous rebaseline

.github/workflows/oswec-passive-yaw-continuous-step.yml

@@ -0,0 +1,77 @@
+name: OSWEC continuous-heading source step gate
+
+on:
+  pull_request:
+    paths:
+      - '.github/workflows/oswec-passive-yaw-continuous-step.yml'
+      - 'tests/matlab/oswec_passive_yaw_continuous_step_refinement.m'
+      - 'tools/oswec_continuous_yaw_step_audit.py'
+      - 'tests/test_oswec_continuous_yaw_refinement.py'
+  workflow_dispatch:
+
+jobs:
+  refine:
+    strategy:
+      fail-fast: false
+      matrix:
+        phase_seed: [1, 2, 3]
+    runs-on: ubuntu-latest
+    timeout-minutes: 120
+    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
+      - uses: actions/checkout@v7
+        with:
+          repository: cmccomb/WEC-Sim_Applications
+          ref: 34ca45db2b82bb1d2f344175b2b3c51be72266db
+          path: applications
+          filter: blob:none
+          lfs: true
+          sparse-checkout: |
+            _Common_Input_Files/OSWEC
+            Passive_Yaw/PassiveYawRegression
+      - 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'); oswec_passive_yaw_continuous_step_refinement(${{ matrix.phase_seed }})
+      - uses: actions/setup-python@v7
+        with:
+          python-version: '3.12'
+      - run: python -m pip install -e . pytest
+      - id: source_gate
+        continue-on-error: true
+        run: >-
+          python tools/oswec_continuous_yaw_step_audit.py
+          --source-dir matlab-oswec-passive-yaw-continuous-refinement
+          --seed ${{ matrix.phase_seed }}
+          --output matlab-oswec-passive-yaw-continuous-refinement/seed-${{ matrix.phase_seed }}-summary.json
+          --enforce
+      - id: physical_gate
+        continue-on-error: true
+        run: python -m pytest -q -s tests/test_oswec_continuous_yaw_refinement.py
+        env:
+          WEC_SIM_APPLICATIONS_DIR: applications
+          WEC_SIM_CONTINUOUS_REFINEMENT_DIR: matlab-oswec-passive-yaw-continuous-refinement
+          WEC_SIM_PHASE_SEED: ${{ matrix.phase_seed }}
+      - name: Require both source stability and physical parity
+        if: always()
+        env:
+          SOURCE_GATE: ${{ steps.source_gate.outcome }}
+          PHYSICAL_GATE: ${{ steps.physical_gate.outcome }}
+        run: test "$SOURCE_GATE" = success && test "$PHYSICAL_GATE" = success
+      - uses: actions/upload-artifact@v4
+        if: always()
+        with:
+          name: matlab-oswec-continuous-yaw-step-seed-${{ matrix.phase_seed }}
+          path: |
+            matlab-oswec-passive-yaw-continuous-refinement/
+            applications/_Common_Input_Files/OSWEC/hydroData/oswec.h5

.github/workflows/oswec-passive-yaw-step-audit.yml

@@ -5,14 +5,16 @@ on:
     paths:
       - '.github/workflows/oswec-passive-yaw-step-audit.yml'
       - 'tests/matlab/oswec_passive_yaw_step_refinement.m'
+      - 'tools/oswec_held_yaw_step_audit.py'
+      - 'tools/oswec_held_yaw_physical_run.py'
   workflow_dispatch:
 
 jobs:
   refine:
     strategy:
       fail-fast: false
       matrix:
-        phase_seed: [1, 3]
+        phase_seed: [1, 2, 3]
     runs-on: ubuntu-latest
     timeout-minutes: 90
     steps:
@@ -41,6 +43,23 @@ jobs:
       - uses: matlab-actions/run-command@v3
         with:
           command: addpath('tests/matlab'); oswec_passive_yaw_step_refinement(${{ matrix.phase_seed }})
+      - uses: actions/setup-python@v7
+        with:
+          python-version: '3.12'
+      - run: python -m pip install -e .
+      - run: >-
+          python tools/oswec_held_yaw_physical_run.py
+          --hydro applications/_Common_Input_Files/OSWEC/hydroData/oswec.h5
+          --seed ${{ matrix.phase_seed }}
+          --dt 0.0025
+          --output matlab-oswec-passive-yaw-refinement/seed-${{ matrix.phase_seed }}-physical-dt-0.0025.npz
+      - run: >-
+          python tools/oswec_held_yaw_step_audit.py
+          --source-dir matlab-oswec-passive-yaw-refinement
+          --hydro applications/_Common_Input_Files/OSWEC/hydroData/oswec.h5
+          --seed ${{ matrix.phase_seed }}
+          --physical matlab-oswec-passive-yaw-refinement/seed-${{ matrix.phase_seed }}-physical-dt-0.0025.npz
+          --output matlab-oswec-passive-yaw-refinement/seed-${{ matrix.phase_seed }}-summary.json
       - uses: actions/upload-artifact@v4
         if: always()
         with:

PARITY.md

@@ -461,6 +461,46 @@ stable control here: across seeds 1 and 3, halving its step from `0.01` to
 `0.005` s changes PTO work by at most `0.0033%` and yaw by at most
 `4.5e-5` rad, consistent with its separately paired 250 s MATLAB cases.
 
+The proposed [Applications correction](https://github.com/WEC-Sim/WEC-Sim_Applications/pull/120)
+sets both yaw thresholds to zero and rebaselines its regression from the
+archived continuous-heading source. The new full-duration gate pins that
+candidate application commit, checks MATLAB 0.01/0.005 s self-change for
+each seed, then independently advances the public physical Python path at
+0.005 s against the refined source. Its predeclared yaw, speed, PTO torque,
+wave, and work limits remain unpromoted until all three seeds pass in CI.
+
+The remaining seed-2 source refinement uses the same pinned R2025b model,
+HDF5, and phase seed with only `simu.dt` changed from `0.01` to `0.005` to
+`0.0025` s. The local published-step rerun reproduces the archived seed-2
+trajectory within `4.1e-15` rad yaw, `8.7e-16` rad/s speed, and
+`2.4e-10` N m torque. Wave elevation is identical at matched samples across
+all three steps. MATLAB self-selected heading updates number `52`, `51`,
+and `39`; PTO work is `1990.936`, `2045.774`, and `1521.927` J. The
+half/quarter-step source self-change reaches `0.125657` rad yaw,
+`0.024061` rad/s speed, and `2887.33` N m torque, while work changes
+`-25.606%`. An independently advanced public physical Python run at
+`0.0025` s selects `57` events and absorbs `2681.039` J, differing from
+that source step by `0.128375` rad yaw, `0.030576` rad/s speed,
+`3669.08` N m torque, and `76.161%` work. Its own `0.00125` s refinement
+selects `53` events and absorbs `1746.504` J, a `-34.857%` work change
+and `0.322703` rad maximum yaw change, with identical matched wave samples.
+The independently advanced Python run at `0.0025` s also completes the other
+two seas. The full-duration quarter-step comparisons are:
+
+| Phase seed | MATLAB / Python heading updates | Maximum yaw gap (rad) | Maximum speed gap (rad/s) | Maximum PTO torque gap (N m) | Python versus MATLAB PTO work |
+| --- | ---: | ---: | ---: | ---: | ---: |
+| 1 | 121 / 134 | 0.192312 | 0.025656 | 3078.70 | +31.780% |
+| 2 | 39 / 57 | 0.128375 | 0.030576 | 3669.08 | +76.161% |
+| 3 | 125 / 104 | 0.097148 | 0.018322 | 2198.59 | -34.823% |
+
+The [three-seed step runner](tests/matlab/oswec_passive_yaw_step_refinement.m),
+[independent physical runner](tools/oswec_held_yaw_physical_run.py), and
+[audit](tools/oswec_held_yaw_step_audit.py) retain source and physical
+trajectories separately and report event times, force, motion, wave, and
+work with limits recorded before the quarter-step result. None of the
+three one-degree held-heading seeds has a source-stable full-duration
+physical baseline at the tested steps; the published case remains qualified.
+
 The pinned OWC Simulink model has a force-path difference:
 its flexible `GBM` state-space block integrates inverse effective mass times
 the hydrodynamic force, while the piston force joins a separate reported

tests/matlab/oswec_passive_yaw_continuous_step_refinement.m

@@ -0,0 +1,92 @@
+function oswec_passive_yaw_continuous_step_refinement(phaseSeed, steps)
+% Check the corrected application's continuous-heading regression at two
+% full-duration steps for each of the same three 250 s seas.
+if nargin < 2
+    steps = [0.01, 0.005];
+end
+assert(ismember(phaseSeed, [1, 2, 3]) && ...
+    all(ismember(steps, [0.01, 0.005])) && numel(unique(steps)) == numel(steps), ...
+    'Unexpected phase seed or time steps');
+originalFolder = pwd;
+repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath'))));
+addpath(genpath(fullfile(repoRoot, 'matlab-ref', 'source')));
+hydroDir = fullfile(repoRoot, 'applications', '_Common_Input_Files', ...
+    'OSWEC', 'hydroData');
+cd(hydroDir);
+if ~isfile('oswec.h5')
+    bemio;
+end
+sourceDir = fullfile(repoRoot, 'applications', 'Passive_Yaw', ...
+    'PassiveYawRegression');
+original = fileread(fullfile(sourceDir, 'wecSimInputFile.m'));
+assert(contains(original, 'simu.dt = 0.01;') && ...
+    contains(original, 'simu.endTime = 250;') && ...
+    contains(original, 'waves.phaseSeed = 1;') && ...
+    contains(original, 'body(1).yaw.threshold = 0;') && ...
+    contains(original, 'body(2).yaw.threshold = 0;'), ...
+    'Pinned corrected passive-yaw input changed');
+caseDir = fullfile(repoRoot, 'applications', 'Passive_Yaw', ...
+    sprintf('PassiveYawContinuousStepSeed%d', phaseSeed));
+assert(~isfolder(caseDir), 'Disposable continuous-yaw case already exists');
+mkdir(caseDir);
+cleanup = onCleanup(@() cleanupCase(caseDir, originalFolder));
+copyfile(fullfile(sourceDir, 'OSWEC.slx'), fullfile(caseDir, 'OSWEC.slx'));
+inputFile = fullfile(caseDir, 'wecSimInputFile.m');
+outDir = fullfile(repoRoot, 'matlab-oswec-passive-yaw-continuous-refinement');
+mkdir(outDir);
+for dt = steps
+    contents = strrep(original, 'simu.dt = 0.01;', ...
+        sprintf('simu.dt = %.6g;', dt));
+    contents = strrep(contents, 'waves.phaseSeed = 1;', ...
+        sprintf('waves.phaseSeed = %d;', phaseSeed));
+    fid = fopen(inputFile, 'w');
+    assert(fid ~= -1, 'Could not write disposable continuous-yaw input');
+    fprintf(fid, '%s', contents);
+    fclose(fid);
+    cd(caseDir);
+    wecSim;
+    assert(simu.dt == dt && simu.endTime == 250 && ...
+        simu.rampTime == 100 && strcmp(simu.solver, 'ode4') && ...
+        strcmp(waves.type, 'irregular') && waves.phaseSeed == phaseSeed && ...
+        waves.height == 2.5 && waves.period == 8 && ...
+        body(1).yaw.threshold == 0 && body(2).yaw.threshold == 0, ...
+        'Continuous-yaw refinement changed another setting');
+    time = output.bodies(1).time;
+    yaw = output.bodies(1).position(:, 6);
+    yawRate = output.bodies(1).velocity(:, 6);
+    ptoPower = output.ptos(1).powerInternalMechanics(:, 5);
+    ptoForce = output.ptos(1).forceInternalMechanics(:, 5);
+    wave = waves.waveAmpTime;
+    work = -trapz(time, ptoPower);
+    assert(numel(time) == round(250 / dt) + 1 && ...
+        isequal(size(wave), [numel(time), 2]) && ...
+        all(isfinite([time, yaw, yawRate, ptoPower, ptoForce, wave]), 'all') && ...
+        isfinite(work) && work > 0, ...
+        'Continuous-yaw refinement has missing or nonfinite output');
+    save(fullfile(outDir, sprintf('seed-%d-dt-%g.mat', phaseSeed, dt)), ...
+        'phaseSeed', 'dt', 'time', 'yaw', 'yawRate', ...
+        'ptoPower', 'ptoForce', 'wave', 'work', '-v7');
+end
+if numel(steps) == 2
+    coarse = load(fullfile(outDir, sprintf('seed-%d-dt-0.01.mat', phaseSeed)));
+    fine = load(fullfile(outDir, sprintf('seed-%d-dt-0.005.mat', phaseSeed)));
+    assert(max(abs(coarse.time - fine.time(1:2:end))) < 1e-10 && ...
+        max(abs(coarse.wave(:, 2) - fine.wave(1:2:end, 2))) < 1e-11, ...
+        'Continuous-yaw runs did not see the same wave');
+    fprintf('Seed %d continuous-yaw: yaw %.6g rad, speed %.6g rad/s, torque %.6g N m, work %.6g%%\n', ...
+        phaseSeed, max(abs(coarse.yaw - fine.yaw(1:2:end))), ...
+        max(abs(coarse.yawRate - fine.yawRate(1:2:end))), ...
+        max(abs(coarse.ptoForce - fine.ptoForce(1:2:end))), ...
+        100 * (fine.work / coarse.work - 1));
+end
+end
+
+function cleanupCase(caseDir, originalFolder)
+if bdIsLoaded('OSWEC')
+    close_system('OSWEC', 0);
+end
+cd(originalFolder);
+if isfolder(caseDir)
+    rmdir(caseDir, 's');
+end
+end

tests/matlab/oswec_passive_yaw_step_refinement.m

@@ -1,6 +1,7 @@
 function oswec_passive_yaw_step_refinement(phaseSeed)
 % Audit the published held-heading model's sensitivity to fixed time step.
-assert(ismember(phaseSeed, [1, 3]), 'Expected published or third audit seed');
+assert(ismember(phaseSeed, [1, 2, 3]), 'Expected published phase seed');
+originalFolder = pwd;
 repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath'))));
 addpath(genpath(fullfile(repoRoot, 'matlab-ref', 'source')));
 hydroDir = fullfile(repoRoot, 'applications', '_Common_Input_Files', ...
@@ -19,12 +20,14 @@ function oswec_passive_yaw_step_refinement(phaseSeed)
     contains(original, 'body(1).yaw.threshold = 1;') && ...
     contains(original, 'body(2).yaw.threshold = 1;'), ...
     'Pinned passive-yaw case input changed');
+restore = onCleanup(@() restoreInput(inputFile, original, originalFolder));
 outDir = fullfile(repoRoot, 'matlab-oswec-passive-yaw-refinement');
 mkdir(outDir);
 steps = [0.01, 0.005, 0.0025];
 for dt = steps
     runOneStep(caseDir, inputFile, original, outDir, phaseSeed, dt);
 end
+
 for pair = 1:(numel(steps) - 1)
     coarse = load(fullfile(outDir, ...
         sprintf('seed-%d-dt-%g.mat', phaseSeed, steps(pair))));
@@ -75,3 +78,14 @@ function runOneStep(caseDir, inputFile, original, outDir, phaseSeed, dt)
 save(fullfile(outDir, filename), 'phaseSeed', 'dt', 'time', 'yaw', ...
     'yawRate', 'ptoPower', 'ptoForce', 'wave', 'work', '-v7');
 end
+
+function restoreInput(inputFile, original, originalFolder)
+if bdIsLoaded('OSWEC')
+    close_system('OSWEC', 0);
+end
+fid = fopen(inputFile, 'w');
+assert(fid ~= -1, 'Could not restore passive-yaw input');
+fprintf(fid, '%s', original);
+fclose(fid);
+cd(originalFolder);
+end

tests/test_oswec_continuous_yaw_refinement.py

@@ -0,0 +1,58 @@
+"""Full-duration physical parity against a source-stable continuous-yaw step."""
+
+import os
+from pathlib import Path
+
+import numpy as np
+import pytest
+from scipy.io import loadmat
+
+from wecsim import PMWave, WEC, WorldPoint
+
+
+APPLICATIONS = os.environ.get("WEC_SIM_APPLICATIONS_DIR")
+REFERENCE = os.environ.get("WEC_SIM_CONTINUOUS_REFINEMENT_DIR")
+PHASE_SEED = os.environ.get("WEC_SIM_PHASE_SEED")
+pytestmark = pytest.mark.skipif(
+    not (APPLICATIONS and REFERENCE and PHASE_SEED),
+    reason="pinned continuous-yaw source refinement not provided",
+)
+
+
+def test_continuous_heading_physical_full_duration():
+    seed = int(PHASE_SEED)
+    assert seed in (1, 2, 3)
+    hydro = (Path(APPLICATIONS)
+             / "_Common_Input_Files/OSWEC/hydroData/oswec.h5")
+    source = loadmat(
+        Path(REFERENCE) / f"seed-{seed}-dt-0.005.mat", squeeze_me=True,
+    )
+    wec = WEC("OSWEC continuous-heading physical parity")
+    flap = wec.body(
+        "flap", hydro, mass=12700, inertia=(1.85e6,) * 3,
+        passive_yaw=True, passive_yaw_threshold=0,
+    )
+    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=seed, phase_generator="matlab"),
+        dt=0.005, end_time=250, ramp_time=100, radiation_memory=40,
+    )
+    time = np.asarray(source["time"]).ravel()
+    source_yaw = np.asarray(source["yaw"]).ravel()
+    source_speed = np.asarray(source["yawRate"]).ravel()
+    source_torque = np.asarray(source["ptoForce"]).ravel()
+    source_power = np.asarray(source["ptoPower"]).ravel()
+    source_wave = np.asarray(source["wave"])[:, 1]
+    np.testing.assert_allclose(result.time, time, rtol=0, atol=1e-10)
+    assert np.max(np.abs(result.wave_elevation - source_wave)) < 1e-11
+    assert np.max(np.abs(result.coordinates["yaw"].position - source_yaw)) < 1e-4
+    assert np.max(np.abs(result.coordinates["yaw"].velocity - source_speed)) < 2e-5
+    assert np.max(np.abs(result.ptos["hinge"].force - source_torque)) < 2.4
+    source_work = -np.trapezoid(source_power, time)
+    physical_work = np.trapezoid(result.ptos["hinge"].absorbed_power, time)
+    assert source_work > 0
+    assert abs(physical_work / source_work - 1) < 0.0002
+    print(f"seed {seed} continuous 250 s: source work {source_work:.6f} J, "
+          f"physical work {physical_work:.6f} J")