Reproduce corrected OWC piston path for full parity gate

.github/workflows/owc-orifice-corrected-source.yml

@@ -0,0 +1,80 @@
+name: OWC corrected piston source
+
+on:
+  pull_request:
+    paths:
+      - '.github/workflows/owc-orifice-corrected-source.yml'
+      - 'tests/matlab/owc_orifice_corrected_reference.m'
+      - 'tests/matlab/owc_orifice_correct_piston_path.m'
+      - 'tests/test_owc_orifice_corrected_full.py'
+      - 'tools/pack_owc_orifice_corrected.py'
+      - 'tests/fixtures/matlab/owc_orifice_corrected_piston/full_source.npz'
+  workflow_dispatch:
+
+jobs:
+  source:
+    runs-on: ubuntu-latest
+    timeout-minutes: 180
+    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:
+          # Exact model proposed in WEC-Sim/WEC-Sim_Applications#119.
+          repository: cmccomb/WEC-Sim_Applications
+          ref: 3c0c05fefce55125d8ea084deedbd2defb103811
+          path: applications
+          filter: blob:none
+          lfs: true
+          sparse-checkout: OWC/OrificeModel
+      - name: Install pinned cleaned hydrodynamics
+        run: cp tests/fixtures/matlab/owc_orifice.h5 applications/OWC/OrificeModel/hydroData/test17a_clean.h5
+      - uses: matlab-actions/setup-matlab@v3
+        with:
+          release: R2025b
+          products: Simulink Simscape Simscape_Multibody Signal_Processing_Toolbox Control_System_Toolbox
+      - uses: matlab-actions/run-command@v3
+        with:
+          command: run('tests/matlab/owc_orifice_corrected_reference.m')
+      - uses: actions/setup-python@v7
+        with:
+          python-version: '3.12'
+          cache: pip
+      - run: python -m pip install -e . pytest
+      - name: Reproduce archived corrected source from published input
+        run: |
+          python tools/pack_owc_orifice_corrected.py \
+            owc-orifice-corrected-output \
+            owc-orifice-corrected-output/fresh_source.npz
+          python - <<'PY'
+          from pathlib import Path
+          import numpy as np
+
+          archived = Path('tests/fixtures/matlab/owc_orifice_corrected_piston/full_source.npz')
+          fresh = Path('owc-orifice-corrected-output/fresh_source.npz')
+          with np.load(archived) as expected, np.load(fresh) as actual:
+              assert set(expected.files) == set(actual.files)
+              for name in expected.files:
+                  np.testing.assert_allclose(
+                      actual[name], expected[name], rtol=0, atol=1e-9,
+                      err_msg=name,
+                  )
+          PY
+      - name: Gate freshly generated source against independent Python motion and air outputs
+        run: python -m pytest -q tests/test_owc_orifice_corrected_full.py
+        env:
+          WEC_SIM_OWC_CORRECTED_SOURCE: owc-orifice-corrected-output/fresh_source.npz
+      - uses: actions/upload-artifact@v4
+        if: always()
+        with:
+          name: owc-orifice-corrected-source
+          path: |
+            owc-orifice-corrected-output/*.csv
+            owc-orifice-corrected-output/fresh_source.npz
+          if-no-files-found: ignore

PARITY.md

@@ -539,6 +539,25 @@ reaches Mach 1.113. This establishes full-duration physical-path agreement
 against an explicitly **derived** model correction; the unmodified published
 application remains qualified, and badge counts stay unchanged.
 
+The original finer corrected-model source logs now also provide a
+[full-output physical gate](tests/test_owc_orifice_corrected_full.py). Its
+pointwise 130 s maximum Python/source differences are 0.357/0.245 mm/s in
+surge/heave speed, 0.165 mrad/s pitch speed, 0.163 mm/s piston speed,
+0.450 N PTO force, 0.00229 kPa pressure, 0.0000320 m³/s flow, and
+0.0000405 kW PTO power. Independently integrated Python PTO work is 865.042 J
+against 864.422 J from the source's full internal solver trace (0.0717%).
+The source orifice signals also reconstruct from its logged piston speed and
+published air law. The [reproduction workflow](.github/workflows/owc-orifice-corrected-source.yml)
+runs the proposed application in [WEC-Sim_Applications PR #119](https://github.com/WEC-Sim/WEC-Sim_Applications/pull/119)
+from a pinned fork commit and pairs fresh outputs with independent Python.
+WEC-Sim initialization resets the model MaxStep from `simu.dt`; the proposed
+published input therefore sets `simu.dt = 0.00015625` s while retaining 0.005 s
+output and force-sampling clocks. A [full 130 s hosted run](https://github.com/cmudrc/wec-sim-python/actions/runs/39877611225)
+of that input passed the independent motion and air-system gate and reproduced
+the earlier corrected archive within `3.96e-11`. A [full 0.005 s solver-step audit](https://github.com/cmudrc/wec-sim-python/actions/runs/39871054213)
+found 1.687 m surge and a 1.62% PTO-work change. The case remains qualified
+until the correction is merged and pinned as the maintained application revision.
+
 A separate [local 30 s OWC clock audit](tests/fixtures/matlab/owc_orifice_dt_30s/README.md)
 sets MATLAB `simu.dt` and `ode23t` MaxStep together to `0.005`, `0.0025`,
 and `0.00125` s. The first arm is the archived published-clock prefix; all

tests/fixtures/matlab/owc_orifice_corrected_piston/README.md

@@ -8,6 +8,8 @@ published model, [this MATLAB script](../../../matlab/owc_orifice_correct_piston
 adds the existing piston signal to the state-space force input. It leaves the
 reported acceleration route, rigid-body coupling, orifice law, HDF5, and wave
 settings alone. This is a derived physical model, not the published application.
+The corresponding correction to the maintained application is under review in
+[WEC-Sim_Applications PR #119](https://github.com/WEC-Sim/WEC-Sim_Applications/pull/119).
 
 The source checkouts are WEC-Sim `0753b2e47f2457c078751dcfe5d251d1767b80ab`
 and Applications `d53d4d4c9eda2581f04204f5d394a6ef84bb099e`. The original
@@ -74,3 +76,53 @@ is 0.155. That flag does not switch force laws. The one-connection edit
 removes the large physical-path discrepancy; it does **not** establish parity
 with the original published Simulink model. `OWC/OrificeModel` remains
 `qualified_motion`, and strict badge counts do not change.
+
+## Full corrected-model output gate
+
+The original finer-arm CSVs above were packed with
+[`tools/pack_owc_orifice_corrected.py`](../../../../tools/pack_owc_orifice_corrected.py)
+into [`full_source.npz`](full_source.npz) (SHA-256
+`aecd4eea595a24f3450a637600871b6de4b8b7e3fe100c2362a2e433ea906c9e`).
+The archive retains all 26,001 output-grid body positions and velocities,
+flexible outputs, source piston/orifice force, pressure, flow, power, wave,
+and PM components. The orifice outputs were interpolated from the original
+832,001 strictly increasing internal samples onto the unchanged 0.005 s grid;
+the source work and maximum Mach number were calculated from the full internal
+record. The unchanged 130 s source work is 864.422 J and maximum Mach is
+0.1550, below the configured 0.3 warning threshold.
+
+The [full-output gate](../../../test_owc_orifice_corrected_full.py) independently
+advances the public Python physical model from the pinned HDF5 and wave phases.
+It first checks the exported source orifice signals against the separately
+saved source piston speed and the published air law, then compares Python
+motion and air outputs pointwise. Its measured maximum gaps are:
+
+| Quantity | Maximum Python/source gap | Declared limit |
+| --- | ---: | ---: |
+| Surge, heave, pitch, flexible displacement | 2.246 mm, 0.116 mm, 0.359 mrad, 0.022 mm | 20 mm, 10 mm, 1 mrad, 2 mm |
+| Surge, heave, pitch, piston speed | 0.357 mm/s, 0.245 mm/s, 0.165 mrad/s, 0.163 mm/s | 20 mm/s, 20 mm/s, 2 mrad/s, 5 mm/s |
+| PTO force, air pressure, flow, PTO power | 0.450 N, 0.00229 kPa, 0.0000320 m³/s, 0.0000405 kW | 3% of source peak for force, pressure, flow; 4% for power |
+| Integrated PTO work | 0.619 J, 0.0717% | 1% |
+
+The source reconstruction and Python gates pass on the local Apple silicon
+environment. [`owc_orifice_corrected_reference.m`](../../../matlab/owc_orifice_corrected_reference.m)
+and its [MATLAB workflow](../../../../.github/workflows/owc-orifice-corrected-source.yml)
+run the exact corrected model proposed in upstream PR #119 at fork commit
+`3c0c05fefce55125d8ea084deedbd2defb103811` and pinned WEC-Sim core on
+CI, export independently logged source output, and pair that fresh output with
+Python. This candidate is saved in the original R2023b format (SHA-256
+`5f9f1da3bf73e85b52224a3ff483991c3ec60774b92856e3c45cda4347927b59`).
+WEC-Sim initialization resets the model MaxStep from `simu.dt`, so the proposed
+published input sets `simu.dt = 0.00015625` s while explicitly preserving
+0.005 s output, radiation-convolution, nonlinear-force, and Morison-force
+clocks. The proposed input file has SHA-256
+`e5f19816c78e5e4316b30c426035d1f06b22c4c313ac5ec373bfeca6b084f10f`.
+A [full 130 s hosted run of that input](https://github.com/cmudrc/wec-sim-python/actions/runs/39877611225)
+passed the independent Python gate. Its fresh body, flexible, orifice, wave,
+component, and summary arrays differ from the archived refined trace by at
+most `3.96e-11`; work is 864.422 J and peak Mach is 0.155. A separate
+[0.005 s solver-step audit](https://github.com/cmudrc/wec-sim-python/actions/runs/39871054213)
+showed 1.687 m surge and 1.62% PTO-work changes. The workflow now requires
+`1e-9` reproduction of the archive from the proposed published input.
+The case remains qualified until the upstream correction is merged and pinned
+as the maintained application revision.

tests/matlab/owc_orifice_corrected_reference.m

@@ -0,0 +1,100 @@
+% Re-run the pinned corrected OWC model proposed upstream in PR #119.
+% Run as a script so Simulink mask and variant initialization can resolve the
+% WEC-Sim variables in the same workspace as initializeWecSim.
+repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath'))));
+coreRoot = fullfile(repoRoot, 'matlab-ref');
+publishedDir = fullfile(repoRoot, 'applications', 'OWC', 'OrificeModel');
+caseDir = fullfile(repoRoot, 'owc-orifice-corrected-case');
+outputDir = fullfile(repoRoot, 'owc-orifice-corrected-output');
+assert(isfile(fullfile(publishedDir, 'OWC_GBM.slx')) && ...
+    isfile(fullfile(publishedDir, 'hydroData', 'test17a_clean.h5')));
+assert(isfolder(fullfile(coreRoot, 'source')));
+assert(~isfolder(caseDir), 'Do not reuse a previously modified model');
+[copied, message] = copyfile(publishedDir, caseDir);
+assert(copied, message);
+mkdir(outputDir);
+
+% Avoid opening the Mechanics Explorer or post-run plots on the CI runner.
+inputFile = fullfile(caseDir, 'wecSimInputFile.m');
+input = fileread(inputFile);
+old = "simu.explorer = 'on';";
+assert(count(input, old) == 1, 'Published input changed');
+input = replace(input, old, "simu.explorer = 'off';");
+fid = fopen(inputFile, 'w');
+assert(fid ~= -1);
+fprintf(fid, '%s', input);
+fclose(fid);
+fid = fopen(fullfile(caseDir, 'userDefinedFunctions.m'), 'w');
+assert(fid ~= -1);
+fprintf(fid, '%% Plotting disabled for the source comparison.\n');
+fclose(fid);
+
+addpath(genpath(fullfile(coreRoot, 'source')));
+modelFile = fullfile(caseDir, 'OWC_GBM.slx');
+load_system(modelFile);
+sub = 'OWC_GBM/Flex Body/FlexibleBody';
+sumBlock = [sub '/Piston input to GBM'];
+assert(strcmp(get_param(sumBlock, 'BlockType'), 'Sum'));
+assert(strcmp(get_param(sumBlock, 'Inputs'), '++'));
+productPorts = get_param([sub '/Product4'], 'PortHandles');
+sumPorts = get_param(sumBlock, 'PortHandles');
+expectedInputs = {[sub '/Excitation'], [sub '/Gain']};
+for i = 1:2
+    line = get_param(sumPorts.Inport(i), 'Line');
+    source = get_param(line, 'SrcBlockHandle');
+    assert(strcmp(getfullname(source), expectedInputs{i}));
+end
+line = get_param(productPorts.Inport(2), 'Line');
+source = get_param(line, 'SrcBlockHandle');
+assert(strcmp(getfullname(source), sumBlock));
+close_system('OWC_GBM', 0);
+previousDir = pwd;
+cd(caseDir);
+run('initializeWecSim');
+assert(simu.endTime == 130 && simu.dt == .00015625 && ...
+    simu.dtOut == .005 && simu.cicDt == .005 && ...
+    simu.nonlinearDt == .005 && simu.morisonDt == .005 && ...
+    strcmp(simu.solver, 'ode23t') && ...
+    waves.phaseSeed == 128 && waves.height == 1 && waves.period == 4, ...
+    'Published OWC settings changed');
+[~, model, ~] = fileparts(simu.simMechanicsFile);
+assert(abs(str2double(get_param(model, 'MaxStep')) - .00015625) < 1e-12, ...
+    'Published OWC input must resolve the refined MaxStep');
+try
+    % matlab-actions runs this script in the base workspace; the model's
+    % masked blocks resolve the initialized body and variant values there.
+    sim(simu.simMechanicsFile);
+catch error
+    fprintf(2, '%s\n', getReport(error, 'extended', 'hyperlinks', 'off'));
+    rethrow(error);
+end
+run('stopWecSim');
+
+time = output.bodies(1).time(:);
+assert(numel(time) == 26001 && numel(Flex_out.time) == numel(time));
+assert(max(abs(Flex_out.time(:) - time)) < 1e-10);
+assert(numel(orificeData.time) > 800000 && ...
+    all(diff(orificeData.time(:)) > 0));
+assert(size(orificeData.signals.values, 2) == 5);
+orifice = interp1(orificeData.time(:), ...
+    orificeData.signals.values, time, 'linear');
+assert(all(isfinite(orifice), 'all'));
+assert(~any(orificeData.signals.values(:, 2)), ...
+    'Corrected source exceeded the configured Mach limit');
+
+moving = output.bodies(1);
+writematrix([time, moving.position, moving.velocity], ...
+    fullfile(outputDir, 'body.csv'));
+writematrix([time, Flex_out.signals.values], ...
+    fullfile(outputDir, 'flex.csv'));
+writematrix([time, orifice], fullfile(outputDir, 'orifice.csv'));
+writematrix([output.wave.time(:), output.wave.elevation(:)], ...
+    fullfile(outputDir, 'wave.csv'));
+writematrix([waves.omega(:), waves.amplitude(:), waves.dOmega(:), ...
+    waves.phase(:)], fullfile(outputDir, 'components.csv'));
+writematrix([trapz(orificeData.time(:), ...
+    orificeData.signals.values(:, 5)) * 1000, ...
+    max(abs(orificeData.signals.values(:, 4) / A2) / 343.2)], ...
+    fullfile(outputDir, 'summary.csv'));
+close_system(model, 0);
+cd(previousDir);

tests/matlab/owc_orifice_saved_step_reference.m

@@ -0,0 +1,98 @@
+% Audit the corrected upstream OWC model at its unchanged saved MaxStep.
+% Run as a script so Simulink mask and variant initialization can resolve the
+% WEC-Sim variables in the same workspace as initializeWecSim.
+repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath'))));
+coreRoot = fullfile(repoRoot, 'matlab-ref');
+publishedDir = fullfile(repoRoot, 'applications', 'OWC', 'OrificeModel');
+caseDir = fullfile(repoRoot, 'owc-orifice-corrected-case');
+outputDir = fullfile(repoRoot, 'owc-orifice-saved-step-output');
+assert(isfile(fullfile(publishedDir, 'OWC_GBM.slx')) && ...
+    isfile(fullfile(publishedDir, 'hydroData', 'test17a_clean.h5')));
+assert(isfolder(fullfile(coreRoot, 'source')));
+assert(~isfolder(caseDir), 'Do not reuse a previously modified model');
+[copied, message] = copyfile(publishedDir, caseDir);
+assert(copied, message);
+mkdir(outputDir);
+
+% Avoid opening the Mechanics Explorer or post-run plots on the CI runner.
+inputFile = fullfile(caseDir, 'wecSimInputFile.m');
+input = fileread(inputFile);
+old = "simu.explorer = 'on';";
+assert(count(input, old) == 1, 'Published input changed');
+input = replace(input, old, "simu.explorer = 'off';");
+fid = fopen(inputFile, 'w');
+assert(fid ~= -1);
+fprintf(fid, '%s', input);
+fclose(fid);
+fid = fopen(fullfile(caseDir, 'userDefinedFunctions.m'), 'w');
+assert(fid ~= -1);
+fprintf(fid, '%% Plotting disabled for the source comparison.\n');
+fclose(fid);
+
+addpath(genpath(fullfile(coreRoot, 'source')));
+modelFile = fullfile(caseDir, 'OWC_GBM.slx');
+load_system(modelFile);
+sub = 'OWC_GBM/Flex Body/FlexibleBody';
+sumBlock = [sub '/Piston input to GBM'];
+assert(strcmp(get_param(sumBlock, 'BlockType'), 'Sum'));
+assert(strcmp(get_param(sumBlock, 'Inputs'), '++'));
+productPorts = get_param([sub '/Product4'], 'PortHandles');
+sumPorts = get_param(sumBlock, 'PortHandles');
+expectedInputs = {[sub '/Excitation'], [sub '/Gain']};
+for i = 1:2
+    line = get_param(sumPorts.Inport(i), 'Line');
+    source = get_param(line, 'SrcBlockHandle');
+    assert(strcmp(getfullname(source), expectedInputs{i}));
+end
+line = get_param(productPorts.Inport(2), 'Line');
+source = get_param(line, 'SrcBlockHandle');
+assert(strcmp(getfullname(source), sumBlock));
+close_system('OWC_GBM', 0);
+previousDir = pwd;
+cd(caseDir);
+run('initializeWecSim');
+assert(simu.endTime == 130 && simu.dt == .005 && ...
+    simu.dtOut == .005 && strcmp(simu.solver, 'ode23t') && ...
+    waves.phaseSeed == 128 && waves.height == 1 && waves.period == 4, ...
+    'Published OWC settings changed');
+[~, model, ~] = fileparts(simu.simMechanicsFile);
+assert(abs(str2double(get_param(model, 'MaxStep')) - .005) < 1e-12);
+try
+    % matlab-actions runs this script in the base workspace; the model's
+    % masked blocks resolve the initialized body and variant values there.
+    sim(simu.simMechanicsFile);
+catch error
+    fprintf(2, '%s\n', getReport(error, 'extended', 'hyperlinks', 'off'));
+    rethrow(error);
+end
+run('stopWecSim');
+
+time = output.bodies(1).time(:);
+assert(numel(time) == 26001 && numel(Flex_out.time) == numel(time));
+assert(max(abs(Flex_out.time(:) - time)) < 1e-10);
+assert(numel(orificeData.time) > 1000 && ...
+    all(diff(orificeData.time(:)) > 0));
+assert(size(orificeData.signals.values, 2) == 5);
+orifice = interp1(orificeData.time(:), ...
+    orificeData.signals.values, time, 'linear');
+assert(all(isfinite(orifice), 'all'));
+
+moving = output.bodies(1);
+writematrix([time, moving.position, moving.velocity], ...
+    fullfile(outputDir, 'body.csv'));
+writematrix([time, Flex_out.signals.values], ...
+    fullfile(outputDir, 'flex.csv'));
+writematrix([time, orifice], fullfile(outputDir, 'orifice.csv'));
+writematrix([output.wave.time(:), output.wave.elevation(:)], ...
+    fullfile(outputDir, 'wave.csv'));
+writematrix([waves.omega(:), waves.amplitude(:), waves.dOmega(:), ...
+    waves.phase(:)], fullfile(outputDir, 'components.csv'));
+writematrix([trapz(orificeData.time(:), ...
+    orificeData.signals.values(:, 5)) * 1000, ...
+    max(abs(orificeData.signals.values(:, 4) / A2) / 343.2), ...
+    numel(orificeData.time), ...
+    mean(orificeData.signals.values(:, 2) ~= 0), ...
+    max(abs(orificeData.signals.values(:, 3)))], ...
+    fullfile(outputDir, 'summary.csv'));
+close_system(model, 0);
+cd(previousDir);