Screen finest ellipsoid source step with coarser output

.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-output.yml

@@ -0,0 +1,45 @@
+name: Ellipsoid ode45 finest physical step with coarser output
+
+on:
+  pull_request:
+    paths:
+      - '.github/workflows/ellipsoid-ode45-finest-output.yml'
+      - 'tests/matlab/ellipsoid_ode45_finest_output.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_finest_output('Regular')
+      - uses: actions/upload-artifact@v4
+        if: always()
+        with:
+          name: ellipsoid-ode45-finest-output-Regular
+          path: |
+            matlab-ellipsoid-ode45-finest-output/Regular/
+            applications/Nonlinear_Hydro/hydroData/ellipsoid.h5
+            applications/Nonlinear_Hydro/geometry/elipsoid.stl

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, 0.000625];
+labels = {'step005', 'step0025', 'step00125', 'step000625'};
+
+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/matlab/ellipsoid_ode45_finest_output.m

@@ -0,0 +1,98 @@
+function ellipsoid_ode45_finest_output(waveType)
+% Diagnose the finest physical clock with a coarser output clock.
+% The 0.00125 s A/B established output neutrality only at that step; this
+% 0.000625 s run needs its own source and independent physical gates.
+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-finest-output', waveType);
+mkdir(outDir);
+steps = 0.000625;
+outputStep = 0.005;
+labels = {'step000625'};
+
+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;\nsimu.dtOut = %.6f;', ...
+        step, outputStep));
+    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 - outputStep) < 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 finest mixed-output ellipsoid ode45 settings changed');
+    n = round(simu.endTime / outputStep) + 1;
+    ratio = round(outputStep / step);
+    assert(abs(ratio * step - outputStep) < 1e-12);
+    bodyRecord = output.bodies(1);
+    ptoRecord = output.ptos(1);
+    waveRecord = output.wave;
+    if numel(waveRecord.time) == n
+        waveTime = waveRecord.time;
+        waveElevation = waveRecord.elevation;
+    else
+        assert(numel(waveRecord.time) == round(simu.endTime / step) + 1);
+        waveTime = waveRecord.time(1:ratio:end);
+        waveElevation = waveRecord.elevation(1:ratio:end);
+    end
+    values = [bodyRecord.time(:), bodyRecord.position(:, 3), ...
+        bodyRecord.velocity(:, 3), ptoRecord.position(:, 3), ...
+        ptoRecord.velocity(:, 3), ptoRecord.forceInternalMechanics(:, 3), ...
+        -ptoRecord.forceInternalMechanics(:, 3) .* ptoRecord.velocity(:, 3), ...
+        waveElevation(:), 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 && ...
+        numel(waveTime) == n && ...
+        max(abs(bodyRecord.time(:) - waveTime(:))) < 1e-10, ...
+        'The finest mixed-output ellipsoid history 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