Validate ellipsoid nonlinear heave with radiation convolution

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

@@ -62,7 +62,7 @@ on:
       model:
         description: Run all models or a selected RM3 application
         type: choice
-        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]
+        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]
         default: all
 
 jobs:
@@ -72,7 +72,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        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"]') }}
+        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"]') }}
     name: ${{ matrix.model }}
     steps:
       - uses: actions/checkout@v7
@@ -86,7 +86,7 @@ jobs:
             source
             examples/${{ matrix.model }}
       - uses: actions/checkout@v7
-        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'
+        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'
         with:
           repository: WEC-Sim/WEC-Sim_Applications
           ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e
@@ -142,8 +142,9 @@ 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'
+        if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC'
         env:
+          WEC_SIM_REFERENCE_MODEL: ${{ matrix.model }}
           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

PARITY.md

@@ -49,8 +49,9 @@ 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. |
-| 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. |
+| 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. |
+| 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)
 confirms that MATLAB's applied added-mass matrix and adjusted rigid mass sum

README.md

@@ -439,8 +439,9 @@ 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
+The Python builder supports the published `Nonlinear_Hydro/ode4/Regular` and
+`ode4/RegularCIC` ellipsoid cases 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
@@ -462,12 +463,16 @@ result = wec.run(RegularWave(height=4, period=6), dt=0.05,
                  end_time=150, ramp_time=50, rho=1025)
 ```
 
+For convolution radiation, use `RegularCICWave(height=4, period=6)` and
+`radiation_memory=60` in `wec.run`.
+
 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.
+horizontal origin in zero-direction regular waves, with either constant or
+convolution radiation.
 Other motions and sea states raise an error until their mesh force and dynamics
 checks are paired with MATLAB.
 

tests/matlab/reference_model_baseline.m

@@ -40,16 +40,21 @@ function reference_model_baseline(model)
         end
         cases = ["0m", "1m", "1m-ME", "3m", "5m"];
         caseDirs = fullfile(repoRoot, 'applications', 'Free_Decay', cases);
-    case "ELLIPSOID_NLH_REG"
+    case {"ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"}
         hydroDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', ...
             'hydroData');
         cd(hydroDir);
         if ~isfile('ellipsoid.h5')
             bemio;
         end
-        cases = "ode4_Regular";
+        if string(model) == "ELLIPSOID_NLH_REG"
+            waveDir = "Regular";
+        else
+            waveDir = "RegularCIC";
+        end
+        cases = "ode4_" + waveDir;
         caseDirs = string(fullfile(repoRoot, 'applications', ...
-            'Nonlinear_Hydro', 'ode4', 'Regular'));
+            'Nonlinear_Hydro', 'ode4', waveDir));
     case "SPHERE_MEAN_DRIFT"
         hydroDir = fullfile(repoRoot, 'applications', 'Mean_Drift', 'hydroData');
         cd(hydroDir);
@@ -403,18 +408,23 @@ 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"
+    if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"])
+        if string(model) == "ELLIPSOID_NLH_REG"
+            waveType = 'regular';
+        else
+            waveType = 'regularCIC';
+        end
         assert(simu.dt == 0.05 && simu.endTime == 150 && ...
             simu.rampTime == 50 && simu.rho == 1025 && ...
-            strcmp(waves.type, 'regular') && waves.height == 4 && ...
+            strcmp(waves.type, waveType) && 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'));
+            fullfile(outDir, string(model) + '_wave.csv'));
         writematrix([body(1).mass, body(1).inertia], ...
-            fullfile(outDir, 'ELLIPSOID_NLH_REG_mass.csv'));
+            fullfile(outDir, string(model) + '_mass.csv'));
     end
     if string(model) == "OSWEC_MULTI_WAVE"
         assert(simu.dt == 0.1 && simu.endTime == 100 && ...
@@ -608,7 +618,7 @@ function reference_model_baseline(model)
             values = [values, response.forceRadiationDamping, ...
                 response.forceAddedMass, response.forceRestoring];
         end
-        if string(model) == "ELLIPSOID_NLH_REG"
+        if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"])
             values = [values, response.forceRadiationDamping, ...
                 response.forceAddedMass, response.forceRestoring, ...
                 response.forceMorisonAndViscous, response.acceleration];