Isolate finest ellipsoid ode45 source step

.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: 120
+    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

.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

tests/matlab/ellipsoid_ode45_convergence.m

@@ -0,0 +1,92 @@
+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);
+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, 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);
+    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 = %.6f;', 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(:, 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 && ...
+        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

tests/test_ellipsoid_ode45_convergence.py

@@ -0,0 +1,195 @@
+"""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,
+and 0.000625 s.
+"""
+
+import os
+from pathlib import Path
+
+import numpy as np
+import pytest
+from scipy.integrate import trapezoid
+
+from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint
+from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro
+
+
+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, "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"),
+            (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, finest = (
+        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")
+    _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
+    ):
+        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,
+           "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, finest[::2, 0], 1e-9,
+           "Python/source time")
+    _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, finest[::2, column], limit,
+               f"independent physical {name}")
+    _check(python_one_quarter.ptos["PTO1"].stroke,
+           finest[::2, 3], PYTHON_LIMITS[CASE][0],
+           "independent physical PTO stroke")
+    source_work = float(trapezoid(finest[:, 6], finest[:, 0]))
+    python_work = float(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