Skip to content

Document ellipsoid ode45 source-sampling gap #35

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cmccomb merged 2 commits into from
Oct 8, 2026
Merged

Document ellipsoid ode45 source-sampling gap #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6c402aa
Export pinned ellipsoid ode45 baselines [skip ci]
cmccomb Oct 8, 2026
5dbf7e4
Track ellipsoid ode45 source sampling gap
cmccomb Oct 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions .github/workflows/reference-model-baselines.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ on:
model:
description: Run all models or a selected RM3 application
type: choice
options: [all, ELLIPSOID_NLH_REG, ELLIPSOID_NLH_CIC, OSWEC_PASSIVE_YAW_IRR_CONT, OSWEC_PASSIVE_YAW_IRR, OSWEC_PASSIVE_YAW, OSWEC_FULL_DIR, OSWEC_MULTI_WAVE, RM3_B2B, RM3_MOORING_MATRIX, SPHERE_MEAN_DRIFT, RM3_MCR_SEASTATE, RM3_END_STOPS, RM3_END_STOPS_STEP, RM3_END_STOPS_STEP_FINE, RM3_END_STOPS_STEP_FINER, RM3_END_STOPS_FULL_FINE, RM3_END_STOPS_FULL_FINER, RM3_END_STOPS_FULL_PAIR]
options: [all, ELLIPSOID_NLH_REG, ELLIPSOID_NLH_CIC, ELLIPSOID_NLH_ODE45, OSWEC_PASSIVE_YAW_IRR_CONT, OSWEC_PASSIVE_YAW_IRR, OSWEC_PASSIVE_YAW, OSWEC_FULL_DIR, OSWEC_MULTI_WAVE, RM3_B2B, RM3_MOORING_MATRIX, SPHERE_MEAN_DRIFT, RM3_MCR_SEASTATE, RM3_END_STOPS, RM3_END_STOPS_STEP, RM3_END_STOPS_STEP_FINE, RM3_END_STOPS_STEP_FINER, RM3_END_STOPS_FULL_FINE, RM3_END_STOPS_FULL_FINER, RM3_END_STOPS_FULL_PAIR]
default: all

