Add paired ellipsoid nonlinear heave dynamics

.github/workflows/reference-model-baselines.yml

@@ -29,6 +29,7 @@ on:
       - 'tests/test_sphere_declutching_parity.py'
       - 'tests/test_sphere_latching_parity.py'
       - 'tests/test_case_dynamics_parity.py'
+      - 'tests/test_ellipsoid_nonlinear_hydro_parity.py'
       - 'tests/test_oswec_hinge_dynamics_parity.py'
       - 'tests/test_oswec_nonhydro_parity.py'
       - 'tests/test_oswec_irregular_wave_parity.py'
@@ -51,6 +52,7 @@ on:
       - 'wecsim/controls.py'
       - 'wecsim/ptoConnections.py'
       - 'wecsim/passiveYaw.py'
+      - 'wecsim/nonlinearHydro.py'
       - 'wecsim/api.py'
       - 'wecsim/__init__.py'
       - 'wecsim/cli.py'
@@ -60,7 +62,7 @@ on:
       model:
         description: Run all models or a selected RM3 application
         type: choice
-        options: [all, 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, 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:
@@ -70,7 +72,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        model: ${{ fromJSON(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 == '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"]') }}
     name: ${{ matrix.model }}
     steps:
       - uses: actions/checkout@v7
@@ -84,7 +86,7 @@ jobs:
             source
             examples/${{ matrix.model }}
       - uses: actions/checkout@v7
-        if: 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 == '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
@@ -119,6 +121,7 @@ jobs:
             Multiple_Condition_Runs/RM3_MCROPT3_SeaState
             Mooring/MooringMatrix
             Mean_Drift
+            Nonlinear_Hydro
       - uses: matlab-actions/setup-matlab@v3
         with:
           release: R2025b
@@ -138,6 +141,11 @@ jobs:
           WEC_SIM_APPLICATIONS_DIR: applications
           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'
+        env:
+          WEC_SIM_APPLICATIONS_DIR: applications
+          WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output
       - run: python -m pytest -q tests/test_live_current_h5.py
         if: matrix.model == 'RM3' || matrix.model == 'OSWEC'
         env:
@@ -323,6 +331,8 @@ jobs:
             applications/_Common_Input_Files/RM3/hydroData/rm3.h5
             applications/_Common_Input_Files/OSWEC/hydroData/oswec.h5
             applications/Mean_Drift/hydroData/sphere.h5
+            applications/Nonlinear_Hydro/hydroData/ellipsoid.h5
+            applications/Nonlinear_Hydro/geometry/elipsoid.stl
 
   compare-end-stops-full:
     if: inputs.model == 'RM3_END_STOPS_FULL_PAIR'

PARITY.md

@@ -49,7 +49,8 @@ the production Python code. The live wave comparison passed on 6 October
 | OSWEC irregular passive yaw | Pinned MATLAB Applications `PassiveYawRegression` input, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectory | The public `PMWave` replays the published 500-frequency phase realization over 250 s with 40 s radiation memory. Frequency bins, spectral amplitudes, and widths agree within `6e-15`; elevation agrees within `1e-13` m. MATLAB holds its excitation-heading coefficients for 1° of relative yaw; reconstructing that rule on the saved MATLAB yaw path matches all six logged force components within `1e-7` N or N m. With continuous heading interpolation, Python differs from the published trajectory by up to 0.452 rad in flap yaw, 0.101 rad/s in yaw speed, 12.2 kN m in PTO torque, and 1.06 kW in source-signed power. Replaying MATLAB's logged six-component force through the Python dynamics reduces these maximum differences to 0.00180 rad, 0.0000935 rad/s, 11.3 N m, and 2.02 W. This isolates the large discrepancy to excitation heading updates without changing the Python default. Published-case trajectory parity with Python's native excitation is **not** established. |
 | OSWEC irregular passive yaw, continuous-heading control | Same pinned MATLAB input and phase realization, with only `body(1:2).yaw.threshold` changed from 1° to 0° in a temporary application copy | MATLAB and Python then agree over all 25,001 time samples: maximum differences are `1.1e-7` N or N m across the six excitation components on the MATLAB yaw path, `3.16e-5` rad flap/PTO angle, `7.30e-6` rad/s yaw/PTO speed, `0.875` N m PTO torque, and `0.0123` W source-signed power. Wave elevation differs by less than `1e-13` m and the fixed base remains stationary. This control validates the continuous-heading Python dynamics; it does not erase the published 1° threshold difference. |
 | 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. |
-| 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, and configured Sphere PTO checks agree with MATLAB. Arbitrary Simscape layouts, general moorings, nonlinear hydro, and other application cases remain unsupported. |
+| 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. Only the single-body, pure-heave, zero-direction regular-wave mode is supported so far. |
+| 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 one instantaneous nonlinear-hydro case 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)
 confirms that MATLAB's applied added-mass matrix and adjusted rigid mass sum

