diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index 88fbbdf..c2dbe75 100644 --- a/.github/workflows/reference-model-baselines.yml +++ b/.github/workflows/reference-model-baselines.yml @@ -62,7 +62,7 @@ on: model: description: Run all models or a selected RM3 application type: choice - options: [all, ELLIPSOID_NLH_REG, OSWEC_PASSIVE_YAW_IRR_CONT, OSWEC_PASSIVE_YAW_IRR, OSWEC_PASSIVE_YAW, OSWEC_FULL_DIR, OSWEC_MULTI_WAVE, RM3_B2B, RM3_MOORING_MATRIX, SPHERE_MEAN_DRIFT, RM3_MCR_SEASTATE, RM3_END_STOPS, RM3_END_STOPS_STEP, RM3_END_STOPS_STEP_FINE, RM3_END_STOPS_STEP_FINER, RM3_END_STOPS_FULL_FINE, RM3_END_STOPS_FULL_FINER, RM3_END_STOPS_FULL_PAIR] + options: [all, ELLIPSOID_NLH_REG, ELLIPSOID_NLH_CIC, OSWEC_PASSIVE_YAW_IRR_CONT, OSWEC_PASSIVE_YAW_IRR, OSWEC_PASSIVE_YAW, OSWEC_FULL_DIR, OSWEC_MULTI_WAVE, RM3_B2B, RM3_MOORING_MATRIX, SPHERE_MEAN_DRIFT, RM3_MCR_SEASTATE, RM3_END_STOPS, RM3_END_STOPS_STEP, RM3_END_STOPS_STEP_FINE, RM3_END_STOPS_STEP_FINER, RM3_END_STOPS_FULL_FINE, RM3_END_STOPS_FULL_FINER, RM3_END_STOPS_FULL_PAIR] default: all jobs: @@ -72,7 +72,7 @@ jobs: strategy: fail-fast: false matrix: - model: ${{ fromJSON(inputs.model == 'ELLIPSOID_NLH_REG' && '["ELLIPSOID_NLH_REG"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' && '["OSWEC_PASSIVE_YAW_IRR_CONT"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR' && '["OSWEC_PASSIVE_YAW_IRR"]' || inputs.model == 'OSWEC_PASSIVE_YAW' && '["OSWEC_PASSIVE_YAW"]' || inputs.model == 'OSWEC_FULL_DIR' && '["OSWEC_FULL_DIR"]' || inputs.model == 'OSWEC_MULTI_WAVE' && '["OSWEC_MULTI_WAVE"]' || inputs.model == 'RM3_B2B' && '["RM3_B2B"]' || inputs.model == 'RM3_MOORING_MATRIX' && '["RM3_MOORING_MATRIX"]' || inputs.model == 'SPHERE_MEAN_DRIFT' && '["SPHERE_MEAN_DRIFT"]' || inputs.model == 'RM3_MCR_SEASTATE' && '["RM3_MCR_SEASTATE"]' || inputs.model == 'RM3_END_STOPS' && '["RM3_END_STOPS"]' || inputs.model == 'RM3_END_STOPS_STEP' && '["RM3_END_STOPS_STEP"]' || inputs.model == 'RM3_END_STOPS_STEP_FINE' && '["RM3_END_STOPS_STEP_FINE"]' || inputs.model == 'RM3_END_STOPS_STEP_FINER' && '["RM3_END_STOPS_STEP_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_FINE' && '["RM3_END_STOPS_FULL_FINE"]' || inputs.model == 'RM3_END_STOPS_FULL_FINER' && '["RM3_END_STOPS_FULL_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_PAIR' && '["RM3_END_STOPS_FULL_FINE", "RM3_END_STOPS_FULL_FINER"]' || '["RM3", "OSWEC", "OSWEC_Nonhydro", "OSWEC_FULL_DIR", "OSWEC_MULTI_WAVE", "OSWEC_PASSIVE_YAW", "Sphere", "SPHERE_MEAN_DRIFT", "RM3_B2B", "RM3_END_STOPS", "RM3_PTO_Extension", "RM3_Radiation_Options", "Sphere_Passive", "Sphere_PTO_Config", "Sphere_Reactive_PI", "Sphere_Declutching", "Sphere_Latching", "RM3_MCR", "RM3_MCR_ARRAY", "RM3_MCR_EXCEL", "RM3_MCR_MAT", "RM3_MCR_SEASTATE"]') }} + model: ${{ fromJSON(inputs.model == 'ELLIPSOID_NLH_REG' && '["ELLIPSOID_NLH_REG"]' || inputs.model == 'ELLIPSOID_NLH_CIC' && '["ELLIPSOID_NLH_CIC"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' && '["OSWEC_PASSIVE_YAW_IRR_CONT"]' || inputs.model == 'OSWEC_PASSIVE_YAW_IRR' && '["OSWEC_PASSIVE_YAW_IRR"]' || inputs.model == 'OSWEC_PASSIVE_YAW' && '["OSWEC_PASSIVE_YAW"]' || inputs.model == 'OSWEC_FULL_DIR' && '["OSWEC_FULL_DIR"]' || inputs.model == 'OSWEC_MULTI_WAVE' && '["OSWEC_MULTI_WAVE"]' || inputs.model == 'RM3_B2B' && '["RM3_B2B"]' || inputs.model == 'RM3_MOORING_MATRIX' && '["RM3_MOORING_MATRIX"]' || inputs.model == 'SPHERE_MEAN_DRIFT' && '["SPHERE_MEAN_DRIFT"]' || inputs.model == 'RM3_MCR_SEASTATE' && '["RM3_MCR_SEASTATE"]' || inputs.model == 'RM3_END_STOPS' && '["RM3_END_STOPS"]' || inputs.model == 'RM3_END_STOPS_STEP' && '["RM3_END_STOPS_STEP"]' || inputs.model == 'RM3_END_STOPS_STEP_FINE' && '["RM3_END_STOPS_STEP_FINE"]' || inputs.model == 'RM3_END_STOPS_STEP_FINER' && '["RM3_END_STOPS_STEP_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_FINE' && '["RM3_END_STOPS_FULL_FINE"]' || inputs.model == 'RM3_END_STOPS_FULL_FINER' && '["RM3_END_STOPS_FULL_FINER"]' || inputs.model == 'RM3_END_STOPS_FULL_PAIR' && '["RM3_END_STOPS_FULL_FINE", "RM3_END_STOPS_FULL_FINER"]' || '["RM3", "OSWEC", "OSWEC_Nonhydro", "OSWEC_FULL_DIR", "OSWEC_MULTI_WAVE", "OSWEC_PASSIVE_YAW", "Sphere", "SPHERE_MEAN_DRIFT", "RM3_B2B", "RM3_END_STOPS", "RM3_PTO_Extension", "RM3_Radiation_Options", "Sphere_Passive", "Sphere_PTO_Config", "Sphere_Reactive_PI", "Sphere_Declutching", "Sphere_Latching", "RM3_MCR", "RM3_MCR_ARRAY", "RM3_MCR_EXCEL", "RM3_MCR_MAT", "RM3_MCR_SEASTATE"]') }} name: ${{ matrix.model }} steps: - uses: actions/checkout@v7 @@ -86,7 +86,7 @@ jobs: source examples/${{ matrix.model }} - uses: actions/checkout@v7 - if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'Sphere' || matrix.model == 'SPHERE_MEAN_DRIFT' || matrix.model == 'Sphere_Passive' || matrix.model == 'Sphere_PTO_Config' || matrix.model == 'Sphere_Reactive_PI' || matrix.model == 'Sphere_Declutching' || matrix.model == 'Sphere_Latching' || matrix.model == 'RM3_B2B' || matrix.model == 'RM3_MOORING_MATRIX' || matrix.model == 'RM3_END_STOPS' || matrix.model == 'RM3_END_STOPS_STEP' || matrix.model == 'RM3_END_STOPS_STEP_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINER' || matrix.model == 'RM3_END_STOPS_STEP_FINER' || matrix.model == 'RM3_PTO_Extension' || matrix.model == 'RM3_Radiation_Options' || matrix.model == 'RM3_MCR' || matrix.model == 'RM3_MCR_ARRAY' || matrix.model == 'RM3_MCR_EXCEL' || matrix.model == 'RM3_MCR_MAT' || matrix.model == 'RM3_MCR_SEASTATE' || matrix.model == 'OSWEC_Nonhydro' || matrix.model == 'OSWEC_MULTI_WAVE' || matrix.model == 'OSWEC_FULL_DIR' || matrix.model == 'OSWEC_PASSIVE_YAW' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' + if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'Sphere' || matrix.model == 'SPHERE_MEAN_DRIFT' || matrix.model == 'Sphere_Passive' || matrix.model == 'Sphere_PTO_Config' || matrix.model == 'Sphere_Reactive_PI' || matrix.model == 'Sphere_Declutching' || matrix.model == 'Sphere_Latching' || matrix.model == 'RM3_B2B' || matrix.model == 'RM3_MOORING_MATRIX' || matrix.model == 'RM3_END_STOPS' || matrix.model == 'RM3_END_STOPS_STEP' || matrix.model == 'RM3_END_STOPS_STEP_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINE' || matrix.model == 'RM3_END_STOPS_FULL_FINER' || matrix.model == 'RM3_END_STOPS_STEP_FINER' || matrix.model == 'RM3_PTO_Extension' || matrix.model == 'RM3_Radiation_Options' || matrix.model == 'RM3_MCR' || matrix.model == 'RM3_MCR_ARRAY' || matrix.model == 'RM3_MCR_EXCEL' || matrix.model == 'RM3_MCR_MAT' || matrix.model == 'RM3_MCR_SEASTATE' || matrix.model == 'OSWEC_Nonhydro' || matrix.model == 'OSWEC_MULTI_WAVE' || matrix.model == 'OSWEC_FULL_DIR' || matrix.model == 'OSWEC_PASSIVE_YAW' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR' || matrix.model == 'OSWEC_PASSIVE_YAW_IRR_CONT' with: repository: WEC-Sim/WEC-Sim_Applications ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e @@ -142,8 +142,9 @@ jobs: WEC_SIM_SPHERE_H5: applications/_Common_Input_Files/Sphere/hydroData/sphere.h5 WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - run: python -m pytest -q tests/test_ellipsoid_nonlinear_hydro_parity.py - if: matrix.model == 'ELLIPSOID_NLH_REG' + if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' env: + WEC_SIM_REFERENCE_MODEL: ${{ matrix.model }} WEC_SIM_APPLICATIONS_DIR: applications WEC_SIM_MATLAB_MODEL_OUTPUT_DIR: matlab-reference-model-output - run: python -m pytest -q tests/test_live_current_h5.py diff --git a/PARITY.md b/PARITY.md index 7d591c4..3145cf7 100644 --- a/PARITY.md +++ b/PARITY.md @@ -49,8 +49,9 @@ the production Python code. The live wave comparison passed on 6 October | OSWEC irregular passive yaw | Pinned MATLAB Applications `PassiveYawRegression` input, BEMIO-generated OSWEC HDF5, and fresh R2025b trajectory | The public `PMWave` replays the published 500-frequency phase realization over 250 s with 40 s radiation memory. Frequency bins, spectral amplitudes, and widths agree within `6e-15`; elevation agrees within `1e-13` m. MATLAB holds its excitation-heading coefficients for 1° of relative yaw; reconstructing that rule on the saved MATLAB yaw path matches all six logged force components within `1e-7` N or N m. With continuous heading interpolation, Python differs from the published trajectory by up to 0.452 rad in flap yaw, 0.101 rad/s in yaw speed, 12.2 kN m in PTO torque, and 1.06 kW in source-signed power. Replaying MATLAB's logged six-component force through the Python dynamics reduces these maximum differences to 0.00180 rad, 0.0000935 rad/s, 11.3 N m, and 2.02 W. This isolates the large discrepancy to excitation heading updates without changing the Python default. Published-case trajectory parity with Python's native excitation is **not** established. | | OSWEC irregular passive yaw, continuous-heading control | Same pinned MATLAB input and phase realization, with only `body(1:2).yaw.threshold` changed from 1° to 0° in a temporary application copy | MATLAB and Python then agree over all 25,001 time samples: maximum differences are `1.1e-7` N or N m across the six excitation components on the MATLAB yaw path, `3.16e-5` rad flap/PTO angle, `7.30e-6` rad/s yaw/PTO speed, `0.875` N m PTO torque, and `0.0123` W source-signed power. Wave elevation differs by less than `1e-13` m and the fixed base remains stationary. This control validates the continuous-heading Python dynamics; it does not erase the published 1° threshold difference. | | OSWEC fixed nonhydrodynamic base | Pinned MATLAB Applications `Nonhydro_Body` case and its BEMIO-generated OSWEC HDF5 | The regular-wave solver reports the stationary base and nonlinear flap motion about the PTO hinge. Against a fresh 400 s MATLAB R2025b run (4,001 samples), maximum flap position differences are 20.8 mm surge, 7.9 mm heave, and 0.00444 rad pitch; velocity differences are 15.6 mm/s surge, 6.8 mm/s heave, and 0.00328 rad/s pitch. All six excitation-force components differ by less than 1 N, the base position and velocity agree exactly, and zero PTO torque differs only by MATLAB numerical noise below `3.3e-7` N m. The base's ground-constraint reaction forces are not calculated. | -| Ellipsoid instantaneous nonlinear hydro, `ode4/Regular` | Pinned MATLAB Applications `Nonlinear_Hydro` and BEMIO-generated ellipsoid HDF5/STL | The Python builder combines mesh buoyancy, instantaneous Froude–Krylov correction, quadratic heave drag, BEM linear excitation/radiation, and a configured PTO. On all 3,001 saved MATLAB states, mesh equilibrium mass agrees within `2e-9` kg, buoyancy within `6e-9` N, drag within `7e-11` N, and total heave excitation within `1e-6` N. The independent 150 s Python trajectory differs by at most 5.16 mm heave, 5.86 mm/s velocity, 5.16 mm PTO stroke, 7.03 kN PTO force, and 7.21 kW absorbed power. MATLAB splits the BEM added mass between Simscape mass and an applied force; Python uses the combined effective mass. Only the single-body, pure-heave, zero-direction regular-wave mode is supported so far. | -| Case-driven dynamics runner | Current MATLAB RM3 and OSWEC examples plus Sphere and RM3 Applications cases | One generalized-coordinate engine assembles rigid inertia, hydrodynamic added mass and radiation, hydrostatic restoring, excitation, and linear PTO forces. The heave, fixed-hinge, and floating-joint layouts run the paired cases above, including imported elevation and a joint surge spring for RM3 MooringMatrix. A fourth `linear_subspace` layout maps independent coordinates into arbitrary bodies; its RM3 two-body heave, Sphere free decay, configured Sphere PTO, and one instantaneous nonlinear-hydro case have paired MATLAB checks. Arbitrary Simscape layouts, general moorings, and other nonlinear-hydro/application cases remain unsupported. | +| Ellipsoid instantaneous nonlinear hydro, `ode4/Regular` | Pinned MATLAB Applications `Nonlinear_Hydro` and BEMIO-generated ellipsoid HDF5/STL | The Python builder combines mesh buoyancy, instantaneous Froude–Krylov correction, quadratic heave drag, BEM linear excitation/radiation, and a configured PTO. On all 3,001 saved MATLAB states, mesh equilibrium mass agrees within `2e-9` kg, buoyancy within `6e-9` N, drag within `7e-11` N, and total heave excitation within `1e-6` N. The independent 150 s Python trajectory differs by at most 5.16 mm heave, 5.86 mm/s velocity, 5.16 mm PTO stroke, 7.03 kN PTO force, and 7.21 kW absorbed power. MATLAB splits the BEM added mass between Simscape mass and an applied force; Python uses the combined effective mass. This case validates the single-body, pure-heave, zero-direction constant-radiation mode. | +| Ellipsoid instantaneous nonlinear hydro, `ode4/RegularCIC` | Pinned MATLAB Applications `Nonlinear_Hydro` RegularCIC input and BEMIO-generated ellipsoid HDF5/STL | The same Python mesh model uses a 60 s radiation convolution. On all 3,001 saved MATLAB states, mesh buoyancy and quadratic drag agree within `6e-9` and `7e-11` N, linear plus nonlinear heave excitation within `4e-7` N, and radiation convolution within `3e-10` N. The independent 150 s Python trajectory differs by at most 6.39 mm heave, 7.41 mm/s velocity, 6.39 mm PTO stroke, 8.89 kN PTO force, and 8.98 kW absorbed power. The memory-step fixed-point iteration was allowed more iterations to converge near the moving waterline; no force coefficient was tuned. | +| Case-driven dynamics runner | Current MATLAB RM3 and OSWEC examples plus Sphere and RM3 Applications cases | One generalized-coordinate engine assembles rigid inertia, hydrodynamic added mass and radiation, hydrostatic restoring, excitation, and linear PTO forces. The heave, fixed-hinge, and floating-joint layouts run the paired cases above, including imported elevation and a joint surge spring for RM3 MooringMatrix. A fourth `linear_subspace` layout maps independent coordinates into arbitrary bodies; its RM3 two-body heave, Sphere free decay, configured Sphere PTO, and two instantaneous nonlinear-hydro cases have paired MATLAB checks. Arbitrary Simscape layouts, general moorings, and other nonlinear-hydro/application cases remain unsupported. | The targeted [RM3 sea-state matrix export](https://github.com/cmudrc/wec-sim-python/actions/runs/37682612044) confirms that MATLAB's applied added-mass matrix and adjusted rigid mass sum diff --git a/README.md b/README.md index a36f247..140a5d5 100644 --- a/README.md +++ b/README.md @@ -439,8 +439,9 @@ reactive PI laws, and hydraulic PTO models are not implemented. ### Instantaneous nonlinear hydrodynamics for heave -The Python builder supports the published `Nonlinear_Hydro/ode4/Regular` -ellipsoid through a heave-only mesh mode. Supply the BEM HDF5 file and an STL +The Python builder supports the published `Nonlinear_Hydro/ode4/Regular` and +`ode4/RegularCIC` ellipsoid cases through a heave-only mesh mode. Supply the +BEM HDF5 file and an STL whose triangle coordinates are relative to the body's center of gravity: ```python @@ -462,12 +463,16 @@ result = wec.run(RegularWave(height=4, period=6), dt=0.05, end_time=150, ramp_time=50, rho=1025) ``` +For convolution radiation, use `RegularCICWave(height=4, period=6)` and +`radiation_memory=60` in `wec.run`. + The mode integrates mesh buoyancy, instantaneous free-surface Froude–Krylov correction, heave quadratic drag, BEM diffraction/radiation, and the configured PTO. The STL determines equilibrium mass when `mass="equilibrium"`; this can differ from the HDF5 displaced volume. Current validation covers one pure-heave body with its center of gravity at -horizontal origin in zero-direction regular waves. +horizontal origin in zero-direction regular waves, with either constant or +convolution radiation. Other motions and sea states raise an error until their mesh force and dynamics checks are paired with MATLAB. diff --git a/tests/matlab/reference_model_baseline.m b/tests/matlab/reference_model_baseline.m index 8f64668..7a964d7 100644 --- a/tests/matlab/reference_model_baseline.m +++ b/tests/matlab/reference_model_baseline.m @@ -40,16 +40,21 @@ function reference_model_baseline(model) end cases = ["0m", "1m", "1m-ME", "3m", "5m"]; caseDirs = fullfile(repoRoot, 'applications', 'Free_Decay', cases); - case "ELLIPSOID_NLH_REG" + case {"ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"} hydroDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', ... 'hydroData'); cd(hydroDir); if ~isfile('ellipsoid.h5') bemio; end - cases = "ode4_Regular"; + if string(model) == "ELLIPSOID_NLH_REG" + waveDir = "Regular"; + else + waveDir = "RegularCIC"; + end + cases = "ode4_" + waveDir; caseDirs = string(fullfile(repoRoot, 'applications', ... - 'Nonlinear_Hydro', 'ode4', 'Regular')); + 'Nonlinear_Hydro', 'ode4', waveDir)); case "SPHERE_MEAN_DRIFT" hydroDir = fullfile(repoRoot, 'applications', 'Mean_Drift', 'hydroData'); cd(hydroDir); @@ -403,18 +408,23 @@ function reference_model_baseline(model) fullfile(outDir, 'OSWEC_wave_directions.csv')); writematrix(waves.waveAmpTime, fullfile(outDir, 'OSWEC_wave_elevation.csv')); end - if string(model) == "ELLIPSOID_NLH_REG" + if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"]) + if string(model) == "ELLIPSOID_NLH_REG" + waveType = 'regular'; + else + waveType = 'regularCIC'; + end assert(simu.dt == 0.05 && simu.endTime == 150 && ... simu.rampTime == 50 && simu.rho == 1025 && ... - strcmp(waves.type, 'regular') && waves.height == 4 && ... + strcmp(waves.type, waveType) && waves.height == 4 && ... waves.period == 6 && body(1).nonlinearHydro == 2 && ... pto(1).damping == 1200000 && ... isequal(constraint(1).location, [0 0 -12.5]), ... 'The pinned nonlinear-hydro input changed'); writematrix([output.wave.time(:), output.wave.elevation(:)], ... - fullfile(outDir, 'ELLIPSOID_NLH_REG_wave.csv')); + fullfile(outDir, string(model) + '_wave.csv')); writematrix([body(1).mass, body(1).inertia], ... - fullfile(outDir, 'ELLIPSOID_NLH_REG_mass.csv')); + fullfile(outDir, string(model) + '_mass.csv')); end if string(model) == "OSWEC_MULTI_WAVE" assert(simu.dt == 0.1 && simu.endTime == 100 && ... @@ -608,7 +618,7 @@ function reference_model_baseline(model) values = [values, response.forceRadiationDamping, ... response.forceAddedMass, response.forceRestoring]; end - if string(model) == "ELLIPSOID_NLH_REG" + if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"]) values = [values, response.forceRadiationDamping, ... response.forceAddedMass, response.forceRestoring, ... response.forceMorisonAndViscous, response.acceleration]; diff --git a/tests/test_ellipsoid_nonlinear_hydro_parity.py b/tests/test_ellipsoid_nonlinear_hydro_parity.py index 8afa388..2a2f8bb 100644 --- a/tests/test_ellipsoid_nonlinear_hydro_parity.py +++ b/tests/test_ellipsoid_nonlinear_hydro_parity.py @@ -5,14 +5,18 @@ import numpy as np import pytest +from scipy.signal import fftconvolve -from wecsim.api import RegularWave, WEC, WorldPoint +from wecsim.api import RegularCICWave, RegularWave, WEC, WorldPoint from wecsim.bodyClass import BodyClass from wecsim.nonlinearHydro import HeaveMeshHydro APPLICATIONS = os.environ.get("WEC_SIM_APPLICATIONS_DIR") REFERENCE = os.environ.get("WEC_SIM_MATLAB_MODEL_OUTPUT_DIR") +MODEL = os.environ.get("WEC_SIM_REFERENCE_MODEL", "ELLIPSOID_NLH_REG") +CASE = ("ode4_RegularCIC" if MODEL == "ELLIPSOID_NLH_CIC" + else "ode4_Regular") pytestmark = pytest.mark.skipif( not (APPLICATIONS and REFERENCE), reason="paired MATLAB nonlinear-hydro output and Applications absent", @@ -21,11 +25,11 @@ def _source(): source = Path(REFERENCE) - body = np.loadtxt(source / "ELLIPSOID_NLH_REG_ode4_Regular_body1.csv", + body = np.loadtxt(source / f"{MODEL}_{CASE}_body1.csv", delimiter=",") - pto = np.loadtxt(source / "ELLIPSOID_NLH_REG_ode4_Regular_pto1.csv", + pto = np.loadtxt(source / f"{MODEL}_{CASE}_pto1.csv", delimiter=",") - mass = np.loadtxt(source / "ELLIPSOID_NLH_REG_mass.csv", delimiter=",") + mass = np.loadtxt(source / f"{MODEL}_mass.csv", delimiter=",") assert body.shape == (3001, 55) and pto.shape == (3001, 25) return body, pto, mass @@ -70,9 +74,12 @@ def test_mesh_forces_on_matlab_trajectory(): "characteristicArea": np.zeros(6)} bem.linearDamping = np.zeros((6, 6)) omega = 2 * np.pi / 6 - bem.hydroForcePre(omega, [0], 1, np.array([0.0]), [], .05, - 1025, 9.81, "regular", np.zeros((2, 2)), - 1, 1, 0, 0, 0) + cic = MODEL == "ELLIPSOID_NLH_CIC" + convolution_time = np.arange(1201) * .05 if cic else np.array([0.0]) + bem.hydroForcePre(omega, [0], len(convolution_time), convolution_time, + [], .05, 1025, 9.81, + "regularCIC" if cic else "regular", + np.zeros((2, 2)), 1, 1, 0, 0, 0) re = np.asarray(bem.hydroForce["fExt"]["re"])[2] im = np.asarray(bem.hydroForce["fExt"]["im"])[2] ramp = np.array([mesh.ramp(t) for t in body[:, 0]]) @@ -80,6 +87,17 @@ def test_mesh_forces_on_matlab_trajectory(): - im * np.sin(omega * body[:, 0])) _max_error(linear + forces[:, 1], body[:, 21], 1e-5, "linear plus nonlinear Froude-Krylov excitation") + if cic: + kernel = np.asarray(bem.hydroForce["irkb"])[:, 2, 2] + velocity = body[:, 9] + radiation = .05 * fftconvolve(velocity, kernel)[:len(velocity)] + radiation -= .05 / 2 * kernel[0] * velocity + last_lag = np.minimum(np.arange(len(velocity)), len(kernel) - 1) + radiation -= .05 / 2 * kernel[last_lag] * ( + velocity[np.arange(len(velocity)) - last_lag] + ) + _max_error(radiation, body[:, 27], 1e-6, + "regularCIC radiation convolution") def test_public_python_configuration_matches_matlab_motion_and_pto(): @@ -94,16 +112,20 @@ def test_public_python_configuration_matches_matlab_motion_and_pto(): wec.coordinate("heave", ellipsoid.move("heave")) wec.pto("PTO1", WorldPoint(0, 0, -12.5), ellipsoid.at(0, 0, 0), damping=1_200_000) - result = wec.run(RegularWave(4, 6), dt=.05, end_time=150, - ramp_time=50, rho=1025) + cic = MODEL == "ELLIPSOID_NLH_CIC" + result = wec.run(RegularCICWave(4, 6) if cic else RegularWave(4, 6), + dt=.05, end_time=150, ramp_time=50, + radiation_memory=60 if cic else None, rho=1025) np.testing.assert_allclose(result.time, body[:, 0], rtol=0, atol=1e-10) + limits = ((.0075, .0085, 10_000, 10_000) if cic + else (.006, .0065, 8_000, 8_000)) _max_error(result.bodies["ellipsoid"].position[:, 2], body[:, 3], - .006, "heave position") + limits[0], "heave position") _max_error(result.bodies["ellipsoid"].velocity[:, 2], body[:, 9], - .0065, "heave velocity") - _max_error(result.ptos["PTO1"].stroke, pto[:, 3], .006, + limits[1], "heave velocity") + _max_error(result.ptos["PTO1"].stroke, pto[:, 3], limits[0], "PTO stroke") - _max_error(result.ptos["PTO1"].force, pto[:, 15], 8_000, + _max_error(result.ptos["PTO1"].force, pto[:, 15], limits[2], "PTO force") _max_error(result.ptos["PTO1"].absorbed_power, - -pto[:, 15] * pto[:, 9], 8_000, "PTO absorbed power") + -pto[:, 15] * pto[:, 9], limits[3], "PTO absorbed power") diff --git a/wecsim/caseDynamics.py b/wecsim/caseDynamics.py index 587abc6..de81002 100644 --- a/wecsim/caseDynamics.py +++ b/wecsim/caseDynamics.py @@ -795,13 +795,14 @@ def _run_linear_subspace(case, sim, wave, constraint, bodies, hydro, heave_map[2, 0] = 1 if (len(bodies) != 1 or nonlinear_indices != [0] or n != 1 or not np.array_equal(maps[0], heave_map) - or wave["type"] != "regular" or direction != 0 or b2b + or wave["type"] not in ("regular", "regularCIC") + or direction != 0 or b2b or bodies[0]["nonlinear_hydro"] != "instantaneous" or bodies[0].get("mean_drift", "none") != "none" or bodies[0].get("passive_yaw", False)): raise ValueError( "instantaneous nonlinear hydro currently needs one pure-heave " - "body and a zero-direction regular wave" + "body and a zero-direction regular or regularCIC wave" ) for spec in bodies: if spec.get("nonlinear_hydro") is None and any( diff --git a/wecsim/generalDynamics.py b/wecsim/generalDynamics.py index eea8713..738efc4 100644 --- a/wecsim/generalDynamics.py +++ b/wecsim/generalDynamics.py @@ -442,7 +442,9 @@ def known_radiation(step): self.delayed_body_acceleration = tuple(delayed) known = known_radiation(step) trial_speed = v[step - 1].copy() - for _ in range(12): + # Mesh-pressure forces can make the trapezoidal fixed point + # slower to converge near the instantaneous waterline. + for _ in range(30): trial_coordinate = q[step - 1] + dt * (v[step - 1] + trial_speed) / 2 trial_acceleration = self.acceleration( time[step], trial_coordinate, trial_speed, diff --git a/wecsim/nonlinearHydro.py b/wecsim/nonlinearHydro.py index acf5c37..1305a4e 100644 --- a/wecsim/nonlinearHydro.py +++ b/wecsim/nonlinearHydro.py @@ -1,7 +1,8 @@ """Instantaneous free-surface forces on a triangulated heaving body. This implements WEC-Sim's nonlinearHydro=2 regular-wave force construction -for a body constrained to heave. Mesh triangles are in body coordinates about +for a body constrained to heave, with constant or convolution radiation. +Mesh triangles are in body coordinates about the center of gravity, as required by WEC-Sim's geometry import. """