jobs:
Expand All @@ -72,7 +72,7 @@ jobs:
strategy:
fail-fast: false
matrix:
model: ${{ fromJSON(inputs.model == 'ELLIPSOID_NLH_REG' && '["ELLIPSOID_NLH_REG"]' || inputs.model == 'ELLIPSOID_NLH_CIC' && '["ELLIPSOID_NLH_CIC"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' && '["OSWEC_PASSIVE_YAW_IRR_CONT"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR' && '["OSWEC_PASSIVE_YAW_IRR"]' || inputs.model == 'OSWEC_PASSIVE_YAW' && '["OSWEC_PASSIVE_YAW"]' || inputs.model == 'OSWEC_FULL_DIR' && '["OSWEC_FULL_DIR"]' || inputs.model == 'OSWEC_MULTI_WAVE' && '["OSWEC_MULTI_WAVE"]' || inputs.model == 'RM3_B2B' && '["RM3_B2B"]' || inputs.model == 'RM3_MOORING_MATRIX' && '["RM3_MOORING_MATRIX"]' || inputs.model == 'SPHERE_MEAN_DRIFT' && '["SPHERE_MEAN_DRIFT"]' || inputs.model == 'RM3_MCR_SEASTATE' && '["RM3_MCR_SEASTATE"]' || inputs.model == 'RM3_END_STOPS' && '["RM3_END_STOPS"]' || inputs.model == 'RM3_END_STOPS_STEP' && '["RM3_END_STOPS_STEP"]' || inputs.model == 'RM3_END_STOPS_STEP_FINE' && '["RM3_END_STOPS_STEP_FINE"]' || inputs.model == 'RM3_END_STOPS_STEP_FINER' && '["RM3_END_STOPS_STEP_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_FINE' && '["RM3_END_STOPS_FULL_FINE"]' || inputs.model == 'RM3_END_STOPS_FULL_FINER' && '["RM3_END_STOPS_FULL_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_PAIR' && '["RM3_END_STOPS_FULL_FINE", "RM3_END_STOPS_FULL_FINER"]' || '["RM3", "OSWEC", "OSWEC_Nonhydro", "OSWEC_FULL_DIR", "OSWEC_MULTI_WAVE", "OSWEC_PASSIVE_YAW", "Sphere", "SPHERE_MEAN_DRIFT", "RM3_B2B", "RM3_END_STOPS", "RM3_PTO_Extension", "RM3_Radiation_Options", "Sphere_Passive", "Sphere_PTO_Config", "Sphere_Reactive_PI", "Sphere_Declutching", "Sphere_Latching", "RM3_MCR", "RM3_MCR_ARRAY", "RM3_MCR_EXCEL", "RM3_MCR_MAT", "RM3_MCR_SEASTATE"]') }}
model: ${{ fromJSON(inputs.model == 'ELLIPSOID_NLH_REG' && '["ELLIPSOID_NLH_REG"]' || inputs.model == 'ELLIPSOID_NLH_CIC' && '["ELLIPSOID_NLH_CIC"]' || inputs.model == 'ELLIPSOID_NLH_ODE45' && '["ELLIPSOID_NLH_ODE45"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' && '["OSWEC_PASSIVE_YAW_IRR_CONT"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR' && '["OSWEC_PASSIVE_YAW_IRR"]' || inputs.model == 'OSWEC_PASSIVE_YAW' && '["OSWEC_PASSIVE_YAW"]' || inputs.model == 'OSWEC_FULL_DIR' && '["OSWEC_FULL_DIR"]' || inputs.model == 'OSWEC_MULTI_WAVE' && '["OSWEC_MULTI_WAVE"]' || inputs.model == 'RM3_B2B' && '["RM3_B2B"]' || inputs.model == 'RM3_MOORING_MATRIX' && '["RM3_MOORING_MATRIX"]' || inputs.model == 'SPHERE_MEAN_DRIFT' && '["SPHERE_MEAN_DRIFT"]' || inputs.model == 'RM3_MCR_SEASTATE' && '["RM3_MCR_SEASTATE"]' || inputs.model == 'RM3_END_STOPS' && '["RM3_END_STOPS"]' || inputs.model == 'RM3_END_STOPS_STEP' && '["RM3_END_STOPS_STEP"]' || inputs.model == 'RM3_END_STOPS_STEP_FINE' && '["RM3_END_STOPS_STEP_FINE"]' || inputs.model == 'RM3_END_STOPS_STEP_FINER' && '["RM3_END_STOPS_STEP_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_FINE' && '["RM3_END_STOPS_FULL_FINE"]' || inputs.model == 'RM3_END_STOPS_FULL_FINER' && '["RM3_END_STOPS_FULL_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_PAIR' && '["RM3_END_STOPS_FULL_FINE", "RM3_END_STOPS_FULL_FINER"]' || '["RM3", "OSWEC", "OSWEC_Nonhydro", "OSWEC_FULL_DIR", "OSWEC_MULTI_WAVE", "OSWEC_PASSIVE_YAW", "Sphere", "SPHERE_MEAN_DRIFT", "RM3_B2B", "RM3_END_STOPS", "RM3_PTO_Extension", "RM3_Radiation_Options", "Sphere_Passive", "Sphere_PTO_Config", "Sphere_Reactive_PI", "Sphere_Declutching", "Sphere_Latching", "RM3_MCR", "RM3_MCR_ARRAY", "RM3_MCR_EXCEL", "RM3_MCR_MAT", "RM3_MCR_SEASTATE", "ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", "ELLIPSOID_NLH_ODE45"]') }}
name: ${{ matrix.model }}
steps:
- uses: actions/checkout@v7
Expand All @@ -86,7 +86,7 @@ jobs:
source
examples/${{ matrix.model }}
- uses: actions/checkout@v7
if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'Sphere' || matrix.model == 'SPHERE_MEAN_DRIFT' || matrix.model == 'Sphere_Passive' || matrix.model == 'Sphere_PTO_Config' || matrix.model == 'Sphere_Reactive_PI' || matrix.model == 'Sphere_Declutching' || matrix.model == 'Sphere_Latching' || matrix.model == 'RM3_B2B' || matrix.model == 'RM3_MOORING_MATRIX' || matrix.model == 'RM3_END_STOPS' || matrix.model == 'RM3_END_STOPS_STEP' || matrix.model == 'RM3_END_STOPS_STEP_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINER' || matrix.model == 'RM3_END_STOPS_STEP_FINER' || matrix.model == 'RM3_PTO_Extension' || matrix.model == 'RM3_Radiation_Options' || matrix.model == 'RM3_MCR' || matrix.model == 'RM3_MCR_ARRAY' || matrix.model == 'RM3_MCR_EXCEL' || matrix.model == 'RM3_MCR_MAT' || matrix.model == 'RM3_MCR_SEASTATE' || matrix.model == 'OSWEC_Nonhydro' || matrix.model == 'OSWEC_MULTI_WAVE' || matrix.model == 'OSWEC_FULL_DIR' || matrix.model == 'OSWEC_PASSIVE_YAW' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR_CONT'
if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45' || matrix.model == 'Sphere' || matrix.model == 'SPHERE_MEAN_DRIFT' || matrix.model == 'Sphere_Passive' || matrix.model == 'Sphere_PTO_Config' || matrix.model == 'Sphere_Reactive_PI' || matrix.model == 'Sphere_Declutching' || matrix.model == 'Sphere_Latching' || matrix.model == 'RM3_B2B' || matrix.model == 'RM3_MOORING_MATRIX' || matrix.model == 'RM3_END_STOPS' || matrix.model == 'RM3_END_STOPS_STEP' || matrix.model == 'RM3_END_STOPS_STEP_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINER' || matrix.model == 'RM3_END_STOPS_STEP_FINER' || matrix.model == 'RM3_PTO_Extension' || matrix.model == 'RM3_Radiation_Options' || matrix.model == 'RM3_MCR' || matrix.model == 'RM3_MCR_ARRAY' || matrix.model == 'RM3_MCR_EXCEL' || matrix.model == 'RM3_MCR_MAT' || matrix.model == 'RM3_MCR_SEASTATE' || matrix.model == 'OSWEC_Nonhydro' || matrix.model == 'OSWEC_MULTI_WAVE' || matrix.model == 'OSWEC_FULL_DIR' || matrix.model == 'OSWEC_PASSIVE_YAW' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR_CONT'
with:
repository: WEC-Sim/WEC-Sim_Applications
ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e
Expand Down Expand Up @@ -142,7 +142,7 @@ jobs:
WEC_SIM_SPHERE_H5: applications/_Common_Input_Files/Sphere/hydroData/sphere.h5
WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output
- run: python -m pytest -q tests/test_ellipsoid_nonlinear_hydro_parity.py
if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC'
if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45'
env:
WEC_SIM_REFERENCE_MODEL: ${{ matrix.model }}
WEC_SIM_APPLICATIONS_DIR: applications
Expand Down
1 change: 1 addition & 0 deletions PARITY.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ the production Python code. The live wave comparison passed on 6 October
| OSWEC fixed nonhydrodynamic base | Pinned MATLAB Applications `Nonhydro_Body` case and its BEMIO-generated OSWEC HDF5 | The regular-wave solver reports the stationary base and nonlinear flap motion about the PTO hinge. Against a fresh 400 s MATLAB R2025b run (4,001 samples), maximum flap position differences are 20.8 mm surge, 7.9 mm heave, and 0.00444 rad pitch; velocity differences are 15.6 mm/s surge, 6.8 mm/s heave, and 0.00328 rad/s pitch. All six excitation-force components differ by less than 1 N, the base position and velocity agree exactly, and zero PTO torque differs only by MATLAB numerical noise below `3.3e-7` N m. The base's ground-constraint reaction forces are not calculated. |
| Ellipsoid instantaneous nonlinear hydro, `ode4/Regular` | Pinned MATLAB Applications `Nonlinear_Hydro` and BEMIO-generated ellipsoid HDF5/STL | The Python builder combines mesh buoyancy, instantaneous Froude–Krylov correction, quadratic heave drag, BEM linear excitation/radiation, and a configured PTO. On all 3,001 saved MATLAB states, mesh equilibrium mass agrees within `2e-9` kg, buoyancy within `6e-9` N, drag within `7e-11` N, and total heave excitation within `1e-6` N. The independent 150 s Python trajectory differs by at most 5.16 mm heave, 5.86 mm/s velocity, 5.16 mm PTO stroke, 7.03 kN PTO force, and 7.21 kW absorbed power. MATLAB splits the BEM added mass between Simscape mass and an applied force; Python uses the combined effective mass. This case validates the single-body, pure-heave, zero-direction constant-radiation mode. |
| Ellipsoid instantaneous nonlinear hydro, `ode4/RegularCIC` | Pinned MATLAB Applications `Nonlinear_Hydro` RegularCIC input and BEMIO-generated ellipsoid HDF5/STL | The same Python mesh model uses a 60 s radiation convolution. On all 3,001 saved MATLAB states, mesh buoyancy and quadratic drag agree within `6e-9` and `7e-11` N, linear plus nonlinear heave excitation within `4e-7` N, and radiation convolution within `3e-10` N. The independent 150 s Python trajectory differs by at most 6.39 mm heave, 7.41 mm/s velocity, 6.39 mm PTO stroke, 8.89 kN PTO force, and 8.98 kW absorbed power. The memory-step fixed-point iteration was allowed more iterations to converge near the moving waterline; no force coefficient was tuned. |
| Ellipsoid `ode45/Regular` and `ode45/RegularCIC` diagnostics | Pinned MATLAB Applications `Nonlinear_Hydro` ode45 inputs and BEMIO-generated ellipsoid HDF5/STL | The physical inputs match the two ode4 cases, but MATLAB ode45 applies nonlinear buoyancy from the preceding 0.05 s output sample: its restoring log matches the mesh law on the prior body state within `6e-9` N, while the current-state mismatch exceeds 54 kN. Its logged total force and adjusted Simscape mass reconstruct the logged acceleration within `2e-9` N, confirming that this is an applied source force. Drag and wave excitation remain current-state forces; the RegularCIC radiation log differs from convolution of output-sampled velocity by at most 266 N. MATLAB ode4 versus ode45 differs by up to 15.9 mm heave and 35.4 kN PTO force with identical physical settings. Python's instantaneous-force trajectories differ from MATLAB ode45 by at most 19.94 mm heave, 30.05 mm/s velocity, 36.05 kN PTO force, and 28.78 kW absorbed power. These are solver-envelope diagnostics, **not established ode45 numerical parity**; Python does not add a one-sample buoyancy delay to reproduce the source solver artifact. |
| Case-driven dynamics runner | Current MATLAB RM3 and OSWEC examples plus Sphere and RM3 Applications cases | One generalized-coordinate engine assembles rigid inertia, hydrodynamic added mass and radiation, hydrostatic restoring, excitation, and linear PTO forces. The heave, fixed-hinge, and floating-joint layouts run the paired cases above, including imported elevation and a joint surge spring for RM3 MooringMatrix. A fourth `linear_subspace` layout maps independent coordinates into arbitrary bodies; its RM3 two-body heave, Sphere free decay, configured Sphere PTO, and two instantaneous nonlinear-hydro cases have paired MATLAB checks. Arbitrary Simscape layouts, general moorings, and other nonlinear-hydro/application cases remain unsupported. |

