From 5b2a25721a36270b0beb5ace32026f76a3b3135d Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Sat, 24 Oct 2026 09:47:35 -0400 Subject: [PATCH 1/4] Screen ellipsoid ode45 source convergence and independent motion --- .../workflows/ellipsoid-ode45-convergence.yml | 63 +++++++ tests/matlab/ellipsoid_ode45_convergence.m | 84 +++++++++ tests/test_ellipsoid_ode45_convergence.py | 166 ++++++++++++++++++ 3 files changed, 313 insertions(+) create mode 100644 .github/workflows/ellipsoid-ode45-convergence.yml create mode 100644 tests/matlab/ellipsoid_ode45_convergence.m create mode 100644 tests/test_ellipsoid_ode45_convergence.py diff --git a/.github/workflows/ellipsoid-ode45-convergence.yml b/.github/workflows/ellipsoid-ode45-convergence.yml new file mode 100644 index 00000000..af2f08ed --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-convergence.yml @@ -0,0 +1,63 @@ +name: Ellipsoid ode45 full-duration source convergence + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-convergence.yml' + - 'tests/matlab/ellipsoid_ode45_convergence.m' + - 'tests/test_ellipsoid_ode45_convergence.py' + - 'tests/fixtures/matlab/ellipsoid_ode45_convergence/**' + - 'wecsim/hydrodynamics/nonlinearHydro.py' + - 'wecsim/dynamics/caseDynamics.py' + workflow_dispatch: + +jobs: + compare: + strategy: + fail-fast: false + matrix: + wave: [Regular, RegularCIC] + runs-on: ubuntu-latest + timeout-minutes: 90 + name: ${{ matrix.wave }} + steps: + - uses: actions/checkout@v7 + - uses: actions/checkout@v7 + with: + repository: WEC-Sim/WEC-Sim + ref: 0753b2e47f2457c078751dcfe5d251d1767b80ab + path: matlab-ref + filter: blob:none + sparse-checkout: source + - uses: actions/checkout@v7 + with: + repository: WEC-Sim/WEC-Sim_Applications + ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e + path: applications + filter: blob:none + lfs: true + sparse-checkout: Nonlinear_Hydro + - uses: matlab-actions/setup-matlab@v3 + with: + release: R2025b + products: Simulink Simscape Simscape_Multibody + - uses: matlab-actions/run-command@v3 + with: + command: addpath('tests/matlab'); ellipsoid_ode45_convergence('${{ matrix.wave }}') + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e . pytest + - run: python -m pytest -q -s tests/test_ellipsoid_ode45_convergence.py + env: + WEC_SIM_APPLICATIONS_DIR: applications + WEC_SIM_ELLIPSOID_CONVERGENCE_DIR: matlab-ellipsoid-ode45-convergence + WEC_SIM_ELLIPSOID_CASE: ${{ matrix.wave }} + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-convergence-${{ matrix.wave }} + path: | + matlab-ellipsoid-ode45-convergence/${{ matrix.wave }}/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + applications/Nonlinear_Hydro/geometry/elipsoid.stl diff --git a/tests/matlab/ellipsoid_ode45_convergence.m b/tests/matlab/ellipsoid_ode45_convergence.m new file mode 100644 index 00000000..e77816c3 --- /dev/null +++ b/tests/matlab/ellipsoid_ode45_convergence.m @@ -0,0 +1,84 @@ +function ellipsoid_ode45_convergence(waveType) +% Screen the pinned ellipsoid ode45 application at successively finer clocks. +% Work only in copied application folders; never edit the pinned source input. +waveType = string(waveType); +assert(any(waveType == ["Regular", "RegularCIC"]), ... + 'Choose a pinned ellipsoid ode45 wave case'); +repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath')))); +addpath(genpath(fullfile(repoRoot, 'matlab-ref', 'source'))); +hydroDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', 'hydroData'); +cd(hydroDir); +if ~isfile('ellipsoid.h5') + bemio; +end +caseDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro', ... + 'ode45', waveType); +outDir = fullfile(repoRoot, 'matlab-ellipsoid-ode45-convergence', waveType); +mkdir(outDir); +steps = [0.005, 0.0025, 0.00125]; +labels = {'step005', 'step0025', 'step00125'}; + +for iStep = 1:numel(steps) + step = steps(iStep); + label = labels{iStep}; + targetDir = fullfile(fileparts(caseDir), waveType + "_" + label); + mkdir(targetDir); + copyfile(fullfile(caseDir, 'wecSimInputFile.m'), ... + fullfile(targetDir, 'wecSimInputFile.m')); + copyfile(fullfile(caseDir, 'ellipsoid.slx'), ... + fullfile(targetDir, 'ellipsoid.slx')); + inputFile = fullfile(targetDir, 'wecSimInputFile.m'); + contents = fileread(inputFile); + assert(contains(contents, 'simu.dt = 0.05;') && ... + ~contains(contents, 'simu.dtOut') && ... + ~contains(contents, 'simu.nonlinearDt') && ... + ~contains(contents, 'simu.cicDt'), ... + 'The pinned ellipsoid ode45 sample settings changed'); + contents = strrep(contents, 'simu.dt = 0.05;', ... + sprintf('simu.dt = %.5f;', step)); + fid = fopen(inputFile, 'w'); + assert(fid > 0); + cleanup = onCleanup(@() fclose(fid)); + fwrite(fid, contents); + clear cleanup; + + cd(targetDir); + wecSim; + maxStep = str2double(get_param('ellipsoid', 'MaxStep')); + assert(strcmp(simu.solver, 'ode45') && ... + abs(simu.dt - step) < 1e-12 && ... + abs(simu.dtOut - step) < 1e-12 && ... + abs(simu.nonlinearDt - step) < 1e-12 && ... + abs(simu.cicDt - step) < 1e-12 && ... + abs(maxStep - step) < 1e-12 && ... + simu.endTime == 150 && simu.rampTime == 50 && ... + strcmp(waves.type, lowerFirst(waveType)) && ... + waves.height == 4 && waves.period == 6 && ... + body(1).nonlinearHydro == 2 && pto(1).damping == 1200000, ... + 'The refined ellipsoid ode45 settings changed'); + n = round(simu.endTime / step) + 1; + bodyRecord = output.bodies(1); + ptoRecord = output.ptos(1); + waveRecord = output.wave; + values = [bodyRecord.time(:), bodyRecord.position(:, 3), ... + bodyRecord.velocity(:, 3), ptoRecord.position(:, 2), ... + ptoRecord.velocity(:, 2), ptoRecord.forceInternalMechanics(:, 2), ... + -ptoRecord.forceInternalMechanics(:, 2) .* ptoRecord.velocity(:, 2), ... + waveRecord.elevation(:), bodyRecord.forceRestoring(:, 3), ... + bodyRecord.forceExcitation(:, 3), bodyRecord.acceleration(:, 3)]; + assert(size(values, 1) == n && size(values, 2) == 11 && ... + all(isfinite(values), 'all') && ... + max(abs(bodyRecord.time(:) - ptoRecord.time(:))) < 1e-10 && ... + max(abs(bodyRecord.time(:) - waveRecord.time(:))) < 1e-10, ... + 'The refined ellipsoid output is incomplete'); + writematrix(values, fullfile(outDir, label + "_history.csv")); + writematrix([simu.dt, simu.dtOut, simu.nonlinearDt, ... + simu.cicDt, maxStep], fullfile(outDir, label + "_settings.csv")); + close_system('ellipsoid', 0); +end +end + +function value = lowerFirst(value) +value = char(value); +value(1) = lower(value(1)); +end diff --git a/tests/test_ellipsoid_ode45_convergence.py b/tests/test_ellipsoid_ode45_convergence.py new file mode 100644 index 00000000..28d81f3d --- /dev/null +++ b/tests/test_ellipsoid_ode45_convergence.py @@ -0,0 +1,166 @@ +"""Full-duration independent motion against a source-step-screened ode45 input. + +The unchanged 0.05 s published cases retain their separate diagnostic gates. +This test concerns copied, derived source inputs at 0.005/0.0025/0.00125 s. +""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint + + +ROOT = Path(__file__).resolve().parents[1] +APPLICATIONS = os.environ.get("WEC_SIM_APPLICATIONS_DIR") +RECORDS = os.environ.get("WEC_SIM_ELLIPSOID_CONVERGENCE_DIR") +CASE = os.environ.get("WEC_SIM_ELLIPSOID_CASE") +ARCHIVED_005 = (ROOT / "tests/fixtures/matlab/ellipsoid_ode45_maxstep_only" + / "refined_005_common_grid.npz") +pytestmark = pytest.mark.skipif( + not (APPLICATIONS and RECORDS and CASE), + reason="fresh pinned R2025b ode45 convergence output absent", +) + +# These original 150 s physical limits come from the established refined-step +# gate. They are deliberately unchanged; the source screen adds an independent +# requirement that the final source change be below one third of these limits. +PYTHON_LIMITS = { + "Regular": (.003, .004, 5000, 4000), + "RegularCIC": (.0035, .0045, 6000, 5000), +} +SOURCE_STABILITY_LIMITS = { + "Regular": (.001, .0013, 1600, 1300), + "RegularCIC": (.0011, .0015, 2000, 1600), +} +PYTHON_STEP_LIMITS = (.0001, .0001, 100, 100) +STEPS = {"step005": .005, "step0025": .0025, + "step00125": .00125} +# Columns: time, heave, speed, PTO stroke, stroke speed, PTO force, +# absorbed power, wave, restoring, excitation, acceleration. +CHANNELS = ((1, "heave position"), (2, "heave speed"), + (5, "PTO force"), (6, "absorbed power")) + + +def _error(a, b): + assert a.shape == b.shape + assert np.isfinite(a).all() and np.isfinite(b).all() + return float(np.max(np.abs(a - b))) + + +def _check(a, b, limit, label): + error = _error(a, b) + print(f"{label}: {error:.8g} (gate {limit:.8g})") + assert error < limit, f"{label}: {error:.8g} >= {limit:.8g}" + return error + + +def _read_source(case): + folder = Path(RECORDS) / case + result = {} + for label, step in STEPS.items(): + history = np.loadtxt(folder / f"{label}_history.csv", + delimiter=",", ndmin=2) + settings = np.loadtxt(folder / f"{label}_settings.csv", + delimiter=",", ndmin=2) + assert history.shape == (round(150 / step) + 1, 11), label + assert np.isfinite(history).all(), label + np.testing.assert_allclose(settings, np.full((1, 5), step), + rtol=0, atol=1e-12) + np.testing.assert_allclose(history[:, 0], + np.arange(len(history)) * step, + rtol=0, atol=1e-9) + result[label] = history + return result + + +def _python(case, step, app): + wec = WEC("source-step-screened ellipsoid ode45") + body = wec.body( + "ellipsoid", app / "hydroData/ellipsoid.h5", + mass="equilibrium", inertia=(1.375264e6, 1.375264e6, 1.341721e6), + geometry_file=app / "geometry/elipsoid.stl", + nonlinear_hydro="instantaneous", drag_coefficient=1, + drag_area=np.pi * 25, + ) + wec.coordinate("heave", body.move("heave")) + wec.pto("PTO1", WorldPoint(0, 0, -12.5), body.at(0, 0, 0), + damping=1_200_000) + wave = RegularCICWave(4, 6) if case == "RegularCIC" else RegularWave(4, 6) + return wec.run(wave, dt=step, end_time=150, ramp_time=50, + radiation_memory=60 if case == "RegularCIC" else None, + rho=1025) + + +def _python_channels(result): + body = result.bodies["ellipsoid"] + pto = result.ptos["PTO1"] + return (body.position[:, 2], body.velocity[:, 2], + pto.force, pto.absorbed_power) + + +def test_pinned_source_step_screen_and_full_duration_independent_motion(): + assert CASE in PYTHON_LIMITS + source = _read_source(CASE) + five, two_half, one_quarter = ( + source[label] for label in STEPS + ) + key = "regular_cic" if CASE == "RegularCIC" else "regular" + with np.load(ARCHIVED_005, allow_pickle=False) as archived: + for column, suffix, limit in ((0, "time", 1e-10), + (1, "heave", 1e-7), + (2, "speed", 1e-7), + (5, "pto_force", .1), + (7, "wave", 1e-10), + (8, "restoring", .1)): + _check(five[::10, column], archived[f"{key}_{suffix}"], + limit, f"regenerated 0.005 s {suffix}") + + _check(five[:, 7], two_half[::2, 7], 1e-10, + "source 0.005/0.0025 s wave identity") + _check(two_half[:, 7], one_quarter[::2, 7], 1e-10, + "source 0.0025/0.00125 s wave identity") + for (column, channel), stability in zip( + CHANNELS, SOURCE_STABILITY_LIMITS[CASE], strict=True + ): + previous = _error(five[:, column], two_half[::2, column]) + final = _check(two_half[:, column], one_quarter[::2, column], + stability, f"source final-step {channel}") + assert final < .7 * previous, ( + f"source {channel} did not contract: {previous:.8g} -> {final:.8g}" + ) + + app = Path(APPLICATIONS) / "Nonlinear_Hydro" + python_two_half = _python(CASE, .0025, app) + python_one_quarter = _python(CASE, .00125, app) + _check(python_two_half.time, python_one_quarter.time[::2], 1e-9, + "Python common times") + for (_, name), coarse, fine, limit in zip( + CHANNELS, _python_channels(python_two_half), + _python_channels(python_one_quarter), PYTHON_STEP_LIMITS, strict=True + ): + _check(coarse, fine[::2], limit, f"Python final-step {name}") + _check(python_one_quarter.time, one_quarter[:, 0], 1e-9, + "Python/source time") + _check(python_one_quarter.wave_elevation, one_quarter[:, 7], 1e-10, + "Python/source wave") + for (column, name), actual, limit in zip( + CHANNELS, _python_channels(python_one_quarter), + PYTHON_LIMITS[CASE], strict=True + ): + _check(actual, one_quarter[:, column], limit, + f"independent physical {name}") + _check(python_one_quarter.ptos["PTO1"].stroke, + one_quarter[:, 3], PYTHON_LIMITS[CASE][0], + "independent physical PTO stroke") + source_work = float(np.trapezoid(one_quarter[:, 6], one_quarter[:, 0])) + python_work = float(np.trapezoid( + python_one_quarter.ptos["PTO1"].absorbed_power, + python_one_quarter.time, + )) + assert source_work > 1e3 + work_error = abs(python_work / source_work - 1) + print(f"independent signed PTO work: {work_error:.6%} (gate 0.1%)") + assert work_error < .001 From 55973c00b560e500f8cb6ed0459a0f263c8ad1f9 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Sat, 24 Oct 2026 09:54:32 -0400 Subject: [PATCH 2/4] Extract the vertical PTO axis in ellipsoid source screen --- tests/matlab/ellipsoid_ode45_convergence.m | 6 +++--- tests/test_ellipsoid_ode45_convergence.py | 5 +++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/matlab/ellipsoid_ode45_convergence.m b/tests/matlab/ellipsoid_ode45_convergence.m index e77816c3..03d4df56 100644 --- a/tests/matlab/ellipsoid_ode45_convergence.m +++ b/tests/matlab/ellipsoid_ode45_convergence.m @@ -61,9 +61,9 @@ function ellipsoid_ode45_convergence(waveType) ptoRecord = output.ptos(1); waveRecord = output.wave; values = [bodyRecord.time(:), bodyRecord.position(:, 3), ... - bodyRecord.velocity(:, 3), ptoRecord.position(:, 2), ... - ptoRecord.velocity(:, 2), ptoRecord.forceInternalMechanics(:, 2), ... - -ptoRecord.forceInternalMechanics(:, 2) .* ptoRecord.velocity(:, 2), ... + bodyRecord.velocity(:, 3), ptoRecord.position(:, 3), ... + ptoRecord.velocity(:, 3), ptoRecord.forceInternalMechanics(:, 3), ... + -ptoRecord.forceInternalMechanics(:, 3) .* ptoRecord.velocity(:, 3), ... waveRecord.elevation(:), bodyRecord.forceRestoring(:, 3), ... bodyRecord.forceExcitation(:, 3), bodyRecord.acceleration(:, 3)]; assert(size(values, 1) == n && size(values, 2) == 11 && ... diff --git a/tests/test_ellipsoid_ode45_convergence.py b/tests/test_ellipsoid_ode45_convergence.py index 28d81f3d..8a1a6a8c 100644 --- a/tests/test_ellipsoid_ode45_convergence.py +++ b/tests/test_ellipsoid_ode45_convergence.py @@ -9,6 +9,7 @@ import numpy as np import pytest +from scipy.integrate import trapezoid from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint @@ -155,8 +156,8 @@ def test_pinned_source_step_screen_and_full_duration_independent_motion(): _check(python_one_quarter.ptos["PTO1"].stroke, one_quarter[:, 3], PYTHON_LIMITS[CASE][0], "independent physical PTO stroke") - source_work = float(np.trapezoid(one_quarter[:, 6], one_quarter[:, 0])) - python_work = float(np.trapezoid( + source_work = float(trapezoid(one_quarter[:, 6], one_quarter[:, 0])) + python_work = float(trapezoid( python_one_quarter.ptos["PTO1"].absorbed_power, python_one_quarter.time, )) From dc7c2fe2e73cf31c132272deee78ab77920233d6 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Sat, 24 Oct 2026 09:58:32 -0400 Subject: [PATCH 3/4] Advance ellipsoid source screen to the energy-stable step --- .../workflows/ellipsoid-ode45-convergence.yml | 2 +- tests/matlab/ellipsoid_ode45_convergence.m | 6 +- tests/test_ellipsoid_ode45_convergence.py | 48 ++++-- tools/archive_ellipsoid_ode45_convergence.py | 137 ++++++++++++++++++ 4 files changed, 179 insertions(+), 14 deletions(-) create mode 100644 tools/archive_ellipsoid_ode45_convergence.py diff --git a/.github/workflows/ellipsoid-ode45-convergence.yml b/.github/workflows/ellipsoid-ode45-convergence.yml index af2f08ed..182677c0 100644 --- a/.github/workflows/ellipsoid-ode45-convergence.yml +++ b/.github/workflows/ellipsoid-ode45-convergence.yml @@ -18,7 +18,7 @@ jobs: matrix: wave: [Regular, RegularCIC] runs-on: ubuntu-latest - timeout-minutes: 90 + timeout-minutes: 120 name: ${{ matrix.wave }} steps: - uses: actions/checkout@v7 diff --git a/tests/matlab/ellipsoid_ode45_convergence.m b/tests/matlab/ellipsoid_ode45_convergence.m index 03d4df56..84e01738 100644 --- a/tests/matlab/ellipsoid_ode45_convergence.m +++ b/tests/matlab/ellipsoid_ode45_convergence.m @@ -15,8 +15,8 @@ function ellipsoid_ode45_convergence(waveType) 'ode45', waveType); outDir = fullfile(repoRoot, 'matlab-ellipsoid-ode45-convergence', waveType); mkdir(outDir); -steps = [0.005, 0.0025, 0.00125]; -labels = {'step005', 'step0025', 'step00125'}; +steps = [0.005, 0.0025, 0.00125, 0.000625]; +labels = {'step005', 'step0025', 'step00125', 'step000625'}; for iStep = 1:numel(steps) step = steps(iStep); @@ -35,7 +35,7 @@ function ellipsoid_ode45_convergence(waveType) ~contains(contents, 'simu.cicDt'), ... 'The pinned ellipsoid ode45 sample settings changed'); contents = strrep(contents, 'simu.dt = 0.05;', ... - sprintf('simu.dt = %.5f;', step)); + sprintf('simu.dt = %.6f;', step)); fid = fopen(inputFile, 'w'); assert(fid > 0); cleanup = onCleanup(@() fclose(fid)); diff --git a/tests/test_ellipsoid_ode45_convergence.py b/tests/test_ellipsoid_ode45_convergence.py index 8a1a6a8c..52c76230 100644 --- a/tests/test_ellipsoid_ode45_convergence.py +++ b/tests/test_ellipsoid_ode45_convergence.py @@ -1,7 +1,8 @@ """Full-duration independent motion against a source-step-screened ode45 input. The unchanged 0.05 s published cases retain their separate diagnostic gates. -This test concerns copied, derived source inputs at 0.005/0.0025/0.00125 s. +This test concerns copied, derived source inputs at 0.005, 0.0025, 0.00125, +and 0.000625 s. """ import os @@ -12,6 +13,7 @@ from scipy.integrate import trapezoid from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro ROOT = Path(__file__).resolve().parents[1] @@ -38,7 +40,7 @@ } PYTHON_STEP_LIMITS = (.0001, .0001, 100, 100) STEPS = {"step005": .005, "step0025": .0025, - "step00125": .00125} + "step00125": .00125, "step000625": .000625} # Columns: time, heave, speed, PTO stroke, stroke speed, PTO force, # absorbed power, wave, restoring, excitation, acceleration. CHANNELS = ((1, "heave position"), (2, "heave speed"), @@ -105,7 +107,7 @@ def _python_channels(result): def test_pinned_source_step_screen_and_full_duration_independent_motion(): assert CASE in PYTHON_LIMITS source = _read_source(CASE) - five, two_half, one_quarter = ( + five, two_half, one_quarter, finest = ( source[label] for label in STEPS ) key = "regular_cic" if CASE == "RegularCIC" else "regular" @@ -123,17 +125,43 @@ def test_pinned_source_step_screen_and_full_duration_independent_motion(): "source 0.005/0.0025 s wave identity") _check(two_half[:, 7], one_quarter[::2, 7], 1e-10, "source 0.0025/0.00125 s wave identity") + _check(one_quarter[:, 7], finest[::2, 7], 1e-10, + "source 0.00125/0.000625 s wave identity") for (column, channel), stability in zip( CHANNELS, SOURCE_STABILITY_LIMITS[CASE], strict=True ): - previous = _error(five[:, column], two_half[::2, column]) - final = _check(two_half[:, column], one_quarter[::2, column], + first = _error(five[:, column], two_half[::2, column]) + previous = _error(two_half[:, column], one_quarter[::2, column]) + final = _check(one_quarter[:, column], finest[::2, column], stability, f"source final-step {channel}") + assert previous < .7 * first, ( + f"source {channel} did not contract: {first:.8g} -> {previous:.8g}" + ) assert final < .7 * previous, ( f"source {channel} did not contract: {previous:.8g} -> {final:.8g}" ) app = Path(APPLICATIONS) / "Nonlinear_Hydro" + mesh = HeaveMeshHydro.from_stl( + app / "geometry/elipsoid.stl", center_z=-2, rho=1025, + gravity=9.81, depth=70, period=6, height=4, ramp_time=50, + mass=None, drag_coefficient=1, drag_area=np.pi * 25, + ) + samples = np.linspace(1, len(finest) - 1, 1001, dtype=int) + previous = np.array([ + mesh.forces(finest[i - 1, 0], finest[i - 1, 1] + 2, + finest[i - 1, 2])[0] + for i in samples + ]) + current = np.array([ + mesh.forces(finest[i, 0], finest[i, 1] + 2, + finest[i, 2])[0] + for i in samples + ]) + _check(previous, -finest[samples, 8], 1e-6, + "finest source preceding-state restoring law") + assert _error(current, -finest[samples, 8]) > 500 + python_two_half = _python(CASE, .0025, app) python_one_quarter = _python(CASE, .00125, app) _check(python_two_half.time, python_one_quarter.time[::2], 1e-9, @@ -143,20 +171,20 @@ def test_pinned_source_step_screen_and_full_duration_independent_motion(): _python_channels(python_one_quarter), PYTHON_STEP_LIMITS, strict=True ): _check(coarse, fine[::2], limit, f"Python final-step {name}") - _check(python_one_quarter.time, one_quarter[:, 0], 1e-9, + _check(python_one_quarter.time, finest[::2, 0], 1e-9, "Python/source time") - _check(python_one_quarter.wave_elevation, one_quarter[:, 7], 1e-10, + _check(python_one_quarter.wave_elevation, finest[::2, 7], 1e-10, "Python/source wave") for (column, name), actual, limit in zip( CHANNELS, _python_channels(python_one_quarter), PYTHON_LIMITS[CASE], strict=True ): - _check(actual, one_quarter[:, column], limit, + _check(actual, finest[::2, column], limit, f"independent physical {name}") _check(python_one_quarter.ptos["PTO1"].stroke, - one_quarter[:, 3], PYTHON_LIMITS[CASE][0], + finest[::2, 3], PYTHON_LIMITS[CASE][0], "independent physical PTO stroke") - source_work = float(trapezoid(one_quarter[:, 6], one_quarter[:, 0])) + source_work = float(trapezoid(finest[:, 6], finest[:, 0])) python_work = float(trapezoid( python_one_quarter.ptos["PTO1"].absorbed_power, python_one_quarter.time, diff --git a/tools/archive_ellipsoid_ode45_convergence.py b/tools/archive_ellipsoid_ode45_convergence.py new file mode 100644 index 00000000..d3f63d43 --- /dev/null +++ b/tools/archive_ellipsoid_ode45_convergence.py @@ -0,0 +1,137 @@ +"""Archive full 150 s ellipsoid ode45 source-screen artifacts for portable CI. + +Example: + gh run download RUN_ID --dir /tmp/ellipsoid-source + python tools/archive_ellipsoid_ode45_convergence.py \ + /tmp/ellipsoid-source tests/fixtures/matlab/ellipsoid_ode45_convergence +""" + +import argparse +import hashlib +from pathlib import Path +import shutil + +import h5py +import numpy as np + + +ROOT = Path(__file__).resolve().parents[1] +ARCHIVED_005 = (ROOT / "tests/fixtures/matlab/ellipsoid_ode45_maxstep_only" + / "refined_005_common_grid.npz") +STEPS = {"step005": .005, "step0025": .0025, + "step00125": .00125, "step000625": .000625} +STABILITY_LIMITS = { + "Regular": (.001, .0013, 1600, 1300), + "RegularCIC": (.0011, .0015, 2000, 1600), +} +CHANNELS = (1, 2, 5, 6) + + +def _sha(path): + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def _maximum(a, b): + assert a.shape == b.shape + assert np.isfinite(a).all() and np.isfinite(b).all() + return float(np.max(np.abs(a - b))) + + +def _datasets(path): + result = {} + with h5py.File(path) as source: + def collect(name, obj): + if isinstance(obj, h5py.Dataset): + result[name] = obj[()] + source.visititems(collect) + return result + + +def _read_case(artifacts, case): + base = artifacts / f"ellipsoid-ode45-convergence-{case}" + output = base / "matlab-ellipsoid-ode45-convergence" / case + histories = {} + for label, step in STEPS.items(): + history = np.loadtxt(output / f"{label}_history.csv", + delimiter=",", ndmin=2) + settings = np.loadtxt(output / f"{label}_settings.csv", + delimiter=",", ndmin=2) + assert history.shape == (round(150 / step) + 1, 11) + assert np.isfinite(history).all() + np.testing.assert_allclose(settings, np.full((1, 5), step), + rtol=0, atol=1e-12) + np.testing.assert_allclose(history[:, 0], + np.arange(len(history)) * step, + rtol=0, atol=1e-9) + histories[label] = history + return base, histories + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("artifacts", type=Path) + parser.add_argument("output", type=Path) + args = parser.parse_args() + records = {} + inputs = {} + with np.load(ARCHIVED_005, allow_pickle=False) as archived: + for case in ("Regular", "RegularCIC"): + base, histories = _read_case(args.artifacts, case) + key = "regular_cic" if case == "RegularCIC" else "regular" + five = histories["step005"] + for column, suffix, limit in ((0, "time", 1e-10), + (1, "heave", 1e-7), + (2, "speed", 1e-7), + (5, "pto_force", .1), + (7, "wave", 1e-10), + (8, "restoring", .1)): + assert _maximum(five[::10, column], + archived[f"{key}_{suffix}"]) < limit + half = histories["step0025"] + quarter = histories["step00125"] + finest = histories["step000625"] + assert _maximum(five[:, 7], half[::2, 7]) < 1e-10 + assert _maximum(half[:, 7], quarter[::2, 7]) < 1e-10 + assert _maximum(quarter[:, 7], finest[::2, 7]) < 1e-10 + for column, limit in zip(CHANNELS, STABILITY_LIMITS[case], + strict=True): + first = _maximum(five[:, column], half[::2, column]) + previous = _maximum(half[:, column], quarter[::2, column]) + final = _maximum(quarter[:, column], finest[::2, column]) + print(case, column, "source changes", first, previous, final) + assert (final < limit and previous < .7 * first + and final < .7 * previous) + for label, history in histories.items(): + records[f"{key}_{label}"] = history + inputs[case] = { + "hydro": base / "applications/Nonlinear_Hydro/hydroData/ellipsoid.h5", + "geometry": base / "applications/Nonlinear_Hydro/geometry/elipsoid.stl", + } + + assert _sha(inputs["Regular"]["geometry"]) == _sha( + inputs["RegularCIC"]["geometry"]) + regular_hydro = _datasets(inputs["Regular"]["hydro"]) + cic_hydro = _datasets(inputs["RegularCIC"]["hydro"]) + assert regular_hydro.keys() == cic_hydro.keys() + for name, values in regular_hydro.items(): + other = cic_hydro[name] + if np.issubdtype(np.asarray(values).dtype, np.number): + np.testing.assert_allclose(values, other, rtol=0, atol=0, + equal_nan=True, err_msg=name) + else: + np.testing.assert_array_equal(values, other, err_msg=name) + + args.output.mkdir(parents=True, exist_ok=True) + destination = args.output / "source.npz" + np.savez_compressed(destination, **records) + shutil.copyfile(inputs["Regular"]["hydro"], + args.output / "ellipsoid.h5") + shutil.copyfile(inputs["Regular"]["geometry"], + args.output / "elipsoid.stl") + for file in (destination, args.output / "ellipsoid.h5", + args.output / "elipsoid.stl"): + print(file, _sha(file), file.stat().st_size) + + +if __name__ == "__main__": + main() From c8a4e3f41ad5a30bd77545944caf9817f4efe129 Mon Sep 17 00:00:00 2001 From: Chris McComb Date: Sat, 24 Oct 2026 11:08:10 -0400 Subject: [PATCH 4/4] Run finest ellipsoid source step in isolated MATLAB process --- .../ellipsoid-ode45-finest-isolated.yml | 45 +++++++++++++++++++ tests/matlab/ellipsoid_ode45_convergence.m | 10 ++++- 2 files changed, 54 insertions(+), 1 deletion(-) create mode 100644 .github/workflows/ellipsoid-ode45-finest-isolated.yml diff --git a/.github/workflows/ellipsoid-ode45-finest-isolated.yml b/.github/workflows/ellipsoid-ode45-finest-isolated.yml new file mode 100644 index 00000000..288ecfb3 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-finest-isolated.yml @@ -0,0 +1,45 @@ +name: Ellipsoid ode45 isolated finest source step + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-finest-isolated.yml' + - 'tests/matlab/ellipsoid_ode45_convergence.m' + workflow_dispatch: + +jobs: + regular: + runs-on: ubuntu-latest + timeout-minutes: 120 + steps: + - uses: actions/checkout@v7 + - uses: actions/checkout@v7 + with: + repository: WEC-Sim/WEC-Sim + ref: 0753b2e47f2457c078751dcfe5d251d1767b80ab + path: matlab-ref + filter: blob:none + sparse-checkout: source + - uses: actions/checkout@v7 + with: + repository: WEC-Sim/WEC-Sim_Applications + ref: d53d4d4c9eda2581f04204f5d394a6ef84bb099e + path: applications + filter: blob:none + lfs: true + sparse-checkout: Nonlinear_Hydro + - uses: matlab-actions/setup-matlab@v3 + with: + release: R2025b + products: Simulink Simscape Simscape_Multibody + - uses: matlab-actions/run-command@v3 + with: + command: addpath('tests/matlab'); ellipsoid_ode45_convergence('Regular', 'step000625') + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-finest-isolated-Regular + path: | + matlab-ellipsoid-ode45-convergence/Regular/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + applications/Nonlinear_Hydro/geometry/elipsoid.stl diff --git a/tests/matlab/ellipsoid_ode45_convergence.m b/tests/matlab/ellipsoid_ode45_convergence.m index 84e01738..ced2b8a6 100644 --- a/tests/matlab/ellipsoid_ode45_convergence.m +++ b/tests/matlab/ellipsoid_ode45_convergence.m @@ -1,4 +1,4 @@ -function ellipsoid_ode45_convergence(waveType) +function ellipsoid_ode45_convergence(waveType, selectedLabel) % Screen the pinned ellipsoid ode45 application at successively finer clocks. % Work only in copied application folders; never edit the pinned source input. waveType = string(waveType); @@ -17,6 +17,14 @@ function ellipsoid_ode45_convergence(waveType) mkdir(outDir); steps = [0.005, 0.0025, 0.00125, 0.000625]; labels = {'step005', 'step0025', 'step00125', 'step000625'}; +if nargin >= 2 + % Run a single step in a fresh MATLAB process to isolate memory use. + selectedLabel = string(selectedLabel); + selection = strcmp(labels, selectedLabel); + assert(nnz(selection) == 1, 'Choose a supported source step'); + steps = steps(selection); + labels = labels(selection); +end for iStep = 1:numel(steps) step = steps(iStep);