From f525e98aac3f2e7000668cffd49cfde6ade5112b Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 22 Oct 2026 18:15:09 -0400 Subject: [PATCH 1/6] Audit all OSWEC held-yaw seeds under source step refinement --- .../oswec-passive-yaw-step-audit.yml | 13 +- .../oswec_passive_yaw_step_refinement.m | 13 +- tools/oswec_held_yaw_step_audit.py | 166 ++++++++++++++++++ 3 files changed, 190 insertions(+), 2 deletions(-) create mode 100644 tools/oswec_held_yaw_step_audit.py diff --git a/.github/workflows/oswec-passive-yaw-step-audit.yml b/.github/workflows/oswec-passive-yaw-step-audit.yml index 85ce0e43..74707f4c 100644 --- a/.github/workflows/oswec-passive-yaw-step-audit.yml +++ b/.github/workflows/oswec-passive-yaw-step-audit.yml @@ -5,6 +5,7 @@ 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' workflow_dispatch: jobs: @@ -12,7 +13,7 @@ jobs: strategy: fail-fast: false matrix: - phase_seed: [1, 3] + phase_seed: [1, 2, 3] runs-on: ubuntu-latest timeout-minutes: 90 steps: @@ -41,6 +42,16 @@ 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_step_audit.py + --source-dir matlab-oswec-passive-yaw-refinement + --hydro applications/_Common_Input_Files/OSWEC/hydroData/oswec.h5 + --seed ${{ matrix.phase_seed }} + --output matlab-oswec-passive-yaw-refinement/seed-${{ matrix.phase_seed }}-summary.json - uses: actions/upload-artifact@v4 if: always() with: diff --git a/tests/matlab/oswec_passive_yaw_step_refinement.m b/tests/matlab/oswec_passive_yaw_step_refinement.m index cbfc850d..26abbb9b 100644 --- a/tests/matlab/oswec_passive_yaw_step_refinement.m +++ b/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,11 @@ 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) +fid = fopen(inputFile, 'w'); +assert(fid ~= -1, 'Could not restore passive-yaw input'); +fprintf(fid, '%s', original); +fclose(fid); +cd(originalFolder); +end diff --git a/tools/oswec_held_yaw_step_audit.py b/tools/oswec_held_yaw_step_audit.py new file mode 100644 index 00000000..57a4cc16 --- /dev/null +++ b/tools/oswec_held_yaw_step_audit.py @@ -0,0 +1,166 @@ +"""Summarize independent held-yaw step refinements from pinned MATLAB runs. + +The source files come from ``oswec_passive_yaw_step_refinement.m``. Optional +Python input must be a self-advanced physical result; this tool never feeds +source motion, forces, or heading events into that run. +""" + +import argparse +import hashlib +import json +from pathlib import Path + +import h5py +import numpy as np +from scipy.io import loadmat + + +STEPS = (0.01, 0.005, 0.0025) +# The earlier prescribed-heading diagnostic gates yaw, speed, and work at +# 0.003 rad, 0.0002 rad/s, and 0.1%. The PTO is a 120,000 N m s/rad damper, +# so its speed limit implies 24 N m of torque. Require the same event count +# and event times within one published 0.01 s sample before using a source +# step as a physical reference. +LIMITS = { + "yaw_rad": 0.003, + "yaw_speed_rad_s": 0.0002, + "pto_torque_nm": 24, + "pto_work_percent": 0.1, + "heading_event_s": 0.01, + "wave_m": 1e-11, +} + + +def _heading_events(time, yaw, directions, incident_direction=10, threshold=1): + """Replay the published coefficient-selection rule on accepted yaw states.""" + previous = None + headings = np.empty(len(yaw)) + for index, angle in enumerate(yaw): + relative = incident_direction - np.degrees(angle) + distance = (relative - directions + 180) % 360 - 180 + nearest = int(np.argmin(np.abs(distance))) + if abs(distance[nearest]) <= threshold: + previous = relative - distance[nearest] + elif previous is None or abs(relative - previous) > threshold: + previous = relative + headings[index] = previous + changes = np.flatnonzero(np.abs(np.diff(headings)) > 1e-5) + 1 + return time[changes] + + +def _load_source(path, dt): + data = loadmat(path, squeeze_me=True) + run = {key: np.asarray(data[key], dtype=float).ravel() + for key in ("time", "yaw", "yawRate", "ptoForce", "ptoPower")} + run["wave"] = np.asarray(data["wave"], dtype=float)[:, 1] + size = round(250 / dt) + 1 + assert all(len(value) == size for value in run.values()), path + assert np.max(np.abs(run["time"] - np.arange(size) * dt)) < 1e-10, path + assert all(np.isfinite(value).all() for value in run.values()), path + work = -np.trapezoid(run["ptoPower"], run["time"]) + assert work > 0 and abs(work - float(data["work"])) < 1e-8, path + return run, work + + +def _event_comparison(first, second): + same_count = len(first) == len(second) + return { + "same_heading_event_count": same_count, + "max_heading_event_time_difference_s": ( + float(np.max(np.abs(first - second))) if same_count and len(first) + else 0.0 if same_count else None + ), + } + + +def _within_limits(metrics): + event_gap = metrics["max_heading_event_time_difference_s"] + return (metrics["max_wave_difference"] < LIMITS["wave_m"] + and metrics["max_yaw_difference"] < LIMITS["yaw_rad"] + and metrics["max_yawRate_difference"] < LIMITS["yaw_speed_rad_s"] + and metrics["max_ptoForce_difference"] < LIMITS["pto_torque_nm"] + and abs(metrics["pto_work_difference_percent"]) < LIMITS["pto_work_percent"] + and metrics["same_heading_event_count"] + and event_gap is not None + and event_gap <= LIMITS["heading_event_s"]) + + +def audit(source_dir, hydro, seed, physical=None): + with h5py.File(hydro) as file: + directions = np.asarray(file["simulation_parameters/wave_dir"]).ravel() + runs = {} + summary = {"seed": seed, "limits": LIMITS, + "matlab": {}, "source_step_pairs": {}} + for dt in STEPS: + path = source_dir / f"seed-{seed}-dt-{dt:g}.mat" + run, work = _load_source(path, dt) + events = _heading_events(run["time"], run["yaw"], directions) + runs[dt] = (run, work) + summary["matlab"][f"{dt:g}"] = { + "sha256": hashlib.sha256(path.read_bytes()).hexdigest(), + "pto_work_j": work, + "heading_events_s": events.tolist(), + } + for coarse_dt, fine_dt in zip(STEPS, STEPS[1:], strict=False): + coarse, coarse_work = runs[coarse_dt] + fine, fine_work = runs[fine_dt] + pair = {} + for key in ("time", "wave", "yaw", "yawRate", "ptoForce"): + pair[f"max_{key}_difference"] = float(np.max( + np.abs(coarse[key] - fine[key][::2]))) + assert pair["max_time_difference"] < 1e-10 + assert pair["max_wave_difference"] < 1e-11 + pair["pto_work_difference_percent"] = 100 * (fine_work / coarse_work - 1) + pair.update(_event_comparison( + np.asarray(summary["matlab"][f"{coarse_dt:g}"]["heading_events_s"]), + np.asarray(summary["matlab"][f"{fine_dt:g}"]["heading_events_s"]), + )) + pair["within_limits"] = _within_limits(pair) + summary["source_step_pairs"][f"{coarse_dt:g}/{fine_dt:g}"] = pair + if physical is not None: + with np.load(physical) as archive: + result = {key: archive[key] for key in archive.files} + fine, fine_work = runs[STEPS[-1]] + assert np.max(np.abs(result["time"] - fine["time"])) < 1e-10 + assert np.max(np.abs(result["wave"] - fine["wave"])) < 1e-11 + physical_events = _heading_events( + result["time"], result["yaw"], directions, + ) + physical_work = np.trapezoid(result["power"], result["time"]) + comparison = { + "sha256": hashlib.sha256(physical.read_bytes()).hexdigest(), + "pto_work_j": physical_work, + "heading_events_s": physical_events.tolist(), + "max_wave_difference": float(np.max(np.abs(result["wave"] - fine["wave"]))), + "max_yaw_difference": float(np.max(np.abs(result["yaw"] - fine["yaw"]))), + "max_yawRate_difference": float(np.max(np.abs(result["speed"] - fine["yawRate"]))), + "max_ptoForce_difference": float(np.max(np.abs(result["torque"] - fine["ptoForce"]))), + "pto_work_difference_percent": 100 * (physical_work / fine_work - 1), + } + comparison.update(_event_comparison( + np.asarray(summary["matlab"][f"{STEPS[-1]:g}"]["heading_events_s"]), + physical_events, + )) + comparison["within_limits"] = _within_limits(comparison) + summary["physical_python_at_0.0025"] = comparison + return summary + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--source-dir", required=True, type=Path) + parser.add_argument("--hydro", required=True, type=Path) + parser.add_argument("--seed", required=True, type=int, choices=(1, 2, 3)) + parser.add_argument("--physical", type=Path) + parser.add_argument("--output", type=Path) + args = parser.parse_args() + summary = audit(args.source_dir, args.hydro, args.seed, args.physical) + contents = json.dumps(summary, indent=2, sort_keys=True) + "\n" + if args.output: + args.output.write_text(contents) + else: + print(contents, end="") + + +if __name__ == "__main__": + main() From 4c586c3042972c84dd9c77f2435ce107c5155b47 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 22 Oct 2026 18:30:11 -0400 Subject: [PATCH 2/6] Record three-seed held-yaw instability and predeclare continuous source gate --- .../oswec-passive-yaw-continuous-step.yml | 61 ++++++++++++ .../oswec-passive-yaw-step-audit.yml | 8 ++ PARITY.md | 32 ++++++ ...c_passive_yaw_continuous_step_refinement.m | 97 +++++++++++++++++++ .../oswec_passive_yaw_step_refinement.m | 3 + tools/oswec_continuous_yaw_step_audit.py | 86 ++++++++++++++++ tools/oswec_held_yaw_physical_run.py | 59 +++++++++++ 7 files changed, 346 insertions(+) create mode 100644 .github/workflows/oswec-passive-yaw-continuous-step.yml create mode 100644 tests/matlab/oswec_passive_yaw_continuous_step_refinement.m create mode 100644 tools/oswec_continuous_yaw_step_audit.py create mode 100644 tools/oswec_held_yaw_physical_run.py diff --git a/.github/workflows/oswec-passive-yaw-continuous-step.yml b/.github/workflows/oswec-passive-yaw-continuous-step.yml new file mode 100644 index 00000000..3146a496 --- /dev/null +++ b/.github/workflows/oswec-passive-yaw-continuous-step.yml @@ -0,0 +1,61 @@ +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' + 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: WEC-Sim/WEC-Sim_Applications + ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e + 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 . + - 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 + - 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 diff --git a/.github/workflows/oswec-passive-yaw-step-audit.yml b/.github/workflows/oswec-passive-yaw-step-audit.yml index 74707f4c..37e2d288 100644 --- a/.github/workflows/oswec-passive-yaw-step-audit.yml +++ b/.github/workflows/oswec-passive-yaw-step-audit.yml @@ -6,6 +6,7 @@ on: - '.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: @@ -46,11 +47,18 @@ jobs: 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() diff --git a/PARITY.md b/PARITY.md index d2c0b387..2ac524e7 100644 --- a/PARITY.md +++ b/PARITY.md @@ -461,6 +461,38 @@ 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 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 diff --git a/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m b/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m new file mode 100644 index 00000000..bc78eb06 --- /dev/null +++ b/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m @@ -0,0 +1,97 @@ +function oswec_passive_yaw_continuous_step_refinement(phaseSeed, steps) +% Check whether disabling the published 1-degree coefficient hold stabilizes +% the same three 250 s seas. The pinned application input is never edited. +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 = 1;') && ... + contains(original, 'body(2).yaw.threshold = 1;'), ... + 'Pinned 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)); + for iBody = 1:2 + contents = strrep(contents, ... + sprintf('body(%d).yaw.threshold = 1;', iBody), ... + sprintf('body(%d).yaw.threshold = 0;', iBody)); + end + 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 diff --git a/tests/matlab/oswec_passive_yaw_step_refinement.m b/tests/matlab/oswec_passive_yaw_step_refinement.m index 26abbb9b..a9637021 100644 --- a/tests/matlab/oswec_passive_yaw_step_refinement.m +++ b/tests/matlab/oswec_passive_yaw_step_refinement.m @@ -80,6 +80,9 @@ function runOneStep(caseDir, inputFile, original, outDir, phaseSeed, dt) 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); diff --git a/tools/oswec_continuous_yaw_step_audit.py b/tools/oswec_continuous_yaw_step_audit.py new file mode 100644 index 00000000..ebfc7971 --- /dev/null +++ b/tools/oswec_continuous_yaw_step_audit.py @@ -0,0 +1,86 @@ +"""Gate the corrected continuous-heading source on two full-duration steps.""" + +import argparse +import hashlib +import json +from pathlib import Path + +import numpy as np +from scipy.io import loadmat + + +# Chosen before evaluating the half-step source run. These limits exceed the +# prior three-seed source/Python 0.01 s gaps but still require small full-run +# changes after halving the source step. +LIMITS = { + "wave_m": 1e-11, + "yaw_rad": 1e-4, + "yaw_speed_rad_s": 2e-5, + "pto_torque_nm": 2.4, + "pto_work_percent": 0.02, +} + + +def _load(path, dt): + data = loadmat(path, squeeze_me=True) + arrays = {key: np.asarray(data[key], dtype=float).ravel() + for key in ("time", "yaw", "yawRate", "ptoForce", "ptoPower")} + arrays["wave"] = np.asarray(data["wave"], dtype=float)[:, 1] + length = round(250 / dt) + 1 + assert all(len(value) == length for value in arrays.values()), path + assert np.max(np.abs(arrays["time"] - np.arange(length) * dt)) < 1e-10, path + assert all(np.isfinite(value).all() for value in arrays.values()), path + work = -np.trapezoid(arrays["ptoPower"], arrays["time"]) + assert work > 0 and abs(work - float(data["work"])) < 1e-8, path + return arrays, work + + +def audit(source_dir, seed): + files = {dt: source_dir / f"seed-{seed}-dt-{dt:g}.mat" + for dt in (0.01, 0.005)} + coarse, coarse_work = _load(files[0.01], 0.01) + fine, fine_work = _load(files[0.005], 0.005) + differences = { + key: float(np.max(np.abs(coarse[key] - fine[key][::2]))) + for key in ("time", "wave", "yaw", "yawRate", "ptoForce") + } + assert differences["time"] < 1e-10 + work_change = 100 * (fine_work / coarse_work - 1) + within_limits = ( + differences["wave"] < LIMITS["wave_m"] + and differences["yaw"] < LIMITS["yaw_rad"] + and differences["yawRate"] < LIMITS["yaw_speed_rad_s"] + and differences["ptoForce"] < LIMITS["pto_torque_nm"] + and abs(work_change) < LIMITS["pto_work_percent"] + ) + return { + "seed": seed, + "limits": LIMITS, + "source_sha256": {f"{dt:g}": hashlib.sha256(path.read_bytes()).hexdigest() + for dt, path in files.items()}, + "pto_work_j": {"0.01": coarse_work, "0.005": fine_work}, + "max_common_grid_difference": differences, + "pto_work_change_percent": work_change, + "within_limits": within_limits, + } + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--source-dir", required=True, type=Path) + parser.add_argument("--seed", required=True, type=int, choices=(1, 2, 3)) + parser.add_argument("--output", type=Path) + parser.add_argument("--enforce", action="store_true") + args = parser.parse_args() + report = audit(args.source_dir, args.seed) + contents = json.dumps(report, indent=2, sort_keys=True) + "\n" + if args.output: + args.output.write_text(contents) + else: + print(contents, end="") + if args.enforce and not report["within_limits"]: + raise SystemExit("Continuous-heading source did not meet step limits") + + +if __name__ == "__main__": + main() diff --git a/tools/oswec_held_yaw_physical_run.py b/tools/oswec_held_yaw_physical_run.py new file mode 100644 index 00000000..85af6aae --- /dev/null +++ b/tools/oswec_held_yaw_physical_run.py @@ -0,0 +1,59 @@ +"""Advance a published OSWEC held-yaw sea through the public physical solver. + +Only the pinned HDF5 and MATLAB-compatible generated wave phases enter the +run. Source trajectories, forces, and heading events are never inputs. +""" + +import argparse +from pathlib import Path + +import numpy as np + +from wecsim import PMWave, WEC, WorldPoint + + +def run(hydro, seed, dt): + wec = WEC("OSWEC published held yaw, physical implicit path") + flap = wec.body( + "flap", hydro, mass=12700, inertia=(1.85e6,) * 3, + passive_yaw=True, passive_yaw_threshold=1, + ) + 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=dt, end_time=250, ramp_time=100, radiation_memory=40, + ) + data = { + "time": result.time, + "wave": result.wave_elevation, + "yaw": result.coordinates["yaw"].position, + "speed": result.coordinates["yaw"].velocity, + "torque": result.ptos["hinge"].force, + "power": result.ptos["hinge"].absorbed_power, + } + expected_size = round(250 / dt) + 1 + assert all(len(value) == expected_size for value in data.values()) + assert all(np.isfinite(value).all() for value in data.values()) + np.testing.assert_allclose(data["torque"], -120000 * data["speed"], + rtol=0, atol=1e-7) + assert np.trapezoid(data["power"], data["time"]) > 0 + return data + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--hydro", required=True, type=Path) + parser.add_argument("--seed", required=True, type=int, choices=(1, 2, 3)) + parser.add_argument("--dt", required=True, type=float) + parser.add_argument("--output", required=True, type=Path) + args = parser.parse_args() + data = run(args.hydro, args.seed, args.dt) + np.savez_compressed(args.output, **data) + print(f"Seed {args.seed}, dt={args.dt:g}: " + f"physical PTO work {np.trapezoid(data['power'], data['time']):.6f} J") + + +if __name__ == "__main__": + main() From 4872199e7bf6d9adfb2a8e036895c3c1633dce33 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 22 Oct 2026 18:34:22 -0400 Subject: [PATCH 3/6] Predeclare full-duration continuous-heading physical parity gate --- .../oswec-passive-yaw-continuous-step.yml | 8 ++- tests/test_oswec_continuous_yaw_refinement.py | 58 +++++++++++++++++++ 2 files changed, 65 insertions(+), 1 deletion(-) create mode 100644 tests/test_oswec_continuous_yaw_refinement.py diff --git a/.github/workflows/oswec-passive-yaw-continuous-step.yml b/.github/workflows/oswec-passive-yaw-continuous-step.yml index 3146a496..27a3d347 100644 --- a/.github/workflows/oswec-passive-yaw-continuous-step.yml +++ b/.github/workflows/oswec-passive-yaw-continuous-step.yml @@ -6,6 +6,7 @@ on: - '.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: @@ -45,13 +46,18 @@ jobs: - uses: actions/setup-python@v7 with: python-version: '3.12' - - run: python -m pip install -e . + - run: python -m pip install -e . pytest - 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 + - 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 }} - uses: actions/upload-artifact@v4 if: always() with: diff --git a/tests/test_oswec_continuous_yaw_refinement.py b/tests/test_oswec_continuous_yaw_refinement.py new file mode 100644 index 00000000..b17cd75d --- /dev/null +++ b/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") From 340fcf8cb2c4d1aa1988a848f4ec802d5666a8b9 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 22 Oct 2026 18:42:33 -0400 Subject: [PATCH 4/6] Collect both continuous source and physical gate outcomes --- .../oswec-passive-yaw-continuous-step.yml | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/.github/workflows/oswec-passive-yaw-continuous-step.yml b/.github/workflows/oswec-passive-yaw-continuous-step.yml index 27a3d347..f35b7056 100644 --- a/.github/workflows/oswec-passive-yaw-continuous-step.yml +++ b/.github/workflows/oswec-passive-yaw-continuous-step.yml @@ -47,17 +47,27 @@ jobs: with: python-version: '3.12' - run: python -m pip install -e . pytest - - run: >- + - 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 - - run: python -m pytest -q -s tests/test_oswec_continuous_yaw_refinement.py + - 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: From 61e988c853b6a6a810efffb0b6bc65c42ba29b3e Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 22 Oct 2026 18:48:15 -0400 Subject: [PATCH 5/6] Pin corrected application input for continuous physical gate --- .../oswec-passive-yaw-continuous-step.yml | 4 ++-- PARITY.md | 8 ++++++++ ...oswec_passive_yaw_continuous_step_refinement.m | 15 +++++---------- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/.github/workflows/oswec-passive-yaw-continuous-step.yml b/.github/workflows/oswec-passive-yaw-continuous-step.yml index f35b7056..8af1c4b9 100644 --- a/.github/workflows/oswec-passive-yaw-continuous-step.yml +++ b/.github/workflows/oswec-passive-yaw-continuous-step.yml @@ -28,8 +28,8 @@ jobs: sparse-checkout: source - uses: actions/checkout@v7 with: - repository: WEC-Sim/WEC-Sim_Applications - ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e + repository: cmccomb/WEC-Sim_Applications + ref: 34ca45db2b82bb1d2f344175b2b3c51be72266db path: applications filter: blob:none lfs: true diff --git a/PARITY.md b/PARITY.md index 2ac524e7..0761b8cb 100644 --- a/PARITY.md +++ b/PARITY.md @@ -461,6 +461,14 @@ 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 diff --git a/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m b/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m index bc78eb06..414e2a2f 100644 --- a/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m +++ b/tests/matlab/oswec_passive_yaw_continuous_step_refinement.m @@ -1,6 +1,6 @@ function oswec_passive_yaw_continuous_step_refinement(phaseSeed, steps) -% Check whether disabling the published 1-degree coefficient hold stabilizes -% the same three 250 s seas. The pinned application input is never edited. +% 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 @@ -22,9 +22,9 @@ function oswec_passive_yaw_continuous_step_refinement(phaseSeed, steps) assert(contains(original, 'simu.dt = 0.01;') && ... contains(original, 'simu.endTime = 250;') && ... contains(original, 'waves.phaseSeed = 1;') && ... - contains(original, 'body(1).yaw.threshold = 1;') && ... - contains(original, 'body(2).yaw.threshold = 1;'), ... - 'Pinned passive-yaw input changed'); + 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'); @@ -39,11 +39,6 @@ function oswec_passive_yaw_continuous_step_refinement(phaseSeed, steps) sprintf('simu.dt = %.6g;', dt)); contents = strrep(contents, 'waves.phaseSeed = 1;', ... sprintf('waves.phaseSeed = %d;', phaseSeed)); - for iBody = 1:2 - contents = strrep(contents, ... - sprintf('body(%d).yaw.threshold = 1;', iBody), ... - sprintf('body(%d).yaw.threshold = 0;', iBody)); - end fid = fopen(inputFile, 'w'); assert(fid ~= -1, 'Could not write disposable continuous-yaw input'); fprintf(fid, '%s', contents); From fd5a8f320b732cdd0ec2815a898b4bf6f25b165a Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 22 Oct 2026 19:39:09 -0400 Subject: [PATCH 6/6] Serialize continuous source gate result as native boolean --- tools/oswec_continuous_yaw_step_audit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/oswec_continuous_yaw_step_audit.py b/tools/oswec_continuous_yaw_step_audit.py index ebfc7971..22e1c5bd 100644 --- a/tools/oswec_continuous_yaw_step_audit.py +++ b/tools/oswec_continuous_yaw_step_audit.py @@ -46,7 +46,7 @@ def audit(source_dir, seed): } assert differences["time"] < 1e-10 work_change = 100 * (fine_work / coarse_work - 1) - within_limits = ( + within_limits = bool( differences["wave"] < LIMITS["wave_m"] and differences["yaw"] < LIMITS["yaw_rad"] and differences["yawRate"] < LIMITS["yaw_speed_rad_s"]