README.md

@@ -437,6 +437,40 @@ rotational moment arm.
 Force limits, hard stops, controllers other than the published declutching and
 reactive PI laws, and hydraulic PTO models are not implemented.
 
+### Instantaneous nonlinear hydrodynamics for heave
+
+The Python builder supports the published `Nonlinear_Hydro/ode4/Regular`
+ellipsoid through a heave-only mesh mode. Supply the BEM HDF5 file and an STL
+whose triangle coordinates are relative to the body's center of gravity:
+
+```python
+import numpy as np
+from wecsim import RegularWave, WEC, WorldPoint
+
+wec = WEC("ellipsoid")
+body = wec.body(
+    "ellipsoid", "hydroData/ellipsoid.h5",
+    geometry_file="geometry/elipsoid.stl",
+    nonlinear_hydro="instantaneous", mass="equilibrium",
+    inertia=(1_375_264, 1_375_264, 1_341_721),
+    drag_coefficient=1, drag_area=np.pi * 5**2,
+)
+wec.coordinate("heave", body.move("heave"))
+wec.pto("PTO1", WorldPoint(0, 0, -12.5), body.at(0, 0, 0),
+        damping=1_200_000)
+result = wec.run(RegularWave(height=4, period=6), dt=0.05,
+                 end_time=150, ramp_time=50, rho=1025)
+```
+
+The mode integrates mesh buoyancy, instantaneous free-surface
+Froude–Krylov correction, heave quadratic drag, BEM diffraction/radiation,
+and the configured PTO. The STL determines equilibrium mass when
+`mass="equilibrium"`; this can differ from the HDF5 displaced volume.
+Current validation covers one pure-heave body with its center of gravity at
+horizontal origin in zero-direction regular waves.
+Other motions and sea states raise an error until their mesh force and dynamics
+checks are paired with MATLAB.
+
 The published Sphere free-decay cases can also be calculated with the focused solver:
 
 ```python

tests/matlab/reference_model_baseline.m

@@ -40,6 +40,16 @@ function reference_model_baseline(model)
         end
         cases = ["0m", "1m", "1m-ME", "3m", "5m"];
         caseDirs = fullfile(repoRoot, 'applications', 'Free_Decay', cases);
+    case "ELLIPSOID_NLH_REG"
+        hydroDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', ...
+            'hydroData');
+        cd(hydroDir);
+        if ~isfile('ellipsoid.h5')
+            bemio;
+        end
+        cases = "ode4_Regular";
+        caseDirs = string(fullfile(repoRoot, 'applications', ...
+            'Nonlinear_Hydro', 'ode4', 'Regular'));
     case "SPHERE_MEAN_DRIFT"
         hydroDir = fullfile(repoRoot, 'applications', 'Mean_Drift', 'hydroData');
         cd(hydroDir);
@@ -393,6 +403,19 @@ function reference_model_baseline(model)
             fullfile(outDir, 'OSWEC_wave_directions.csv'));
         writematrix(waves.waveAmpTime, fullfile(outDir, 'OSWEC_wave_elevation.csv'));
     end
+    if string(model) == "ELLIPSOID_NLH_REG"
+        assert(simu.dt == 0.05 && simu.endTime == 150 && ...
+            simu.rampTime == 50 && simu.rho == 1025 && ...
+            strcmp(waves.type, 'regular') && waves.height == 4 && ...
+            waves.period == 6 && body(1).nonlinearHydro == 2 && ...
+            pto(1).damping == 1200000 && ...
+            isequal(constraint(1).location, [0 0 -12.5]), ...
+            'The pinned nonlinear-hydro input changed');
+        writematrix([output.wave.time(:), output.wave.elevation(:)], ...
+            fullfile(outDir, 'ELLIPSOID_NLH_REG_wave.csv'));
+        writematrix([body(1).mass, body(1).inertia], ...
+            fullfile(outDir, 'ELLIPSOID_NLH_REG_mass.csv'));
+    end
     if string(model) == "OSWEC_MULTI_WAVE"
         assert(simu.dt == 0.1 && simu.endTime == 100 && ...
             simu.rampTime == 0.1 && simu.cicEndTime == 40 && ...
@@ -585,6 +608,11 @@ function reference_model_baseline(model)
             values = [values, response.forceRadiationDamping, ...
                 response.forceAddedMass, response.forceRestoring];
         end
+        if string(model) == "ELLIPSOID_NLH_REG"
+            values = [values, response.forceRadiationDamping, ...
+                response.forceAddedMass, response.forceRestoring, ...
+                response.forceMorisonAndViscous, response.acceleration];
+        end
         assert(all(isfinite(values), 'all'), 'The MATLAB response contains nonfinite values');
         filename = sprintf('%s_%s_body%d.csv', model, cases(iCase), iBody);
         writematrix(values, fullfile(outDir, filename));