diff --git a/.github/workflows/reference-model-baselines.yml b/.github/workflows/reference-model-baselines.yml index c2dbe75..7ae6c97 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, 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] + options: [all, 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: @@ -72,7 +72,7 @@ jobs: strategy: fail-fast: false matrix: - 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"]') }} + model: ${{ fromJSON(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"]' || '["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", "ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", "ELLIPSOID_NLH_ODE45"]') }} 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 == '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' + if: 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 @@ -142,7 +142,7 @@ 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' || matrix.model == 'ELLIPSOID_NLH_CIC' + if: matrix.model == 'ELLIPSOID_NLH_REG' || matrix.model == 'ELLIPSOID_NLH_CIC' || matrix.model == 'ELLIPSOID_NLH_ODE45' env: WEC_SIM_REFERENCE_MODEL: ${{ matrix.model }} WEC_SIM_APPLICATIONS_DIR: applications diff --git a/PARITY.md b/PARITY.md index 3145cf7..13ebeb7 100644 --- a/PARITY.md +++ b/PARITY.md @@ -51,6 +51,7 @@ the production Python code. The live wave comparison passed on 6 October | 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. 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. | | 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) diff --git a/README.md b/README.md index 140a5d5..00e13e4 100644 --- a/README.md +++ b/README.md @@ -473,6 +473,10 @@ and the configured PTO. The STL determines equilibrium mass when Current validation covers one pure-heave body with its center of gravity at horizontal origin in zero-direction regular waves, with either constant or convolution radiation. +The published ode45 variants are tracked separately: MATLAB applies mesh +buoyancy from the preceding 0.05 s sample while the Python mode evaluates it +at the current state. Their motion comparisons are recorded as solver +diagnostics in [PARITY.md](PARITY.md), not as ode45 numerical parity. 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 7a964d7..712134e 100644 --- a/tests/matlab/reference_model_baseline.m +++ b/tests/matlab/reference_model_baseline.m @@ -40,21 +40,28 @@ function reference_model_baseline(model) end cases = ["0m", "1m", "1m-ME", "3m", "5m"]; caseDirs = fullfile(repoRoot, 'applications', 'Free_Decay', cases); - case {"ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"} + case {"ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", ... + "ELLIPSOID_NLH_ODE45"} hydroDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', ... 'hydroData'); cd(hydroDir); if ~isfile('ellipsoid.h5') bemio; end - if string(model) == "ELLIPSOID_NLH_REG" - waveDir = "Regular"; + if string(model) == "ELLIPSOID_NLH_ODE45" + cases = ["ode45_Regular", "ode45_RegularCIC"]; + caseDirs = fullfile(repoRoot, 'applications', ... + 'Nonlinear_Hydro', 'ode45', ["Regular", "RegularCIC"]); else - waveDir = "RegularCIC"; + if string(model) == "ELLIPSOID_NLH_REG" + waveDir = "Regular"; + else + waveDir = "RegularCIC"; + end + cases = "ode4_" + waveDir; + caseDirs = string(fullfile(repoRoot, 'applications', ... + 'Nonlinear_Hydro', 'ode4', waveDir)); end - cases = "ode4_" + waveDir; - caseDirs = string(fullfile(repoRoot, 'applications', ... - 'Nonlinear_Hydro', 'ode4', waveDir)); case "SPHERE_MEAN_DRIFT" hydroDir = fullfile(repoRoot, 'applications', 'Mean_Drift', 'hydroData'); cd(hydroDir); @@ -408,12 +415,21 @@ function reference_model_baseline(model) fullfile(outDir, 'OSWEC_wave_directions.csv')); writematrix(waves.waveAmpTime, fullfile(outDir, 'OSWEC_wave_elevation.csv')); end - if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"]) - if string(model) == "ELLIPSOID_NLH_REG" + if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", ... + "ELLIPSOID_NLH_ODE45"]) + if string(cases(iCase)) == "ode4_Regular" || ... + string(cases(iCase)) == "ode45_Regular" waveType = 'regular'; else waveType = 'regularCIC'; end + if string(model) == "ELLIPSOID_NLH_ODE45" + assert(strcmp(simu.solver, 'ode45'), ... + 'The pinned nonlinear-hydro solver changed'); + prefix = string(model) + "_" + cases(iCase); + else + prefix = string(model); + end assert(simu.dt == 0.05 && simu.endTime == 150 && ... simu.rampTime == 50 && simu.rho == 1025 && ... strcmp(waves.type, waveType) && waves.height == 4 && ... @@ -422,9 +438,9 @@ function reference_model_baseline(model) isequal(constraint(1).location, [0 0 -12.5]), ... 'The pinned nonlinear-hydro input changed'); writematrix([output.wave.time(:), output.wave.elevation(:)], ... - fullfile(outDir, string(model) + '_wave.csv')); + fullfile(outDir, prefix + '_wave.csv')); writematrix([body(1).mass, body(1).inertia], ... - fullfile(outDir, string(model) + '_mass.csv')); + fullfile(outDir, prefix + '_mass.csv')); end if string(model) == "OSWEC_MULTI_WAVE" assert(simu.dt == 0.1 && simu.endTime == 100 && ... @@ -618,7 +634,8 @@ function reference_model_baseline(model) values = [values, response.forceRadiationDamping, ... response.forceAddedMass, response.forceRestoring]; end - if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC"]) + if any(string(model) == ["ELLIPSOID_NLH_REG", "ELLIPSOID_NLH_CIC", ... + "ELLIPSOID_NLH_ODE45"]) 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 2a2f8bb..ff91631 100644 --- a/tests/test_ellipsoid_nonlinear_hydro_parity.py +++ b/tests/test_ellipsoid_nonlinear_hydro_parity.py @@ -1,4 +1,4 @@ -"""Paired pinned Nonlinear_Hydro/ode4/Regular ellipsoid validation.""" +"""Paired pinned Nonlinear_Hydro heaving-ellipsoid validation.""" import os from pathlib import Path @@ -15,21 +15,25 @@ 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") +CASES = ({"ELLIPSOID_NLH_REG": ("ode4_Regular",), + "ELLIPSOID_NLH_CIC": ("ode4_RegularCIC",), + "ELLIPSOID_NLH_ODE45": ("ode45_Regular", "ode45_RegularCIC")} + .get(MODEL, ())) pytestmark = pytest.mark.skipif( not (APPLICATIONS and REFERENCE), reason="paired MATLAB nonlinear-hydro output and Applications absent", ) -def _source(): +def _source(case_name): source = Path(REFERENCE) - body = np.loadtxt(source / f"{MODEL}_{CASE}_body1.csv", + body = np.loadtxt(source / f"{MODEL}_{case_name}_body1.csv", delimiter=",") - pto = np.loadtxt(source / f"{MODEL}_{CASE}_pto1.csv", + pto = np.loadtxt(source / f"{MODEL}_{case_name}_pto1.csv", delimiter=",") - mass = np.loadtxt(source / f"{MODEL}_mass.csv", delimiter=",") + mass_prefix = (f"{MODEL}_{case_name}" if MODEL == "ELLIPSOID_NLH_ODE45" + else MODEL) + mass = np.loadtxt(source / f"{mass_prefix}_mass.csv", delimiter=",") assert body.shape == (3001, 55) and pto.shape == (3001, 25) return body, pto, mass @@ -46,8 +50,9 @@ def _max_error(actual, expected, limit, label): assert error < limit, f"{label}: {error:.6g} exceeds {limit}" -def test_mesh_forces_on_matlab_trajectory(): - body, _, mass = _source() +@pytest.mark.parametrize("case_name", CASES) +def test_mesh_forces_on_matlab_trajectory(case_name): + body, pto, mass = _source(case_name) hydro_file, geometry_file = _files() mesh = HeaveMeshHydro.from_stl( geometry_file, center_z=-2, rho=1025, gravity=9.81, @@ -59,8 +64,15 @@ def test_mesh_forces_on_matlab_trajectory(): mesh.forces(t, z + 2, v) for t, z, v in zip(body[:, 0], body[:, 3], body[:, 9]) ]) - _max_error(forces[:, 0], -body[:, 39], 1e-6, - "mesh buoyancy minus weight") + if MODEL == "ELLIPSOID_NLH_ODE45": + # The source's ode45 nonlinear buoyancy block holds the preceding + # 0.05 s sample while motion and the other force logs advance. + _max_error(forces[:-1, 0], -body[1:, 39], 1e-6, + "sampled mesh buoyancy minus weight") + assert np.max(np.abs(forces[:, 0] + body[:, 39])) > 40_000 + else: + _max_error(forces[:, 0], -body[:, 39], 1e-6, + "mesh buoyancy minus weight") _max_error(forces[:, 2], -body[:, 45], 1e-6, "quadratic drag") # WEC-Sim logs its linear BEM excitation plus the nonlinear FK correction. @@ -74,7 +86,7 @@ def test_mesh_forces_on_matlab_trajectory(): "characteristicArea": np.zeros(6)} bem.linearDamping = np.zeros((6, 6)) omega = 2 * np.pi / 6 - cic = MODEL == "ELLIPSOID_NLH_CIC" + cic = case_name.endswith("RegularCIC") 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, @@ -87,6 +99,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 MODEL == "ELLIPSOID_NLH_ODE45": + # The lagged restoring force is in the applied source force sum, + # not merely a reporting offset in the body output. + total = (body[:, 21] - body[:, 27] - body[:, 33] + - body[:, 39] - body[:, 45]) + _max_error(total, body[:, 15], 1e-6, "source force assembly") + adjusted_mass = (mesh.mass + 2 * np.trace( + np.asarray(bem.hydroForce["fAddedMass"])[:3, :3] + )) + _max_error(adjusted_mass * body[:, 51], body[:, 15] + pto[:, 15], + 1e-6, "source acceleration and adjusted mass") if cic: kernel = np.asarray(bem.hydroForce["irkb"])[:, 2, 2] velocity = body[:, 9] @@ -96,12 +119,14 @@ def test_mesh_forces_on_matlab_trajectory(): radiation -= .05 / 2 * kernel[last_lag] * ( velocity[np.arange(len(velocity)) - last_lag] ) - _max_error(radiation, body[:, 27], 1e-6, + _max_error(radiation, body[:, 27], + 300 if MODEL == "ELLIPSOID_NLH_ODE45" else 1e-6, "regularCIC radiation convolution") -def test_public_python_configuration_matches_matlab_motion_and_pto(): - body, pto, mass = _source() +@pytest.mark.parametrize("case_name", CASES) +def test_public_python_configuration_against_matlab_motion_and_pto(case_name): + body, pto, mass = _source(case_name) hydro_file, geometry_file = _files() wec = WEC("ellipsoid") ellipsoid = wec.body( @@ -112,13 +137,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) - cic = MODEL == "ELLIPSOID_NLH_CIC" + cic = case_name.endswith("RegularCIC") 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)) + if MODEL == "ELLIPSOID_NLH_ODE45": + # The published MATLAB ode4 and ode45 runs themselves differ by up + # to 16 mm and 35 kN with identical physical settings. This is a + # solver-envelope check, not a claim of matching ode45 internals. + limits = ((.023, .034, 40_000, 33_000) if cic + else (.022, .032, 38_000, 30_000)) + else: + 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], limits[0], "heave position") _max_error(result.bodies["ellipsoid"].velocity[:, 2], body[:, 9],