diff --git a/.github/workflows/owc-orifice-corrected-source.yml b/.github/workflows/owc-orifice-corrected-source.yml new file mode 100644 index 00000000..310b4916 --- /dev/null +++ b/.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 diff --git a/PARITY.md b/PARITY.md index c97bacf7..241fd7d7 100644 --- a/PARITY.md +++ b/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 diff --git a/tests/fixtures/matlab/owc_orifice_corrected_piston/README.md b/tests/fixtures/matlab/owc_orifice_corrected_piston/README.md index 01a382d2..53c0330d 100644 --- a/tests/fixtures/matlab/owc_orifice_corrected_piston/README.md +++ b/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. diff --git a/tests/fixtures/matlab/owc_orifice_corrected_piston/full_source.npz b/tests/fixtures/matlab/owc_orifice_corrected_piston/full_source.npz new file mode 100644 index 00000000..52b7d1bd Binary files /dev/null and b/tests/fixtures/matlab/owc_orifice_corrected_piston/full_source.npz differ diff --git a/tests/matlab/owc_orifice_corrected_reference.m b/tests/matlab/owc_orifice_corrected_reference.m new file mode 100644 index 00000000..a9a43f11 --- /dev/null +++ b/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); diff --git a/tests/matlab/owc_orifice_saved_step_reference.m b/tests/matlab/owc_orifice_saved_step_reference.m new file mode 100644 index 00000000..0c9bc2a3 --- /dev/null +++ b/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); diff --git a/tests/test_owc_orifice_corrected_full.py b/tests/test_owc_orifice_corrected_full.py new file mode 100644 index 00000000..b1c16b4a --- /dev/null +++ b/tests/test_owc_orifice_corrected_full.py @@ -0,0 +1,128 @@ +"""Gate the independently advanced OWC against the corrected source model. + +The corrected model adds piston force to the flexible state input. Its source +trace is derived from the published application, so this is not a published +input badge gate. +""" + +import os +from pathlib import Path + +import numpy as np + +from wecsim import OrificePTO, PMWave, WEC + + +FIXTURES = Path(__file__).parent / "fixtures/matlab" +SOURCE = Path(os.environ.get( + "WEC_SIM_OWC_CORRECTED_SOURCE", + FIXTURES / "owc_orifice_corrected_piston/full_source.npz", +)) +ORIGINAL = FIXTURES / "owc_orifice_corrected_piston/source.npz" +PUBLISHED = FIXTURES / "owc_orifice_130s/source.npz" +HYDRO = FIXTURES / "owc_orifice.h5" + + +def _peak(actual, expected, limit, name): + assert actual.shape == expected.shape, name + assert np.isfinite(actual).all() and np.isfinite(expected).all(), name + error = float(np.max(np.abs(actual - expected))) + assert error <= limit, f"{name}: {error:.8g} > {limit:.8g}" + + +def test_corrected_owc_full_motion_and_air_system(tmp_path): + with np.load(SOURCE, allow_pickle=False) as saved: + source = {name: saved[name] for name in saved.files} + with np.load(ORIGINAL, allow_pickle=False) as saved: + existing = {name: saved[name] for name in saved.files} + with np.load(PUBLISHED, allow_pickle=False) as saved: + phases = saved["phases"].copy() + parameters = saved["orifice_parameters"].copy() + + body = source["body"] + flex = source["flex"] + air = source["orifice"] + time = body[:, 0] + _peak(time, existing["time"], 1e-10, "output time") + _peak(source["wave"][:, 1], existing["wave"], 1e-10, + "incident wave") + _peak(source["components"], existing["components"], 1e-12, + "PM components and phases") + _peak(source["components"][:, 3], phases, 1e-12, "PM phase input") + assert not np.any(air[:, 2]) + assert source["summary"][1] < parameters[4] + + # The exported fluid signals must agree with the source piston speed and + # orifice law before comparing independently advanced Python outputs. + orifice = OrificePTO(*parameters) + source_air = orifice.evaluate(flex[:, 2]) + _peak(source_air.force, air[:, 1], 1e-5, "source PTO force law") + _peak(source_air.pressure_drop / 1000, air[:, 3], 1e-7, + "source pressure law") + _peak(source_air.flow_rate, air[:, 4], 1e-10, "source flow law") + _peak(source_air.absorbed_power / 1000, air[:, 5], 1e-7, + "source power law") + + phase_file = tmp_path / "phases.csv" + np.savetxt(phase_file, phases, delimiter=",") + wec = WEC("OWC corrected piston-path gate") + floating = wec.body("OWC", HYDRO, inertia=(99.28, 11.04, 99.2)) + wec.floating_gbm( + floating, orifice=orifice, + heave_linear_damping=100, mode_linear_damping=100, + heave_drag_cd=1.2, heave_drag_area=8, + pitch_drag_cd=1.2, pitch_drag_area=8, + ) + result = wec.run( + PMWave(height=1, period=4, phase_file=phase_file), + dt=.005, end_time=130, ramp_time=10, radiation_memory=15, + ) + _peak(result.time, time, 1e-10, "Python time") + _peak(result.wave_elevation, source["wave"][:, 1], 1e-10, + "Python incident wave") + + pose = result.bodies["OWC"].position + velocity = result.bodies["OWC"].velocity + flexible = result.flexible_modes["OWC"] + source_pitch = np.r_[0, np.cumsum( + (body[1:, 11] + body[:-1, 11]) * .005 / 2, + )] + source_motion = np.column_stack(( + body[:, 1], body[:, 3], source_pitch, flex[:, 1], + )) + source_reproduction = np.max(np.abs( + source_motion - existing["fine_motion"]), axis=0) + assert np.all(source_reproduction <= [.020, .010, .001, .002]), ( + source_reproduction + ) + assert abs(float(source["summary"][0]) + - float(existing["fine_work_internal"])) / float( + existing["fine_work_internal"]) <= .005 + _peak(pose[:, 0], body[:, 1], .020, "surge") + _peak(pose[:, 2], body[:, 3], .010, "heave") + _peak(result.coordinates["pitch_unwrapped"].position, + source_pitch, .001, "pitch") + _peak(flexible.position[:, 0], flex[:, 1], .002, + "flexible displacement") + _peak(velocity[:, 0], body[:, 7], .020, "surge speed") + _peak(velocity[:, 2], body[:, 9], .020, "heave speed") + _peak(velocity[:, 4], body[:, 11], .002, "pitch speed") + _peak(flexible.velocity[:, 0], flex[:, 2], .005, "piston speed") + + pto = result.ptos["orifice"] + python_air = orifice.evaluate(pto.velocity) + force_scale = np.max(np.abs(air[:, 1])) + pressure_scale = np.max(np.abs(air[:, 3])) + flow_scale = np.max(np.abs(air[:, 4])) + power_scale = np.max(np.abs(air[:, 5])) + assert min(force_scale, pressure_scale, flow_scale, power_scale) > 0 + _peak(pto.force, air[:, 1], .03 * force_scale, "PTO force") + _peak(python_air.pressure_drop / 1000, air[:, 3], + .03 * pressure_scale, "air pressure") + _peak(python_air.flow_rate, air[:, 4], .03 * flow_scale, "air flow") + _peak(pto.absorbed_power / 1000, air[:, 5], + .04 * power_scale, "PTO power") + assert not np.any(python_air.compressibility_flag) + source_work = float(source["summary"][0]) + python_work = float(np.trapezoid(pto.absorbed_power, time)) + assert abs(python_work - source_work) / source_work <= .01 diff --git a/tools/pack_owc_orifice_corrected.py b/tools/pack_owc_orifice_corrected.py new file mode 100644 index 00000000..79b09eac --- /dev/null +++ b/tools/pack_owc_orifice_corrected.py @@ -0,0 +1,63 @@ +"""Pack a corrected OWC MATLAB CI artifact into the checked source fixture.""" + +import argparse +from pathlib import Path + +import numpy as np + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("artifact_dir", type=Path) + parser.add_argument("output", type=Path) + args = parser.parse_args() + + names = ("body", "flex", "orifice", "wave", "components") + source = { + name: np.loadtxt(args.artifact_dir / f"{name}.csv", delimiter=",") + for name in names + } + source["body"] = source["body"][:, :13] + time = np.arange(26_001) * .005 + internal_orifice = source["orifice"] + if internal_orifice.shape[0] > len(time): + assert internal_orifice.shape[1] == 6 + assert np.all(np.diff(internal_orifice[:, 0]) > 0) + with np.load( + Path(__file__).parents[1] + / "tests/fixtures/matlab/owc_orifice_130s/source.npz", + allow_pickle=False, + ) as published: + orifice_area = float(published["orifice_parameters"][1]) + source["summary"] = np.array([ + np.trapezoid(internal_orifice[:, 5] * 1000, + internal_orifice[:, 0]), + np.max(np.abs(internal_orifice[:, 4] / orifice_area) / 343.2), + ]) + source["orifice"] = np.column_stack(( + time, + *[np.interp(time, internal_orifice[:, 0], internal_orifice[:, i]) + for i in range(1, 6)], + )) + else: + source["summary"] = np.loadtxt( + args.artifact_dir / "summary.csv", delimiter=",", + ) + assert source["body"].shape == (26_001, 13) + assert source["flex"].shape == (26_001, 11) + assert source["orifice"].shape == (26_001, 6) + assert source["wave"].shape == (26_001, 2) + assert source["components"].shape == (500, 4) + assert source["summary"].shape == (2,) + for name in ("body", "flex", "orifice", "wave"): + np.testing.assert_allclose(source[name][:, 0], time, rtol=0, atol=1e-10) + assert np.isfinite(source[name]).all() + assert np.isfinite(source["components"]).all() + assert source["summary"][0] > 0 + assert 0 <= source["summary"][1] < .3 + args.output.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(args.output, **source) + + +if __name__ == "__main__": + main()