diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index b2ec0ac..c50b067 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -19,10 +19,12 @@ on: - 'tests/test_rm3_radiation_options_parity.py' - 'tests/test_rm3_pto_extension_parity.py' - 'tests/test_rm3_mcr_parity.py' + - 'tests/test_rm3_public_floating_joint_parity.py' - 'tests/test_rm3_mcr_orchestration_parity.py' - 'tests/test_rm3_mcr_seastate_parity.py' - 'tests/test_rm3_mooring_source.py' - 'tests/test_rm3_mooring_parity.py' + - 'tests/test_rm3_public_mooring_parity.py' - 'tests/test_sphere_mean_drift_source.py' - 'tests/test_sphere_mean_drift_parity.py' - 'tests/test_sphere_pto_config_parity.py' @@ -82,7 +84,7 @@ on: model: description: Run all models or a selected reference application type: choice - options: [all, SPHERE_ELEVATION_IMPORT, MORISON_FIXED, SPHERE_MOVING_MORISON, SPHERE_MOVING_MORISON_WAVE, RM3_DD_PTO, SPHERE_MPC, GBM_BARGE, SPHERE_REACTIVE_DDPTO, SPHERE_VARIABLE_MASS, 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] + options: [all, RM3_MCR_ARRAY, SPHERE_ELEVATION_IMPORT, MORISON_FIXED, SPHERE_MOVING_MORISON, SPHERE_MOVING_MORISON_WAVE, RM3_DD_PTO, SPHERE_MPC, GBM_BARGE, SPHERE_REACTIVE_DDPTO, SPHERE_VARIABLE_MASS, 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: @@ -92,7 +94,7 @@ jobs: strategy: fail-fast: false matrix: - model: ${{ fromJSON(inputs.model == 'SPHERE_ELEVATION_IMPORT' && '["SPHERE_ELEVATION_IMPORT"]' || inputs.model == 'MORISON_FIXED' && '["MORISON_FIXED"]' || inputs.model == 'SPHERE_MOVING_MORISON' && '["SPHERE_MOVING_MORISON"]' || inputs.model == 'SPHERE_MOVING_MORISON_WAVE' && '["SPHERE_MOVING_MORISON_WAVE"]' || inputs.model == 'RM3_DD_PTO' && '["RM3_DD_PTO"]' || inputs.model == 'SPHERE_MPC' && '["SPHERE_MPC"]' || inputs.model == 'GBM_BARGE' && '["GBM_BARGE"]' || inputs.model == 'SPHERE_REACTIVE_DDPTO' && '["SPHERE_REACTIVE_DDPTO"]' || inputs.model == 'SPHERE_VARIABLE_MASS' && '["SPHERE_VARIABLE_MASS"]' || 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"]' || '["SPHERE_ELEVATION_IMPORT", "MORISON_FIXED", "SPHERE_MOVING_MORISON", "SPHERE_MOVING_MORISON_WAVE", "RM3_DD_PTO", "RM3", "OSWEC", "OSWEC_Nonhydro", "GBM_BARGE", "SPHERE_VARIABLE_MASS", "SPHERE_REACTIVE_DDPTO", "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"]') }} + model: ${{ fromJSON(inputs.model == 'RM3_MCR_ARRAY' && '["RM3_MCR_ARRAY"]' || inputs.model == 'SPHERE_ELEVATION_IMPORT' && '["SPHERE_ELEVATION_IMPORT"]' || inputs.model == 'MORISON_FIXED' && '["MORISON_FIXED"]' || inputs.model == 'SPHERE_MOVING_MORISON' && '["SPHERE_MOVING_MORISON"]' || inputs.model == 'SPHERE_MOVING_MORISON_WAVE' && '["SPHERE_MOVING_MORISON_WAVE"]' || inputs.model == 'RM3_DD_PTO' && '["RM3_DD_PTO"]' || inputs.model == 'SPHERE_MPC' && '["SPHERE_MPC"]' || inputs.model == 'GBM_BARGE' && '["GBM_BARGE"]' || inputs.model == 'SPHERE_REACTIVE_DDPTO' && '["SPHERE_REACTIVE_DDPTO"]' || inputs.model == 'SPHERE_VARIABLE_MASS' && '["SPHERE_VARIABLE_MASS"]' || 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"]' || '["SPHERE_ELEVATION_IMPORT", "MORISON_FIXED", "SPHERE_MOVING_MORISON", "SPHERE_MOVING_MORISON_WAVE", "RM3_DD_PTO", "RM3", "OSWEC", "OSWEC_Nonhydro", "GBM_BARGE", "SPHERE_VARIABLE_MASS", "SPHERE_REACTIVE_DDPTO", "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 @@ -294,6 +296,11 @@ jobs: WEC_SIM_MCR_LABEL: ${{ matrix.model }} WEC_SIM_APPLICATIONS_DIR: applications WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output + - run: python -m pytest -q tests/test_rm3_public_floating_joint_parity.py + if: matrix.model == 'RM3_MCR_ARRAY' + env: + WEC_SIM_APPLICATIONS_DIR: applications + WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - run: python -m pytest -q tests/test_rm3_mcr_seastate_parity.py if: matrix.model == 'RM3_MCR_SEASTATE' env: @@ -309,6 +316,11 @@ jobs: env: WEC_SIM_APPLICATIONS_DIR: applications WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output + - run: python -m pytest -q tests/test_rm3_public_mooring_parity.py + if: matrix.model == 'RM3_MOORING_MATRIX' + env: + WEC_SIM_APPLICATIONS_DIR: applications + WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - run: python -m pytest -q tests/test_sphere_mean_drift_source.py if: matrix.model == 'SPHERE_MEAN_DRIFT' env: diff --git a/PARITY.md b/PARITY.md index a469c20..7ef188b 100644 --- a/PARITY.md +++ b/PARITY.md @@ -44,6 +44,7 @@ the production Python code. The live wave comparison passed on 6 October | RM3 PTO extension float and spar free decays | Pinned MATLAB Applications `RM3_PTO_Extension` inputs and BEMIO-generated RM3 HDF5 | The no-wave floating-joint solver initializes the float at +5 m or the spar at −5 m, reproducing each published case's +5 m PTO stroke. Over 30 s, maximum active-body heave differences are 47.2 mm for float and 6.85 mm for spar; velocity differences are 86.1 mm/s and 1.87 mm/s. Both initial poses and the passive body's heave agree within numerical precision; PTO force and power are zero in Python and within `1.7e-8` in MATLAB output. | | RM3 Multiple Condition Runs Option 1 physical sweep | Pinned MATLAB Applications RM3 input, expanded to eight scalar height × period × PTO-damping combinations | The Python `regularCIC` floating-joint solver agrees across all eight 400 s conditions within 75.1 mm surge, 0.387 mm heave, 0.000183 rad pitch, 0.413 mm PTO stroke, 0.231 mm/s PTO speed, 0.479 kN PTO force, and 0.487 kW PTO power. These scalar runs validate each physical condition; all three MCR drivers are separately checked below. | | RM3 Multiple Condition Runs Options 1–3 orchestration | Pinned MATLAB Applications arrays, Option 2 Excel grid, Option 3 MAT file, and three actual MATLAB R2025b `wecSimMCR` runs | The Python loaders produce the same ordered eight-case table for all three published inputs. Each driver has a separate paired gate for every body's 4,001 position and velocity samples, every PTO stroke, speed, force, and power trace, the eight mean powers, and two 2×2 power matrices with period rows and height columns. Python mean absorbed powers differ from MATLAB by at most 0.155 kW or 0.059% across eight conditions in the MAT driver. The Option 2 Excel and Option 3 MAT driver outputs agree numerically exactly; Option 1 array and Option 3 MAT body positions differ by at most `1.9e-12` m. MATLAB's PTO power sign is opposite Python's positive absorbed-power convention. Phase-seed sweeps, multiple PTOs, and other postprocessing are not yet paired. | +| Public RM3 floating-joint configuration | Pinned MATLAB MCR Option 1 array condition 1 and published MooringMatrix application; fresh paired [MCR](https://github.com/cmudrc/wec-sim-python/actions/runs/37787926107) and [mooring](https://github.com/cmudrc/wec-sim-python/actions/runs/37787936962) runs | `WEC.floating_joint` routes a Python-configured two-body device through the already paired pitched-slider solver and returns named body, four-coordinate, and PTO histories. In the 400 s MCR case, maximum differences are 35.1 mm surge, 0.228 mm heave, `3.46e-6` rad pitch, 0.243 mm PTO stroke, 163 N PTO force, and 118 W absorbed power; explicit gates cover positions, speeds, stroke, force, and power at all 4,001 samples. A second public configuration reproduces the 400 s imported-elevation MooringMatrix run, including its initial spar offset and joint surge spring, under the existing 40,001-sample body/PTO/mooring gates. Both fresh MATLAB jobs pass. The builder's implicit added-mass default remains physical; source-compatible delayed feedback is explicit for the MCR comparison. This layout does not represent arbitrary Simscape joints or off-axis PTO endpoints. | | RM3 imported-spectrum sea-state MCR | Pinned MATLAB Applications `RM3_MCROPT3_SeaState` and three actual MATLAB MCR runs | All three imported component tables agree to `4.9e-15`, wave elevations to `1.0e-13` m, and active excitation forces and moments to `3.8e-8` N or N m. The public `WEC.run(ImportedSpectrumWave(...))` path now replays all 4,001 incident-wave samples in each 400 s case within `1.0e-13` m; its general WEC dynamics are not claimed to reproduce the published four-coordinate RM3 motion. On MATLAB body velocities, the 60 s radiation convolution matches every logged force component within `5.4e-9` N or N m. The applied added-mass force reconstructs from the exported matrix and delayed acceleration within `7e-7` N or N m. The exact pitched-slider geometry and source-compatible mass feedback bring all three 400 s paired trajectories within 14.4 mm surge, 1.03 mm heave, and 0.000638 rad pitch; velocity differences are below 3.26 mm/s surge, 0.883 mm/s heave, and 0.000256 rad/s pitch. PTO stroke, speed, force, and absorbed power differ by at most 2.44 mm, 1.41 mm/s, 1.69 kN, and 1.19 kW. Mean absorbed power differs by at most 29 W. Explicit paired gates passed for the specialized MCR solver; PR #18 was merged. | | RM3 PTO-Sim direct linear generator | Pinned MATLAB Applications `PTO-Sim/RM3/RM3_DD_PTO`, its BEMIO-generated RM3 HDF5, and a fresh 400 s `ode4` run at 0.0005 s | Both the focused Python two-heave runner and a public `WEC.pto(linear_generator=...)` configuration integrate the source block's d/q flux states and electrical angle with float and spar motion. The public configuration uses the existing body-local endpoint and coordinate maps. Source-only checks verify PTO relative motion, friction, load voltage, and mechanical/electrical power identities. Across all 40,001 saved samples, both paired gates require body heaves/velocities and PTO stroke/speed within 10 µm or µm/s, generator force and powers within 0.005 N or W, and all three phase currents/voltages within 0.0001 A/0.01 V. The [fresh source run](https://github.com/cmudrc/wec-sim-python/actions/runs/37749845161) and local paired tests pass. This validates the published vertical two-body regular-wave layout, not general PTO-Sim electrical networks, radiation-memory coupling, or arbitrary attachment geometry. | | Sphere passive controller | Pinned MATLAB Applications `Controls/Passive (P)` and MATLAB-generated Sphere HDF5 | The Python heave model represents the published proportional controller as an 860,870 N s/m damper. Maximum differences over 400 s are 0.063 mm position, 0.065 mm/s velocity, 56 N controller force, and 32 W controller power after accounting for MATLAB's opposite force and power signs. | @@ -245,8 +246,10 @@ that body only heaves, the shifted point does not affect its trajectory; attachment-location dynamics remain unpaired. `wecsim.WEC` provides a Python builder for this same validated path and returns named NumPy body, coordinate, and PTO histories. The JSON case runner -remains available for saved cases; the Python builder currently covers the -`linear_subspace` layout only. +remains available for saved cases; the Python builder covers the mapped +`linear_subspace`, one-body `floating_gbm`, and two-body RM3-style +`floating_joint` layouts. The latter uses a relative-heave PTO and does not +accept arbitrary attachment points. `python -m wecsim CASE.json --output motion.csv` runs a supported dynamics configuration. The case declares wave, body, constraint, PTO, and time settings; the result includes all six body position and diff --git a/README.md b/README.md index c0f0180..27d9158 100644 --- a/README.md +++ b/README.md @@ -71,10 +71,36 @@ two-body device with a shared pitch pivot and body-local PTO points. Run it with `python -m examples.configurable_rm3_pto`. Relative HDF5 paths in the Python API resolve from the current directory unless `base_dir` is passed to `wec.run`. Body order must match the HDF5 hydrodynamic body order. The -Python builder covers `linear_subspace` with supported waves and the -single-body `floating_gbm` regular-wave layout. Its named coordinates use -small-motion kinematics and fixed-axis PTOs. Other validated reference -layouts remain available through the case runner and focused solver functions. +Python builder covers `linear_subspace`, the single-body `floating_gbm` +regular-wave layout, and the paired two-body `floating_joint` layout. Mapped +coordinates use small-motion kinematics and fixed-axis PTOs. The floating +joint uses its own pitched-slider geometry and a relative-heave PTO; arbitrary +PTO attachment points are not part of that reduced layout. + +For the published RM3 floating joint, configure the two bodies in HDF5 order: + +```python +from wecsim import RegularCICWave, WEC + +wec = WEC("RM3 floating joint") +float_body = wec.body("float", "path/to/rm3.h5", + inertia=(0, 21_306_090.66, 0)) +spar = wec.body("spar", "path/to/rm3.h5", + inertia=(0, 94_407_091.24, 0)) +wec.floating_joint(float_body, spar, damping=1_200_000) +result = wec.run(RegularCICWave(2.5, 8), dt=0.1, end_time=400, + ramp_time=100, radiation_memory=60) +stroke = result.ptos["relative_heave"].stroke +``` + +The default uses implicit added mass. For comparisons with the pinned MATLAB +Simulink block, `added_mass_scheme="simulink_delay"` selects its numerical +feedback setting explicitly. `floating_joint` also accepts a joint surge +spring, PTO stiffness and equilibrium, hard stops, and convolution/FIR +radiation where the underlying paired solver supports them. The returned +`absorbed_power` counts the linear damper; spring energy exchange can be +computed from `-force * velocity`. The returned stroke is float heave minus +spar heave, so an initial body offset appears in the first stroke sample. For the published generalized-body-mode barge, generate its HDF5 with the Applications `Generalized_Body_Modes/hydroData/bemio.m`, then run: diff --git a/tests/test_python_api.py b/tests/test_python_api.py index 7abb799..9b90040 100644 --- a/tests/test_python_api.py +++ b/tests/test_python_api.py @@ -164,6 +164,32 @@ def test_python_builder_runs_imported_elevation_with_relative_mat_path(): assert dict(response.raw.extra_outputs)["body1_excitation_force"].shape == (3, 6) +def test_python_builder_configures_paired_floating_joint_and_pto(): + wec = WEC("RM3 floating joint") + float_body = wec.body("float", HYDRO, inertia=(0, 21_306_090.66, 0)) + spar = wec.body("spar", HYDRO, inertia=(0, 94_407_091.24, 0)) + wec.floating_joint( + float_body, spar, damping=1_200_000, stiffness=100_000, + equilibrium_position=0.1, pto_name="main", + ) + result = wec.run(RegularWave(2.5, 8), dt=0.1, end_time=0.2, + ramp_time=1) + assert set(result.coordinates) == {"surge", "float_heave", + "spar_heave", "pitch"} + pto = result.ptos["main"] + np.testing.assert_allclose( + pto.force, + -1_200_000 * pto.velocity - 100_000 * (pto.stroke - 0.1), + ) + np.testing.assert_allclose(pto.absorbed_power, + 1_200_000 * pto.velocity**2) + np.testing.assert_allclose( + pto.stroke, + result.coordinates["float_heave"].position + - result.coordinates["spar_heave"].position, + ) + + def test_python_builder_rejects_foreign_body_attachment(): first = WEC("first") second = WEC("second") diff --git a/tests/test_rm3_public_floating_joint_parity.py b/tests/test_rm3_public_floating_joint_parity.py new file mode 100644 index 0000000..10a469d --- /dev/null +++ b/tests/test_rm3_public_floating_joint_parity.py @@ -0,0 +1,80 @@ +"""Pair the public Python floating-joint builder with published RM3 MCR.""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +from wecsim import RegularCICWave, WEC + + +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 RM3 MCR output not provided", +) + + +def _max_error(actual, expected, limit, label): + actual = np.asarray(actual) + expected = np.asarray(expected) + assert actual.shape == expected.shape, label + assert np.isfinite(actual).all() and np.isfinite(expected).all(), label + error = float(np.max(np.abs(actual - expected))) + assert error < limit, f"{label}: {error:.6g} exceeds {limit}" + + +def test_public_floating_joint_tracks_published_mcr_condition(): + root = Path(APPLICATIONS) + reference = Path(REFERENCE) + hydro = "_Common_Input_Files/RM3/hydroData/rm3.h5" + wec = WEC("RM3 MCR condition 1") + float_body = wec.body("float", hydro, + inertia=(0, 21_306_090.66, 0)) + spar = wec.body("spar", hydro, + inertia=(0, 94_407_091.24, 0)) + wec.floating_joint( + float_body, spar, damping=1_200_000, + radiation_method="convolution", + added_mass_scheme="simulink_delay", + ) + result = wec.run( + RegularCICWave(1.5, 6), dt=0.1, end_time=400, + ramp_time=100, radiation_memory=60, base_dir=root, + ) + assert result.case["constraint"]["kind"] == "floating_joint" + assert result.time.shape == (4001,) + for body_index, name in enumerate(("float", "spar"), start=1): + expected = np.loadtxt( + reference / f"RM3_MCR_ARRAY_case1_body{body_index}.csv", + delimiter=",", + ) + _max_error(result.time, expected[:, 0], 1e-10, "time") + for dof, position_limit, speed_limit in ( + (0, 0.06, 0.004), + (2, 0.002, 0.001), + (4, 0.0001, 0.00005), + ): + _max_error(result.bodies[name].position[:, dof], + expected[:, 1 + dof], position_limit, + f"{name} position DOF {dof}") + _max_error(result.bodies[name].velocity[:, dof], + expected[:, 7 + dof], speed_limit, + f"{name} velocity DOF {dof}") + expected_pto = np.loadtxt( + reference / "RM3_MCR_ARRAY_case1_pto1.csv", delimiter=",", + ) + pto = result.ptos["relative_heave"] + _max_error(pto.stroke, expected_pto[:, 3], 0.002, "PTO stroke") + _max_error(pto.velocity, expected_pto[:, 9], 0.001, "PTO speed") + _max_error(pto.force, expected_pto[:, 15], 500, "PTO force") + _max_error(pto.absorbed_power, -expected_pto[:, 21], 500, + "PTO absorbed power") + np.testing.assert_allclose( + pto.stroke, + result.coordinates["float_heave"].position + - result.coordinates["spar_heave"].position, + rtol=0, atol=1e-12, + ) diff --git a/tests/test_rm3_public_mooring_parity.py b/tests/test_rm3_public_mooring_parity.py new file mode 100644 index 0000000..758447d --- /dev/null +++ b/tests/test_rm3_public_mooring_parity.py @@ -0,0 +1,85 @@ +"""Pair a configured Python RM3 joint with the published MooringMatrix case.""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +from wecsim import ImportedElevationWave, WEC + + +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 RM3 MooringMatrix output not provided", +) + + +def _max_error(actual, expected, limit, label): + actual = np.asarray(actual) + expected = np.asarray(expected) + assert actual.shape == expected.shape, label + assert np.isfinite(actual).all() and np.isfinite(expected).all(), label + error = float(np.max(np.abs(actual - expected))) + assert error < limit, f"{label}: {error:.6g} exceeds {limit}" + + +def test_public_joint_with_imported_elevation_and_surge_mooring(): + root = Path(APPLICATIONS) + reference = Path(REFERENCE) + hydro = "_Common_Input_Files/RM3/hydroData/rm3.h5" + wec = WEC("RM3 MooringMatrix") + float_body = wec.body("float", hydro, + inertia=(0, 21_306_090.66, 0)) + spar = wec.body("spar", hydro, + inertia=(0, 94_407_091.24, 0)) + wec.floating_joint( + float_body, spar, damping=1_200_000, + mooring_surge_stiffness=100_000, + ) + result = wec.run( + ImportedElevationWave( + "Mooring/MooringMatrix/etaData.mat", + reapply_force_ramp=True, + ), + dt=0.01, end_time=400, ramp_time=40, + radiation_memory=60, + initial_coordinate={"spar_heave": -0.21}, + base_dir=root, + ) + wave = np.loadtxt(reference / "RM3_MOORING_MATRIX_wave.csv", + delimiter=",") + _max_error(result.wave_elevation, wave[:, 1], 1e-10, "wave elevation") + for index, name in enumerate(("float", "spar"), start=1): + saved = np.loadtxt( + reference / f"RM3_MOORING_MATRIX_MooringMatrix_body{index}.csv", + delimiter=",") + _max_error(result.time, saved[:, 0], 1e-8, "time") + for dof, position_limit, speed_limit in ( + (0, 0.015, 0.006), + (2, 0.004, 0.0006), + (4, 0.0009, 0.0004), + ): + _max_error(result.bodies[name].position[:, dof], + saved[:, 1 + dof], position_limit, + f"{name} position DOF {dof}") + _max_error(result.bodies[name].velocity[:, dof], + saved[:, 7 + dof], speed_limit, + f"{name} velocity DOF {dof}") + source_mooring = np.loadtxt( + reference / "RM3_MOORING_MATRIX_mooring1.csv", delimiter=",") + outputs = dict(result.raw.extra_outputs) + _max_error(outputs["mooring_surge_force"], source_mooring[:, 13], + 1_500, "mooring surge force") + source_pto = np.loadtxt( + reference / "RM3_MOORING_MATRIX_MooringMatrix_pto1.csv", + delimiter=",") + pto = result.ptos["relative_heave"] + _max_error(pto.stroke - pto.stroke[0], source_pto[:, 3], + 0.004, "PTO stroke from initial datum") + _max_error(pto.velocity, source_pto[:, 9], 0.0005, "PTO velocity") + _max_error(pto.force, source_pto[:, 15], 500, "PTO force") + _max_error(pto.absorbed_power, -source_pto[:, 21], 500, + "PTO absorbed power") diff --git a/wecsim/api.py b/wecsim/api.py index 7c74415..c445d63 100644 --- a/wecsim/api.py +++ b/wecsim/api.py @@ -7,7 +7,7 @@ from __future__ import annotations -from dataclasses import dataclass, field +from dataclasses import asdict, dataclass, field from numbers import Real from pathlib import Path from typing import Mapping, Sequence @@ -17,6 +17,7 @@ from .caseDynamics import CaseResponse, run_case from .controls import DeclutchingControl, LatchingControl from .directLinearGenerator import DirectLinearGenerator +from .hardStops import LinearHardStops from .morison import MorisonElement @@ -99,6 +100,21 @@ def move(self, dof: str, *, scale: float = 1.0, return Motion(self, dof, scale, pivot) +@dataclass(frozen=True) +class _FloatingJoint: + float_body: Body + spar_body: Body + location: WorldPoint + pto_name: str + damping: float + stiffness: float + equilibrium_position: float + mooring_surge_stiffness: float + hard_stops: LinearHardStops | None + radiation_method: str | None + added_mass_scheme: str + + @dataclass(frozen=True) class Coordinate: name: str @@ -350,6 +366,7 @@ def __init__(self, name: str = "WEC", *, body_to_body: bool = False): self.ptos: list[LinearPTO] = [] self.rotational_ptos: list[RotationalPTO] = [] self._floating_gbm_body: Body | None = None + self._floating_joint: _FloatingJoint | None = None self._morison_elements: list[tuple[Body, MorisonElement]] = [] def body(self, name: str, hydro_file: str | Path, *, @@ -447,6 +464,41 @@ def floating_gbm(self, body: Body) -> None: raise ValueError("a floating GBM body is already selected") self._floating_gbm_body = body + def floating_joint(self, float_body: Body, spar_body: Body, *, + location: WorldPoint = WorldPoint(0, 0, 0), + pto_name: str = "relative_heave", + damping: float = 0.0, stiffness: float = 0.0, + equilibrium_position: float = 0.0, + mooring_surge_stiffness: float = 0.0, + hard_stops: LinearHardStops | None = None, + radiation_method: str | None = None, + added_mass_scheme: str = "implicit") -> None: + """Select the paired two-body surge/heave/pitch slider joint. + + The PTO acts on float heave minus spar heave. The reduced joint does + not model off-axis PTO endpoints or arbitrary Simscape constraints. + """ + if (len(self.bodies) != 2 or self.bodies[0] is not float_body + or self.bodies[1] is not spar_body or self._floating_joint is not None + or self._floating_gbm_body is not None or self.coordinates + or self.ptos or self.rotational_ptos or self._morison_elements): + raise ValueError("floating_joint needs exactly two ordered bodies and no other layout or PTO") + if not isinstance(location, WorldPoint) or location.x != 0 or location.y != 0: + raise ValueError("floating_joint location needs a world point on the z axis") + if not isinstance(pto_name, str) or not pto_name: + raise ValueError("floating_joint PTO name must be nonempty") + if hard_stops is not None and not isinstance(hard_stops, LinearHardStops): + raise TypeError("hard_stops must be LinearHardStops") + if radiation_method not in (None, "constant", "convolution", "fir"): + raise ValueError("unsupported floating_joint radiation method") + if added_mass_scheme not in ("implicit", "simulink_delay"): + raise ValueError("unsupported floating_joint added-mass scheme") + self._floating_joint = _FloatingJoint( + float_body, spar_body, location, pto_name, damping, stiffness, + equilibrium_position, mooring_surge_stiffness, hard_stops, + radiation_method, added_mass_scheme, + ) + def coordinate(self, name: str, *motions: Motion) -> Coordinate: if any(existing.name == name for existing in self.coordinates): raise ValueError(f"coordinate name already exists: {name}") @@ -544,6 +596,10 @@ def to_case( ): if value is not None: simulation[key] = value + if self._floating_joint is not None: + return self._floating_joint_case( + wave, simulation, initial_coordinate, initial_speed, + ) bodies = [] for index, body in enumerate(self.bodies, start=1): if body.fixed: @@ -657,6 +713,63 @@ def to_case( for pto in self.rotational_ptos]) return case + def _floating_joint_case(self, wave, simulation, + initial_coordinate, initial_speed) -> dict: + joint = self._floating_joint + if (len(self.bodies) != 2 or self.bodies[0] is not joint.float_body + or self.bodies[1] is not joint.spar_body or self.coordinates + or self.ptos or self.rotational_ptos or self._morison_elements + or self._floating_gbm_body is not None): + raise ValueError("floating_joint cannot combine with other bodies, coordinates, or PTOs") + if not isinstance(wave, (RegularWave, RegularCICWave, + ImportedElevationWave, NoWave)): + raise ValueError("floating_joint supports regular, regularCIC, imported elevation, or no waves") + if (not np.isfinite(joint.location.coordinates()).all() + or not np.isfinite([joint.damping, joint.stiffness, + joint.equilibrium_position, + joint.mooring_surge_stiffness]).all() + or joint.mooring_surge_stiffness < 0): + raise ValueError("floating_joint location, PTO, and mooring settings must be finite") + bodies = [] + for index, body in enumerate((joint.float_body, joint.spar_body), start=1): + if (body.fixed or body.hydro_file is None or body.mass != "equilibrium" + or body.mean_drift != "none" or body.passive_yaw + or body.geometry_file is not None or body.nonlinear_hydro is not None + or body.variable_hydro is not None or body.drag_coefficient + or body.drag_area or body.passive_yaw_threshold): + raise ValueError("floating_joint needs equilibrium-mass hydrodynamic bodies without extra force models") + inertia = np.asarray(body.inertia, dtype=float) + if inertia.shape != (3,) or not np.isfinite(inertia).all() or inertia[1] <= 0: + raise ValueError("floating_joint needs a positive pitch inertia on each body") + bodies.append({ + "hydro_file": str(body.hydro_file), + "hydro_body": body.hydro_body if body.hydro_body is not None else index, + "mass": "equilibrium", "pitch_inertia": float(inertia[1]), + }) + simulation["added_mass_scheme"] = joint.added_mass_scheme + if joint.radiation_method is not None: + simulation["radiation_method"] = joint.radiation_method + constraint = {"kind": "floating_joint", + "location": joint.location.coordinates()} + if initial_coordinate is not None: + constraint["initial_coordinate"] = _state(initial_coordinate) + if initial_speed is not None: + constraint["initial_speed"] = _state(initial_speed) + pto = {"kind": "relative_heave", "damping": joint.damping, + "stiffness": joint.stiffness, + "equilibrium_position": joint.equilibrium_position} + if joint.hard_stops is not None: + pto["hard_stops"] = asdict(joint.hard_stops) + case = {"name": self.name, "simulation": simulation, + "wave": wave.as_case(), "bodies": bodies, + "constraint": constraint, "pto": pto} + if self.body_to_body: + case["body_to_body"] = True + if joint.mooring_surge_stiffness: + case["mooring"] = {"kind": "joint_surge_spring", + "stiffness": joint.mooring_surge_stiffness} + return case + def run( self, wave: RegularWave | RegularCICWave | PMWave | JONSWAPWave | ImportedSpectrumWave | ImportedElevationWave | NoWave, *, dt: float, end_time: float, @@ -682,6 +795,27 @@ def run( ) for index, body in enumerate(self.bodies) } + if self._floating_joint is not None: + coordinates = { + name: MotionHistory( + response.coordinate_position[:, index], + response.coordinate_velocity[:, index], + ) + for index, name in enumerate( + ("surge", "float_heave", "spar_heave", "pitch") + ) + } + pto = PTOHistory( + response.pto_stroke, + response.pto_velocity, + response.pto_force, + response.pto_absorbed_power, + ) + return WECResult( + response.time, bodies, coordinates, + {self._floating_joint.pto_name: pto}, + response.wave_elevation, case, response, + ) coordinates = { coordinate.name: MotionHistory( extras[f"coordinate_{coordinate.name}_position"], diff --git a/wecsim/caseDynamics.py b/wecsim/caseDynamics.py index b6329f9..e8d1372 100644 --- a/wecsim/caseDynamics.py +++ b/wecsim/caseDynamics.py @@ -61,6 +61,11 @@ class CaseResponse: auxiliary_files: tuple[Path, ...] = () pto_generalized_force: np.ndarray | None = None extra_outputs: tuple[tuple[str, np.ndarray], ...] = () + coordinate_position: np.ndarray | None = None + coordinate_velocity: np.ndarray | None = None + pto_stroke: np.ndarray | None = None + pto_velocity: np.ndarray | None = None + pto_absorbed_power: np.ndarray | None = None def _section(value, name, required, allowed): @@ -757,6 +762,11 @@ def run_case(case: Mapping, *, base_dir: str | Path = ".") -> CaseResponse: wave_elevation=elevation, auxiliary_files=auxiliary_files, extra_outputs=tuple(extra_outputs), + coordinate_position=solved.coordinate_position, + coordinate_velocity=solved.coordinate_velocity, + pto_stroke=solved.pto_stroke, + pto_velocity=solved.pto_velocity, + pto_absorbed_power=solved.pto_dissipated_power, ) raise ValueError(f"unsupported constraint layout: {kind}") diff --git a/wecsim/rm3Regular.py b/wecsim/rm3Regular.py index 2345bc4..541536d 100644 --- a/wecsim/rm3Regular.py +++ b/wecsim/rm3Regular.py @@ -25,6 +25,8 @@ class RM3RegularResponse: time: np.ndarray body_position: np.ndarray body_velocity: np.ndarray + coordinate_position: np.ndarray + coordinate_velocity: np.ndarray pto_force: np.ndarray pto_stroke: np.ndarray pto_velocity: np.ndarray @@ -323,6 +325,7 @@ def stop_generalized_force(coordinate, speed): return RM3RegularResponse( time=solved.time, body_position=solved.body_position, body_velocity=solved.body_velocity, + coordinate_position=q, coordinate_velocity=v, pto_force=pto_force, pto_stroke=pto_stroke, pto_velocity=pto_velocity, pto_mechanical_power=-pto_force * pto_velocity,