diff --git a/.github/workflows/ellipsoid-ode45-mass-correction.yml b/.github/workflows/ellipsoid-ode45-mass-correction.yml new file mode 100644 index 00000000..7534c157 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-mass-correction.yml @@ -0,0 +1,74 @@ +name: Ellipsoid ode45 derived heave mass correction + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-mass-correction.yml' + - 'tests/matlab/ellipsoid_ode45_mass_correction.m' + - 'tests/test_ellipsoid_ode45_mass_correction.py' + - 'tests/fixtures/matlab/ellipsoid_ode45_mass_correction/**' + - 'tools/pack_ellipsoid_ode45_mass_correction.py' + workflow_dispatch: + +jobs: + compare: + strategy: + fail-fast: false + matrix: + wave: [Regular, RegularCIC] + runs-on: ubuntu-latest + timeout-minutes: 150 + name: ${{ matrix.wave }} + 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: Nonlinear_Hydro + - 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'); waveCase='${{ matrix.wave }}'; mode='published_clock'; ellipsoid_ode45_mass_correction + - uses: matlab-actions/run-command@v3 + with: + command: addpath('tests/matlab'); waveCase='${{ matrix.wave }}'; mode='refined_clock'; ellipsoid_ode45_mass_correction + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e ".[test]" + - run: python -m pytest -q tests/test_ellipsoid_ode45_delay_prefix.py + - name: Repack and compare the archived source traces + run: >- + python tools/pack_ellipsoid_ode45_mass_correction.py + "${{ matrix.wave }}" + "matlab-ellipsoid-ode45-mass-correction/${{ matrix.wave }}" + "matlab-ellipsoid-ode45-mass-correction/${{ matrix.wave }}/packed.npz" + --expect "tests/fixtures/matlab/ellipsoid_ode45_mass_correction/${{ matrix.wave == 'Regular' && 'regular.npz' || 'regular_cic.npz' }}" + - name: Check derived mass identity and independent 150 s motion + run: python -m pytest -q -s tests/test_ellipsoid_ode45_mass_correction.py + env: + WEC_SIM_APPLICATIONS_DIR: applications + WEC_SIM_ELLIPSOID_MASS_CORRECTION_DIR: matlab-ellipsoid-ode45-mass-correction + WEC_SIM_ELLIPSOID_CASE: ${{ matrix.wave }} + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-mass-correction-${{ matrix.wave }} + path: | + matlab-ellipsoid-ode45-mass-correction/${{ matrix.wave }}/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + if-no-files-found: ignore diff --git a/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/README.md b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/README.md new file mode 100644 index 00000000..c5ad9a93 --- /dev/null +++ b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/README.md @@ -0,0 +1,28 @@ +# Derived ellipsoid ode45 heave mass correction + +This is an on-demand investigation of [issue #599](https://github.com/cmudrc/wec-sim-python/issues/599), not a replacement for either published `Nonlinear_Hydro/ode45` input. The source is WEC-Sim core `0753b2e47f2457c078751dcfe5d251d1767b80ab` and Applications `d53d4d4c9eda2581f04204f5d394a6ef84bb099e`. The [source script](../../../matlab/ellipsoid_ode45_mass_correction.m) copies each application model and input into a disposable sibling directory. It does not write to the pinned WEC-Sim or Applications repositories. + +The core `bodyClass.adjustMassMatrix` shifts `D = adjMassFactor × trace(A_xyz)` into the rigid translational mass and leaves `A_zz − D` in the acceleration force sent through a 100 ns Transport Delay. The copied input sets `adjMassFactor = A_zz / trace(A_xyz)` from the same HDF5 and wave mode. For these **single active heave-coordinate** cases, the corrected model therefore has effective heave mass `m + A_zz` and zero active delayed coefficient. It retains the linked core block and the original wave, geometry, damping, PTO, and 150 s duration. The scalar shift also changes inactive surge and sway masses, so this is not a general multidimensional core fix. + +The [R2025b reproduction run](https://github.com/cmudrc/wec-sim-python/actions/runs/40102870958) at Python-repository head `963f97442da8258578c1fde94c9f6d3bcb78050e` completed both corrected full-duration source arms for each wave case. The first uses the published 0.05 s solver, output, nonlinear-force, and CIC clocks. The second changes solver, nonlinear-force, and CIC clocks to 0.005 s while retaining the 0.05 s output grid. The original published input is unchanged. Both arms have 3,001 output samples and identical wave elevation at shared times. The initialized heave delayed coefficient is zero, the actual applied delayed heave force stays below `6.0e-22` N, and all inactive coordinate accelerations stay below `1.3e-21` in their respective units. + +| Maximum over 150 s | Regular | RegularCIC | +| --- | ---: | ---: | +| Derived mass factor | 0.655990137167 | 0.783566018956 | +| Corrected source 0.05→0.005 s heave change | 21.268 mm | 21.414 mm | +| Unchanged source 0.05→0.005 s heave change | 17.090 mm | 17.919 mm | +| Corrected 0.05 s versus independent Python heave | 23.525 mm | 23.770 mm | +| Corrected 0.05 s versus Python speed | 29.450 mm/s | 31.858 mm/s | +| Corrected 0.05 s versus Python PTO force | 35.340 kN | 38.230 kN | +| Corrected 0.05 s versus Python absorbed power | 28.745 kW | 30.983 kW | +| Corrected 0.005 s versus independent Python heave | 2.265 mm | 2.378 mm | +| Corrected 0.005 s versus Python speed | 2.886 mm/s | 2.995 mm/s | +| Corrected 0.005 s versus Python PTO force | 3.464 kN | 3.594 kN | +| Corrected 0.005 s versus Python absorbed power | 2.804 kW | 2.983 kW | +| Corrected 0.005 s versus unchanged 0.005 s heave | 0.651 mm | 0.531 mm | + +The [portable Python gate](../../../test_ellipsoid_ode45_mass_correction.py) advances motion and PTO independently from the archived source-generated `ellipsoid.h5`, the byte-identical pinned STL, and the published wave settings. Its refined 0.005 s comparisons pass the existing pointwise physical limits: 3/3.5 mm heave, 4/4.5 mm/s speed, 5/6 kN force, and 4/5 kW power for Regular/RegularCIC. Strict expected failures apply those same physical limits to the corrected published clock and retain the source-step stability limits on the corrected copy; an unexpected pass fails the test. The separate gate against each **unchanged original** published input remains unchanged and unpaired. The corrected published-clock heave misses both the 3/3.5 mm physical limits and the broader 22/23 mm numerical solver envelope used by the unchanged-input test. The corrected source's 21 mm self-change exceeds the established 3/3.75 mm heave stability screens. Removing the heave added-mass feedback alone therefore does not close the published-clock gap. Case status and badges remain unchanged. + +`regular.npz` and `regular_cic.npz` archive the corrected source's time, heave, speed, PTO force/speed, wave, applied delayed heave force, inactive speeds/accelerations, and mass identity for both clocks. They were packed from the two uploaded MATLAB artifacts with [`tools/pack_ellipsoid_ode45_mass_correction.py`](../../../../tools/pack_ellipsoid_ode45_mass_correction.py); that script checks their full output grids before packing. The 882 KiB `ellipsoid.h5` here is from the Regular source artifact (SHA-256 `0092f6dd4a00c1ce6c6836df1a06bb3047ffa0e3f32038d96bcad4a408b5664c`). The two newly generated HDF5 files have different binary hashes but all 38 HDF5 datasets agree numerically, treating matching NaNs as equal. The older bundled test HDF5 has some additional datasets and minor differences in displaced volume and stored impulse-response arrays, so this portable gate uses the archived source HDF5. The STL in Applications and the bundled test geometry are byte-identical (SHA-256 `a225213457fb9273a327a141273b9d68266e650cf973a98f84e84e8607cbb7ec`). + +The [workflow](../../../../.github/workflows/ellipsoid-ode45-mass-correction.yml) regenerates source evidence on request or when its diagnostic code and fixtures change, and checks fresh traces against the archived fixtures. The compact fixtures make the independent refined physical gate and both strict expected original-clock failures run in ordinary Python CI without MATLAB. diff --git a/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/ellipsoid.h5 b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/ellipsoid.h5 new file mode 100644 index 00000000..eb52154a Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/ellipsoid.h5 differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/regular.npz b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/regular.npz new file mode 100644 index 00000000..52a9cfc0 Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/regular.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/regular_cic.npz b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/regular_cic.npz new file mode 100644 index 00000000..292bb56a Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_mass_correction/regular_cic.npz differ diff --git a/tests/matlab/ellipsoid_ode45_mass_correction.m b/tests/matlab/ellipsoid_ode45_mass_correction.m new file mode 100644 index 00000000..eb1a2482 --- /dev/null +++ b/tests/matlab/ellipsoid_ode45_mass_correction.m @@ -0,0 +1,139 @@ +% Test a heave-only, exact added-mass loop break on a disposable source input. +% Invoke as a script with waveCase and mode in the base workspace. +assert(any(strcmp(waveCase, {'Regular', 'RegularCIC'}))); +assert(any(strcmp(mode, {'published_clock', 'refined_clock'}))); +repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath')))); +sourceDir = fullfile(repoRoot, 'matlab-ref', 'source'); +appDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro'); +hydroDir = fullfile(appDir, 'hydroData'); +originalDir = fullfile(appDir, 'ode45', waveCase); +addpath(genpath(sourceDir)); +cd(hydroDir); +if ~isfile('ellipsoid.h5') + bemio; +end + +% The published bodyClass adjustment uses m + 2*trace(A_xyz) as its rigid +% translational mass and feeds the remaining A_zz - 2*trace(A_xyz) through +% the 100 ns Transport Delay. This application has one active heave DOF. +% Moving exactly A_zz into the rigid mass makes the active delay coefficient +% zero while keeping the physical effective heave mass m + A_zz unchanged. +hydro = readBEMIOH5(fullfile(hydroDir, 'ellipsoid.h5'), 1, 0); +if strcmp(waveCase, 'Regular') + omega = 2*pi/6; + frequencies = hydro.simulation_parameters.w; + diagonal = zeros(3, 1); + for axis = 1:3 + diagonal(axis) = interp1(frequencies, ... + squeeze(hydro.hydro_coeffs.added_mass.all(axis, axis, :)), ... + omega, 'spline'); + end +else + matrix = hydro.hydro_coeffs.added_mass.inf_freq; + diagonal = diag(matrix(1:3, 1:3)); +end +factor = diagonal(3)/sum(diagonal); +assert(isfinite(factor) && factor > 0 && factor < 1, ... + 'The pinned heave-only added-mass coefficient is invalid'); + +original = fileread(fullfile(originalDir, 'wecSimInputFile.m')); +assert(contains(original, 'simu.dt = 0.05;') && ... + ~contains(original, 'simu.dtOut') && ... + ~contains(original, 'simu.nonlinearDt') && ... + ~contains(original, 'simu.cicDt') && ... + contains(original, 'body(1).nonlinearHydro = 2;') && ... + ~contains(original, 'body(1).adjMassFactor'), ... + 'The pinned published input has changed'); +caseDir = fullfile(fileparts(originalDir), ... + [waveCase '_mass_correction_' mode]); +assert(~isfolder(caseDir), 'Disposable corrected case already exists'); +mkdir(caseDir); +copyfile(fullfile(originalDir, 'ellipsoid.slx'), ... + fullfile(caseDir, 'ellipsoid.slx')); +corrected = strrep(original, 'body(1).nonlinearHydro = 2;', sprintf( ... + ['body(1).nonlinearHydro = 2;\n' ... + 'body(1).adjMassFactor = %.17g; %% heave-only exact mass shift'], factor)); +if strcmp(mode, 'refined_clock') + corrected = strrep(corrected, 'simu.dt = 0.05;', sprintf( ... + ['simu.dt = 0.005;\n' ... + 'simu.dtOut = 0.05;\n' ... + 'simu.nonlinearDt = 0.005;\n' ... + 'simu.cicDt = 0.005;'])); +end +fid = fopen(fullfile(caseDir, 'wecSimInputFile.m'), 'w'); +assert(fid > 0, 'Could not write disposable corrected input'); +fprintf(fid, '%s', corrected); +fclose(fid); + +cd(caseDir); +run('initializeWecSim'); +expectedStep = 0.05; +if strcmp(mode, 'refined_clock') + expectedStep = 0.005; +end +assert(strcmp(simu.solver, 'ode45') && ... + abs(simu.dt - expectedStep) < 1e-12 && ... + abs(simu.dtOut - 0.05) < 1e-12 && ... + abs(simu.nonlinearDt - expectedStep) < 1e-12 && ... + abs(simu.cicDt - expectedStep) < 1e-12 && ... + simu.endTime == 150 && simu.rampTime == 50, ... + 'An unrelated source setting changed'); +model = 'ellipsoid'; +step = str2double(get_param(model, 'MaxStep')); +assert(abs(step - expectedStep) < 1e-12); +delayBlock = [model '/Rigid Body/Hydrodynamic Body/Transport Delay']; +assert(abs(str2double(get_param(delayBlock, 'DelayTime')) - 1e-7) < 1e-16, ... + 'The copied model has an unexpected core delay'); +force = body(1).hydroForce.hf1; +originalMass = force.storage.fAddedMass; +assert(abs(originalMass(3,3) - simu.rho*diagonal(3)) < 1e-7, ... + 'The derived HDF5 and initialized heave added mass differ'); +assert(abs(sum(diag(originalMass(1:3,1:3)))*factor - originalMass(3,3)) < 1e-7); +assert(abs(force.fAddedMass(3,3)) < 1e-7, ... + 'Active heave acceleration still passes through the delay'); +assert(abs(force.adjustedMass - (force.mass + originalMass(3,3))) < 1e-7, ... + 'The corrected heave mass is not the physical effective mass'); +inactive = [1 2 4 5 6]; +massRecord = [factor, sum(diag(originalMass(1:3,1:3))), ... + originalMass(3,3), force.mass, force.adjustedMass, ... + force.fAddedMass(3,3), step, force.fAddedMass(3,inactive)]; +fprintf(['MASS_CORRECTION case=%s mode=%s factor=%.17g ' ... + 'activeAddedMass=%.9g residual=%.9g effectiveMass=%.9g\n'], ... + waveCase, mode, factor, originalMass(3,3), ... + force.fAddedMass(3,3), force.adjustedMass); + +sim(simu.simMechanicsFile, [], simset('SrcWorkspace', 'parent')); +run('stopWecSim'); +record = output.bodies(1); +% postProcessWecSim recomputes the reported added-mass force. Its storage +% field retains the force actually sent through the delayed Simulink input, +% which must vanish in active heave here. +appliedValues = [record.time(:), ... + body(1).hydroForce.hf1.storage.output_forceAddedMass]; +bodyValues = [record.time(:), record.position, record.velocity, ... + record.forceTotal, record.forceExcitation, ... + record.forceRadiationDamping, record.forceAddedMass, ... + record.forceRestoring, record.forceMorisonAndViscous, ... + record.acceleration]; +pto = output.ptos(1); +ptoValues = [pto.time(:), pto.position, pto.velocity, ... + pto.forceInternalMechanics, pto.powerInternalMechanics]; +waveValues = [output.wave.time(:), output.wave.elevation(:)]; +assert(isequal(size(bodyValues), [3001, 55]) && ... + isequal(size(appliedValues), [3001, 7]) && ... + isequal(size(ptoValues), [3001, 25]) && ... + isequal(size(waveValues), [3001, 2]) && ... + all(isfinite(bodyValues), 'all') && ... + all(isfinite(appliedValues), 'all') && ... + all(isfinite(ptoValues), 'all') && ... + all(isfinite(waveValues), 'all'), 'Incomplete corrected source output'); +outDir = fullfile(repoRoot, 'matlab-ellipsoid-ode45-mass-correction', waveCase); +if ~isfolder(outDir) + mkdir(outDir); +end +writematrix(bodyValues, fullfile(outDir, [mode '_body.csv'])); +writematrix(appliedValues, fullfile(outDir, [mode '_applied_added.csv'])); +writematrix(ptoValues, fullfile(outDir, [mode '_pto.csv'])); +writematrix(waveValues, fullfile(outDir, [mode '_wave.csv'])); +writematrix(massRecord, fullfile(outDir, [mode '_mass.csv'])); +close_system(model, 0); diff --git a/tests/test_ellipsoid_ode45_mass_correction.py b/tests/test_ellipsoid_ode45_mass_correction.py new file mode 100644 index 00000000..00e4252f --- /dev/null +++ b/tests/test_ellipsoid_ode45_mass_correction.py @@ -0,0 +1,262 @@ +"""Portable 150 s gate for a derived, heave-only source mass-loop correction.""" + +import json +import os +from pathlib import Path + +import numpy as np +import pytest + +from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint + + +ROOT = Path(__file__).resolve().parents[1] +APPLICATIONS = os.environ.get("WEC_SIM_APPLICATIONS_DIR") +RECORDS = os.environ.get("WEC_SIM_ELLIPSOID_MASS_CORRECTION_DIR") +SELECTED_CASE = os.environ.get("WEC_SIM_ELLIPSOID_CASE") +ARCHIVE = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_mass_correction" +ORIGINAL = (ROOT / "tests/fixtures/matlab/ellipsoid_ode45_nonlinear_clock" + / "source_maxstep_common_grid.npz") +REFINED = (ROOT / "tests/fixtures/matlab/ellipsoid_ode45_maxstep_only" + / "refined_005_common_grid.npz") + +# The 0.005 s limits are the established independent refined-step gates. +REFINED_PHYSICAL_LIMITS = { + "Regular": (.003, .004, 5000, 4000), + "RegularCIC": (.0035, .0045, 6000, 5000), +} +# The unchanged published-input test treats these as ode45 solver-envelope +# limits, not as a claim of physical parity with ode45 internals. +PUBLISHED_SOLVER_ENVELOPE_LIMITS = { + "Regular": (.022, .032, 38_000, 30_000), + "RegularCIC": (.023, .034, 40_000, 33_000), +} +# Established source step-stability screens from the prior refinement gate. +SOURCE_STEP_LIMITS = { + "Regular": (.003, .00325, 4000, 4000), + "RegularCIC": (.00375, .00425, 5000, 5000), +} + + +class CorrectedPublishedPhysicalGap(AssertionError): + """The corrected 0.05 s source still misses physical parity limits.""" + + +class CorrectedSourceStepGap(AssertionError): + """The corrected source still changes beyond its step-stability screen.""" + + +def _csv(folder, mode, name, shape): + values = np.loadtxt(folder / f"{mode}_{name}.csv", + delimiter=",", ndmin=2) + assert values.shape == shape and np.isfinite(values).all(), (mode, name) + return values + + +def _read(wave_case, mode): + if RECORDS: + folder = Path(RECORDS) / wave_case + body = _csv(folder, mode, "body", (3001, 55)) + pto = _csv(folder, mode, "pto", (3001, 25)) + wave = _csv(folder, mode, "wave", (3001, 2)) + applied = _csv(folder, mode, "applied_added", (3001, 7)) + mass = _csv(folder, mode, "mass", (1, 12))[0] + record = { + "time": body[:, 0], "heave": body[:, 3], + "speed": body[:, 9], "pto_force": pto[:, 15], + "pto_speed": pto[:, 9], "wave": wave[:, 1], + "applied_heave": applied[:, 3], + "inactive_speed": body[:, [7, 8, 10, 11, 12]], + "inactive_acceleration": body[:, [49, 50, 52, 53, 54]], + "mass": mass, + } + for name, values in (("pto", pto), ("wave", wave), + ("applied", applied)): + np.testing.assert_allclose(values[:, 0], record["time"], + rtol=0, atol=1e-10, err_msg=name) + else: + filename = "regular.npz" if wave_case == "Regular" else "regular_cic.npz" + with np.load(ARCHIVE / filename, allow_pickle=False) as source: + assert str(source["case"]) == wave_case + record = {name: source[f"{mode}_{name}"] for name in ( + "time", "heave", "speed", "pto_force", "pto_speed", + "wave", "applied_heave", "inactive_speed", + "inactive_acceleration", "mass", + )} + time = np.arange(3001) * .05 + np.testing.assert_allclose(record["time"], time, rtol=0, atol=1e-10) + for name in ("heave", "speed", "pto_force", "pto_speed", "wave", + "applied_heave"): + assert record[name].shape == (3001,) and np.isfinite(record[name]).all() + for name in ("inactive_speed", "inactive_acceleration"): + assert record[name].shape == (3001, 5) + assert np.isfinite(record[name]).all() + factor, trace, a33, rigid, adjusted, residual, step = record["mass"][:7] + assert record["mass"].shape == (12,) + np.testing.assert_allclose(factor, a33 / trace, rtol=0, atol=1e-13) + np.testing.assert_allclose(adjusted, rigid + a33, rtol=0, atol=1e-7) + assert abs(residual) < 1e-7 + assert abs(step - (.05 if mode == "published_clock" else .005)) < 1e-12 + assert np.max(np.abs(record["inactive_speed"])) < 1e-5 + assert np.max(np.abs(record["inactive_acceleration"])) < 1e-5 + cross_force = record["inactive_acceleration"] @ record["mass"][7:] + assert np.max(np.abs(cross_force)) < 1e-5 + assert np.max(np.abs(record["applied_heave"])) < 1e-5 + return record + + +def _error(left, right): + assert left.shape == right.shape + return float(np.max(np.abs(left - right))) + + +def _source_power(record): + return -record["pto_force"] * record["pto_speed"] + + +def _gaps(position, speed, force, power, source): + return ( + _error(position, source["heave"]), + _error(speed, source["speed"]), + _error(force, source["pto_force"]), + _error(power, _source_power(source)), + ) + + +def _check_gaps(gaps, limits, error_type, label): + names = ("heave position", "heave speed", "PTO force", "absorbed power") + failures = [f"{name}: {actual:.8g} >= {limit:.8g}" + for name, actual, limit in zip(names, gaps, limits) + if not np.isfinite(actual) or actual >= limit] + if failures: + raise error_type(f"{label}: " + "; ".join(failures)) + + +@pytest.fixture(scope="module", params=( + [SELECTED_CASE] if SELECTED_CASE else ["Regular", "RegularCIC"] +)) +def paired_case(request): + wave_case = request.param + assert wave_case in REFINED_PHYSICAL_LIMITS + published = _read(wave_case, "published_clock") + refined = _read(wave_case, "refined_clock") + key = "regular" if wave_case == "Regular" else "regular_cic" + with (np.load(ORIGINAL, allow_pickle=False) as original, + np.load(REFINED, allow_pickle=False) as original_refined): + unchanged_heave = original[f"{key}_published_heave"] + unchanged_speed = original[f"{key}_published_speed"] + unchanged_wave = original[f"{key}_published_wave"] + fine_heave = original_refined[f"{key}_heave"] + fine_speed = original_refined[f"{key}_speed"] + fine_wave = original_refined[f"{key}_wave"] + for source in (published, refined): + assert _error(source["wave"], unchanged_wave) < 1e-10 + assert _error(source["wave"], fine_wave) < 1e-10 + if APPLICATIONS: + app = Path(APPLICATIONS) / "Nonlinear_Hydro" + hydro = app / "hydroData/ellipsoid.h5" + geometry = app / "geometry/elipsoid.stl" + else: + hydro = ARCHIVE / "ellipsoid.h5" + geometry = (ROOT / "tests/test_objects/test_bodyclass/testData" + / "body_10_test/elipsoid.stl") + wec = WEC("derived ellipsoid heave mass correction") + body = wec.body( + "ellipsoid", hydro, + mass="equilibrium", inertia=(1.375264e6, 1.375264e6, 1.341721e6), + geometry_file=geometry, nonlinear_hydro="instantaneous", + drag_coefficient=1, drag_area=np.pi * 25, + ) + wec.coordinate("heave", body.move("heave")) + wec.pto("PTO1", WorldPoint(0, 0, -12.5), body.at(0, 0, 0), + damping=1_200_000) + wave = RegularWave(4, 6) if wave_case == "Regular" else RegularCICWave(4, 6) + python = wec.run( + wave, dt=.05, end_time=150, ramp_time=50, rho=1025, + radiation_memory=60 if wave_case == "RegularCIC" else None, + ) + return { + "case": wave_case, "published": published, "refined": refined, + "unchanged_heave": unchanged_heave, "unchanged_speed": unchanged_speed, + "fine_heave": fine_heave, "fine_speed": fine_speed, + "position": python.bodies["ellipsoid"].position[:, 2], + "speed": python.bodies["ellipsoid"].velocity[:, 2], + "force": python.ptos["PTO1"].force, + "power": python.ptos["PTO1"].absorbed_power, + } + + +def test_derived_heave_mass_shift_and_refined_physical_gate(paired_case): + d = paired_case + case, published, refined = d["case"], d["published"], d["refined"] + refined_gaps = _gaps(d["position"], d["speed"], d["force"], + d["power"], refined) + published_gaps = _gaps(d["position"], d["speed"], d["force"], + d["power"], published) + published_solver_envelope_limits = PUBLISHED_SOLVER_ENVELOPE_LIMITS[case] + source_step_gaps = _gaps(published["heave"], published["speed"], + published["pto_force"], _source_power(published), + refined) + metrics = { + "case": case, + "corrected_mass_factor": float(published["mass"][0]), + "corrected_published_to_refined_heave_m": source_step_gaps[0], + "corrected_published_to_refined_speed_mps": source_step_gaps[1], + "corrected_published_to_refined_pto_force_n": source_step_gaps[2], + "corrected_published_to_refined_power_w": source_step_gaps[3], + "unchanged_published_to_refined_heave_m": _error( + d["unchanged_heave"], d["fine_heave"]), + "corrected_published_vs_unchanged_heave_m": _error( + published["heave"], d["unchanged_heave"]), + "corrected_refined_vs_unchanged_heave_m": _error( + refined["heave"], d["fine_heave"]), + "corrected_published_vs_unchanged_speed_mps": _error( + published["speed"], d["unchanged_speed"]), + "corrected_refined_vs_unchanged_speed_mps": _error( + refined["speed"], d["fine_speed"]), + "python_vs_corrected_refined_heave_m": refined_gaps[0], + "python_vs_corrected_refined_speed_mps": refined_gaps[1], + "python_vs_corrected_refined_pto_force_n": refined_gaps[2], + "python_vs_corrected_refined_power_w": refined_gaps[3], + "python_vs_corrected_published_heave_m": published_gaps[0], + "python_vs_corrected_published_speed_mps": published_gaps[1], + "python_vs_corrected_published_pto_force_n": published_gaps[2], + "python_vs_corrected_published_power_w": published_gaps[3], + "corrected_published_exceeds_original_solver_envelope": any( + gap >= limit for gap, limit in zip( + published_gaps, published_solver_envelope_limits)), + "corrected_nonheave_speed_max_mps": max( + float(np.max(np.abs(source["inactive_speed"]))) + for source in (published, refined)), + "corrected_nonheave_acceleration_max": max( + float(np.max(np.abs(source["inactive_acceleration"]))) + for source in (published, refined)), + } + print("derived ellipsoid mass correction:", json.dumps(metrics, sort_keys=True)) + _check_gaps(refined_gaps, REFINED_PHYSICAL_LIMITS[case], AssertionError, + "derived refined source against independent Python") + + +@pytest.mark.xfail( + strict=True, raises=CorrectedPublishedPhysicalGap, + reason="derived 0.05 s source still misses refined physical limits", +) +def test_derived_corrected_published_clock_physical_gate(paired_case): + d = paired_case + gaps = _gaps(d["position"], d["speed"], d["force"], d["power"], + d["published"]) + _check_gaps(gaps, REFINED_PHYSICAL_LIMITS[d["case"]], + CorrectedPublishedPhysicalGap, "derived published clock") + + +@pytest.mark.xfail( + strict=True, raises=CorrectedSourceStepGap, + reason="derived 0.05 s source has not passed its step-stability screen", +) +def test_derived_corrected_source_step_stability(paired_case): + d = paired_case + published, refined = d["published"], d["refined"] + gaps = _gaps(published["heave"], published["speed"], + published["pto_force"], _source_power(published), refined) + _check_gaps(gaps, SOURCE_STEP_LIMITS[d["case"]], + CorrectedSourceStepGap, "derived source 0.05 to 0.005 s") diff --git a/tools/pack_ellipsoid_ode45_mass_correction.py b/tools/pack_ellipsoid_ode45_mass_correction.py new file mode 100644 index 00000000..1b6656ee --- /dev/null +++ b/tools/pack_ellipsoid_ode45_mass_correction.py @@ -0,0 +1,82 @@ +"""Pack one pinned source job's corrected ellipsoid traces for portable CI. + +Usage: + python tools/pack_ellipsoid_ode45_mass_correction.py CASE CSV_DIR OUTPUT.npz + [--expect ARCHIVED.npz] + +CASE is Regular or RegularCIC. CSV_DIR is the case folder uploaded by the +Ellipsoid ode45 derived heave mass correction workflow. +""" + +import argparse +from pathlib import Path + +import numpy as np + + +SOURCE_SHAPES = {"body": (3001, 55), "pto": (3001, 25), + "wave": (3001, 2), "applied_added": (3001, 7), + "mass": (1, 12)} + + +def _read(folder: Path, mode: str, name: str) -> np.ndarray: + values = np.loadtxt(folder / f"{mode}_{name}.csv", delimiter=",", ndmin=2) + assert values.shape == SOURCE_SHAPES[name], (mode, name, values.shape) + assert np.isfinite(values).all(), (mode, name) + return values + + +def pack(case: str, folder: Path, output: Path, + expected: Path | None = None) -> None: + if case not in {"Regular", "RegularCIC"}: + raise ValueError(f"Unsupported wave case: {case}") + packed = {"case": np.array(case)} + common_time = np.arange(3001) * .05 + for mode, expected_step in (("published_clock", .05), + ("refined_clock", .005)): + body, pto, wave, applied, mass = ( + _read(folder, mode, name) for name in SOURCE_SHAPES + ) + for name, record in (("body", body), ("pto", pto), + ("wave", wave), ("applied_added", applied)): + np.testing.assert_allclose(record[:, 0], common_time, + rtol=0, atol=1e-10, err_msg=name) + assert abs(mass[0, 6] - expected_step) < 1e-12 + prefix = f"{mode}_" + packed.update({ + prefix + "time": common_time, + prefix + "heave": body[:, 3], + prefix + "speed": body[:, 9], + prefix + "pto_force": pto[:, 15], + prefix + "pto_speed": pto[:, 9], + prefix + "wave": wave[:, 1], + prefix + "applied_heave": applied[:, 3], + prefix + "inactive_speed": body[:, [7, 8, 10, 11, 12]], + prefix + "inactive_acceleration": body[:, [49, 50, 52, 53, 54]], + prefix + "mass": mass[0], + }) + output.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(output, **packed) + print(f"Packed {case} corrected source to {output}") + if expected is not None: + with np.load(expected, allow_pickle=False) as archive: + assert set(archive.files) == set(packed), expected + for name, values in packed.items(): + if name == "case": + assert str(archive[name]) == str(values) + else: + np.testing.assert_allclose( + values, archive[name], rtol=1e-9, atol=1e-9, + err_msg=f"{case} {name} differs from {expected}", + ) + print(f"Matched archived source traces in {expected}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("case", choices=("Regular", "RegularCIC")) + parser.add_argument("csv_dir", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--expect", type=Path) + args = parser.parse_args() + pack(args.case, args.csv_dir, args.output, args.expect)