diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index d9a9cb6..b2ec0ac 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -54,6 +54,7 @@ on: - 'tests/test_sphere_wave_morison_parity.py' - 'tests/test_passive_yaw_configuration.py' - 'tests/test_sphere_dynamics_parity.py' + - 'tests/test_sphere_elevation_import_parity.py' - 'wecsim/linearHeave.py' - 'wecsim/hingePitch.py' - 'wecsim/irregularWave.py' @@ -81,7 +82,7 @@ on: model: description: Run all models or a selected reference application type: choice - options: [all, 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, 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: @@ -91,7 +92,7 @@ jobs: strategy: fail-fast: false matrix: - model: ${{ fromJSON(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"]' || '["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 == '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 @@ -105,7 +106,7 @@ jobs: source examples/${{ matrix.model }} - uses: actions/checkout@v7 - if: matrix.model == 'MORISON_FIXED' || matrix.model == 'SPHERE_MOVING_MORISON' || matrix.model == 'SPHERE_MOVING_MORISON_WAVE' || matrix.model == 'RM3_DD_PTO' || matrix.model == 'SPHERE_MPC' || matrix.model == 'GBM_BARGE' || matrix.model == 'SPHERE_REACTIVE_DDPTO' || matrix.model == 'SPHERE_VARIABLE_MASS' || matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45' || 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 == 'SPHERE_ELEVATION_IMPORT' || matrix.model == 'MORISON_FIXED' || matrix.model == 'SPHERE_MOVING_MORISON' || matrix.model == 'SPHERE_MOVING_MORISON_WAVE' || matrix.model == 'RM3_DD_PTO' || matrix.model == 'SPHERE_MPC' || matrix.model == 'GBM_BARGE' || matrix.model == 'SPHERE_REACTIVE_DDPTO' || matrix.model == 'SPHERE_VARIABLE_MASS' || matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45' || 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 @@ -379,6 +380,11 @@ jobs: env: 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_sphere_elevation_import_parity.py + if: matrix.model == 'SPHERE_ELEVATION_IMPORT' + env: + WEC_SIM_APPLICATIONS_DIR: applications + WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - run: python -m pytest -q tests/test_sphere_passive_parity.py if: matrix.model == 'Sphere_Passive' env: @@ -411,6 +417,7 @@ jobs: path: | matlab-reference-model-output/ applications/_Common_Input_Files/Sphere/hydroData/sphere.h5 + applications/Free_Decay/0m/etaData.mat applications/_Common_Input_Files/RM3/hydroData/rm3.h5 applications/_Common_Input_Files/OSWEC/hydroData/oswec.h5 applications/Mean_Drift/hydroData/sphere.h5 diff --git a/PARITY.md b/PARITY.md index b398218..a469c20 100644 --- a/PARITY.md +++ b/PARITY.md @@ -26,6 +26,7 @@ the production Python code. The live wave comparison passed on 6 October | RM3 body interaction force preprocessing | Original MATLAB-generated `body_4_test` through `body_9_test` files | All six regular/regularCIC variants, including body-to-body coupling on/off and state-space radiation on/off, agree for both RM3 bodies on restoring stiffness, added mass, excitation, radiation IRF, and state-space matrices where present. | | OSWEC directional irregular force preprocessing | Original MATLAB-generated `body_3_test` files | Production `BodyClass` agrees for restoring stiffness, added mass, radiation IRF, and three-direction excitation after replacing removed SciPy `interp2d`. | | Sphere `noWaveCIC` heave free decay (0 m, 1 m, 1 m with Morison elements, 3 m, 5 m) | Current MATLAB WEC-Sim and MATLAB-generated Sphere HDF5 | The focused Python linear heave solver agrees over 40 s to maximum differences of 0.45 mm position, 0.64 mm/s velocity, and 363 N total force in the 5 m case. The Morison element has only x-direction coefficients in the published 1 m case, so it does not affect heave. | +| Sphere sampled-elevation heave | Pinned `Free_Decay/0m` Sphere model with only its no-wave input replaced by a reproducible two-frequency `elevationImport` record; [fresh R2025b paired run](https://github.com/cmudrc/wec-sim-python/actions/runs/37785332545) | The public `WEC.run(ImportedElevationWave(...))` path reuses the production imported-elevation convolution and runs all 4,001 samples over 40 s. With the source body block's second force ramp explicitly selected, maximum MATLAB/Python differences are `6.7e-15` m elevation, `3.7e-9` N or N m across six excitation components, 0.0243 mm heave, and 0.0323 mm/s speed. Paired gates are `1e-12` m, `1e-6` N or N m, 0.1 mm, and 0.1 mm/s. This derived case validates one-body heave dynamics with a sampled sea; the published two-body RM3 MooringMatrix imported-elevation motion remains separately paired through its floating-joint solver. | | Fixed monopile Cartesian Morison force | Pinned MATLAB Applications `Morison_Element/morisonElement` and R2025b run, with no HDF5 bodies | The public `WEC.fixed_body` and `WEC.morison_element` configuration replays the published 400 s, three-heading PM sea and its six-component stationary-body force. All 500 equal-energy frequencies, spectral amplitudes, widths, and finite-depth wavenumbers agree within `1e-12`; elevation agrees within `1.9e-13` m. On all 40,001 samples, the largest force/moment difference is `7.3e-6` N m against a `1e-3` paired gate, and both body positions and velocities agree exactly. The source logs the negative of physical Morison force. Its `irregWaveMorison.m` uses the first random-phase column for every force heading even though wave elevation uses heading-specific phases; Python's explicit `phase_mode="matlab_shared"` reproduces that source behavior while the default uses heading-specific phases. This validates the published fixed Cartesian element, not moving elements, current profiles, or normal/tangential coefficient mode. | | Moving Cartesian Morison source force in regular waves | Pinned `regWaveMorison.m` option 1, called at eight prescribed six-DOF states during the [MORISON_FIXED MATLAB baseline](https://github.com/cmudrc/wec-sim-python/actions/runs/37768083149) | The separate `regular_morison_source_force` diagnostic matches all six MATLAB force and moment components within `5.5e-12` N or N m, including nonzero body velocity and acceleration, angular motion, wave ramp, and an emerged zero-force state. These are prescribed-state force checks, **not coupled moving-body trajectory parity**. The pinned source rotation matrix is nonorthogonal at the tested nonzero attitudes and its angular kinematics use the unrotated local point; the diagnostic reproduces those source rules without changing the default WEC dynamics. General moving Morison layouts remain unsupported in the public `WEC` runner. | | Sphere moving Morison heave free decay | The pinned `Free_Decay/1m-ME` application with its originally surge-only element changed to axial heave coefficients (`Cd=Ca=1`, area 100 m², volume 20 m³); [fresh R2025b run](https://github.com/cmudrc/wec-sim-python/actions/runs/37770321922) | The public `WEC.morison_element` couples body-local axial drag and acceleration-dependent added mass to a hydrodynamic heave body in still water. Over all 4,001 samples, independent Python motion differs by at most 0.325 mm heave and 0.952 mm/s speed. After the first second, physical Morison force differs by at most 43.2 N; the 40 s signed force-impulse difference is 198 N s. MATLAB logs its acceleration-feedback Morison force as zero at the first two samples and has a 41.2 kN pointwise difference from Python during startup; reconstructing MATLAB force from its *saved* acceleration only agrees within 14 N after 0.5 s. Python keeps implicit added mass rather than reproducing that source feedback transient. This derived case validates coupled single-body heave in still water. | diff --git a/README.md b/README.md index 3e8a86f..c0f0180 100644 --- a/README.md +++ b/README.md @@ -134,7 +134,7 @@ Supported combinations are: | `floating_joint` | `regular` or `regularCIC` | Two equilibrium-mass bodies, pitch inertias, relative-heave PTO; optional `body_to_body` | Constant-frequency radiation, impulse-response convolution, or sampled FIR radiation | | `floating_joint` | `elevationImport` | Two equilibrium-mass bodies, relative-heave PTO, optional joint surge spring | Imported MAT elevation and radiation convolution | | `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`, 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 | +| `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 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 @@ -326,6 +326,20 @@ result = wec.run(ImportedSpectrumWave("spectrumData1.mat"), The file path resolves from `base_dir`. This selects an incident sea for the configured Python WEC; the specialized RM3 floating-joint MCR runner remains the paired solver for the published four-coordinate RM3 sea-state motion. +For a sampled time/elevation MAT record, use the same builder with +`ImportedElevationWave`: + +```python +from wecsim import ImportedElevationWave + +result = wec.run(ImportedElevationWave("etaData.mat"), + dt=0.01, end_time=40, ramp_time=10, + radiation_memory=15, base_dir="path/to/inputs") +``` + +The named MAT variable defaults to `etaData` and must contain increasing time +and elevation columns. `reapply_force_ramp=True` explicitly reproduces the +second force ramp used by the pinned MATLAB body block in paired comparisons. The inherited `WaveClass` now uses the pinned MATLAB PM and JONSWAP spectrum definitions, including height-dependent PM energy and JONSWAP's inferred `gamma` when it is unspecified. Its seeded phases use a local NumPy generator; diff --git a/tests/matlab/reference_model_baseline.m b/tests/matlab/reference_model_baseline.m index 626331b..b100e57 100644 --- a/tests/matlab/reference_model_baseline.m +++ b/tests/matlab/reference_model_baseline.m @@ -41,6 +41,33 @@ function reference_model_baseline(model) end cases = ["0m", "1m", "1m-ME", "3m", "5m"]; caseDirs = fullfile(repoRoot, 'applications', 'Free_Decay', cases); + case "SPHERE_ELEVATION_IMPORT" + hydroDir = fullfile(repoRoot, 'applications', '_Common_Input_Files', ... + 'Sphere', 'hydroData'); + cd(hydroDir); + if ~isfile('sphere.h5') + bemio; + end + cases = "0m"; + caseDirs = string(fullfile(repoRoot, 'applications', 'Free_Decay', cases)); + inputFile = fullfile(caseDirs, 'wecSimInputFile.m'); + contents = fileread(inputFile); + oldWave = "waves = waveClass('noWaveCIC');"; + assert(contains(contents, oldWave), ... + 'The pinned Sphere free-decay wave input changed'); + newWave = sprintf(['waves = waveClass(''elevationImport'');\n' ... + 'waves.elevationFile = ''etaData.mat'';\n' ... + 'simu.rampTime = 10;\n' ... + 'simu.solver = ''ode4'';']); + contents = strrep(contents, char(oldWave), char(newWave)); + fid = fopen(inputFile, 'w'); + assert(fid ~= -1, 'Could not edit the Sphere elevation input'); + fprintf(fid, '%s', contents); + fclose(fid); + sampleTime = (0:0.01:40)'; + etaData = [sampleTime, 0.6*cos(2*pi*sampleTime/8) + ... + 0.2*cos(2*pi*sampleTime/3)]; + save(fullfile(caseDirs, 'etaData.mat'), 'etaData'); case {"SPHERE_MOVING_MORISON", "SPHERE_MOVING_MORISON_WAVE"} hydroDir = fullfile(repoRoot, 'applications', '_Common_Input_Files', ... 'Sphere', 'hydroData'); @@ -936,6 +963,17 @@ function reference_model_baseline(model) writematrix(stiffness, fullfile(outDir, ... 'RM3_MOORING_MATRIX_stiffness.csv')); end + if string(model) == "SPHERE_ELEVATION_IMPORT" + assert(simu.dt == 0.01 && simu.endTime == 40 && ... + simu.rampTime == 10 && simu.cicEndTime == 15 && ... + strcmp(simu.solver, 'ode4') && ... + strcmp(waves.type, 'elevationImport'), ... + 'The derived Sphere imported-elevation settings changed'); + writematrix([output.wave.time(:), output.wave.elevation(:)], ... + fullfile(outDir, 'SPHERE_ELEVATION_IMPORT_wave.csv')); + writematrix(etaData, ... + fullfile(outDir, 'SPHERE_ELEVATION_IMPORT_input.csv')); + end if string(model) == "SPHERE_MEAN_DRIFT" assert(simu.dt == 0.01 && simu.endTime == 100 && ... simu.rampTime == 20 && strcmp(waves.type, 'regularCIC') && ... diff --git a/tests/test_python_api.py b/tests/test_python_api.py index e882b06..7abb799 100644 --- a/tests/test_python_api.py +++ b/tests/test_python_api.py @@ -5,10 +5,11 @@ import numpy as np import pytest +from scipy.io import loadmat from examples.configurable_rm3_pto import HYDRO, build_wec from wecsim.caseDynamics import run_case -from wecsim import (ImportedSpectrumWave, JONSWAPWave, LatchingControl, NoWave, PMWave, +from wecsim import (ImportedElevationWave, ImportedSpectrumWave, JONSWAPWave, LatchingControl, NoWave, PMWave, RegularWave, WEC, WorldPoint) from wecsim.irregularWave import (imported_spectrum_components, jonswap_equal_energy_components, @@ -143,6 +144,26 @@ def test_python_builder_runs_imported_spectrum_with_relative_mat_path(): assert np.isfinite(response.ptos["main"].force).all() +def test_python_builder_runs_imported_elevation_with_relative_mat_path(): + source = ("tests/test_objects/test_waveclass/testData/" + "etaImport_1_test/etaData.mat") + wec = build_wec() + response = wec.run( + ImportedElevationWave(source), dt=0.1, end_time=0.2, + ramp_time=1, radiation_memory=0.2, base_dir=ROOT, + ) + raw = loadmat(ROOT / source)["etaData"] + expected = np.interp(response.time, raw[:, 0], raw[:, 1]) + ramp = (1 - np.cos(np.pi * response.time)) / 2 + np.testing.assert_allclose(response.wave_elevation, expected * ramp, + rtol=0, atol=1e-12) + assert response.raw.auxiliary_files == ((ROOT / source).resolve(),) + assert np.isfinite(response.bodies["float"].position).all() + assert np.isfinite(response.bodies["spar"].position).all() + assert np.isfinite(response.ptos["main"].force).all() + assert dict(response.raw.extra_outputs)["body1_excitation_force"].shape == (3, 6) + + def test_python_builder_rejects_foreign_body_attachment(): first = WEC("first") second = WEC("second") diff --git a/tests/test_sphere_elevation_import_parity.py b/tests/test_sphere_elevation_import_parity.py new file mode 100644 index 0000000..df13a0c --- /dev/null +++ b/tests/test_sphere_elevation_import_parity.py @@ -0,0 +1,63 @@ +"""Pair a sampled-elevation Sphere heave WEC with pinned MATLAB WEC-Sim.""" + +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="derived MATLAB Sphere imported-elevation 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}: max error {error:.6g} exceeds {limit}" + + +def test_public_sphere_imported_elevation_motion_and_force(): + root = Path(APPLICATIONS) + reference = Path(REFERENCE) + source = np.loadtxt( + reference / "SPHERE_ELEVATION_IMPORT_0m_body1.csv", delimiter=",", + ) + source_wave = np.loadtxt( + reference / "SPHERE_ELEVATION_IMPORT_wave.csv", delimiter=",", + ) + wec = WEC("Sphere with imported elevation") + sphere = wec.body( + "sphere", "_Common_Input_Files/Sphere/hydroData/sphere.h5", + inertia=(20_907_301, 21_306_090.66, 37_085_481.11), + ) + wec.coordinate("heave", sphere.move("heave")) + result = wec.run( + ImportedElevationWave( + "Free_Decay/0m/etaData.mat", reapply_force_ramp=True, + ), + dt=0.01, end_time=40, ramp_time=10, radiation_memory=15, + base_dir=root, + ) + assert source.shape == (4001, 25) + _max_error(result.time, source[:, 0], 1e-10, "time") + _max_error(result.wave_elevation, source_wave[:, 1], 1e-12, + "wave elevation") + force = dict(result.raw.extra_outputs)["body1_excitation_force"] + _max_error(force, source[:, 19:25], 1e-6, "six excitation forces") + _max_error(result.bodies["sphere"].position[:, 2], source[:, 3], 1e-4, + "heave position") + _max_error(result.bodies["sphere"].velocity[:, 2], source[:, 9], 1e-4, + "heave velocity") + assert result.raw.auxiliary_files == ( + (root / "Free_Decay/0m/etaData.mat").resolve(), + ) diff --git a/wecsim/__init__.py b/wecsim/__init__.py index 7a088a1..2d1c2d4 100644 --- a/wecsim/__init__.py +++ b/wecsim/__init__.py @@ -1,7 +1,7 @@ """Public Python interface for supported WEC-Sim device dynamics.""" from .api import ( - Body, BodyPoint, Coordinate, DirectDriveHistory, FlexibleModeHistory, HydroState, ImportedSpectrumWave, JONSWAPWave, LinearGeneratorHistory, LinearPTO, Motion, MotionHistory, NoWave, PMWave, + Body, BodyPoint, Coordinate, DirectDriveHistory, FlexibleModeHistory, HydroState, ImportedElevationWave, ImportedSpectrumWave, JONSWAPWave, LinearGeneratorHistory, LinearPTO, Motion, MotionHistory, NoWave, PMWave, RotationalPTO, PTOHistory, RegularCICWave, RegularWave, SimpleDirectDrive, VariableHydro, WEC, WECResult, WorldPoint, ) @@ -27,7 +27,7 @@ __all__ = [ "Body", "BodyPoint", "CaseResponse", "Coordinate", "DeclutchingControl", "DirectDriveHistory", "DirectLinearGenerator", "DirectLinearGeneratorSignals", "FullDirectionalComponents", "IrregularResponse", "LatchingControl", "LinearHardStops", "LinearPTO", "Motion", - "MotionHistory", "FlexibleModeHistory", "ImportedSpectrumWave", "JONSWAPWave", "MorisonElement", "NoWave", "PMWave", "PTOHistory", "LinearGeneratorHistory", "RotationalPTO", "RegularCICWave", + "MotionHistory", "FlexibleModeHistory", "ImportedElevationWave", "ImportedSpectrumWave", "JONSWAPWave", "MorisonElement", "NoWave", "PMWave", "PTOHistory", "LinearGeneratorHistory", "RotationalPTO", "RegularCICWave", "RegularWave", "SimpleDirectDrive", "VariableHydro", "HydroState", "WEC", "WECResult", "WorldPoint", "MCRCondition", "MCRPowerMatrix", "MCRResult", "MCRSeaStateResult", diff --git a/wecsim/api.py b/wecsim/api.py index 3a2993f..7c74415 100644 --- a/wecsim/api.py +++ b/wecsim/api.py @@ -243,6 +243,24 @@ def as_case(self) -> dict: return {"type": "spectrumImport", "file": str(self.file)} +@dataclass(frozen=True) +class ImportedElevationWave: + """Sampled time/elevation MAT record in WEC-Sim's ``elevationImport`` mode. + + ``file`` resolves from ``WEC.run(base_dir=...)``. The named MAT variable + contains increasing time and elevation columns in seconds and metres. + """ + + file: str | Path + variable: str = "etaData" + reapply_force_ramp: bool = False + + def as_case(self) -> dict: + return {"type": "elevationImport", "file": str(self.file), + "variable": self.variable, + "reapply_force_ramp": self.reapply_force_ramp} + + @dataclass(frozen=True) class NoWave: def as_case(self) -> dict: @@ -506,7 +524,7 @@ def rotational_pto(self, name: str, coordinate: Coordinate, *, return pto def to_case( - self, wave: RegularWave | RegularCICWave | PMWave | JONSWAPWave | ImportedSpectrumWave | NoWave, *, + self, wave: RegularWave | RegularCICWave | PMWave | JONSWAPWave | ImportedSpectrumWave | ImportedElevationWave | NoWave, *, dt: float, end_time: float, ramp_time: float | None = None, radiation_memory: float | None = None, @@ -516,8 +534,9 @@ def to_case( ) -> dict: """Return the case mapping used by the validated dynamics runner.""" if not isinstance(wave, (RegularWave, RegularCICWave, PMWave, - JONSWAPWave, ImportedSpectrumWave, NoWave)): - raise TypeError("wave must be RegularWave, RegularCICWave, PMWave, JONSWAPWave, ImportedSpectrumWave, or NoWave") + JONSWAPWave, ImportedSpectrumWave, + ImportedElevationWave, NoWave)): + raise TypeError("wave must be RegularWave, RegularCICWave, PMWave, JONSWAPWave, ImportedSpectrumWave, ImportedElevationWave, or NoWave") simulation = {"dt": dt, "end_time": end_time} for key, value in ( ("ramp_time", ramp_time), ("radiation_memory", radiation_memory), @@ -639,7 +658,7 @@ def to_case( return case def run( - self, wave: RegularWave | RegularCICWave | PMWave | JONSWAPWave | ImportedSpectrumWave | NoWave, *, + self, wave: RegularWave | RegularCICWave | PMWave | JONSWAPWave | ImportedSpectrumWave | ImportedElevationWave | NoWave, *, dt: float, end_time: float, ramp_time: float | None = None, radiation_memory: float | None = None, diff --git a/wecsim/caseDynamics.py b/wecsim/caseDynamics.py index ca1111d..b6329f9 100644 --- a/wecsim/caseDynamics.py +++ b/wecsim/caseDynamics.py @@ -193,8 +193,8 @@ def _hard_stops(spec): return LinearHardStops(**values) -def _imported_elevation(wave, base, hydro_file, time, dt, ramp_time, rho, g): - """Build RM3 body forces from one sampled MATLAB elevation record.""" +def _sampled_elevation(wave, base, time, ramp_time): + """Load and ramp a sampled MATLAB elevation record.""" _section(wave, "imported wave", {"type", "file"}, {"type", "file", "variable", "direction", "reapply_force_ramp"}) raw = wave["file"] @@ -208,7 +208,7 @@ def _imported_elevation(wave, base, hydro_file, time, dt, ramp_time, rho, g): raise ValueError("wave.variable must be a nonempty MAT variable name") direction = _number(wave.get("direction", 0), "wave.direction") if direction != 0: - raise ValueError("floating-joint imported elevation currently supports 0-degree waves") + raise ValueError("imported elevation currently supports 0-degree waves") second_ramp = wave.get("reapply_force_ramp", False) if not isinstance(second_ramp, bool): raise ValueError("wave.reapply_force_ramp must be a boolean") @@ -230,6 +230,13 @@ def _imported_elevation(wave, base, hydro_file, time, dt, ramp_time, rho, g): early = time < ramp_time ramp[early] = (1 - np.cos(np.pi * time[early] / ramp_time)) / 2 elevation = np.interp(time, samples[:, 0], samples[:, 1]) * ramp + return elevation, ramp, path + + +def _imported_elevation(wave, base, hydro_file, time, dt, ramp_time, rho, g): + """Build RM3 body forces from one sampled MATLAB elevation record.""" + elevation, ramp, path = _sampled_elevation(wave, base, time, ramp_time) + direction = _number(wave.get("direction", 0), "wave.direction") force = np.zeros((len(time), 2, 6)) for number in (1, 2): body = BodyClass(str(hydro_file)) @@ -239,7 +246,7 @@ def _imported_elevation(wave, base, hydro_file, time, dt, ramp_time, rho, g): body.hydroForce["userDefinedFe"] = np.zeros((len(time), 6)) body.userDefinedExcitation(np.vstack((time, elevation)), dt, [direction], rho, g) force[:, number - 1] = body.hydroForce["userDefinedFe"] - if second_ramp: + if wave.get("reapply_force_ramp", False): # The pinned MATLAB body block ramps force after waveClass ramped elevation. force *= ramp[:, None, None] return elevation, force, path @@ -1107,18 +1114,18 @@ def _run_fixed_morison(case, sim, wave, constraint, bodies, hydro, def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, b2b, dt, end_time, ramp_time, rho, g, base_dir): - """Run mapped coordinates with regular, PM, or no incident waves.""" + """Run mapped coordinates with supported regular, sampled, or no waves.""" if set(constraint) - {"kind", "initial_coordinate", "initial_speed", "coordinates"}: raise ValueError("linear_subspace uses coordinate maps, not joint locations") if wave["type"] not in ("regular", "regularCIC", "pm", "jonswap", - "spectrumImport", "none"): - raise ValueError("linear_subspace supports regular, regularCIC, PM, JONSWAP, spectrumImport, or no waves") + "spectrumImport", "elevationImport", "none"): + raise ValueError("linear_subspace supports regular, regularCIC, PM, JONSWAP, spectrumImport, elevationImport, or no waves") if b2b and len(set(hydro)) != 1: raise ValueError("body-to-body hydrodynamics need one shared HDF5 file") if wave["type"] == "none" and set(wave) != {"type"}: raise ValueError("no-wave cases have no wave height or period") - if (wave["type"] in ("pm", "jonswap", "spectrumImport") + if (wave["type"] in ("pm", "jonswap", "spectrumImport", "elevationImport") and any(body.get("mean_drift", "none") != "none" for body in bodies)): raise ValueError("irregular linear-subspace mean-drift forcing is not supported") components = None @@ -1174,6 +1181,8 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, spectrum_file = (base_dir / wave["file"]).expanduser().resolve(strict=True) components = imported_spectrum_components(hydro[0], spectrum_file) auxiliary_files.append(spectrum_file) + elif wave["type"] == "elevationImport": + pass # The sampled record is loaded after the simulation grid is built. else: if "ramp_time" in sim: raise ValueError("ramp_time is inapplicable to no-wave dynamics") @@ -1189,6 +1198,13 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, time = np.arange(round(end_time / dt) + 1) * dt if not np.isclose(time[-1], end_time, atol=1e-10): raise ValueError("end_time must be an integer multiple of dt") + imported_elevation = None + imported_ramp = None + if wave["type"] == "elevationImport": + imported_elevation, imported_ramp, wave_path = _sampled_elevation( + wave, base_dir, time, ramp_time, + ) + auxiliary_files.append(wave_path) body_names = [] for index, body_spec in enumerate(bodies, start=1): _body_number(body_spec, index) @@ -1440,7 +1456,10 @@ def moving_terms(at_time, coordinate, speed): "characteristicArea": np.zeros(6), } body.linearDamping = np.zeros((6, 6)) - wave_amp = np.vstack((time, np.zeros_like(time))) + wave_amp = np.vstack(( + time, imported_elevation if imported_elevation is not None + else np.zeros_like(time), + )) if wave["type"] == "regular": body.hydroForcePre( frequency, [direction], 1, np.array([0.0]), [], dt, rho, g, @@ -1462,6 +1481,7 @@ def moving_terms(at_time, coordinate, speed): len(components.omega) if irregular else [], dt, rho, g, ("regularCIC" if regular_memory else "spectrumImport" if wave["type"] == "spectrumImport" else + "elevationImport" if wave["type"] == "elevationImport" else "irregular" if irregular else "noWaveCIC"), wave_amp, index, len(bodies), 0, 0, int(b2b), @@ -1512,6 +1532,19 @@ def excitation(at_time, *, real=real, imaginary=imaginary, real * np.cos(frequency * at_time) - imaginary * np.sin(frequency * at_time) ) + (height / 2)**2 * ramp * drift) + elif wave["type"] == "elevationImport": + sampled_force = np.asarray(hydro_force["userDefinedFe"]) + if wave.get("reapply_force_ramp", False): + sampled_force = sampled_force * imported_ramp[:, None] + + def excitation(at_time, *, sampled_force=sampled_force): + if len(time) == 1: + return sampled_force[0] + sample = min(max(at_time / dt, 0.0), float(len(time) - 1)) + left = min(int(sample), len(time) - 2) + fraction = sample - left + return ((1 - fraction) * sampled_force[left] + + fraction * sampled_force[left + 1]) elif wave["type"] in ("pm", "jonswap", "spectrumImport"): if body_spec.get("passive_yaw", False): def excitation(at_time): @@ -1644,7 +1677,8 @@ def state_inertia(at_time, coordinate, speed, *, terms=moving_terms): ramp[early] = (1 - np.cos(np.pi * response.time[early] / ramp_time)) / 2 elevation = height / 2 * ramp * np.cos(frequency * response.time) else: - elevation = pm_elevation + elevation = (imported_elevation if imported_elevation is not None + else pm_elevation) coordinate_outputs = tuple( output for index, name in enumerate(coordinate_names) for output in ( @@ -1715,6 +1749,13 @@ def state_inertia(at_time, coordinate, speed, *, terms=moving_terms): np.stack([dynamic.excitation(t) for t in response.time])), ) ) if wave["type"] in ("regular", "regularCIC") else () + imported_outputs = tuple( + (f"body{index}_excitation_force", + (np.asarray(body.hydroForce["userDefinedFe"]) + * (imported_ramp[:, None] if wave.get("reapply_force_ramp", False) + else 1))) + for index, body in enumerate(loaded_bodies, start=1) + ) if wave["type"] == "elevationImport" else () passive_outputs = () if passive_model is not None: body_index = passive_indices[0] + 1 @@ -1761,5 +1802,6 @@ def state_inertia(at_time, coordinate, speed, *, terms=moving_terms): generalized_pto if "pto" in case or "ptos" in case else None ), extra_outputs=(coordinate_outputs + pto_outputs + tuple(generator_outputs) + drift_outputs + + imported_outputs + passive_outputs + nonlinear_outputs + morison_outputs), )