Investigate heave-only added-mass loop correction for ellipsoid ode45

.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

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.

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);