The targeted [RM3 sea-state matrix export](https://github.com/cmudrc/wec-sim-python/actions/runs/37682612044)
Expand Down
4 changes: 4 additions & 0 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -473,6 +473,10 @@ and the configured PTO. The STL determines equilibrium mass when
Current validation covers one pure-heave body with its center of gravity at
horizontal origin in zero-direction regular waves, with either constant or
convolution radiation.
The published ode45 variants are tracked separately: MATLAB applies mesh
buoyancy from the preceding 0.05 s sample while the Python mode evaluates it
at the current state. Their motion comparisons are recorded as solver
diagnostics in [PARITY.md](PARITY.md), not as ode45 numerical parity.
Other motions and sea states raise an error until their mesh force and dynamics
checks are paired with MATLAB.

Expand Down
41 changes: 29 additions & 12 deletions tests/matlab/reference_model_baseline.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,21 +40,28 @@ function reference_model_baseline(model)
end
cases = ["0m", "1m", "1m-ME", "3m", "5m"];
caseDirs = fullfile(repoRoot, 'applications', 'Free_Decay', cases);
case {"ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"}
case {"ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", ...
"ELLIPSOID_NLH_ODE45"}
hydroDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', ...
'hydroData');
cd(hydroDir);
if ~isfile('ellipsoid.h5')
bemio;
end
if string(model) == "ELLIPSOID_NLH_REG"
waveDir = "Regular";
if string(model) == "ELLIPSOID_NLH_ODE45"
cases = ["ode45_Regular", "ode45_RegularCIC"];
caseDirs = fullfile(repoRoot, 'applications', ...
'Nonlinear_Hydro', 'ode45', ["Regular", "RegularCIC"]);
else
waveDir = "RegularCIC";
if string(model) == "ELLIPSOID_NLH_REG"
waveDir = "Regular";
else
waveDir = "RegularCIC";
end
cases = "ode4_" + waveDir;
caseDirs = string(fullfile(repoRoot, 'applications', ...
'Nonlinear_Hydro', 'ode4', waveDir));
end
cases = "ode4_" + waveDir;
caseDirs = string(fullfile(repoRoot, 'applications', ...
'Nonlinear_Hydro', 'ode4', waveDir));
case "SPHERE_MEAN_DRIFT"
hydroDir = fullfile(repoRoot, 'applications', 'Mean_Drift', 'hydroData');
cd(hydroDir);
Expand Down Expand Up @@ -408,12 +415,21 @@ function reference_model_baseline(model)
fullfile(outDir, 'OSWEC_wave_directions.csv'));
writematrix(waves.waveAmpTime, fullfile(outDir, 'OSWEC_wave_elevation.csv'));
end
if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"])
if string(model) == "ELLIPSOID_NLH_REG"
if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", ...
"ELLIPSOID_NLH_ODE45"])
if string(cases(iCase)) == "ode4_Regular" || ...
string(cases(iCase)) == "ode45_Regular"
waveType = 'regular';
else
waveType = 'regularCIC';
end
if string(model) == "ELLIPSOID_NLH_ODE45"
assert(strcmp(simu.solver, 'ode45'), ...
'The pinned nonlinear-hydro solver changed');
prefix = string(model) + "_" + cases(iCase);
else
prefix = string(model);
end
assert(simu.dt == 0.05 && simu.endTime == 150 && ...
simu.rampTime == 50 && simu.rho == 1025 && ...
strcmp(waves.type, waveType) && waves.height == 4 && ...
Expand All @@ -422,9 +438,9 @@ function reference_model_baseline(model)
isequal(constraint(1).location, [0 0 -12.5]), ...
'The pinned nonlinear-hydro input changed');
writematrix([output.wave.time(:), output.wave.elevation(:)], ...
fullfile(outDir, string(model) + '_wave.csv'));
fullfile(outDir, prefix + '_wave.csv'));
writematrix([body(1).mass, body(1).inertia], ...
fullfile(outDir, string(model) + '_mass.csv'));
fullfile(outDir, prefix + '_mass.csv'));
end
if string(model) == "OSWEC_MULTI_WAVE"
assert(simu.dt == 0.1 && simu.endTime == 100 && ...
Expand Down Expand Up @@ -618,7 +634,8 @@ function reference_model_baseline(model)
values = [values, response.forceRadiationDamping, ...
response.forceAddedMass, response.forceRestoring];
end
if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"])
if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", ...
"ELLIPSOID_NLH_ODE45"])
values = [values, response.forceRadiationDamping, ...
response.forceAddedMass, response.forceRestoring, ...
response.forceMorisonAndViscous, response.acceleration];
Expand Down
Loading
Loading