From 01b76d4d7104f9e344c6ad851ecd2878db399708 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 00:09:07 -0400 Subject: [PATCH 1/3] Export pinned ellipsoid nonlinear-hydro baseline [skip ci] --- .../workflows/reference-model-baselines.yml | 9 ++++-- tests/matlab/reference_model_baseline.m | 28 +++++++++++++++++++ 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index 2597254..9336ab6 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -60,7 +60,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 +70,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 +84,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 +119,7 @@ jobs: Multiple_Condition_Runs/RM3_MCROPT3_SeaState Mooring/MooringMatrix Mean_Drift + Nonlinear_Hydro - uses: matlab-actions/setup-matlab@v3 with: release: R2025b @@ -323,6 +324,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' diff --git a/tests/matlab/reference_model_baseline.m b/tests/matlab/reference_model_baseline.m index 9f9d090..8f64668 100644 --- a/tests/matlab/reference_model_baseline.m +++ b/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)); From f6099a813edb8351b54d0a14f5d8936e136ffd00 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 00:23:58 -0400 Subject: [PATCH 2/3] Add paired ellipsoid instantaneous nonlinear heave dynamics --- .../workflows/reference-model-baselines.yml | 7 ++ PARITY.md | 3 +- README.md | 33 ++++++ .../test_ellipsoid_nonlinear_hydro_parity.py | 109 ++++++++++++++++++ wecsim/api.py | 20 +++- wecsim/caseDynamics.py | 88 ++++++++++++-- wecsim/nonlinearHydro.py | 97 ++++++++++++++++ 7 files changed, 347 insertions(+), 10 deletions(-) create mode 100644 tests/test_ellipsoid_nonlinear_hydro_parity.py create mode 100644 wecsim/nonlinearHydro.py diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index 9336ab6..88fbbdf 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.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' @@ -139,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: diff --git a/PARITY.md b/PARITY.md index 8a94025..7d591c4 100644 --- a/PARITY.md +++ b/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 diff --git a/README.md b/README.md index a433a14..d0d0023 100644 --- a/README.md +++ b/README.md @@ -437,6 +437,39 @@ 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 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 diff --git a/tests/test_ellipsoid_nonlinear_hydro_parity.py b/tests/test_ellipsoid_nonlinear_hydro_parity.py new file mode 100644 index 0000000..8afa388 --- /dev/null +++ b/tests/test_ellipsoid_nonlinear_hydro_parity.py @@ -0,0 +1,109 @@ +"""Paired pinned Nonlinear_Hydro/ode4/Regular ellipsoid validation.""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +from wecsim.api import RegularWave, WEC, WorldPoint +from wecsim.bodyClass import BodyClass +from wecsim.nonlinearHydro import HeaveMeshHydro + + +APPLICATIONS = os.environ.get("WEC_SIM_APPLICATIONS_DIR") +REFERENCE = os.environ.get("WEC_SIM_MATLAB_MODEL_OUTPUT_DIR") +pytestmark = pytest.mark.skipif( + not (APPLICATIONS and REFERENCE), + reason="paired MATLAB nonlinear-hydro output and Applications absent", +) + + +def _source(): + source = Path(REFERENCE) + body = np.loadtxt(source / "ELLIPSOID_NLH_REG_ode4_Regular_body1.csv", + delimiter=",") + pto = np.loadtxt(source / "ELLIPSOID_NLH_REG_ode4_Regular_pto1.csv", + delimiter=",") + mass = np.loadtxt(source / "ELLIPSOID_NLH_REG_mass.csv", delimiter=",") + assert body.shape == (3001, 55) and pto.shape == (3001, 25) + return body, pto, mass + + +def _files(): + app = Path(APPLICATIONS) / "Nonlinear_Hydro" + return app / "hydroData/ellipsoid.h5", app / "geometry/elipsoid.stl" + + +def _max_error(actual, expected, limit, label): + assert actual.shape == expected.shape, label + assert np.isfinite(actual).all() and np.isfinite(expected).all(), label + error = np.max(np.abs(actual - expected)) + assert error < limit, f"{label}: {error:.6g} exceeds {limit}" + + +def test_mesh_forces_on_matlab_trajectory(): + body, _, mass = _source() + hydro_file, geometry_file = _files() + mesh = HeaveMeshHydro.from_stl( + geometry_file, center_z=-2, rho=1025, gravity=9.81, + depth=70, period=6, height=4, ramp_time=50, + mass=None, drag_coefficient=1, drag_area=np.pi * 25, + ) + _max_error(np.array([mesh.mass]), mass[:1], 1e-6, "mesh equilibrium mass") + forces = np.array([ + mesh.forces(t, z + 2, v) + for t, z, v in zip(body[:, 0], body[:, 3], body[:, 9]) + ]) + _max_error(forces[:, 0], -body[:, 39], 1e-6, + "mesh buoyancy minus weight") + _max_error(forces[:, 2], -body[:, 45], 1e-6, "quadratic drag") + + # WEC-Sim logs its linear BEM excitation plus the nonlinear FK correction. + bem = BodyClass(str(hydro_file)) + bem.bodyNumber = bem.bodyTotal = 1 + bem.readH5file() + bem.mass = mesh.mass + bem.inertia = mass[1:].tolist() + bem.hydroStiffness = np.zeros((6, 6)) + bem.viscDrag = {"Drag": np.zeros((6, 6)), "cd": np.zeros(6), + "characteristicArea": np.zeros(6)} + bem.linearDamping = np.zeros((6, 6)) + omega = 2 * np.pi / 6 + bem.hydroForcePre(omega, [0], 1, np.array([0.0]), [], .05, + 1025, 9.81, "regular", np.zeros((2, 2)), + 1, 1, 0, 0, 0) + re = np.asarray(bem.hydroForce["fExt"]["re"])[2] + im = np.asarray(bem.hydroForce["fExt"]["im"])[2] + ramp = np.array([mesh.ramp(t) for t in body[:, 0]]) + linear = 2 * ramp * (re * np.cos(omega * body[:, 0]) + - im * np.sin(omega * body[:, 0])) + _max_error(linear + forces[:, 1], body[:, 21], 1e-5, + "linear plus nonlinear Froude-Krylov excitation") + + +def test_public_python_configuration_matches_matlab_motion_and_pto(): + body, pto, mass = _source() + hydro_file, geometry_file = _files() + wec = WEC("ellipsoid") + ellipsoid = wec.body( + "ellipsoid", hydro_file, mass="equilibrium", inertia=mass[1:], + geometry_file=geometry_file, nonlinear_hydro="instantaneous", + drag_coefficient=1, drag_area=np.pi * 25, + ) + wec.coordinate("heave", ellipsoid.move("heave")) + wec.pto("PTO1", WorldPoint(0, 0, -12.5), ellipsoid.at(0, 0, 0), + damping=1_200_000) + result = wec.run(RegularWave(4, 6), dt=.05, end_time=150, + ramp_time=50, rho=1025) + np.testing.assert_allclose(result.time, body[:, 0], rtol=0, atol=1e-10) + _max_error(result.bodies["ellipsoid"].position[:, 2], body[:, 3], + .006, "heave position") + _max_error(result.bodies["ellipsoid"].velocity[:, 2], body[:, 9], + .0065, "heave velocity") + _max_error(result.ptos["PTO1"].stroke, pto[:, 3], .006, + "PTO stroke") + _max_error(result.ptos["PTO1"].force, pto[:, 15], 8_000, + "PTO force") + _max_error(result.ptos["PTO1"].absorbed_power, + -pto[:, 15] * pto[:, 9], 8_000, "PTO absorbed power") diff --git a/wecsim/api.py b/wecsim/api.py index 7292a9a..0de485d 100644 --- a/wecsim/api.py +++ b/wecsim/api.py @@ -59,6 +59,10 @@ class Body: hydro_body: int | None = None mean_drift: str = "none" passive_yaw: bool = False + geometry_file: str | Path | None = None + nonlinear_hydro: str | None = None + drag_coefficient: float = 0.0 + drag_area: float = 0.0 def at(self, x: float, y: float, z: float) -> BodyPoint: """Locate a PTO endpoint or rotation pivot relative to this body's CG.""" @@ -212,11 +216,16 @@ def body(self, name: str, hydro_file: str | Path, *, inertia: Sequence[float] = (0, 0, 0), hydro_body: int | None = None, mean_drift: str = "none", - passive_yaw: bool = False) -> Body: + passive_yaw: bool = False, + geometry_file: str | Path | None = None, + nonlinear_hydro: str | None = None, + drag_coefficient: float = 0.0, + drag_area: float = 0.0) -> Body: if any(existing.name == name for existing in self.bodies): raise ValueError(f"body name already exists: {name}") body = Body(name, hydro_file, mass, tuple(inertia), hydro_body, - mean_drift, passive_yaw) + mean_drift, passive_yaw, geometry_file, nonlinear_hydro, + drag_coefficient, drag_area) self.bodies.append(body) return body @@ -312,6 +321,13 @@ def to_case( body_case["mean_drift"] = body.mean_drift if body.passive_yaw: body_case["passive_yaw"] = True + if body.geometry_file is not None: + body_case["geometry_file"] = str(body.geometry_file) + if body.nonlinear_hydro is not None: + body_case["nonlinear_hydro"] = body.nonlinear_hydro + if body.drag_coefficient or body.drag_area: + body_case["drag_coefficient"] = body.drag_coefficient + body_case["drag_area"] = body.drag_area bodies.append(body_case) constraint = {"kind": "linear_subspace", "coordinates": []} for coordinate in self.coordinates: diff --git a/wecsim/caseDynamics.py b/wecsim/caseDynamics.py index 98d69ff..bcc99c9 100644 --- a/wecsim/caseDynamics.py +++ b/wecsim/caseDynamics.py @@ -26,6 +26,7 @@ ) from .linearCoordinates import build_coordinate_maps, initial_coordinate from .linearHeave import solve_heave_free_decay +from .nonlinearHydro import HeaveMeshHydro from .passiveYaw import PassiveYawExcitation, SampledPassiveYawExcitation from .ptoConnections import build_linear_ptos from .rm3Regular import solve_rm3_regular @@ -92,7 +93,8 @@ def _hydro_file(body, base_dir): _section(body, "body", {"hydro_file"}, {"hydro_file", "hydro_body", "mass", "pitch_inertia", "inertia", "coordinate_map", "name", "mean_drift", "fixed", - "passive_yaw"}) + "passive_yaw", "geometry_file", "nonlinear_hydro", + "drag_coefficient", "drag_area"}) raw = body["hydro_file"] if not isinstance(raw, str) or not raw: raise ValueError("body.hydro_file must be a file path") @@ -784,6 +786,28 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, "passive yaw currently needs one pure-yaw body, stationary others, " "regular or PM waves, and independent radiation" ) + nonlinear_indices = [ + index for index, spec in enumerate(bodies) + if spec.get("nonlinear_hydro") is not None + ] + if nonlinear_indices: + heave_map = np.zeros((6, 1)) + heave_map[2, 0] = 1 + if (len(bodies) != 1 or nonlinear_indices != [0] or n != 1 + or not np.array_equal(maps[0], heave_map) + or wave["type"] != "regular" or direction != 0 or b2b + or bodies[0]["nonlinear_hydro"] != "instantaneous" + or bodies[0].get("mean_drift", "none") != "none" + or bodies[0].get("passive_yaw", False)): + raise ValueError( + "instantaneous nonlinear hydro currently needs one pure-heave " + "body and a zero-direction regular wave" + ) + for spec in bodies: + if spec.get("nonlinear_hydro") is None and any( + key in spec for key in ("geometry_file", "drag_coefficient", "drag_area") + ): + raise ValueError("mesh geometry and drag settings require nonlinear_hydro") initial_q = initial_coordinate( constraint.get("initial_coordinate", [0] * n), coordinate_names, "constraint.initial_coordinate", @@ -819,6 +843,7 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, connections = () dynamic_bodies = [] + nonlinear_models = [] passive_model = None pm_elevation = None for index, (body_spec, body, mapping) in enumerate( @@ -826,6 +851,32 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, mass_setting = body_spec.get("mass", "equilibrium") if mass_setting != "equilibrium": mass_setting = _number(mass_setting, "body.mass", positive=True) + mesh_model = None + if index - 1 in nonlinear_indices: + geometry_file = body_spec.get("geometry_file") + if not isinstance(geometry_file, str) or not geometry_file: + raise ValueError("instantaneous nonlinear hydro needs geometry_file") + geometry_path = (base_dir / geometry_file).expanduser().resolve(strict=True) + cd = _number(body_spec.get("drag_coefficient", 0), + "body.drag_coefficient", nonnegative=True) + drag_area = _number(body_spec.get("drag_area", 0), + "body.drag_area", nonnegative=True) + depth = np.asarray( + body.hydroData["simulation_parameters"]["water_depth"] + ).ravel() + if depth.size != 1 or not np.isfinite(depth[0]) or depth[0] <= 0: + raise ValueError("nonlinear hydro needs positive HDF5 water depth") + mesh_model = HeaveMeshHydro.from_stl( + geometry_path, center_z=centers[index - 1][2], rho=rho, + gravity=g, depth=float(depth[0]), period=period, + height=height, ramp_time=ramp_time, + mass=None if mass_setting == "equilibrium" else mass_setting, + drag_coefficient=cd, drag_area=drag_area, + ) + if mass_setting == "equilibrium": + mass_setting = mesh_model.mass + auxiliary_files.append(geometry_path) + nonlinear_models.append(mesh_model) body.mass = mass_setting inertia = np.asarray(body_spec.get("inertia", [0, 0, 0]), dtype=float) if (inertia.shape != (3,) or not np.isfinite(inertia).all() @@ -949,16 +1000,27 @@ def excitation(at_time): def state_excitation(at_time, coordinate, speed, *, model=passive_model): return model.force(at_time, coordinate[0]) + if mesh_model is not None: + def state_excitation(at_time, coordinate, speed, *, + model=mesh_model, linear=excitation): + buoyancy, fk, drag = model.forces( + at_time, coordinate[0], speed[0] + ) + result = linear(at_time).copy() + result[2] += buoyancy + fk + drag + return result dynamic_bodies.append(DynamicBody( rigid_mass=rigid_mass, added_mass=tuple(added_mass), damping=tuple(radiation_damping), - restoring=np.asarray(hydro_force["linearHydroRestCoef"]), - static_force=np.array([ - 0, 0, (rho * float(np.asarray(body.dispVol).item()) - physical_mass) * g, - 0, 0, 0, - ]), + restoring=(np.zeros((6, 6)) if mesh_model is not None else + np.asarray(hydro_force["linearHydroRestCoef"])), + static_force=(np.zeros(6) if mesh_model is not None else + np.array([ + 0, 0, (rho * float(np.asarray(body.dispVol).item()) + - physical_mass) * g, 0, 0, 0, + ])), reference_position=np.r_[center, np.zeros(3)], motion=motion, excitation=excitation, @@ -1039,6 +1101,18 @@ def state_excitation(at_time, coordinate, speed, *, for t, angle in zip(response.time, response.coordinate[:, 0]) ]), ),) + nonlinear_outputs = () + if nonlinear_indices: + model = nonlinear_models[0] + components = np.array([ + model.forces(t, q[0], v[0]) + for t, q, v in zip(response.time, response.coordinate, response.speed) + ]) + nonlinear_outputs = ( + ("body1_buoyancy_minus_weight", components[:, 0]), + ("body1_nonlinear_fk_correction", components[:, 1]), + ("body1_quadratic_drag_force", components[:, 2]), + ) return CaseResponse( response.time, response.body_position, response.body_velocity, hydro, wave_elevation=elevation, @@ -1047,5 +1121,5 @@ def state_excitation(at_time, coordinate, speed, *, generalized_pto if "pto" in case or "ptos" in case else None ), extra_outputs=(coordinate_outputs + pto_outputs + drift_outputs - + passive_outputs), + + passive_outputs + nonlinear_outputs), ) diff --git a/wecsim/nonlinearHydro.py b/wecsim/nonlinearHydro.py new file mode 100644 index 0000000..acf5c37 --- /dev/null +++ b/wecsim/nonlinearHydro.py @@ -0,0 +1,97 @@ +"""Instantaneous free-surface forces on a triangulated heaving body. + +This implements WEC-Sim's nonlinearHydro=2 regular-wave force construction +for a body constrained to heave. Mesh triangles are in body coordinates about +the center of gravity, as required by WEC-Sim's geometry import. +""" + +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +from scipy.optimize import brentq +import trimesh + + +@dataclass(frozen=True) +class HeaveMeshHydro: + centers: np.ndarray + area_vectors: np.ndarray + center_z: float + rho: float + gravity: float + depth: float + omega: float + wave_number: float + amplitude: float + ramp_time: float + mass: float + drag_coefficient: float + drag_area: float + + @classmethod + def from_stl(cls, path: str | Path, *, center_z: float, rho: float, + gravity: float, depth: float, period: float, height: float, + ramp_time: float, mass: float | None, + drag_coefficient: float = 0, drag_area: float = 0): + mesh = trimesh.load_mesh(path, process=False) + if (not isinstance(mesh, trimesh.Trimesh) or len(mesh.faces) == 0 + or not np.isfinite(mesh.vertices).all() + or not np.isfinite(mesh.area_faces).all()): + raise ValueError("geometry_file must contain a finite triangle mesh") + centers = mesh.triangles_center + areas = mesh.face_normals * mesh.area_faces[:, None] + displaced_mass = rho * np.sum( + np.minimum(centers[:, 2] + center_z, 0) * areas[:, 2] + ) + if not np.isfinite(displaced_mass) or displaced_mass <= 0: + raise ValueError("geometry_file must produce positive equilibrium mass") + omega = 2 * np.pi / period + wave_number = brentq( + lambda k: gravity * k * np.tanh(k * depth) - omega**2, + np.finfo(float).eps, max(1.0, 2 * omega**2 / gravity), + ) + return cls(centers, areas, center_z, rho, gravity, depth, + omega, wave_number, height / 2, ramp_time, + displaced_mass if mass is None else mass, + drag_coefficient, drag_area) + + def ramp(self, at_time: float) -> float: + return (1.0 if self.ramp_time == 0 or at_time >= self.ramp_time + else (1 - np.cos(np.pi * at_time / self.ramp_time)) / 2) + + def forces(self, at_time: float, heave: float, speed: float): + """Return buoyancy-minus-weight, FK correction, and applied drag.""" + ramp = self.ramp(at_time) + eta = (self.amplitude * ramp + * np.cos(self.wave_number * self.centers[:, 0] + - self.omega * at_time)) + mean_z = self.centers[:, 2] + self.center_z + moved_z = mean_z + heave + + # WEC-Sim clips pressure above the local incident free surface. + submerged_z = np.where(moved_z > eta, 0, moved_z) + buoyancy = self.rho * self.gravity * np.sum( + submerged_z * self.area_vectors[:, 2] + ) - self.mass * self.gravity + + mean_pressure = self.rho * self.gravity * eta * ( + np.cosh(self.wave_number * (mean_z + self.depth)) + / np.cosh(self.wave_number * self.depth) + ) + mean_pressure = np.where(mean_z > 0, 0, mean_pressure) + stretched_z = (moved_z - eta) * self.depth / (self.depth + eta) + moved_pressure = self.rho * self.gravity * eta * ( + np.cosh(self.wave_number * (stretched_z + self.depth)) + / np.cosh(self.wave_number * self.depth) + ) + moved_pressure = np.where(stretched_z > 0, 0, moved_pressure) + fk_correction = -np.sum( + (moved_pressure - mean_pressure) * self.area_vectors[:, 2] + ) + # The source body block ramps its nonlinear FK output after waveClass + # has already ramped the incident elevation. + fk_correction *= ramp + drag = (-0.5 * self.rho * self.drag_coefficient + * self.drag_area * speed * abs(speed)) + return buoyancy, fk_correction, drag From 6171e0b9f963dd0459142e61ae5bce4eb5165b26 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Thu, 8 Oct 2026 00:25:10 -0400 Subject: [PATCH 3/3] Guard unpaired horizontal mesh origins --- README.md | 3 ++- wecsim/caseDynamics.py | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index d0d0023..a36f247 100644 --- a/README.md +++ b/README.md @@ -466,7 +466,8 @@ 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 in zero-direction regular waves. +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. diff --git a/wecsim/caseDynamics.py b/wecsim/caseDynamics.py index bcc99c9..587abc6 100644 --- a/wecsim/caseDynamics.py +++ b/wecsim/caseDynamics.py @@ -866,6 +866,8 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, ).ravel() if depth.size != 1 or not np.isfinite(depth[0]) or depth[0] <= 0: raise ValueError("nonlinear hydro needs positive HDF5 water depth") + if not np.allclose(centers[index - 1][:2], 0, rtol=0, atol=1e-10): + raise ValueError("heave mesh hydro currently needs CG at x=y=0") mesh_model = HeaveMeshHydro.from_stl( geometry_path, center_z=centers[index - 1][2], rho=rho, gravity=g, depth=float(depth[0]), period=period,