diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index e4d45a3..4bdc913 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -84,7 +84,7 @@ on: model: description: Run all models or a selected reference application type: choice - 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_VARIABLE_YAW_2DEG, 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, OSWEC, OSWEC_Nonhydro, 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_VARIABLE_YAW_2DEG, 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: @@ -94,7 +94,7 @@ jobs: strategy: fail-fast: false matrix: - 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_VARIABLE_YAW_2DEG' && '["OSWEC_VARIABLE_YAW_2DEG"]' || 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 == 'OSWEC' && '["OSWEC"]' || inputs.model == 'OSWEC_Nonhydro' && '["OSWEC_Nonhydro"]' || 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_VARIABLE_YAW_2DEG' && '["OSWEC_VARIABLE_YAW_2DEG"]' || 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 diff --git a/PARITY.md b/PARITY.md index f718b21..82b69b0 100644 --- a/PARITY.md +++ b/PARITY.md @@ -53,14 +53,14 @@ the production Python code. The live wave comparison passed on 6 October | Sphere declutching controller | Pinned MATLAB Applications `Controls/Declutching` and MATLAB-generated Sphere HDF5 | A configurable body-local PTO endpoint uses the published 232,020 N s/m engaged damping and switches off for 0.8 s after velocity reversal. In 40,001 paired samples over 400 s, heave differs by at most 1.99 mm and velocity by 4.89 mm/s. Both runs have 83 disengagement events, each 0.8 s long; event times differ by at most one 0.01 s sample, affecting four force samples. Where both runs are engaged, force differs by at most 1.14 kN and power by 1.06 kW; both report zero force while disengaged. Absorbed energy differs by 10.5 kJ over a 22.4 MJ MATLAB total. A full-sample force maximum is 106 kN because one-step mode differences occur at two transitions, so force parity is gated by event timing and within-mode values. | | Sphere latching controller | Pinned MATLAB Applications `Controls/Latching` and MATLAB-generated Sphere HDF5 | The public `LatchingControl` applies the published 49,181 N s/m normal damping and 37,308,296 N s/m latch damping for 2.4 s after velocity reversal. Both force laws dissipate energy. In 40,001 paired samples over 400 s, MATLAB and Python each have 83 latch starts; start times differ by at most 0.05 s, and each completed latch lasts 2.41 s in the sampled source implementation. Maximum heave and velocity differences are 0.379 m and 0.573 m/s; the heave extrema differ by at most 0.033 m. Absorbed energy differs by 2.43 MJ against MATLAB's 299.98 MJ (0.81%). This switched high-gain case has 553 samples in different modes, so pointwise controller force is not a useful single parity metric; the paired gate checks the exact force and power laws, event timing, motion, envelope, and integrated energy. The published linear model does not establish a feasible physical stroke or rigid latch. | | Configured Sphere spring, damper, and PTO attachment | Derived from the pinned Sphere passive case with 50,000 N/m PTO stiffness, 100,000 N s/m native PTO damping, and the PTO moved to x = 1 m | A fresh MATLAB run and the Python `WEC` builder agree within 0.051 mm position and PTO stroke, 0.054 mm/s velocity and PTO speed, 52 N combined force, and 27 W combined power. The x offset has no effect on heave-only motion, so this validates the force settings but not attachment-location dynamics. | -| OSWEC PM equal-energy waves, directional excitation, pitch, and PTO | Current MATLAB OSWEC example and current OSWEC HDF5 | Python recreates all 500 PM equal-energy bins from the HDF5 range, then synthesizes wave elevation and six-component excitation using the saved MATLAB phase matrix. Local component differences are below `6e-15`, elevation below `2e-13` m, and excitation below `1e-7` N. The fixed-hinge solver's pitch differs by at most 0.0021 rad over 400 s; PTO torque from the paired-force check differs by at most 37 N m. A Python seed produces a separate reproducible realization. This model covers pitch about a fixed hinge, not a general six-DOF device. | +| OSWEC PM equal-energy waves, directional excitation, pitch, and PTO | Current MATLAB OSWEC example and current OSWEC HDF5 | Python recreates all 500 PM equal-energy bins from the HDF5 range, then synthesizes wave elevation and six-component excitation using the saved MATLAB phase matrix. Local component differences are below `6e-15`, elevation below `2e-13` m, and excitation below `1e-7` N. The public two-body `WEC.fixed_hinge` call passed a [fresh pinned R2025b paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37795546236) with source phases: maximum flap differences over 4,001 samples were 10.0 mm surge, 1.66 mm heave, 0.00201 rad pitch, 0.00297 rad/s pitch speed, and 35.6 N m PTO torque; the fixed hydrodynamic base agreed exactly. A Python seed produces a separate reproducible realization. This model covers pitch about a fixed hinge, not a general six-DOF device. | | OSWEC `Multiple_Wave_Spectra` with two independent PM seas | Pinned MATLAB Applications `Multiple_Wave_Spectra` input and BEMIO-generated OSWEC HDF5 | The public `run_case` interface sums separate 2 m, 0° and 1 m, 90° seas with saved independent phase realizations over 100 s. Combined elevation agrees within `1e-11` m and both bodies' six excitation-force components within `1e-6` N or N m. The fixed hydrodynamic base remains stationary. With the ordinary implicit added-mass default, maximum flap differences are 106 mm surge, 10.8 mm heave, and 0.0213 rad pitch; PTO torque differs by at most 7.25 kN m. The explicit source-delay comparison reduces these to 10.3 mm, 1.20 mm, 0.00207 rad, and 818 N m. Both dynamics paths have separate paired gates. The source enables passive-yaw preprocessing, but its hinge has zero yaw motion; the source's spline frequency interpolation is selected explicitly for this case. Moving-yaw dynamics and fixed-base reactions remain unsupported. | | OSWEC `Full_Directional_Waves` imported spectrum | Pinned MATLAB Applications `Full_Directional_Waves` input and BEMIO-generated OSWEC HDF5 | The imported 47-frequency, 180-heading spectrum and realized phases reproduce all 8,001 wave-elevation samples within `9e-14` m. The source force block omits heading-bin width although its wave-elevation block includes it: logged excitation is about `1/sqrt(2° in radians) = 5.352` times the physically integrated force. Python keeps heading integration as its default. Explicit `matlab_omitted` quadrature plus `spline_frequency` interpolation matches all six logged excitation components for both bodies within `4e-7` N or N m. With that source forcing and ordinary implicit added mass, the 400 s flap differs by at most 12.8 mm surge, 3.52 mm heave, 0.00258 rad pitch, and 56.6 N m PTO torque. The separate opt-in Simulink-delay comparison reduces these to 1.53 mm, 0.293 mm, 0.000306 rad, and 5.54 N m. Source-compatible agreement does not establish that the omitted-width force is physically correct. The fixed base remains stationary; its ground reactions are not computed. | | OSWEC passive yaw off and on | Pinned MATLAB Applications `PassiveYawOFF` and `PassiveYawON` inputs, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectories | The public `WEC` interface defines a yawing flap, fixed hydrodynamic base, and 120,000 N m s/rad rotational PTO. Both published 600 s cases use 2.5 m, 8 s waves at 10°. With passive yaw off, Python and MATLAB flap yaw agree within `4e-14` rad. With passive yaw on, excitation follows the body-relative wave heading and turns the flap toward the 10° incident heading; yaw and angular-speed differences stay below 0.00149 rad and 0.000193 rad/s. PTO torque and absorbed power differ by at most 23.2 N m and 0.128 W, accounting for MATLAB's negative absorbed-power sign. On the saved MATLAB yaw trajectory, all six reconstructed excitation components differ by at most 100 N or N m except yaw moment, which differs by 467 N m. The source holds heading coefficients until yaw changes by 0.01°, while Python interpolates continuously. The fixed base is stationary; full three-dimensional joint mechanics remain unsupported. | | OSWEC variable-hydrodynamics passive yaw, 2° bank | Pinned MATLAB Applications `Variable_Hydro/Passive_Yaw` with its published regular-wave, 600 s, 0.01 s case and the first published `-40:2:40` direction bank; fresh R2025b [paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37792909630) | The Python `WEC.body(..., passive_yaw=True, yaw_heading_bank=range(-40, 41, 2))` selects the nearest body-relative BEM heading using the source rule. The selected heading matches MATLAB at all 60,001 saved MATLAB states; yaw excitation evaluated on that path differs by at most 867 N m. Before the first heading-event split at 64.57 s, yaw, yaw speed, and PTO torque differ by at most `3.89e-5` rad, `3.32e-5` rad/s, and 3.98 N m. At 64.57 s the paths differ by only `3.7e-5` rad but select different headings for one sample. Discrete events amplify that offset; over the full run the maximum yaw, speed, and PTO torque differences are 0.164 rad, 0.0221 rad/s, and 2.65 kN m. Full-trajectory parity is **unestablished**. This mode covers the published pure-yaw/fixed-base regular-wave bank, not general variable hydrodynamics. | | 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. The public `passive_yaw_threshold=1` selects the published coefficient hold and nearby BEM-heading snap. On MATLAB's saved yaw path, that opt-in force law matches all six logged excitation components within `1e-7` N or N m. The source's yaw radiation convolution reconstructs from its saved speed and pinned HDF5 within `5.4e-9` N m; added mass, hydrodynamic force balance, and rigid inertia balance close within `1.4e-8` N m. In a native Python run a heading update occurs one 0.01 s sample early at 53.45 s: MATLAB's relative heading is 10.99980° and Python's is 11.00004°, on opposite sides of the 11° update boundary. The yaw difference there is `4.3e-6` rad, and replaying MATLAB's logged force produces nearly the same pre-event drift (`4.29e-6` rad). Later event differences accumulate to 0.203 rad maximum yaw and 0.0277 rad/s speed differences over 250 s. With continuous heading interpolation, Python differs from the published trajectory by up to 0.452 rad yaw and 0.101 rad/s yaw speed. Replaying MATLAB's logged six-component force through the Python dynamics reduces these maximum differences to 0.00180 rad and 0.0000935 rad/s. The source force laws and balances are paired, but published-case trajectory parity with native excitation remains **unestablished** because the threshold amplifies small integration differences. | | 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. | +| OSWEC fixed nonhydrodynamic base | Pinned MATLAB Applications `Nonhydro_Body` case and its BEMIO-generated OSWEC HDF5 | The public `WEC.fixed_hinge(flap, base=wec.fixed_body(...))` call passed a [fresh pinned R2025b paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37795143695) and reports the stationary base and nonlinear flap motion about the PTO hinge. Across 4,001 samples over 400 s, maximum flap position differences are 20.8 mm surge, 7.9 mm heave, and 0.00444 rad pitch; velocity differences are 15.6 mm/s surge, 6.8 mm/s heave, and 0.00328 rad/s pitch. All six excitation-force components differ by less than 1 N, the base position and velocity agree exactly, and zero PTO torque differs only by MATLAB numerical noise below `3.3e-7` N m. The base's ground-constraint reaction forces are not calculated. | | Ellipsoid instantaneous nonlinear hydro, `ode4/Regular` | Pinned MATLAB Applications `Nonlinear_Hydro` and BEMIO-generated ellipsoid HDF5/STL | The Python builder combines mesh buoyancy, instantaneous Froude–Krylov correction, quadratic heave drag, BEM linear excitation/radiation, and a configured PTO. On all 3,001 saved MATLAB states, mesh equilibrium mass agrees within `2e-9` kg, buoyancy within `6e-9` N, drag within `7e-11` N, and total heave excitation within `1e-6` N. The independent 150 s Python trajectory differs by at most 5.16 mm heave, 5.86 mm/s velocity, 5.16 mm PTO stroke, 7.03 kN PTO force, and 7.21 kW absorbed power. MATLAB splits the BEM added mass between Simscape mass and an applied force; Python uses the combined effective mass. This case validates the single-body, pure-heave, zero-direction constant-radiation mode. | | Ellipsoid instantaneous nonlinear hydro, `ode4/RegularCIC` | Pinned MATLAB Applications `Nonlinear_Hydro` RegularCIC input and BEMIO-generated ellipsoid HDF5/STL | The same Python mesh model uses a 60 s radiation convolution. On all 3,001 saved MATLAB states, mesh buoyancy and quadratic drag agree within `6e-9` and `7e-11` N, linear plus nonlinear heave excitation within `4e-7` N, and radiation convolution within `3e-10` N. The independent 150 s Python trajectory differs by at most 6.39 mm heave, 7.41 mm/s velocity, 6.39 mm PTO stroke, 8.89 kN PTO force, and 8.98 kW absorbed power. The memory-step fixed-point iteration was allowed more iterations to converge near the moving waterline; no force coefficient was tuned. | | Ellipsoid `ode45/Regular` and `ode45/RegularCIC` diagnostics | Pinned MATLAB Applications `Nonlinear_Hydro` ode45 inputs and BEMIO-generated ellipsoid HDF5/STL | The physical inputs match the two ode4 cases, but MATLAB ode45 applies nonlinear buoyancy from the preceding 0.05 s output sample: its restoring log matches the mesh law on the prior body state within `6e-9` N, while the current-state mismatch exceeds 54 kN. Its logged total force and adjusted Simscape mass reconstruct the logged acceleration within `2e-9` N, confirming that this is an applied source force. Drag and wave excitation remain current-state forces; the RegularCIC radiation log differs from convolution of output-sampled velocity by at most 266 N. MATLAB ode4 versus ode45 differs by up to 15.9 mm heave and 35.4 kN PTO force with identical physical settings. Python's instantaneous-force trajectories differ from MATLAB ode45 by at most 19.94 mm heave, 30.05 mm/s velocity, 36.05 kN PTO force, and 28.78 kW absorbed power. These are solver-envelope diagnostics, **not established ode45 numerical parity**; Python does not add a one-sample buoyancy delay to reproduce the source solver artifact. | @@ -248,9 +248,11 @@ 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 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. +`linear_subspace`, one-body `floating_gbm`, two-body RM3-style +`floating_joint`, and one- or two-body OSWEC `fixed_hinge` layouts. The +floating joint uses a relative-heave PTO; the fixed hinge accepts a torsional +PTO at the published world-axis locations. Neither accepts arbitrary +attachment geometry. `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 6af8778..c64440a 100644 --- a/README.md +++ b/README.md @@ -47,7 +47,7 @@ the supported device dynamics have paired MATLAB checks. The published generalized-body-mode barge now has a coupled rigid and flexible regular-wave runner for its floating three-DOF joint. Other GBM layouts and wave options remain unsupported; see [PARITY.md](PARITY.md). The Python API is -the primary way to configure a linearized device. It constructs bodies, +the primary way to configure a supported device. It constructs bodies, named motions, attachment points, PTOs, and waves as Python objects, then returns NumPy arrays directly. Import it as `wecsim` after installation: @@ -72,7 +72,8 @@ 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`, the single-body `floating_gbm` -regular-wave layout, and the paired two-body `floating_joint` layout. Mapped +regular-wave layout, the paired two-body `floating_joint`, and the OSWEC +`fixed_hinge` 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. @@ -162,6 +163,29 @@ Supported combinations are: | `floating_joint` | `none` | Two equilibrium-mass bodies, named initial coordinates and speeds, relative-heave PTO | Radiation convolution for paired free decay; sampled FIR is also available | | `linear_subspace` | `regular`, `regularCIC`, `pm`, `jonswap`, `spectrumImport`, `elevationImport`, or `none` | Any number of six-DOF hydrodynamic bodies; named motions or 6-by-N maps; optional linear or rotational PTOs; selected mean-drift coefficients for regular waves | Constant-frequency radiation or convolution | +The canonical OSWEC pitch hinge also has a direct Python configuration: + +```python +from wecsim import PMWave, WEC, WorldPoint + +wec = WEC("OSWEC pitch flap") +flap = wec.body("flap", "oswec.h5", mass=127_000, + inertia=(1.85e6,) * 3) +base = wec.body("base", "oswec.h5", mass=999, inertia=(999,) * 3) +wec.fixed_hinge(flap, base, location=WorldPoint(0, 0, -10), + pto_location=WorldPoint(0, 0, -8.9), damping=12_000) +sea = PMWave(2.5, 8, directions=(0, 30, 90), + spreading=(0.1, 0.2, 0.7), seed=1) +result = wec.run(sea, dt=0.1, end_time=400, ramp_time=100, + radiation_memory=30) +pitch = result.coordinates["pitch"].position +torque = result.ptos["hinge"].force +``` + +The hinge can also take a `fixed_body(...)` base. Both hinge and PTO points +currently lie on the world z axis. PTO angle and speed are in radians and +radians per second; `absorbed_power` reports positive damping dissipation. + The published OSWEC passive-yaw cases use one yaw coordinate and a torsional PTO. Set `passive_yaw=True` on the moving body to interpolate excitation at the wave heading relative to its current yaw angle: diff --git a/tests/test_case_dynamics_parity.py b/tests/test_case_dynamics_parity.py index 47f4d90..fc9b808 100644 --- a/tests/test_case_dynamics_parity.py +++ b/tests/test_case_dynamics_parity.py @@ -9,7 +9,7 @@ import numpy as np import pytest -from wecsim import NoWave, RegularWave, WEC +from wecsim import NoWave, PMWave, RegularWave, WEC, WorldPoint ROOT = Path(__file__).resolve().parents[1] MODEL = os.environ.get("WEC_SIM_REFERENCE_MODEL") @@ -261,6 +261,47 @@ def test_oswec_general_runner(tmp_path): _max_error(torque * columns["body1_pitch_velocity"], pto[:, 23], max(25, 0.06 * np.max(np.abs(pto[:, 23]))), "OSWEC PTO power") + configured = WEC("Published OSWEC") + flap = configured.body("flap", hydro, mass=127_000, + inertia=(1.85e6,) * 3) + base = configured.body("base", hydro, mass=999, + inertia=(999,) * 3) + configured.fixed_hinge( + flap, base, location=WorldPoint(0, 0, -10), + pto_location=WorldPoint(0, 0, -8.9), damping=12_000, + ) + public = configured.run( + PMWave(2.5, 8, directions=tuple(directions[:, 0]), + spreading=tuple(directions[:, 1]), phase_file="phase.csv"), + dt=0.1, end_time=400, ramp_time=100, + radiation_memory=30, base_dir=tmp_path, + ) + for dof, label, position_limit, velocity_limit in ( + (0, "surge", 0.02, 0.03), + (2, "heave", 0.005, 0.005), + (4, "pitch", 0.004, 0.005), + ): + _max_error(public.bodies["flap"].position[:, dof], + expected[:, 1 + dof], position_limit, + f"configured OSWEC {label} position") + _max_error(public.bodies["flap"].velocity[:, dof], + expected[:, 7 + dof], velocity_limit, + f"configured OSWEC {label} velocity") + _max_error(public.ptos["hinge"].force, pto[:, 17], + 60, "configured OSWEC PTO torque") + _max_error(-public.ptos["hinge"].absorbed_power, pto[:, 23], + max(25, 0.06 * np.max(np.abs(pto[:, 23]))), + "configured OSWEC PTO power") + expected_base = np.loadtxt( + reference / "OSWEC_OSWEC_body2.csv", delimiter=",", + ) + np.testing.assert_allclose(public.bodies["base"].position, + expected_base[:, 1:7], rtol=0, atol=1e-10) + np.testing.assert_allclose(public.bodies["base"].velocity, + expected_base[:, 7:13], rtol=0, atol=1e-10) + np.testing.assert_allclose(public.coordinates["pitch"].position, + public.ptos["hinge"].stroke, rtol=0, atol=0) + @pytest.mark.skipif(MODEL != "Sphere" or not (SPHERE_H5 and REFERENCE), reason="paired MATLAB Sphere baseline not provided") diff --git a/tests/test_oswec_nonhydro_parity.py b/tests/test_oswec_nonhydro_parity.py index 6e66a2a..21be2ae 100644 --- a/tests/test_oswec_nonhydro_parity.py +++ b/tests/test_oswec_nonhydro_parity.py @@ -6,6 +6,7 @@ import numpy as np import pytest +from wecsim import RegularWave, WEC, WorldPoint from wecsim.caseDynamics import run_case from wecsim.hingePitch import solve_hinged_pitch_regular @@ -85,3 +86,30 @@ def test_published_nonhydrodynamic_base_against_matlab(): ) _max_error(direct.excitation_force, flap[:, 19:25], 1.0, "six-component flap excitation force") + + configured = WEC("OSWEC nonhydrodynamic base") + moving = configured.body("flap", hydro, mass=127_000, + inertia=(0, 1.85e6, 0)) + fixed = configured.fixed_body( + "base", center_gravity=(0, 0, -10.9), mass=999, + inertia=(1, 1, 1), + ) + configured.fixed_hinge( + moving, fixed, location=WorldPoint(0, 0, -10), + pto_location=WorldPoint(0, 0, -8.9), damping=0, + ) + public = configured.run( + RegularWave(2.5, 8), dt=0.1, end_time=400, ramp_time=100, + ) + np.testing.assert_allclose(public.raw.body_position, + response.body_position, rtol=0, atol=1e-12) + np.testing.assert_allclose(public.raw.body_velocity, + response.body_velocity, rtol=0, atol=1e-12) + _max_error(public.bodies["flap"].position[:, 4], flap[:, 5], + 0.007, "configured flap pitch") + np.testing.assert_allclose(public.bodies["base"].position, + base[:, 1:7], rtol=0, atol=1e-10) + np.testing.assert_allclose(public.ptos["hinge"].force, + pto[:, 17], rtol=0, atol=1e-6) + np.testing.assert_allclose(public.coordinates["pitch"].position, + public.ptos["hinge"].stroke, rtol=0, atol=0) diff --git a/tests/test_python_api.py b/tests/test_python_api.py index 9b90040..26cb36b 100644 --- a/tests/test_python_api.py +++ b/tests/test_python_api.py @@ -46,6 +46,38 @@ def test_python_builder_matches_equivalent_case(): json.dumps(case_mapping) +def test_public_fixed_hinge_reports_named_motion_and_dissipated_power(): + hydro = (ROOT / "tests/test_objects/test_bodyclass/testData/hydroData/oswec.h5") + wec = WEC("Hinged OSWEC") + flap = wec.body("flap", hydro, mass=127_000, + inertia=(0, 1.85e6, 0)) + base = wec.fixed_body("base", center_gravity=(0, 0, -10.9), + mass=999, inertia=(1, 1, 1)) + wec.fixed_hinge(flap, base, damping=12_000, stiffness=500, + equilibrium_angle=0.1) + result = wec.run(RegularWave(2.5, 8), dt=0.1, end_time=1, + ramp_time=0) + assert result.case["constraint"]["kind"] == "fixed_hinge" + assert result.case["pto"]["equilibrium_position"] == 0.1 + np.testing.assert_array_equal(result.coordinates["pitch"].position, + result.ptos["hinge"].stroke) + np.testing.assert_array_equal(result.coordinates["pitch"].velocity, + result.ptos["hinge"].velocity) + np.testing.assert_allclose( + result.ptos["hinge"].force, + -12_000 * result.ptos["hinge"].velocity + - 500 * (result.ptos["hinge"].stroke - 0.1), + rtol=0, atol=1e-10, + ) + np.testing.assert_allclose( + result.ptos["hinge"].absorbed_power, + 12_000 * result.ptos["hinge"].velocity**2, + rtol=0, atol=1e-10, + ) + np.testing.assert_array_equal(result.bodies["base"].position[:, 2], + np.full(11, -10.9)) + + def test_python_builder_supports_ground_anchor_and_named_initial_state(): wec = WEC("Anchored float") body = wec.body("float", HYDRO) diff --git a/wecsim/api.py b/wecsim/api.py index 5396410..d7a5eb0 100644 --- a/wecsim/api.py +++ b/wecsim/api.py @@ -1,4 +1,4 @@ -"""Python interface for configuring supported linearized WEC devices. +"""Python interface for configuring supported WEC devices. The objects in this module assemble the same validated case used by the JSON runner. Users can define bodies, named motions, attachment points, PTOs, and @@ -116,6 +116,19 @@ class _FloatingJoint: added_mass_scheme: str +@dataclass(frozen=True) +class _FixedHinge: + flap: Body + base: Body | None + location: WorldPoint + pto_location: WorldPoint + pto_name: str + damping: float + stiffness: float + equilibrium_angle: float + added_mass_scheme: str + + @dataclass(frozen=True) class Coordinate: name: str @@ -350,9 +363,10 @@ class WECResult: class WEC: """Build supported WEC motion/PTO layouts or fixed Morison bodies. - Rotations and PTO stroke are linearized about the reference pose. The - PTO axis stays fixed in world coordinates during a run. Hydrodynamic body - order must match the body order in the supplied HDF5 file. + Mapped-coordinate rotations and PTO stroke are linearized about the + reference pose. The fixed-hinge and floating-joint layouts use their + dedicated kinematics. Hydrodynamic body order must match the body order + in the supplied HDF5 file. """ def __init__(self, name: str = "WEC", *, body_to_body: bool = False): @@ -368,6 +382,7 @@ def __init__(self, name: str = "WEC", *, body_to_body: bool = False): self.rotational_ptos: list[RotationalPTO] = [] self._floating_gbm_body: Body | None = None self._floating_joint: _FloatingJoint | None = None + self._fixed_hinge: _FixedHinge | None = None self._morison_elements: list[tuple[Body, MorisonElement]] = [] def body(self, name: str, hydro_file: str | Path, *, @@ -492,7 +507,8 @@ def floating_joint(self, float_body: Body, spar_body: Body, *, """ 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._floating_gbm_body is not None or self._fixed_hinge 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: @@ -511,6 +527,45 @@ def floating_joint(self, float_body: Body, spar_body: Body, *, radiation_method, added_mass_scheme, ) + def fixed_hinge(self, flap: Body, base: Body | None = None, *, + location: WorldPoint = WorldPoint(0, 0, -10), + pto_location: WorldPoint = WorldPoint(0, 0, -8.9), + pto_name: str = "hinge", damping: float = 0.0, + stiffness: float = 0.0, equilibrium_angle: float = 0.0, + added_mass_scheme: str = "implicit") -> None: + """Select the published OSWEC pitch hinge and torsional PTO. + + The flap moves in surge, heave, and pitch about a fixed base. Both + joint and PTO points currently lie on the world z axis. + """ + if (len(self.bodies) != (1 if base is None else 2) + or self.bodies[0] is not flap + or base is not None and self.bodies[1] is not base + or flap.fixed or self.body_to_body + or self._fixed_hinge is not None 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("fixed_hinge needs one flap, optional ordered base, and no other layout or PTO") + if (not isinstance(location, WorldPoint) + or not isinstance(pto_location, WorldPoint) + or not np.isfinite(location.coordinates()).all() + or not np.isfinite(pto_location.coordinates()).all() + or location.x != 0 or location.y != 0 + or pto_location.x != 0 or pto_location.y != 0): + raise ValueError("fixed_hinge points need finite world locations on the z axis") + if not isinstance(pto_name, str) or not pto_name: + raise ValueError("fixed_hinge PTO name must be nonempty") + if (not np.isfinite([damping, stiffness, equilibrium_angle]).all() + or damping < 0 or stiffness < 0 + or equilibrium_angle and not stiffness): + raise ValueError("fixed_hinge PTO needs finite nonnegative damping and stiffness") + if added_mass_scheme not in ("implicit", "simulink_delay"): + raise ValueError("unsupported fixed_hinge added-mass scheme") + self._fixed_hinge = _FixedHinge( + flap, base, location, pto_location, pto_name, + damping, stiffness, equilibrium_angle, 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}") @@ -612,6 +667,10 @@ def to_case( return self._floating_joint_case( wave, simulation, initial_coordinate, initial_speed, ) + if self._fixed_hinge is not None: + return self._fixed_hinge_case( + wave, simulation, initial_coordinate, initial_speed, + ) bodies = [] for index, body in enumerate(self.bodies, start=1): if body.fixed: @@ -785,6 +844,70 @@ def _floating_joint_case(self, wave, simulation, "stiffness": joint.mooring_surge_stiffness} return case + def _fixed_hinge_case(self, wave, simulation, + initial_coordinate, initial_speed) -> dict: + hinge = self._fixed_hinge + if (len(self.bodies) != (1 if hinge.base is None else 2) + or self.bodies[0] is not hinge.flap + or hinge.base is not None and self.bodies[1] is not hinge.base + or self.coordinates or self.ptos or self.rotational_ptos + or self._morison_elements or self._floating_gbm_body is not None + or self._floating_joint is not None or self.body_to_body + or initial_coordinate is not None or initial_speed is not None): + raise ValueError("fixed_hinge needs its selected bodies and no custom coordinates or PTOs") + if not isinstance(wave, (RegularWave, PMWave)): + raise ValueError("fixed_hinge currently supports regular or PM waves") + if (isinstance(wave, RegularWave) + and (wave.direction != 0 + or hinge.added_mass_scheme != "implicit")): + raise ValueError("regular fixed_hinge needs zero-heading waves and implicit added mass") + if hinge.base is not None and hinge.base.fixed: + if (hinge.base.center_gravity is None + or not np.isfinite(hinge.base.center_gravity).all()): + raise ValueError("fixed nonhydrodynamic base needs a finite center of gravity") + base_case = { + "name": hinge.base.name, "nonhydro": True, "fixed": True, + "center_gravity": list(hinge.base.center_gravity), + "mass": hinge.base.mass, "inertia": list(hinge.base.inertia), + "volume": hinge.base.volume, + } + elif hinge.base is not None: + base_case = {"name": hinge.base.name, + "hydro_file": str(hinge.base.hydro_file), + "hydro_body": hinge.base.hydro_body or 2, + "mass": hinge.base.mass, + "pitch_inertia": float(hinge.base.inertia[1])} + else: + base_case = None + flap = hinge.flap + if (not isinstance(flap.mass, Real) or isinstance(flap.mass, bool) + or not np.isfinite(flap.mass) or flap.mass <= 0 + or len(flap.inertia) != 3 + or not np.isfinite(flap.inertia).all() + or flap.inertia[1] <= 0 + or flap.passive_yaw or flap.yaw_heading_bank is not None + or flap.variable_hydro is not None or flap.nonlinear_hydro is not None + or flap.mean_drift != "none" or flap.drag_coefficient + or flap.drag_area): + raise ValueError("fixed_hinge needs a positive flap mass and pitch inertia without extra force models") + bodies = [{"name": flap.name, "hydro_file": str(flap.hydro_file), + "hydro_body": flap.hydro_body or 1, + "mass": float(flap.mass), + "pitch_inertia": float(flap.inertia[1])}] + if base_case is not None: + bodies.append(base_case) + simulation["added_mass_scheme"] = hinge.added_mass_scheme + return { + "name": self.name, "simulation": simulation, + "wave": wave.as_case(), "bodies": bodies, + "constraint": {"kind": "fixed_hinge", + "location": hinge.location.coordinates()}, + "pto": {"kind": "pitch", "damping": hinge.damping, + "stiffness": hinge.stiffness, + "equilibrium_position": hinge.equilibrium_angle, + "location": hinge.pto_location.coordinates()}, + } + def run( self, wave: RegularWave | RegularCICWave | PMWave | JONSWAPWave | ImportedSpectrumWave | ImportedElevationWave | NoWave, *, dt: float, end_time: float, @@ -831,6 +954,19 @@ def run( {self._floating_joint.pto_name: pto}, response.wave_elevation, case, response, ) + if self._fixed_hinge is not None: + angle = response.body_position[:, 0, 4] + speed = response.body_velocity[:, 0, 4] + pto = PTOHistory( + angle, speed, response.pto_force, + self._fixed_hinge.damping * speed**2, + ) + return WECResult( + response.time, bodies, + {"pitch": MotionHistory(angle, speed)}, + {self._fixed_hinge.pto_name: pto}, + response.wave_elevation, case, response, + ) coordinates = { coordinate.name: MotionHistory( extras[f"coordinate_{coordinate.name}_position"],