diff --git a/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-portable.yml b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-portable.yml new file mode 100644 index 00000000..b26d5a3b --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-portable.yml @@ -0,0 +1,28 @@ +name: Ellipsoid ode45 continuous nonlinear mass-control portable gate + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-portable.yml' + - 'tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py' + - 'tests/test_ellipsoid_ode45_convergence.py' + - 'tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_convergence/**' + - 'tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py' + - 'wecsim/hydrodynamics/nonlinearHydro.py' + - 'wecsim/dynamics/caseDynamics.py' + - 'pyproject.toml' + workflow_dispatch: + +jobs: + compare: + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e ".[test]" + - run: python -m pytest -q -s tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py diff --git a/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step-portable.yml b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step-portable.yml new file mode 100644 index 00000000..43dc536e --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step-portable.yml @@ -0,0 +1,33 @@ +name: Ellipsoid ode45 no-delay MaxStep portable diagnostic + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step-portable.yml' + - '.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step.yml' + - 'tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py' + - 'tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py' + - 'tests/test_ellipsoid_ode45_convergence.py' + - 'tests/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m' + - 'tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_restoring_step/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_convergence/**' + - 'tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py' + - 'wecsim/hydrodynamics/nonlinearHydro.py' + - 'wecsim/dynamics/caseDynamics.py' + - 'pyproject.toml' + workflow_dispatch: + +jobs: + compare: + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e ".[test]" + - run: python -m pytest -q -s tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py diff --git a/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step.yml b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step.yml new file mode 100644 index 00000000..9dd8cdf1 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step.yml @@ -0,0 +1,51 @@ +name: Ellipsoid ode45 continuous nonlinear mass-control MaxStep diagnosis + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step.yml' + - 'tests/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m' + workflow_dispatch: + +jobs: + source-step: + strategy: + fail-fast: false + matrix: + wave: [Regular, RegularCIC] + max_step: ['0.0125', '0.00625'] + runs-on: ubuntu-latest + timeout-minutes: 150 + name: ${{ matrix.wave }} MaxStep ${{ matrix.max_step }} + 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'); waveCase='${{ matrix.wave }}'; maxStep=${{ matrix.max_step }}; ellipsoid_ode45_continuous_nonlinear_mass_control + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-continuous-nonlinear-mass-control-${{ matrix.wave }}-maxstep-${{ matrix.max_step }} + path: | + matlab-ellipsoid-ode45-continuous-nonlinear-mass-control-step/${{ matrix.wave }}/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + if-no-files-found: ignore diff --git a/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control.yml b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control.yml new file mode 100644 index 00000000..cd490068 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control.yml @@ -0,0 +1,50 @@ +name: Ellipsoid ode45 continuous nonlinear mass control + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control.yml' + - 'tests/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m' + workflow_dispatch: + +jobs: + source-control: + strategy: + fail-fast: false + matrix: + wave: [Regular, RegularCIC] + runs-on: ubuntu-latest + timeout-minutes: 150 + 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'); waveCase='${{ matrix.wave }}'; ellipsoid_ode45_continuous_nonlinear_mass_control + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-continuous-nonlinear-mass-control-${{ matrix.wave }} + path: | + matlab-ellipsoid-ode45-continuous-nonlinear-mass-control/${{ matrix.wave }}/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + if-no-files-found: ignore diff --git a/.github/workflows/ellipsoid-ode45-restoring-step-portable.yml b/.github/workflows/ellipsoid-ode45-restoring-step-portable.yml new file mode 100644 index 00000000..de71e7ce --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-restoring-step-portable.yml @@ -0,0 +1,30 @@ +name: Ellipsoid ode45 continuous nonlinear MaxStep portable diagnostic + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-restoring-step-portable.yml' + - 'tests/test_ellipsoid_ode45_restoring_step_portable.py' + - 'tests/test_ellipsoid_ode45_convergence.py' + - 'tools/check_ellipsoid_ode45_restoring_step.py' + - 'tools/pack_ellipsoid_ode45_restoring_step.py' + - 'tools/pack_ellipsoid_ode45_restoring_timing.py' + - 'tests/fixtures/matlab/ellipsoid_ode45_restoring_step/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_convergence/**' + - 'wecsim/hydrodynamics/nonlinearHydro.py' + - 'wecsim/dynamics/caseDynamics.py' + - 'pyproject.toml' + workflow_dispatch: + +jobs: + compare: + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e ".[test]" + - run: python -m pytest -q -s tests/test_ellipsoid_ode45_restoring_step_portable.py diff --git a/.github/workflows/ellipsoid-ode45-restoring-step-screen.yml b/.github/workflows/ellipsoid-ode45-restoring-step-screen.yml new file mode 100644 index 00000000..b90bd7e6 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-restoring-step-screen.yml @@ -0,0 +1,52 @@ +name: Ellipsoid ode45 continuous nonlinear MaxStep screen + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-restoring-step-screen.yml' + - 'tests/matlab/ellipsoid_ode45_restoring_timing_smoke.m' + workflow_dispatch: + +jobs: + source-step: + strategy: + fail-fast: false + matrix: + wave: [Regular, RegularCIC] + # The two finer values continue diagnosis after the predeclared screen. + max_step: ['0.025', '0.0125', '0.00625', '0.003125'] + runs-on: ubuntu-latest + timeout-minutes: 150 + name: ${{ matrix.wave }} MaxStep ${{ matrix.max_step }} + 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'); waveCase='${{ matrix.wave }}'; mode='continuous_nonlinear'; duration='full'; maxStep=${{ matrix.max_step }}; ellipsoid_ode45_restoring_timing_smoke + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-restoring-${{ matrix.wave }}-maxstep-${{ matrix.max_step }} + path: | + matlab-ellipsoid-ode45-restoring-timing-step/${{ matrix.wave }}/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + if-no-files-found: ignore diff --git a/.github/workflows/ellipsoid-ode45-restoring-timing-portable.yml b/.github/workflows/ellipsoid-ode45-restoring-timing-portable.yml new file mode 100644 index 00000000..f890e907 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-restoring-timing-portable.yml @@ -0,0 +1,26 @@ +name: Ellipsoid ode45 force-timing portable gate + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-restoring-timing-portable.yml' + - 'tests/test_ellipsoid_ode45_restoring_timing_portable.py' + - 'tests/test_ellipsoid_ode45_convergence.py' + - 'tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_convergence/**' + - 'tests/fixtures/matlab/ellipsoid_ode45_nonlinear_clock/**' + - 'wecsim/hydrodynamics/nonlinearHydro.py' + - 'wecsim/dynamics/caseDynamics.py' + workflow_dispatch: + +jobs: + compare: + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e ".[test]" + - run: python -m pytest -q -s tests/test_ellipsoid_ode45_restoring_timing_portable.py diff --git a/.github/workflows/ellipsoid-ode45-restoring-timing-smoke.yml b/.github/workflows/ellipsoid-ode45-restoring-timing-smoke.yml new file mode 100644 index 00000000..54904872 --- /dev/null +++ b/.github/workflows/ellipsoid-ode45-restoring-timing-smoke.yml @@ -0,0 +1,77 @@ +name: Ellipsoid ode45 restoring timing diagnostic + +on: + pull_request: + paths: + - '.github/workflows/ellipsoid-ode45-restoring-timing-smoke.yml' + - 'tests/matlab/ellipsoid_ode45_restoring_timing_smoke.m' + - 'tools/check_ellipsoid_ode45_restoring_timing.py' + workflow_dispatch: + inputs: + duration: + description: Source duration (2 s smoke or 150 s full) + required: true + default: smoke + type: choice + options: + - smoke + - full + +jobs: + compare: + strategy: + fail-fast: false + matrix: + wave: [Regular, RegularCIC] + runs-on: ubuntu-latest + timeout-minutes: 150 + 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'); waveCase='${{ matrix.wave }}'; mode='baseline'; duration='${{ inputs.duration || 'smoke' }}'; ellipsoid_ode45_restoring_timing_smoke + - uses: matlab-actions/run-command@v3 + with: + command: addpath('tests/matlab'); waveCase='${{ matrix.wave }}'; mode='continuous_restoring'; duration='${{ inputs.duration || 'smoke' }}'; ellipsoid_ode45_restoring_timing_smoke + - uses: matlab-actions/run-command@v3 + with: + command: addpath('tests/matlab'); waveCase='${{ matrix.wave }}'; mode='continuous_nonlinear'; duration='${{ inputs.duration || 'smoke' }}'; ellipsoid_ode45_restoring_timing_smoke + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: python -m pip install -e ".[test]" + - name: Check restoring timing and source comparison + run: >- + python tools/check_ellipsoid_ode45_restoring_timing.py + "${{ matrix.wave }}" + "matlab-ellipsoid-ode45-restoring-timing/${{ matrix.wave }}" + --duration "${{ inputs.duration || 'smoke' }}" + --applications applications + - uses: actions/upload-artifact@v4 + if: always() + with: + name: ellipsoid-ode45-restoring-timing-${{ matrix.wave }}-${{ inputs.duration || 'smoke' }} + path: | + matlab-ellipsoid-ode45-restoring-timing/${{ matrix.wave }}/ + applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 + if-no-files-found: ignore diff --git a/PARITY.md b/PARITY.md index 1a05bdc8..dde6488b 100644 --- a/PARITY.md +++ b/PARITY.md @@ -119,6 +119,7 @@ saves its MATLAB outputs as an artifact. | Ellipsoid `ode45` MaxStep-only isolation | Pinned 150 s `Regular` and `RegularCIC` source inputs, a local R2025b run of each at the published 0.05 s settings and again with only Simulink `MaxStep=0.005`, and the [archived jointly refined 0.005 s source](https://github.com/cmudrc/wec-sim-python/actions/runs/38110743157) | Each pair has identical wave histories and retains `dtOut=nonlinearDt=cicDt=0.05` s. MaxStep-only changes published heave by 6.19 mm (Regular) and 5.09 mm (RegularCIC), but its gaps to the jointly refined source are 21.24 and 21.39 mm, versus 17.09 and 17.92 mm for the published source. An independent public Python run at 0.05 s is only 1.85 and 2.02 mm from the jointly refined source, but 18.94 and 19.93 mm from published MATLAB and 23.09 and 23.40 mm from MaxStep-only MATLAB. Both source restoring histories still reconstruct from the preceding 0.05 s output state within `1e-6` N. [Fresh paired gates](tests/test_ellipsoid_ode45_maxstep_only.py) verify source settings, wave identity, restoring timing, and independent Python motion. Reducing solver MaxStep alone does not remove the published trajectory gap; this experiment does not separately isolate output, nonlinear-force, or CIC sample timing and does not promote either published case. | | Ellipsoid `ode45` nonlinear-force and CIC clock isolation | Pinned 150 s `Regular` and `RegularCIC` inputs, the archived published and MaxStep-only source traces above, and new local R2025b arms with `dt=MaxStep=0.005` and `dtOut=0.05` s | With wave and solver unchanged, changing only `nonlinearDt` from 0.005 to the published 0.05 s moves source heave by 21.24 mm (Regular) and 21.41 mm (RegularCIC), heave speed by 26.72 and 28.96 mm/s, and restoring force by 43.81 and 46.11 kN. Regular's coarse-force arm exactly reproduces the earlier MaxStep-only heave, speed, restoring, excitation, and PTO histories at all 3,001 saved samples. For RegularCIC, independently restoring `cicDt` from 0.005 to 0.05 s shifts heave a further 0.112 mm and then reproduces MaxStep-only heave and speed within `1.1e-14` and `1.1e-15` in their units; force differences are below `1.1e-9` N. Independent public Python remains 1.85 and 2.02 mm from the fine-clock source, versus 23.09 and 23.42 mm from its coarse-force arm. [Paired clock gates](tests/test_ellipsoid_ode45_nonlinear_clock.py) check clock settings, wave identity, source reconstruction, and independent motion. These controlled derived arms attribute the MaxStep-only-to-refined gap to source sample clocks for the pinned inputs; neither published 0.05 s trajectory is promoted. | | Ellipsoid `ode45/Regular` added-mass feedback | The pinned published input, an independent local R2025b 2 s log of the active `10e-8` s Transport Delay, and the [full 150 s published force output](https://github.com/cmudrc/wec-sim-python/actions/runs/39054809363) | The logged delay input equals source heave acceleration at all 41 prefix samples; its output follows two-sample extrapolation to `t−1e-7` s within `1.2e-17` m/s². Across all 3,001 published samples, the logged added-mass force follows that same delay law within `2.1e-8` N. Independently rebuilding the complete right-hand side from Python's bundled mesh and BEM coefficients, source motion, and logged added-mass force reproduces source acceleration within `5e-13` m/s². The [portable gate](tests/test_ellipsoid_ode45_delay_prefix.py) and [replay script](tests/matlab/ellipsoid_ode45_delay_prefix.m) pin this force law. Holding mesh forces on Python's 0.05 s output clock alone reduced the published heave gap from 18.94 to 15.29 mm in a local prototype, so that option was discarded. This source-state force gate is not an independent Python trajectory gate; the published ode45 case remains qualified. | +| Ellipsoid `ode45` copied-model force timing and heave mass controls | Pinned 150 s `Regular` and `RegularCIC` MATLAB Applications inputs; [archived timing arms](tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/README.md), [copied-model MaxStep continuation](tests/fixtures/matlab/ellipsoid_ode45_restoring_step/README.md), and [matched no-active-delay controls](tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/README.md) | An unchanged copied baseline reproduces the published traces. Making only the restoring block continuous changes its logged force from the preceding to the current output state and cuts independent Python heave gaps from 18.940/19.930 to 5.156/6.390 mm. Making the nonlinear restoring, Froude–Krylov, and wave-elevation blocks continuous cuts them to 3.421/4.307 mm at the original 0.05 s MaxStep, but signed PTO work still differs by 0.880/1.114% (inherited limit 0.1%). Saved-sample excitation and incident wave reconstruct from the current source state within `3.5e-7` N and `1.1e-13` m; this does not establish timing at solver minor steps. With original 0.05 s input/output/nonlinear/CIC clocks held fixed, reducing only the copied model MaxStep through 0.025, 0.0125, 0.00625, and 0.003125 s makes successive source changes contract by roughly half. The **predeclared 0.0125 s source-stability screen fails**. At the final finer diagnostic step, source signed-work changes are 0.063/0.071%, independent Python signed-work gaps are 0.059/0.070%, and all inherited pointwise physical limits pass. Separately, an exact heave mass shift leaves less than `6e-22` N of active delayed heave force at 0.05 s MaxStep while preserving continuous nonlinear blocks; it worsens Python heave gaps to 7.978/8.056 mm and signed-work gaps to 2.108/2.169%. Refining that control to 0.00625 s reduces its heave gaps to 0.956/1.122 mm, but work still differs by 0.269/0.262%, worse than the paired delayed arms. The mass shift also changes inactive inertial masses and solver representation, so it does not uniquely identify the remaining numerical mechanism. Every edited arm is a **derived diagnostic**. The original published inputs remain `qualified_motion`, the strict gate and badges are unchanged, and [#624](https://github.com/cmudrc/wec-sim-python/issues/624) stays unchecked for both cases. | | Sphere variable draft and mass in heave | Pinned MATLAB Applications `Variable_Hydro/Variable_Mass`, nine BEMIO-generated draft HDF5 files, and its full 900 s ode4 run | The Python `WEC.variable_body` selects the published state every 100 s and changes equilibrium draft, rigid mass, added mass, radiation damping, restoring stiffness, and excitation together. On MATLAB's 9,001 saved states, state-selected excitation differs by at most `3e-8` N, radiation damping by `3e-10` N, restoring plus weight by `8e-9` N, and assembled hydrodynamic force by `9e-9` N. The independent Python trajectory differs by at most 2.23 mm heave, 4.14 mm/s velocity, and 827 N PTO force over all nine states. Its newly paired PTO speed and absorbed power differ by at most 4.135 mm/s and 1.052 kW; full-run damper work differs by 0.00376%, below the fixed 1% gate. The source reaches −6.63 m by 900 s; agreement with its linear coefficients does not establish physical accuracy at that excursion. Other variable-hydrodynamic motions remain unsupported. | | Barge generalized body modes | Pinned MATLAB Applications `Generalized_Body_Modes`, its BEMIO-generated ten-DOF barge HDF5, and a fresh 400 s R2025b run including the model's separate `Flex_out` signal | `BodyClass` pairs flexible effective mass, stiffness, damping, hydrostatic, added-mass, radiation-damping, and excitation matrices; the largest absolute matrix difference is `1.91e-6` in applied added mass. The source flexible acceleration equation and six signed force components reconstruct from those matrices, with restoring and radiation differences below `1e-6` N, and sampled added-mass feedback below `0.05` N. The public `WEC.floating_gbm` interface runs a coupled three-rigid/four-flexible-coordinate implicit-mass model in the published 2 m, 8 s regular wave. Over all 8,001 samples, maximum rigid position and velocity differences are 1.06 mm and 0.838 mm/s; flexible displacement, velocity, and acceleration differ by at most 0.043 mm, 0.034 mm/s, and `1.26e-4` m/s². This validates one single-body floating 3-DOF joint with four modes, zero-heading regular waves, and no PTO or mooring; other GBM layouts and wave/radiation options remain unsupported. | | OWC `OrificeModel` wave, force, and early coupled motion | Pinned MATLAB Applications `OWC/OrificeModel`, BEMIO-cleaned seven-DOF HDF5, and a fresh 130 s R2025b `ode23t` [source run](https://github.com/cmudrc/wec-sim-python/actions/runs/37799665808) | `wecsim.OrificePTO` reproduces the embedded Simulink orifice block on its 418,828 logged internal solver evaluations. Maximum force, pressure, and source-power differences are `4.8e-10` N, `2.9e-12` kPa, and `4.2e-13` kW; the Mach flag is identical. Interpolating source flow to the independent 26,001-sample flexible-mode output gives a piston-speed difference below `3.6e-13` m/s. The 500 PM bins agree within `5.2e-15`, wave elevation within `2.8e-14` m, and all six rigid plus one flexible excitation channels within `1.5e-10` N. Flexible mass/stiffness/damping, added mass, and hydrostatic matrices agree within `2.3e-13` in their units. With the [rigid-force diagnostic run](https://github.com/cmudrc/wec-sim-python/actions/runs/37804483371), radiation reconstructed from all seven MATLAB velocities differs by at most `0.0074` N in the flexible channel and `0.22` N across rigid channels (peak rigid magnitude `1,182` N). The rigid restoring and linear-damping forces reproduce the source logs within `8.6e-11` N and `1.1e-12` N, and the rigid hydrodynamic force sum closes within `1.7e-10` N. MATLAB's flexible hydrodynamic force sum closes within `1e-12` N; its effective-mass acceleration including the orifice force closes within `3.5e-8` N. The source rigid surge, heave, and pitch balances close within `1.2e-10` in their force units when using its adjusted translational mass and logged piston reaction. `WEC.floating_gbm(..., orifice=...)` now integrates the seven-DOF PM case with memory radiation, quadratic drag, and the coupled piston force. Against the saved source phases, maximum surge, heave, unwrapped-pitch, and flexible-position differences through 6 s are `0.0113` m, `0.0244` m, `0.0098` rad, and `0.0070` m. At 130 s the physical coupled path has heave and flexible-displacement differences of `1.94` m and `0.280` m. The separate opt-in `orifice_force_path="published_owc"` follows the pinned Simulink routing, which leaves the piston out of the flexible state-space input but retains it in reported acceleration and rigid added-mass feedback. Over all 26,001 samples, its maximum surge, heave, unwrapped-pitch, and flexible-position differences are `0.174` m, `0.0367` m, `0.0108` rad, and `0.0054` m; PTO reaction and absorbed-power maxima differ by `1.98` kN and `1.75` kW, and integrated absorbed work by `4.76%`. These are explicit 130 s source-convention gates, **not physical coupled-path parity**. Python reports SI watts while the source's numeric power is kW. The source first exceeds its Mach threshold at 10.13 s; it is above 0.3 for 61.4% of simulated time and peaks at Mach 1.11. Matching its incompressible law after that point does not establish physical accuracy. | diff --git a/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/README.md b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/README.md new file mode 100644 index 00000000..cbcb3c35 --- /dev/null +++ b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/README.md @@ -0,0 +1,55 @@ +# Ellipsoid ode45 continuous nonlinear, no active heave delay control + +This is a diagnostic for [issue #599](https://github.com/cmudrc/wec-sim-python/issues/599), separate from the unchanged published-input parity gate. The [R2025b source run 40130836521](https://github.com/cmudrc/wec-sim-python/actions/runs/40130836521) completed both 150 s wave cases at repository commit `87d27eb62b0d13d84247e113f22bae5d1fd60c03`. It pinned WEC-Sim core `0753b2e47f2457c078751dcfe5d251d1767b80ab` and Applications `d53d4d4c9eda2581f04204f5d394a6ef84bb099e`. The [MATLAB script](../../../matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m) copies each application input and model into a disposable sibling directory. The original upstream inputs and models remain untouched. + +The control retains `ode45`, 150 s duration, 50 s ramp, 0.05 s solver maximum/output/nonlinear/CIC clocks, the 4 m/6 s wave, mesh, hydro data, PTO, and the 100 ns Simulink Transport Delay. Like the [archived continuous-nonlinear source arm](../ellipsoid_ode45_restoring_timing/README.md), its three nonlinear blocks compile at a continuous sample period. The sole additional input change sets `adjMassFactor = A_zz / trace(A_xyz)` from the same pinned HDF5. In this heave-only model, WEC-Sim moves `A_zz` into the rigid inertial mass and leaves zero active heave acceleration coefficient entering the delay. + +| Source mass identity | Regular | RegularCIC | +| --- | ---: | ---: | +| `adjMassFactor` | 0.655990137167053 | 0.783566018955714 | +| Original heave rigid mass | 133,376.067 kg | 133,376.067 kg | +| Heave added mass `A_zz` | 199,883.891 kg | 149,463.860 kg | +| Adjusted inertial mass `m + A_zz` | 333,259.957 kg | 282,839.927 kg | +| Residual active delayed coefficient | 0 kg | 0 kg | +| Maximum actual applied heave added-mass force | `5.99e-22` N | `3.69e-22` N | + +All 3,001 samples share the same incident wave as the archived source arm and independently run Python within `1e-10` m. The active PTO speed equals heave speed and obeys the 1.2 MN s/m damping and absorbed-power laws. Inactive speeds and accelerations remain below `1e-9` in their units; the remaining off-diagonal heave mass couplings therefore apply less than `1e-6` N. The two regenerated HDF5 files agree in all 38 datasets with the [pinned source HDF5](../ellipsoid_ode45_convergence/hydroData/ellipsoid.h5). + +These are peak differences over the entire 150 s between the matched mass control and the archived continuous-nonlinear source run. Signed PTO work rose by 1.217% and 1.043%, respectively, relative to that source run. + +| Control minus archived source, peak magnitude | Regular | RegularCIC | +| --- | ---: | ---: | +| Heave | 4.661 mm | 3.880 mm | +| Speed | 7.538 mm/s | 7.074 mm/s | +| Logged total heave force | 302.729 kN | 173.441 kN | +| Restoring force | 9.776 kN | 9.215 kN | +| Excitation force | 8.303 kN | 8.465 kN | +| Radiation force | 0.884 kN | 0.608 kN | +| Heave acceleration | 0.03321 m/s² | 0.03476 m/s² | +| PTO force | 9.046 kN | 8.489 kN | +| Absorbed power | 7.762 kW | 7.217 kW | + +The total-force change includes the shifted added-mass term; it is not a change in the pressure-force law at a fixed state. On every control state, the original-mass mesh reconstructs the source restoring and drag forces within `6e-9` and `7e-11` N, respectively. The source restoring block still uses the original `133,376.067` kg equilibrium mass. The change also shifts the inactive surge and sway inertial masses and changes the solver's feedback representation, so this control cannot allocate a numerical effect uniquely to the 100 ns delay. + +The [portable Python test](../../../test_ellipsoid_ode45_continuous_nonlinear_mass_control.py) advances independent motion and PTO from the pinned HDF5/STL and published input. It compares all 3,001 source samples with that trajectory and with the archived continuous-nonlinear run: + +| Case and source arm | Heave | Speed | PTO force | Absorbed power | Signed work versus Python | +| --- | ---: | ---: | ---: | ---: | ---: | +| Regular continuous nonlinear | 3.421 mm | 6.021 mm/s | 7.225 kN | 5.901 kW | +0.880% | +| Regular no active delay | 7.978 mm | 13.356 mm/s | 16.027 kN | 13.663 kW | +2.108% | +| RegularCIC continuous nonlinear | 4.307 mm | 8.196 mm/s | 9.835 kN | 8.199 kW | +1.114% | +| RegularCIC no active delay | 8.056 mm | 15.063 mm/s | 18.075 kN | 15.417 kW | +2.169% | + +Removing the active delayed heave feedback **worsens every inherited motion, PTO, and work gap**. Both controls still fail the established 3/3.5 mm heave, 4/4.5 mm/s speed, 5/6 kN PTO-force, 4/5 kW power, and 0.1% work limits. This rules out the active delayed feedback as a sufficient explanation for the remaining continuous-nonlinear residual. It does not establish the underlying cause or promote either published input to physical parity. + +The [portable PR workflow](../../../../.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-portable.yml) runs this diagnostic without MATLAB. A separate [post-screen source workflow](../../../../.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step.yml) requests only copied-model `MaxStep` values of 0.0125 and 0.00625 s for both waves, while leaving the published 0.05 s clocks and continuous force blocks in place. It exports to a distinct directory and artifact per case and step; those future runs are separate from the 0.05 s fixtures here. + +`regular.npz` and `regular_cic.npz` retain time, active and inactive motion, source force terms, actually applied heave added force, PTO channels, wave, nine checked settings, and the twelve-value mass identity. Their SHA-256 hashes are `69018348739432556f08ad9e91bb32d855175b7b89e304174b8d9082309a98c7` and `364ef737b81a9ba3d2fd3b72ca82b0f3a95f68f42aa4aa16629ec6d6e0cfb435`. Regenerate and verify them with: + +```sh +gh run download 40130836521 -n ellipsoid-ode45-continuous-nonlinear-mass-control-Regular -D .evidence/mass_control/Regular +gh run download 40130836521 -n ellipsoid-ode45-continuous-nonlinear-mass-control-RegularCIC -D .evidence/mass_control/RegularCIC +.venv/bin/python tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py Regular .evidence/mass_control/Regular/matlab-ellipsoid-ode45-continuous-nonlinear-mass-control/Regular .evidence/mass_control/regular-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular.npz --source-hydro .evidence/mass_control/Regular/applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 +.venv/bin/python tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py RegularCIC .evidence/mass_control/RegularCIC/matlab-ellipsoid-ode45-continuous-nonlinear-mass-control/RegularCIC .evidence/mass_control/regular-cic-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular_cic.npz --source-hydro .evidence/mass_control/RegularCIC/applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 +.venv/bin/python -m pytest -q tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py +``` diff --git a/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular.npz b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular.npz new file mode 100644 index 00000000..3a4df4de Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular_cic.npz b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular_cic.npz new file mode 100644 index 00000000..f5e3b7b6 Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control/regular_cic.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/README.md b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/README.md new file mode 100644 index 00000000..dae51204 --- /dev/null +++ b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/README.md @@ -0,0 +1,102 @@ +# Ellipsoid no-active-delay MaxStep continuation + +These fixtures archive the four successful 150 s source controls from +[GitHub Actions run 40131806849](https://github.com/cmudrc/wec-sim-python/actions/runs/40131806849) +at repository commit `f96a2056637f9bf0c548851cc45b7b12d1707865`. +The workflow used MATLAB R2025b, pinned WEC-Sim core +`0753b2e47f2457c078751dcfe5d251d1767b80ab`, and pinned Applications +`d53d4d4c9eda2581f04204f5d394a6ef84bb099e`. Its +[runner](../../../matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m) +copied the published input and model, made the restoring, Froude–Krylov, and +wave-elevation blocks continuous, and used the [archived 0.05 s mass +control](../ellipsoid_ode45_continuous_nonlinear_mass_control/README.md). +Only the copied model's `ode45` MaxStep changed to 0.0125 or 0.00625 s +here. Output, nonlinear, and CIC clocks remained 0.05 s; duration remained +150 s. The 100 ns Transport Delay block remained present. + +The mass control sets `adjMassFactor = A_zz / trace(A_xyz)` from the pinned +HDF5. It moves `A_zz` into the rigid inertial mass, leaving zero active +heave acceleration coefficient applied through the delay. Across all three +MaxSteps, the mass values were unchanged: + +| Source mass identity | Regular | RegularCIC | +| --- | ---: | ---: | +| `adjMassFactor` | 0.655990137167053 | 0.783566018955714 | +| Original rigid mass (kg) | 133,376.067 | 133,376.067 | +| Heave added mass `A_zz` (kg) | 199,883.891 | 149,463.860 | +| Adjusted inertial mass (kg) | 333,259.957 | 282,839.927 | +| Residual active heave coefficient (kg) | 0 | 0 | +| Largest applied heave added-mass force (N) | 5.99e-22 | 3.76e-22 | + +Each NPZ stores both refined steps on the same 3,001-point, 0.05 s output +grid: heave, speed, total and component source forces, acceleration, active +and inactive added-mass evidence, PTO stroke/speed/force/absorbed power, +incident wave, nine checked settings, and twelve mass values. The +[packer](../../../../tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py) +compared all 38 datasets of every exported HDF5 against the pinned source +HDF5 by value, including NaN locations. The saved wave agrees with both the +0.05 s control and paired delayed source within `1e-10` m. On every saved +state, the source mesh reconstructs restoring and viscous drag within +`7e-9` and `7e-11` N. PTO channels obey the 1.2 MN s/m damping and +absorbed-power laws. Inactive motion contributes less than `1e-6` N of +cross-axis added-mass force. + +The following peak absolute gaps use **all 3,001 samples**, and work is +signed trapezoidal absorbed energy over 150 s. The independent Python run +uses the pinned HDF5/STL and published wave/PTO inputs, with no source +trajectory supplied to its dynamics. The paired delayed arms are the +[continuous-nonlinear 0.05 s base](../ellipsoid_ode45_restoring_timing/README.md) +and [refined source fixtures](../ellipsoid_ode45_restoring_step/README.md) +at the same MaxStep. + +| Case | MaxStep (s) | No-delay versus Python: heave (mm) | Speed (mm/s) | PTO force (kN) | Power (kW) | Signed work | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | +| Regular | 0.05 | 7.978 | 13.356 | 16.027 | 13.663 | +2.1080% | +| Regular | 0.0125 | 1.916 | 3.356 | 4.027 | 3.156 | +0.5352% | +| Regular | 0.00625 | 0.956 | 1.716 | 2.060 | 1.640 | +0.2688% | +| RegularCIC | 0.05 | 8.056 | 15.063 | 18.075 | 15.417 | +2.1690% | +| RegularCIC | 0.0125 | 2.073 | 3.699 | 4.439 | 3.410 | +0.5131% | +| RegularCIC | 0.00625 | 1.122 | 1.843 | 2.211 | 1.672 | +0.2618% | + +| Case | MaxStep (s) | No-delay versus paired delayed source: heave (mm) | Speed (mm/s) | PTO force (kN) | Power (kW) | Signed work | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | +| Regular | 0.05 | 4.661 | 7.538 | 9.046 | 7.762 | +1.2171% | +| Regular | 0.0125 | 1.100 | 1.853 | 2.223 | 1.764 | +0.2991% | +| Regular | 0.00625 | 0.537 | 0.951 | 1.141 | 0.916 | +0.1459% | +| RegularCIC | 0.05 | 3.880 | 7.074 | 8.489 | 7.217 | +1.0432% | +| RegularCIC | 0.0125 | 0.924 | 1.699 | 2.039 | 1.548 | +0.2300% | +| RegularCIC | 0.00625 | 0.447 | 0.850 | 1.020 | 0.778 | +0.1204% | + +All motion, PTO, and work gaps contract over these three steps. The +0.0125-to-0.00625 no-delay change is still 0.977 mm heave, 1.784 mm/s +speed, 2.141 kN PTO force, 1.574 kW power, and 0.2650% signed work for +Regular; for RegularCIC it is 1.073 mm, 1.982 mm/s, 2.378 kN, 1.738 kW, +and 0.2500%. These finite-step observations do not establish an +asymptotic limit or isolate the 100 ns delay: `adjMassFactor` also shifts +inactive inertial masses and changes the solver's feedback representation. +The source restoring block still uses the original 133,376.067 kg +equilibrium mass. + +At 0.00625 s the no-delay arms meet the inherited pointwise Python +ceilings (Regular: 3 mm, 4 mm/s, 5 kN, 4 kW; RegularCIC: 3.5 mm, 4.5 mm/s, +6 kN, 5 kW), but **both fail the 0.1% signed-work ceiling**. They are +post-screen diagnostics. They do not change the failed prospective source +step screen or the unchanged published-input parity gate. + +The [portable PR +workflow](../../../../.github/workflows/ellipsoid-ode45-continuous-nonlinear-mass-control-step-portable.yml) +runs the [diagnostic test](../../../test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py) +without MATLAB. The two archives are 484,595 and 593,172 bytes, with +SHA-256 hashes `5d9eaefe9f305c0e5e7b41d811a232dd246ade78384e304601dd979ce079d4da` +and `692ffcb3ea207743ff7b77b7b0072a8ce3319bc4c382edf963c306c57cd19282`. +Repack and compare the downloaded artifacts: + +```sh +gh run download 40131806849 -n ellipsoid-ode45-continuous-nonlinear-mass-control-Regular-maxstep-0.0125 -D .evidence/mass_control_step/Regular_0p0125 +gh run download 40131806849 -n ellipsoid-ode45-continuous-nonlinear-mass-control-Regular-maxstep-0.00625 -D .evidence/mass_control_step/Regular_0p00625 +gh run download 40131806849 -n ellipsoid-ode45-continuous-nonlinear-mass-control-RegularCIC-maxstep-0.0125 -D .evidence/mass_control_step/RegularCIC_0p0125 +gh run download 40131806849 -n ellipsoid-ode45-continuous-nonlinear-mass-control-RegularCIC-maxstep-0.00625 -D .evidence/mass_control_step/RegularCIC_0p00625 +.venv/bin/python tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py Regular .evidence/mass_control_step .evidence/regular-mass-step-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular.npz +.venv/bin/python tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py RegularCIC .evidence/mass_control_step .evidence/regular-cic-mass-step-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular_cic.npz +.venv/bin/python -m pytest -q -s tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py +``` diff --git a/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular.npz b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular.npz new file mode 100644 index 00000000..6656a84b Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular_cic.npz b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular_cic.npz new file mode 100644 index 00000000..694bf1b0 Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control_step/regular_cic.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/README.md b/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/README.md new file mode 100644 index 00000000..970e2595 --- /dev/null +++ b/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/README.md @@ -0,0 +1,74 @@ +# Ellipsoid copied-model continuous nonlinear MaxStep continuation + +These fixtures preserve the 150 s Regular and RegularCIC source traces from +[GitHub Actions run 40131187182](https://github.com/cmudrc/wec-sim-python/actions/runs/40131187182), +which succeeded at Python-repository commit +`b910e47d5a6ad2edd5b6af183374d83d74db1e8a`. The workflow pinned WEC-Sim +core `0753b2e47f2457c078751dcfe5d251d1767b80ab` and Applications +`d53d4d4c9eda2581f04204f5d394a6ef84bb099e`, using MATLAB R2025b. The +[runner](../../../matlab/ellipsoid_ode45_restoring_timing_smoke.m) copied each +published application input and model, made the restoring, Froude–Krylov, and +wave-elevation blocks continuous, then varied only the copied model's `ode45` +MaxStep among these arms. The published input stays at 0.05 s. The physical, output, +nonlinear, and CIC clocks stay at 0.05 s, with the original 150 s duration, +50 s ramp, 4 m/6 s sea, and 1.2 MN s/m PTO damper. + +Each NPZ contains 3,001 output samples for four copied-model MaxSteps: +0.025, 0.0125, 0.00625, and 0.003125 s. It stores heave, speed, restoring +force, PTO stroke/speed/force/absorbed power, incident wave, nine source +settings, and the common time grid. The separate +[full force-timing fixture](../ellipsoid_ode45_restoring_timing/README.md) +supplies the 0.05 s continuous nonlinear base. The packer checked all eight +fresh source HDF5 files against the 38 datasets in the archived source HDF5 +by value, including NaN locations. Binary HDF5 hashes need not match. + +The original predeclared screen applied to the 0.025→0.0125 s change. Its +pointwise limits were 1.0 mm heave, 1.3 mm/s speed, 1.6 kN PTO force, and +1.3 kW power for Regular; 1.1 mm, 1.5 mm/s, 2.0 kN, and 1.6 kW for +RegularCIC. Both had a 0.1% signed-work change limit and required contraction +by a factor below 0.7 from the preceding change. Every measure contracted, +but the **predeclared screen failed**: Regular exceeded speed, force, power, +and work limits; RegularCIC exceeded all five limits. + +The following 0.00625 and 0.003125 s runs are a **post-screen diagnostic**. +Their last changes fall below the old numerical ceilings, but this does not +retroactively pass the predeclared 0.0125 s screen. + +| Case | MaxStep pair (s) | Heave (mm) | Speed (mm/s) | PTO force (kN) | Power (kW) | Signed-work change | +| --- | --- | ---: | ---: | ---: | ---: | ---: | +| Regular | 0.05→0.025 | 1.700 | 3.232 | 3.879 | 3.287 | 0.4186% | +| Regular | 0.025→0.0125 | 0.907 | 1.722 | 2.066 | 1.339 | 0.2238% | +| Regular | 0.0125→0.00625 | 0.440 | 0.797 | 0.956 | 0.693 | 0.1125% | +| Regular | 0.00625→0.003125 | 0.242 | 0.395 | 0.474 | 0.367 | 0.0634% | +| RegularCIC | 0.05→0.025 | 2.218 | 4.511 | 5.413 | 4.598 | 0.5494% | +| RegularCIC | 0.025→0.0125 | 1.115 | 2.356 | 2.828 | 1.940 | 0.2785% | +| RegularCIC | 0.0125→0.00625 | 0.556 | 1.082 | 1.299 | 0.968 | 0.1410% | +| RegularCIC | 0.00625→0.003125 | 0.293 | 0.549 | 0.659 | 0.504 | 0.0712% | + +An independent Python run uses the archived source HDF5/STL and configured +wave/PTO input, without source motion or force histories. At the diagnostic +0.003125 s arm, all 3,001 pointwise samples and signed absorbed work fall +within the inherited Python physical limits: + +| Case | Heave gap (mm) | Speed gap (mm/s) | PTO force gap (kN) | Power gap (kW) | Signed-work gap | +| --- | ---: | ---: | ---: | ---: | ---: | +| Regular | 0.220 | 0.383 | 0.460 | 0.371 | 0.0593% | +| RegularCIC | 0.406 | 0.516 | 0.619 | 0.412 | 0.0700% | + +The inherited physical limits are 3 mm/4 mm/s/5 kN/4 kW for Regular and +3.5 mm/4.5 mm/s/6 kN/5 kW for RegularCIC, plus 0.1% signed work for both. +This result concerns only the copied continuous nonlinear model. The +unchanged published-input gate remains red and its status is unchanged. + +`regular.npz` is 533,989 bytes (SHA-256 +`b19a80a8d2451689056d2dea4ba85ecbaa0961e9cc850a641c4c6271df589246`); +`regular_cic.npz` is 633,842 bytes (SHA-256 +`647c1be3f837912442e380417a98af00245bb7ba1b0da3f9d6c7348d1b3e1546`). +Regenerate and compare the archives from the workflow artifacts: + +```sh +gh run download 40131187182 -D .evidence/force_timing_finer +.venv/bin/python tools/pack_ellipsoid_ode45_restoring_step.py Regular .evidence/force_timing_finer .evidence/regular-step-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular.npz +.venv/bin/python tools/pack_ellipsoid_ode45_restoring_step.py RegularCIC .evidence/force_timing_finer .evidence/regular-cic-step-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular_cic.npz +.venv/bin/python -m pytest -q -s tests/test_ellipsoid_ode45_restoring_step_portable.py +``` diff --git a/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular.npz b/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular.npz new file mode 100644 index 00000000..93e71ead Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular_cic.npz b/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular_cic.npz new file mode 100644 index 00000000..9c6ede4b Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_restoring_step/regular_cic.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/README.md b/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/README.md new file mode 100644 index 00000000..06d9605e --- /dev/null +++ b/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/README.md @@ -0,0 +1,84 @@ +# Ellipsoid ode45 nonlinear-force timing, full 150 s + +This fixture is a diagnostic of [issue #599](https://github.com/cmudrc/wec-sim-python/issues/599). +The [R2025b source run 40129984760](https://github.com/cmudrc/wec-sim-python/actions/runs/40129984760) +completed both Regular and RegularCIC jobs at Python-repository commit +`9ed1e2ef2775f5bba02ace7433c3308772f74ddd` on +`codex/ellipsoid-force-timing-599`. Its pinned source checkouts were WEC-Sim +core `0753b2e47f2457c078751dcfe5d251d1767b80ab` and Applications +`d53d4d4c9eda2581f04204f5d394a6ef84bb099e`. The source runner is +[`ellipsoid_ode45_restoring_timing_smoke.m`](../../../matlab/ellipsoid_ode45_restoring_timing_smoke.m), +called with `duration='full'`; it copies the application input and model for +each arm. Nothing is written to either upstream repository. + +Both cases use the original 150 s, 4 m/6 s regular sea, 50 s ramp, ellipsoid +mesh and hydro data, 1.2 MN s/m PTO damper, `ode45`, 0.05 s solver maximum, +and 0.05 s physical/output/nonlinear/CIC clocks. The three arms differ only +in compiled sample time within the copied model: + +| Arm | Nonlinear restoring | Nonlinear Froude–Krylov | Nonlinear wave elevation | +| --- | ---: | ---: | ---: | +| `baseline` | 0.05 s | 0.05 s | 0.05 s | +| `continuous_restoring` | continuous | 0.05 s | 0.05 s | +| `continuous_nonlinear` | continuous | continuous | continuous | + +Each archive retains all 3,001 shared output samples for time, heave position +and speed, restoring and excitation force, acceleration, PTO stroke/speed/force +and absorbed power, wave elevation, plus all nine checked source settings for +each arm. The copied baseline matches the previously archived unchanged +published trace. The shared incident waves match across all arms within +`1e-10` m. The active PTO vertical speed matches body heave speed within +`1e-9` m/s; its force and power satisfy the source damping law. The 38 HDF5 +datasets freshly generated in each source job agree in value with the +[archived source HDF5](../ellipsoid_ode45_convergence/hydroData/ellipsoid.h5), +although binary file hashes differ. The same archived source +[`elipsoid.stl`](../ellipsoid_ode45_convergence/geometry/elipsoid.stl) is used +by the portable Python test. + +The source restoring log has the opposite sign to the mesh buoyancy-minus-weight +force reconstructed in Python. Over all 150 s, the baseline log follows the +**preceding** output state within `6e-9` N; each edited arm follows the +**current** state within `6e-9` N. Applying the opposite timing law leaves +over 51 kN of mismatch. This confirms that the copied-model intervention +changed the intended active force path. It does not change the original +published source input. + +The [portable gate](../../../test_ellipsoid_ode45_restoring_timing_portable.py) +independently advances Python at the published 0.05 s step using only the +archived source-generated HDF5/STL and the configured wave/PTO inputs. It +checks all 3,001 samples and integrated signed absorbed work. These are the +measured peak gaps to that independent physical Python run: + +| Case and source arm | Heave | Speed | PTO force | Absorbed power | Work | +| --- | ---: | ---: | ---: | ---: | ---: | +| Regular baseline | 18.940 mm | 28.341 mm/s | 34.009 kN | 26.270 kW | 5.819% | +| Regular continuous restoring | 5.156 mm | 5.860 mm/s | 7.032 kN | 7.216 kW | 1.433% | +| Regular continuous nonlinear | 3.421 mm | 6.021 mm/s | 7.225 kN | 5.901 kW | 0.880% | +| RegularCIC baseline | 19.930 mm | 30.042 mm/s | 36.050 kN | 28.771 kW | 6.194% | +| RegularCIC continuous restoring | 6.390 mm | 7.516 mm/s | 9.019 kN | 9.098 kW | 1.772% | +| RegularCIC continuous nonlinear | 4.307 mm | 8.196 mm/s | 9.835 kN | 8.199 kW | 1.114% | + +Both edits improve every listed gap by more than half against the copied +baseline, but **both fail every inherited pointwise limit**. Regular's limits +are 3 mm heave/PTO stroke, 4 mm/s speed, 5 kN force, and 4 kW power; +RegularCIC's are 3.5 mm, 4.5 mm/s, 6 kN, and 5 kW. Each also misses the +unchanged 0.1% signed-work limit. The portable test treats the measured +improvement as a regression screen, not as a prospectively declared passing +threshold. The subsequent [predeclared source-step screen](../ellipsoid_ode45_restoring_step/README.md) +failed at MaxStep 0.0125 s for both waves. Later finer steps are diagnostic; +they do not replace the separate derived fine-step reference. The unchanged +published-input gate stays red and both case rows remain `qualified_motion`. + +`regular.npz` is 504,271 bytes (SHA-256 +`0477f8433ff9dc4cb75913b36a726ab2a04478aa01137004f8fc3fbd85995a33`); +`regular_cic.npz` is 612,579 bytes (SHA-256 +`f4d76114d8ba54c37a3256a19b2bec13a8a6ed754a37824c98cbb4626a96ccdd`). +Regenerate or compare each fixture from its downloaded workflow artifact: + +```sh +gh run download 40129984760 -n ellipsoid-ode45-restoring-timing-Regular-full -D .evidence/force_timing_full/ellipsoid-ode45-restoring-timing-Regular-full +.venv/bin/python tools/pack_ellipsoid_ode45_restoring_timing.py Regular .evidence/force_timing_full/ellipsoid-ode45-restoring-timing-Regular-full/matlab-ellipsoid-ode45-restoring-timing/Regular .evidence/regular-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular.npz --source-hydro .evidence/force_timing_full/ellipsoid-ode45-restoring-timing-Regular-full/applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 +gh run download 40129984760 -n ellipsoid-ode45-restoring-timing-RegularCIC-full -D .evidence/force_timing_full/ellipsoid-ode45-restoring-timing-RegularCIC-full +.venv/bin/python tools/pack_ellipsoid_ode45_restoring_timing.py RegularCIC .evidence/force_timing_full/ellipsoid-ode45-restoring-timing-RegularCIC-full/matlab-ellipsoid-ode45-restoring-timing/RegularCIC .evidence/regular-cic-repacked.npz --expect tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular_cic.npz --source-hydro .evidence/force_timing_full/ellipsoid-ode45-restoring-timing-RegularCIC-full/applications/Nonlinear_Hydro/hydroData/ellipsoid.h5 +.venv/bin/python -m pytest -q tests/test_ellipsoid_ode45_restoring_timing_portable.py +``` diff --git a/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular.npz b/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular.npz new file mode 100644 index 00000000..e698d12c Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular.npz differ diff --git a/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular_cic.npz b/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular_cic.npz new file mode 100644 index 00000000..983eed44 Binary files /dev/null and b/tests/fixtures/matlab/ellipsoid_ode45_restoring_timing/regular_cic.npz differ diff --git a/tests/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m b/tests/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m new file mode 100644 index 00000000..38e4bbc7 --- /dev/null +++ b/tests/matlab/ellipsoid_ode45_continuous_nonlinear_mass_control.m @@ -0,0 +1,232 @@ +% Heave-only no-active-delay control for the 150 s continuous-nonlinear arm. +% Invoke as a script with waveCase = 'Regular' or 'RegularCIC' and optional +% maxStep (0.05 by default; 0.0125 and 0.00625 for post-screen diagnosis). +assert(any(strcmp(waveCase, {'Regular', 'RegularCIC'}))); +if ~exist('maxStep', 'var') + maxStep = .05; +end +assert(isnumeric(maxStep) && isscalar(maxStep) && isfinite(maxStep) && ... + any(abs(maxStep - [.05, .0125, .00625]) < 1e-12), ... + 'MaxStep must be a supported copied-model value'); +refineStep = abs(maxStep - .05) >= 1e-12; +stepTag = strrep(num2str(maxStep, '%.5g'), '.', 'p'); +repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath')))); +sourceDir = fullfile(repoRoot, 'matlab-ref', 'source'); +appDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro'); +hydroDir = fullfile(appDir, 'hydroData'); +originalDir = fullfile(appDir, 'ode45', waveCase); +addpath(genpath(sourceDir)); +cd(hydroDir); +if ~isfile('ellipsoid.h5') + bemio; +end + +% Match PR #658's heave mass identity. This shifts exactly A_zz into the +% body's rigid mass, leaving m + A_zz unchanged and the active delayed +% added-mass coefficient zero. The 100 ns Simulink delay remains present. +hydro = readBEMIOH5(fullfile(hydroDir, 'ellipsoid.h5'), 1, 0); +if strcmp(waveCase, 'Regular') + omega = 2*pi/6; + frequencies = hydro.simulation_parameters.w; + diagonal = zeros(3, 1); + for axis = 1:3 + diagonal(axis) = interp1(frequencies, ... + squeeze(hydro.hydro_coeffs.added_mass.all(axis, axis, :)), ... + omega, 'spline'); + end +else + matrix = hydro.hydro_coeffs.added_mass.inf_freq; + diagonal = diag(matrix(1:3, 1:3)); +end +factor = diagonal(3)/sum(diagonal); +assert(isfinite(factor) && factor > 0 && factor < 1, ... + 'The pinned heave-only added-mass coefficient is invalid'); + +original = fileread(fullfile(originalDir, 'wecSimInputFile.m')); +assert(contains(original, 'simu.solver=''ode45'';') && ... + contains(original, 'simu.dt = 0.05;') && ... + contains(original, 'simu.endTime=150;') && ... + contains(original, 'body(1).nonlinearHydro = 2;') && ... + ~contains(original, 'simu.dtOut') && ... + ~contains(original, 'simu.nonlinearDt') && ... + ~contains(original, 'simu.cicDt') && ... + ~contains(original, 'body(1).adjMassFactor'), ... + 'The pinned published input has changed'); +caseName = [waveCase '_continuous_nonlinear_mass_control']; +if refineStep + caseName = [caseName '_maxstep_' stepTag]; +end +caseDir = fullfile(fileparts(originalDir), caseName); +assert(~isfolder(caseDir), 'Disposable mass-control case already exists'); +mkdir(caseDir); +copyfile(fullfile(originalDir, 'ellipsoid.slx'), ... + fullfile(caseDir, 'ellipsoid.slx')); +corrected = strrep(original, 'body(1).nonlinearHydro = 2;', sprintf( ... + ['body(1).nonlinearHydro = 2;\n' ... + 'body(1).adjMassFactor = %.17g; %% heave-only exact mass shift'], factor)); +assert(~strcmp(corrected, original) && ... + contains(corrected, 'body(1).adjMassFactor = ') && ... + contains(corrected, 'simu.dt = 0.05;'), ... + 'The copied input did not receive only the mass factor'); +fid = fopen(fullfile(caseDir, 'wecSimInputFile.m'), 'w'); +assert(fid > 0, 'Could not write disposable mass-control input'); +fprintf(fid, '%s', corrected); +fclose(fid); +assert(strcmp(fileread(fullfile(caseDir, 'wecSimInputFile.m')), corrected)); + +cd(caseDir); +run('initializeWecSim'); +assert(strcmp(simu.solver, 'ode45') && ... + abs(simu.dt - .05) < 1e-12 && ... + abs(simu.dtOut - .05) < 1e-12 && ... + abs(simu.nonlinearDt - .05) < 1e-12 && ... + abs(simu.cicDt - .05) < 1e-12 && ... + simu.endTime == 150 && simu.rampTime == 50 && ... + abs(body(1).adjMassFactor - factor) < 1e-12, ... + 'A published clock, duration, or derived mass factor changed'); +model = 'ellipsoid'; +assert(abs(str2double(get_param(model, 'MaxStep')) - .05) < 1e-12 && ... + abs(str2double(get_param(model, 'StopTime')) - 150) < 1e-12); +delayBlock = [model '/Rigid Body/Hydrodynamic Body/Transport Delay']; +assert(abs(str2double(get_param(delayBlock, 'DelayTime')) - 1e-7) < 1e-16, ... + 'The copied model has an unexpected core delay'); + +force = body(1).hydroForce.hf1; +originalMass = force.storage.fAddedMass; +assert(abs(originalMass(3,3) - simu.rho*diagonal(3)) < 1e-7, ... + 'The HDF5 and initialized heave added mass differ'); +assert(abs(sum(diag(originalMass(1:3,1:3)))*factor - ... + originalMass(3,3)) < 1e-7); +assert(abs(force.fAddedMass(3,3)) < 1e-7, ... + 'Active heave acceleration still passes through the delay'); +assert(abs(force.adjustedMass - (force.mass + originalMass(3,3))) < 1e-7, ... + 'The corrected heave mass is not m + A_zz'); +inactive = [1 2 4 5 6]; +massRecord = [factor, sum(diag(originalMass(1:3,1:3))), ... + originalMass(3,3), force.mass, force.adjustedMass, ... + force.fAddedMass(3,3), maxStep, force.fAddedMass(3,inactive)]; + +% Apply exactly the three continuous nonlinear evaluations in the +% published-clock timing arm, in a disposable copy of the application model. +bodyLink = [model '/Rigid Body']; +if ~strcmp(get_param(bodyLink, 'LinkStatus'), 'none') + set_param(bodyLink, 'LinkStatus', 'none'); +end +hydroLink = [bodyLink '/Hydrodynamic Body']; +if ~strcmp(get_param(hydroLink, 'LinkStatus'), 'none') + set_param(hydroLink, 'LinkStatus', 'none'); +end +restoringBlock = [hydroLink ... + '/Hydrostatic Restoring Force Calculation' ... + '/Linear and Nonlinear Restoring Force Variant Subsystem' ... + '/Nonlinear Hydrostatic Restoring Force/Restoring Forces']; +fkBlock = [hydroLink ... + '/Wave Diffraction and Excitation Force Calculation' ... + '/Nonlinear FK Force Variant Subsystem' ... + '/With Nonlinear FroudeKrylov Force']; +waveBlock = [hydroLink '/Nonlinear Wave Elevation']; +blocks = {restoringBlock, fkBlock, waveBlock}; +for iBlock = 1:numel(blocks) + assert(strcmp(get_param(blocks{iBlock}, 'TreatAsAtomicUnit'), 'on') && ... + strcmp(get_param(blocks{iBlock}, 'SystemSampleTime'), ... + 'simu.nonlinearDt'), ... + 'A pinned nonlinear sample-time block changed'); + set_param(blocks{iBlock}, 'SystemSampleTime', '0'); + assert(strcmp(get_param(blocks{iBlock}, 'SystemSampleTime'), '0')); +end +if refineStep + set_param(model, 'MaxStep', num2str(maxStep, '%.17g')); +end +set_param(model, 'SimulationCommand', 'update'); +compiledPeriod = zeros(1, numel(blocks)); +for iBlock = 1:numel(blocks) + sample = get_param(blocks{iBlock}, 'CompiledSampleTime'); + fprintf('COMPILED_SAMPLE_TIMES case=%s block=%d\n', waveCase, iBlock); + disp(sample); + if ~iscell(sample) + sample = {sample}; + end + rates = zeros(1, numel(sample)); + for iRate = 1:numel(sample) + pair = sample{iRate}; + if ischar(pair) || isstring(pair) + pair = str2num(char(pair)); %#ok compiled Simulink value + end + assert(isnumeric(pair) && numel(pair) == 2 && ... + abs(pair(2)) < 1e-12, ... + 'Nonlinear block has an unparseable compiled sample time'); + rates(iRate) = pair(1); + end + assert(any(abs(rates) < 1e-12), ... + 'Continuous nonlinear sample period is absent from compiled rates'); + compiledPeriod(iBlock) = rates(find(abs(rates) < 1e-12, 1)); +end +assert(abs(str2double(get_param(model, 'MaxStep')) - maxStep) < 1e-12 && ... + abs(str2double(get_param(delayBlock, 'DelayTime')) - 1e-7) < 1e-16 && ... + abs(str2double(get_param(model, 'StopTime')) - 150) < 1e-12, ... + 'The published solver, feedback delay, or duration changed'); +fprintf(['CONTINUOUS_NONLINEAR_MASS_CONTROL case=%s factor=%.17g ' ... + 'activeAddedMass=%.9g residual=%.9g effectiveMass=%.9g ' ... + 'dt=%.9g dtOut=%.9g nonlinearDt=%.9g cicDt=%.9g MaxStep=%.9g\n'], ... + waveCase, factor, originalMass(3,3), force.fAddedMass(3,3), ... + force.adjustedMass, simu.dt, simu.dtOut, simu.nonlinearDt, ... + simu.cicDt, str2double(get_param(model, 'MaxStep'))); + +sim(simu.simMechanicsFile, [], simset('SrcWorkspace', 'parent')); +assert(abs(str2double(get_param(model, 'MaxStep')) - maxStep) < 1e-12 && ... + abs(str2double(get_param(model, 'StopTime')) - 150) < 1e-12, ... + 'Simulink changed the published model settings'); +run('stopWecSim'); +record = output.bodies(1); +% responseClass postprocessing recomputes reported added mass. The storage +% field below retains the force that actually passed the Simulink delay. +appliedValues = [record.time(:), ... + body(1).hydroForce.hf1.storage.output_forceAddedMass]; +bodyValues = [record.time(:), record.position, record.velocity, ... + record.forceTotal, record.forceExcitation, ... + record.forceRadiationDamping, record.forceAddedMass, ... + record.forceRestoring, record.forceMorisonAndViscous, ... + record.acceleration]; +pto = output.ptos(1); +ptoValues = [pto.time(:), pto.position, pto.velocity, ... + pto.forceInternalMechanics, pto.powerInternalMechanics]; +waveValues = [output.wave.time(:), output.wave.elevation(:)]; +inactiveAcceleration = record.acceleration(:, inactive); +crossHeaveForce = inactiveAcceleration * force.fAddedMass(3, inactive)'; +assert(isequal(size(bodyValues), [3001, 55]) && ... + isequal(size(appliedValues), [3001, 7]) && ... + isequal(size(ptoValues), [3001, 25]) && ... + isequal(size(waveValues), [3001, 2]) && ... + max(abs(bodyValues(:, 1) - (0:3000)' * .05)) < 1e-9 && ... + max(abs(appliedValues(:, 1) - bodyValues(:, 1))) < 1e-9 && ... + max(abs(ptoValues(:, 1) - bodyValues(:, 1))) < 1e-9 && ... + max(abs(waveValues(:, 1) - bodyValues(:, 1))) < 1e-9 && ... + all(isfinite(bodyValues), 'all') && ... + all(isfinite(appliedValues), 'all') && ... + all(isfinite(ptoValues), 'all') && ... + all(isfinite(waveValues), 'all') && ... + max(abs(record.velocity(:, inactive)), [], 'all') < 1e-9 && ... + max(abs(inactiveAcceleration), [], 'all') < 1e-9 && ... + max(abs(crossHeaveForce)) < 1e-6 && ... + max(abs(appliedValues(:, 4))) < 1e-6, ... + 'Incomplete control or nonzero inactive/added heave feedback'); +if refineStep + outDir = fullfile(repoRoot, ... + 'matlab-ellipsoid-ode45-continuous-nonlinear-mass-control-step', ... + waveCase, ['maxstep_' stepTag]); +else + outDir = fullfile(repoRoot, ... + 'matlab-ellipsoid-ode45-continuous-nonlinear-mass-control', waveCase); +end +if ~isfolder(outDir) + mkdir(outDir); +end +writematrix(bodyValues, fullfile(outDir, 'body.csv')); +writematrix(appliedValues, fullfile(outDir, 'applied_added.csv')); +writematrix(ptoValues, fullfile(outDir, 'pto.csv')); +writematrix(waveValues, fullfile(outDir, 'wave.csv')); +writematrix([simu.dt, simu.dtOut, simu.nonlinearDt, ... + simu.cicDt, str2double(get_param(model, 'MaxStep')), 150, ... + compiledPeriod], fullfile(outDir, 'settings.csv')); +writematrix(massRecord, fullfile(outDir, 'mass.csv')); +close_system(model, 0); diff --git a/tests/matlab/ellipsoid_ode45_restoring_timing_smoke.m b/tests/matlab/ellipsoid_ode45_restoring_timing_smoke.m new file mode 100644 index 00000000..47436733 --- /dev/null +++ b/tests/matlab/ellipsoid_ode45_restoring_timing_smoke.m @@ -0,0 +1,212 @@ +% Isolate the nonlinear restoring calculation's sample time in the pinned +% ellipsoid ode45 model. Invoke with waveCase, mode, optional duration, and +% optional maxStep (0.05 by default; 0.025 and 0.0125 for the predeclared +% screen; 0.00625 and 0.003125 for post-screen diagnosis). +% mode: baseline, continuous_restoring, or continuous_nonlinear. +% duration: smoke (2 s, default) or full (150 s). +assert(any(strcmp(waveCase, {'Regular', 'RegularCIC'}))); +assert(any(strcmp(mode, {'baseline', 'continuous_restoring', ... + 'continuous_nonlinear'}))); +if ~exist('duration', 'var') + duration = 'smoke'; +end +assert(any(strcmp(duration, {'smoke', 'full'}))); +if ~exist('maxStep', 'var') + maxStep = .05; +end +assert(isnumeric(maxStep) && isscalar(maxStep) && isfinite(maxStep) && ... + any(abs(maxStep - [.05, .025, .0125, .00625, .003125]) < 1e-12), ... + 'MaxStep must be a supported copied-model value'); +if abs(maxStep - .05) >= 1e-12 + assert(strcmp(mode, 'continuous_nonlinear') && ... + strcmp(duration, 'full'), ... + 'The source-step screen is only the 150 s continuous nonlinear arm'); +end +stepTag = strrep(num2str(maxStep, '%.4g'), '.', 'p'); +stopTime = 2; +if strcmp(duration, 'full') + stopTime = 150; +end + +repoRoot = fileparts(fileparts(fileparts(mfilename('fullpath')))); +sourceDir = fullfile(repoRoot, 'matlab-ref', 'source'); +appDir = fullfile(repoRoot, 'applications', 'Nonlinear_Hydro'); +hydroDir = fullfile(appDir, 'hydroData'); +originalDir = fullfile(appDir, 'ode45', waveCase); +addpath(genpath(sourceDir)); +cd(hydroDir); +if ~isfile('ellipsoid.h5') + bemio; +end + +original = fileread(fullfile(originalDir, 'wecSimInputFile.m')); +assert(contains(original, 'simu.solver=''ode45'';') && ... + contains(original, 'simu.dt = 0.05;') && ... + contains(original, 'simu.endTime=150;') && ... + contains(original, 'body(1).nonlinearHydro = 2;') && ... + ~contains(original, 'simu.dtOut') && ... + ~contains(original, 'simu.nonlinearDt') && ... + ~contains(original, 'simu.cicDt') && ... + ~contains(original, 'body(1).adjMassFactor'), ... + 'The pinned published input has changed'); +caseName = [waveCase '_restoring_timing_' mode '_' duration]; +if abs(maxStep - .05) >= 1e-12 + caseName = [caseName '_maxstep_' stepTag]; +end +caseDir = fullfile(fileparts(originalDir), caseName); +assert(~isfolder(caseDir), 'Disposable timing case already exists'); +mkdir(caseDir); +copyfile(fullfile(originalDir, 'ellipsoid.slx'), ... + fullfile(caseDir, 'ellipsoid.slx')); +copyfile(fullfile(originalDir, 'wecSimInputFile.m'), ... + fullfile(caseDir, 'wecSimInputFile.m')); +assert(strcmp(fileread(fullfile(caseDir, 'wecSimInputFile.m')), original)); + +cd(caseDir); +run('initializeWecSim'); +assert(strcmp(simu.solver, 'ode45') && ... + abs(simu.dt - .05) < 1e-12 && ... + abs(simu.dtOut - .05) < 1e-12 && ... + abs(simu.nonlinearDt - .05) < 1e-12 && ... + abs(simu.cicDt - .05) < 1e-12 && ... + simu.endTime == 150 && simu.rampTime == 50, ... + 'A published physical or output clock changed'); +model = 'ellipsoid'; +assert(abs(str2double(get_param(model, 'MaxStep')) - .05) < 1e-12); +delayBlock = [model '/Rigid Body/Hydrodynamic Body/Transport Delay']; +assert(abs(str2double(get_param(delayBlock, 'DelayTime')) - 1e-7) < 1e-16); + +% All arms detach the same copied library link. The source-step screen +% changes only copied-model MaxStep beyond the continuous nonlinear arm. +bodyLink = [model '/Rigid Body']; +if ~strcmp(get_param(bodyLink, 'LinkStatus'), 'none') + set_param(bodyLink, 'LinkStatus', 'none'); +end +hydroLink = [bodyLink '/Hydrodynamic Body']; +if ~strcmp(get_param(hydroLink, 'LinkStatus'), 'none') + set_param(hydroLink, 'LinkStatus', 'none'); +end +restoringBlock = [hydroLink ... + '/Hydrostatic Restoring Force Calculation' ... + '/Linear and Nonlinear Restoring Force Variant Subsystem' ... + '/Nonlinear Hydrostatic Restoring Force/Restoring Forces']; +fkBlock = [hydroLink ... + '/Wave Diffraction and Excitation Force Calculation' ... + '/Nonlinear FK Force Variant Subsystem' ... + '/With Nonlinear FroudeKrylov Force']; +waveBlock = [hydroLink '/Nonlinear Wave Elevation']; +blocks = {restoringBlock, fkBlock, waveBlock}; +for iBlock = 1:numel(blocks) + assert(strcmp(get_param(blocks{iBlock}, 'TreatAsAtomicUnit'), 'on') && ... + strcmp(get_param(blocks{iBlock}, 'SystemSampleTime'), ... + 'simu.nonlinearDt'), ... + 'A pinned nonlinear sample-time block changed'); +end +if ~strcmp(mode, 'baseline') + set_param(restoringBlock, 'SystemSampleTime', '0'); +end +if strcmp(mode, 'continuous_nonlinear') + set_param(fkBlock, 'SystemSampleTime', '0'); + set_param(waveBlock, 'SystemSampleTime', '0'); +end +expectedPeriod = [.05, .05, .05]; +if ~strcmp(mode, 'baseline') + expectedPeriod(1) = 0; +end +if strcmp(mode, 'continuous_nonlinear') + expectedPeriod(2:3) = 0; +end +for iBlock = 1:numel(blocks) + param = get_param(blocks{iBlock}, 'SystemSampleTime'); + if expectedPeriod(iBlock) == 0 + assert(strcmp(param, '0'), 'Continuous block change did not persist'); + else + assert(strcmp(param, 'simu.nonlinearDt'), ... + 'A non-target nonlinear block clock changed'); + end +end + +if abs(maxStep - .05) >= 1e-12 + set_param(model, 'MaxStep', num2str(maxStep, '%.17g')); +end +set_param(model, 'StopTime', num2str(stopTime)); +set_param(model, 'SimulationCommand', 'update'); +compiledPeriod = zeros(1, numel(blocks)); +for iBlock = 1:numel(blocks) + sample = get_param(blocks{iBlock}, 'CompiledSampleTime'); + fprintf('COMPILED_SAMPLE_TIMES case=%s mode=%s block=%d\n', ... + waveCase, mode, iBlock); + disp(sample); + if ~iscell(sample) + sample = {sample}; + end + rates = zeros(1, numel(sample)); + for iRate = 1:numel(sample) + pair = sample{iRate}; + if ischar(pair) || isstring(pair) + pair = str2num(char(pair)); %#ok compiled Simulink value + end + assert(isnumeric(pair) && numel(pair) == 2 && ... + abs(pair(2)) < 1e-12, ... + 'Nonlinear block has an unparseable compiled sample time'); + rates(iRate) = pair(1); + end + assert(any(abs(rates - expectedPeriod(iBlock)) < 1e-12), ... + 'Expected nonlinear sample period is absent from compiled rates'); + compiledPeriod(iBlock) = expectedPeriod(iBlock); +end +assert(abs(str2double(get_param(model, 'MaxStep')) - maxStep) < 1e-12 && ... + abs(str2double(get_param(delayBlock, 'DelayTime')) - 1e-7) < 1e-16 && ... + abs(str2double(get_param(model, 'StopTime')) - stopTime) < 1e-12, ... + 'The model solver, feedback delay, or duration changed'); +fprintf(['RESTORING_TIMING case=%s mode=%s duration=%s dt=%.9g ' ... + 'dtOut=%.9g nonlinearDt=%.9g cicDt=%.9g MaxStep=%.9g ' ... + 'compiledRestoring=%.9g compiledFK=%.9g compiledWave=%.9g\n'], ... + waveCase, mode, duration, simu.dt, simu.dtOut, simu.nonlinearDt, ... + simu.cicDt, str2double(get_param(model, 'MaxStep')), compiledPeriod); + +sim(simu.simMechanicsFile, [], simset('SrcWorkspace', 'parent')); +assert(abs(str2double(get_param(model, 'MaxStep')) - maxStep) < 1e-12 && ... + abs(str2double(get_param(model, 'StopTime')) - stopTime) < 1e-12, ... + 'Simulink changed the prescribed source-step screen settings'); +if abs(maxStep - .05) < 1e-12 + outDir = fullfile(repoRoot, ... + 'matlab-ellipsoid-ode45-restoring-timing', waveCase); +else + outDir = fullfile(repoRoot, ... + 'matlab-ellipsoid-ode45-restoring-timing-step', waveCase, ... + ['maxstep_' stepTag]); +end +if ~isfolder(outDir) + mkdir(outDir); +end +run('stopWecSim'); +record = output.bodies(1); +bodyValues = [record.time(:), record.position, record.velocity, ... + record.forceTotal, record.forceExcitation, ... + record.forceRadiationDamping, record.forceAddedMass, ... + record.forceRestoring, record.forceMorisonAndViscous, ... + record.acceleration]; +pto = output.ptos(1); +ptoValues = [pto.time(:), pto.position, pto.velocity, ... + pto.forceInternalMechanics, pto.powerInternalMechanics]; +waveValues = [output.wave.time(:), output.wave.elevation(:)]; +expectedRows = round(stopTime / .05) + 1; +assert(isequal(size(bodyValues), [expectedRows, 55]) && ... + isequal(size(ptoValues), [expectedRows, 25]) && ... + isequal(size(waveValues), [expectedRows, 2]) && ... + max(abs(bodyValues(:, 1) - (0:expectedRows-1)' * .05)) < 1e-9 && ... + max(abs(ptoValues(:, 1) - bodyValues(:, 1))) < 1e-9 && ... + max(abs(waveValues(:, 1) - bodyValues(:, 1))) < 1e-9 && ... + all(isfinite(bodyValues), 'all') && ... + all(isfinite(ptoValues), 'all') && ... + all(isfinite(waveValues), 'all'), ... + 'Incomplete published-grid source output'); +writematrix(bodyValues, fullfile(outDir, [mode '_body.csv'])); +writematrix(ptoValues, fullfile(outDir, [mode '_pto.csv'])); +writematrix(waveValues, fullfile(outDir, [mode '_wave.csv'])); +writematrix([simu.dt, simu.dtOut, simu.nonlinearDt, ... + simu.cicDt, str2double(get_param(model, 'MaxStep')), ... + stopTime, compiledPeriod], ... + fullfile(outDir, [mode '_settings.csv'])); +close_system(model, 0); diff --git a/tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py b/tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py new file mode 100644 index 00000000..b8d17cfc --- /dev/null +++ b/tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control.py @@ -0,0 +1,156 @@ +"""Portable 150 s diagnostic for continuous nonlinear force plus mass shift. + +The copied MATLAB model keeps the published 0.05 s clocks and evaluates all +three nonlinear blocks continuously. Moving exactly A_zz into the rigid body +removes active heave acceleration from the 100 ns delay. These source edits +remain separate from the unchanged published-input physical gate. +""" + +from pathlib import Path + +import numpy as np +import pytest +from scipy.integrate import trapezoid + +from test_ellipsoid_ode45_convergence import PYTHON_LIMITS, _python +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro + + +ROOT = Path(__file__).resolve().parents[1] +CONTROL = (ROOT / "tests/fixtures/matlab" + / "ellipsoid_ode45_continuous_nonlinear_mass_control") +TIMING = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_timing" +SOURCE_INPUT = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_convergence" +SOURCE_FIELDS = ("time", "heave", "speed", "restoring", "excitation", + "acceleration", "pto_stroke", "pto_speed", "pto_force", + "pto_power", "wave") +PHYSICAL_FIELDS = ("heave", "speed", "pto_stroke", "pto_speed", "pto_force", + "pto_power") +CONTROL_FIELDS = ("time", "heave", "speed", "total", "excitation", + "radiation", "reported_added_mass", "restoring", "viscous", + "acceleration", "applied_heave_added", "pto_stroke", + "pto_speed", "pto_force", "pto_power", "wave") + + +def _peak(left, right): + assert left.shape == right.shape + assert np.isfinite(left).all() and np.isfinite(right).all() + return float(np.max(np.abs(left - right))) + + +@pytest.mark.parametrize("case", ("Regular", "RegularCIC")) +def test_no_active_delay_control_does_not_close_physical_gap(case): + name = "regular_cic" if case == "RegularCIC" else "regular" + with (np.load(CONTROL / f"{name}.npz", allow_pickle=False) as packed, + np.load(TIMING / f"{name}.npz", allow_pickle=False) as archived): + assert set(packed.files) == ({"case", "settings", "mass", + "inactive_speed", "inactive_acceleration"} + | set(CONTROL_FIELDS)) + assert str(packed["case"]) == case + control = {key: packed[key] for key in packed.files if key != "case"} + old = {field: archived[f"continuous_nonlinear_{field}"] + for field in SOURCE_FIELDS} + + time = np.arange(3001) * .05 + for field in CONTROL_FIELDS: + assert control[field].shape == (3001,), field + assert np.isfinite(control[field]).all(), field + np.testing.assert_allclose(control["time"], time, rtol=0, atol=1e-9) + np.testing.assert_allclose(control["settings"], + [.05, .05, .05, .05, .05, 150, 0, 0, 0], + rtol=0, atol=1e-12) + np.testing.assert_allclose(control["time"], old["time"], + rtol=0, atol=1e-9) + assert _peak(control["wave"], old["wave"]) < 1e-10 + + mass = control["mass"] + assert mass.shape == (12,) + factor, trace, added_zz, original, adjusted, active_residual, step = mass[:7] + np.testing.assert_allclose(factor, added_zz / trace, rtol=0, atol=1e-13) + np.testing.assert_allclose(adjusted, original + added_zz, + rtol=0, atol=1e-7) + assert abs(active_residual) < 1e-7 and abs(step - .05) < 1e-12 + for field in ("inactive_speed", "inactive_acceleration"): + assert control[field].shape == (3001, 5) + assert np.isfinite(control[field]).all() + assert np.max(np.abs(control[field])) < 1e-9 + assert np.max(np.abs(control["inactive_acceleration"] @ mass[7:])) < 1e-6 + assert np.max(np.abs(control["applied_heave_added"])) < 1e-6 + assert np.max(np.abs(control["reported_added_mass"])) < 1e-6 + + assert _peak(control["pto_speed"], control["speed"]) < 1e-9 + assert _peak(control["pto_stroke"], control["heave"] + 2) < 1e-9 + assert _peak(control["pto_force"], + -1_200_000 * control["pto_speed"]) < 1e-6 + assert _peak(control["pto_power"], + -control["pto_force"] * control["pto_speed"]) < 1e-6 + + mesh = HeaveMeshHydro.from_stl( + SOURCE_INPUT / "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, + ) + assert abs(mesh.mass - original) < 1e-7 + forces = np.array([ + mesh.forces(t, z + 2, v) + for t, z, v in zip(control["time"], control["heave"], + control["speed"], strict=True) + ]) + # Source logs these forces with the opposite resisting-force sign. + assert _peak(-control["restoring"], forces[:, 0]) < 1e-6 + assert _peak(-control["viscous"], forces[:, 2]) < 1e-6 + + # These full-grid comparisons quantify the effect of the mass shift on + # the previously archived continuous-nonlinear source arm. + source_change = { + field: _peak(control[field], old[field]) + for field in SOURCE_FIELDS if field not in ("time", "wave") + } + assert source_change["heave"] > .003 + assert source_change["speed"] > .007 + assert source_change["restoring"] > 9000 + assert source_change["excitation"] > 8000 + assert source_change["pto_force"] > 8000 + + # Python independently advances its own 150 s dynamics from the pinned + # HDF5/STL and wave/PTO inputs; no source state or force is fed into it. + python = _python(case, .05, SOURCE_INPUT) + python_body, python_pto = (python.bodies["ellipsoid"], + python.ptos["PTO1"]) + physical = { + "heave": python_body.position[:, 2], + "speed": python_body.velocity[:, 2], + "pto_stroke": python_pto.stroke, + "pto_speed": python_pto.velocity, + "pto_force": python_pto.force, + "pto_power": python_pto.absorbed_power, + } + assert _peak(control["time"], python.time) < 1e-9 + assert _peak(control["wave"], python.wave_elevation) < 1e-10 + old_gaps = {field: _peak(old[field], physical[field]) + for field in PHYSICAL_FIELDS} + control_gaps = {field: _peak(control[field], physical[field]) + for field in PHYSICAL_FIELDS} + limits = PYTHON_LIMITS[case] + for field, limit in zip(("heave", "speed", "pto_force", "pto_power"), + limits, strict=True): + assert control_gaps[field] > limit, (case, field, control_gaps[field]) + assert control_gaps["pto_stroke"] > limits[0] + assert control_gaps["pto_speed"] > limits[1] + for field in PHYSICAL_FIELDS: + assert control_gaps[field] > old_gaps[field], ( + case, field, control_gaps[field], old_gaps[field] + ) + + old_work = float(trapezoid(old["pto_power"], time)) + control_work = float(trapezoid(control["pto_power"], time)) + python_work = float(trapezoid(physical["pto_power"], time)) + assert min(old_work, control_work, python_work) > 1e3 + old_work_gap = abs(old_work / python_work - 1) + control_work_gap = abs(control_work / python_work - 1) + assert control_work_gap > old_work_gap > .001 + print(f"{case} mass-control source change: {source_change}") + print(f"{case} archived/Python peak gaps: {old_gaps}") + print(f"{case} control/Python peak gaps: {control_gaps}") + print(f"{case} source/Python signed-work fractions: " + f"{old_work_gap}, {control_work_gap}") diff --git a/tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py b/tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py new file mode 100644 index 00000000..2c2e1b04 --- /dev/null +++ b/tests/test_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py @@ -0,0 +1,183 @@ +"""Portable post-screen MaxStep diagnosis without active heave delay. + +These copied-model controls follow a failed prospective source-step screen. +Their refined behavior cannot qualify either unchanged published input. +""" + +from pathlib import Path + +import numpy as np +import pytest +from scipy.integrate import trapezoid + +from test_ellipsoid_ode45_convergence import PYTHON_LIMITS, _python +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro + + +ROOT = Path(__file__).resolve().parents[1] +ARCHIVE = (ROOT / "tests/fixtures/matlab" + / "ellipsoid_ode45_continuous_nonlinear_mass_control_step") +BASE = (ROOT / "tests/fixtures/matlab" + / "ellipsoid_ode45_continuous_nonlinear_mass_control") +DELAYED = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_step" +DELAYED_BASE = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_timing" +SOURCE_INPUT = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_convergence" +STEPS = (("0p0125", .0125), ("0p00625", .00625)) +TIME = np.arange(3001) * .05 +CHANNELS = ("heave", "speed", "pto_stroke", "pto_speed", "pto_force", + "pto_power") +FIELDS = ("heave", "speed", "total", "excitation", "radiation", + "reported_added_mass", "restoring", "viscous", "acceleration", + "applied_heave_added", "inactive_speed", "inactive_acceleration", + "pto_stroke", "pto_speed", "pto_force", "pto_power", "wave", + "settings", "mass") + + +def _peak(left, right): + assert left.shape == right.shape + assert np.isfinite(left).all() and np.isfinite(right).all() + return float(np.max(np.abs(left - right))) + + +@pytest.mark.parametrize("case", ("Regular", "RegularCIC")) +def test_post_screen_no_active_delay_gaps_contract_but_work_remains_outside_gate(case): + filename = "regular_cic.npz" if case == "RegularCIC" else "regular.npz" + with (np.load(BASE / filename, allow_pickle=False) as a, + np.load(ARCHIVE / filename, allow_pickle=False) as fine, + np.load(DELAYED_BASE / filename, allow_pickle=False) as b, + np.load(DELAYED / filename, allow_pickle=False) as delayed): + assert set(fine.files) == ({"case", "time"} | { + f"maxstep_{tag}_{field}" for tag, _ in STEPS for field in FIELDS + }) + assert str(fine["case"]) == str(a["case"]) == case + np.testing.assert_allclose(fine["time"], TIME, rtol=0, atol=1e-9) + base = {field: a[field] for field in FIELDS} + base["time"] = a["time"] + delayed_base = { + field: b[f"continuous_nonlinear_{field}"] + for field in CHANNELS + ("wave",) + } + controls = { + tag: {field: fine[f"maxstep_{tag}_{field}"] for field in FIELDS} + for tag, _ in STEPS + } + paired = { + tag: { + field: delayed[f"maxstep_{tag}_{field}"] + for field in CHANNELS + ("wave",) + } + for tag, _ in STEPS + } + + np.testing.assert_allclose(base["time"], TIME, rtol=0, atol=1e-9) + np.testing.assert_allclose(base["settings"], + [.05, .05, .05, .05, .05, 150, 0, 0, 0], + rtol=0, atol=1e-12) + mesh = HeaveMeshHydro.from_stl( + SOURCE_INPUT / "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, + ) + assert abs(mesh.mass - base["mass"][3]) < 1e-7 + records = {"0p05": base, **controls} + matches = {"0p05": delayed_base, **paired} + for tag, step in STEPS: + record = records[tag] + np.testing.assert_allclose(record["settings"], + [.05, .05, .05, .05, step, 150, 0, 0, 0], + rtol=0, atol=1e-12) + for field in FIELDS: + shape = (3001, 5) if field.startswith("inactive_") else ( + (9,) if field == "settings" else (12,) if field == "mass" + else (3001,) + ) + assert record[field].shape == shape, (case, tag, field) + assert np.isfinite(record[field]).all(), (case, tag, field) + mass = record["mass"] + np.testing.assert_allclose(np.delete(mass, 6), + np.delete(base["mass"], 6), + rtol=0, atol=1e-7) + assert abs(mass[5]) < 1e-7 and abs(mass[6] - step) < 1e-12 + assert np.max(np.abs(record["inactive_speed"])) < 1e-9 + assert np.max(np.abs(record["inactive_acceleration"])) < 1e-9 + assert np.max(np.abs(record["inactive_acceleration"] @ mass[7:])) < 1e-6 + assert np.max(np.abs(record["applied_heave_added"])) < 1e-6 + assert np.max(np.abs(record["reported_added_mass"])) < 1e-6 + assert _peak(record["wave"], base["wave"]) < 1e-10 + assert _peak(record["wave"], matches[tag]["wave"]) < 1e-10 + assert _peak(record["pto_stroke"], record["heave"] + 2) < 1e-9 + assert _peak(record["pto_speed"], record["speed"]) < 1e-9 + assert _peak(record["pto_force"], + -1_200_000 * record["pto_speed"]) < 1e-6 + assert _peak(record["pto_power"], + -record["pto_force"] * record["pto_speed"]) < 1e-6 + expected = np.array([ + mesh.forces(t, z + 2, v) + for t, z, v in zip(TIME, record["heave"], record["speed"], + strict=True) + ]) + assert _peak(-record["restoring"], expected[:, 0]) < 1e-6 + assert _peak(-record["viscous"], expected[:, 2]) < 1e-6 + + # Python's dynamics use only the pinned HDF5/STL and wave/PTO inputs. + python = _python(case, .05, SOURCE_INPUT) + body, pto = python.bodies["ellipsoid"], python.ptos["PTO1"] + physical = { + "heave": body.position[:, 2], "speed": body.velocity[:, 2], + "pto_stroke": pto.stroke, "pto_speed": pto.velocity, + "pto_force": pto.force, "pto_power": pto.absorbed_power, + } + assert _peak(python.time, TIME) < 1e-9 + assert _peak(python.wave_elevation, base["wave"]) < 1e-10 + control_gaps = { + tag: {field: _peak(record[field], physical[field]) + for field in CHANNELS} + for tag, record in records.items() + } + paired_gaps = { + tag: {field: _peak(records[tag][field], matches[tag][field]) + for field in CHANNELS} + for tag in records + } + control_work = { + tag: float(trapezoid(record["pto_power"], TIME)) + for tag, record in records.items() + } + paired_work = { + tag: float(trapezoid(matches[tag]["pto_power"], TIME)) + for tag in records + } + python_work = float(trapezoid(physical["pto_power"], TIME)) + assert min(*control_work.values(), *paired_work.values(), python_work) > 1e3 + work_gaps = {tag: abs(work/python_work - 1) + for tag, work in control_work.items()} + paired_work_gaps = {tag: abs(control_work[tag]/paired_work[tag] - 1) + for tag in records} + + order = ("0p05", "0p0125", "0p00625") + for field in CHANNELS: + assert all(control_gaps[later][field] < control_gaps[earlier][field] + for earlier, later in zip(order, order[1:])), (case, field) + assert all(paired_gaps[later][field] < paired_gaps[earlier][field] + for earlier, later in zip(order, order[1:])), (case, field) + assert all(work_gaps[later] < work_gaps[earlier] + for earlier, later in zip(order, order[1:])) + assert all(paired_work_gaps[later] < paired_work_gaps[earlier] + for earlier, later in zip(order, order[1:])) + + # The edited post-screen arms satisfy pointwise inherited limits, while + # their signed work still exceeds 0.1%. The original gate is unchanged. + limits = dict(zip(("heave", "speed", "pto_force", "pto_power"), + PYTHON_LIMITS[case], strict=True)) + limits["pto_stroke"] = limits["heave"] + limits["pto_speed"] = limits["speed"] + for tag, _ in STEPS: + for field, limit in limits.items(): + assert control_gaps[tag][field] < limit, ( + case, tag, field, control_gaps[tag][field], limit + ) + assert work_gaps[tag] > .001, (case, tag, work_gaps[tag]) + print(f"{case} no-delay/Python peak gaps: {control_gaps}") + print(f"{case} no-delay/delayed peak gaps: {paired_gaps}") + print(f"{case} no-delay/Python signed-work fractions: {work_gaps}; " + f"no-delay/delayed fractions: {paired_work_gaps}") diff --git a/tests/test_ellipsoid_ode45_restoring_step_portable.py b/tests/test_ellipsoid_ode45_restoring_step_portable.py new file mode 100644 index 00000000..7f33aa49 --- /dev/null +++ b/tests/test_ellipsoid_ode45_restoring_step_portable.py @@ -0,0 +1,168 @@ +"""Portable copied-model MaxStep evidence after the failed 0.0125 s screen. + +The 0.00625 and 0.003125 s arms are later diagnostics. Their comparison with +Python does not qualify the unchanged published 0.05 s source input. +""" + +from pathlib import Path + +import numpy as np +import pytest +from scipy.integrate import trapezoid + +from test_ellipsoid_ode45_convergence import _python +from tools.check_ellipsoid_ode45_restoring_step import ( + PHYSICAL_LIMITS, SOURCE_LIMITS, +) +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro + + +ROOT = Path(__file__).resolve().parents[1] +ARCHIVE = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_step" +BASE = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_timing" +SOURCE_INPUT = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_convergence" +STEPS = (("0p025", .025), ("0p0125", .0125), + ("0p00625", .00625), ("0p003125", .003125)) +FIELDS = ("heave", "speed", "restoring", "pto_stroke", "pto_speed", + "pto_force", "pto_power", "wave", "settings") +CHANNELS = ("heave", "speed", "pto_force", "pto_power") +PREDECLARED_FAILURES = { + "Regular": {"speed", "pto_force", "pto_power", "signed_work"}, + "RegularCIC": {"heave", "speed", "pto_force", "pto_power", + "signed_work"}, +} + + +def _peak(left, right): + assert left.shape == right.shape + assert np.isfinite(left).all() and np.isfinite(right).all() + return float(np.max(np.abs(left - right))) + + +@pytest.mark.parametrize("case", ("Regular", "RegularCIC")) +def test_copied_model_maxstep_continuation(case): + name = "regular_cic" if case == "RegularCIC" else "regular" + with np.load(BASE / f"{name}.npz", allow_pickle=False) as archive: + assert str(archive["case"]) == case + base = {field: archive[f"continuous_nonlinear_{field}"] + for field in FIELDS} + base["time"] = archive["continuous_nonlinear_time"] + with np.load(ARCHIVE / f"{name}.npz", allow_pickle=False) as archive: + expected = {"case", "time"} | { + f"maxstep_{tag}_{field}" for tag, _ in STEPS for field in FIELDS + } + assert set(archive.files) == expected + assert str(archive["case"]) == case + time = archive["time"] + records = { + tag: {field: archive[f"maxstep_{tag}_{field}"] + for field in FIELDS} + for tag, _ in STEPS + } + + assert time.shape == (3001,) + np.testing.assert_allclose(time, np.arange(3001) * .05, + rtol=0, atol=1e-9) + np.testing.assert_allclose(base["time"], time, rtol=0, atol=1e-9) + np.testing.assert_allclose(base["settings"], + [.05, .05, .05, .05, .05, 150, 0, 0, 0], + rtol=0, atol=1e-12) + + mesh = HeaveMeshHydro.from_stl( + SOURCE_INPUT / "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, + ) + for tag, step in STEPS: + record = records[tag] + np.testing.assert_allclose( + record["settings"], [.05, .05, .05, .05, step, 150, 0, 0, 0], + rtol=0, atol=1e-12, + ) + assert record["settings"].shape == (9,) + for field in FIELDS[:-1]: + assert record[field].shape == (3001,), (case, tag, field) + assert np.isfinite(record[field]).all(), (case, tag, field) + assert _peak(record["wave"], base["wave"]) < 1e-10 + assert _peak(record["pto_stroke"], record["heave"] + 2) < 1e-10 + assert _peak(record["pto_speed"], record["speed"]) < 1e-9 + assert _peak(record["pto_force"], + -1_200_000 * record["pto_speed"]) < 1e-6 + assert _peak(record["pto_power"], + -record["pto_force"] * record["pto_speed"]) < 1e-6 + expected_restoring = np.array([ + mesh.forces(t, z + 2, v)[0] + for t, z, v in zip(time, record["heave"], record["speed"], + strict=True) + ]) + assert _peak(-record["restoring"], expected_restoring) < 1e-6 + + histories = [base] + [records[tag] for tag, _ in STEPS] + source_changes = { + field: [_peak(a[field], b[field]) + for a, b in zip(histories, histories[1:])] + for field in CHANNELS + } + work = [float(trapezoid(record["pto_power"], time)) + for record in histories] + work_changes = [abs(later - earlier) / later + for earlier, later in zip(work, work[1:])] + + # The .025 -> .0125 result was the prospective screen. Preserve its + # observed failure, including the signed-work limit, as historical fact. + failures = set() + for field, limit in zip(CHANNELS, SOURCE_LIMITS[case], strict=True): + coarse, screened = source_changes[field][:2] + assert screened < .7 * coarse + if screened >= limit: + failures.add(field) + assert work_changes[1] < .7 * work_changes[0] + if work_changes[1] >= .001: + failures.add("signed_work") + assert failures == PREDECLARED_FAILURES[case] + print(f"{case} predeclared 0.0125 s screen failures: {sorted(failures)}") + print(f"{case} later source changes: {source_changes}; " + f"signed-work fractions: {work_changes}") + + # The later continuation has its own observed result. Checking it here + # keeps that result reproducible without reclassifying the failed screen. + for field, limit in zip(CHANNELS, SOURCE_LIMITS[case], strict=True): + assert all(later < .7 * earlier for earlier, later in zip( + source_changes[field], source_changes[field][1:] + )) + assert source_changes[field][-1] < limit + assert all(later < .7 * earlier for earlier, later in zip( + work_changes, work_changes[1:] + )) + assert work_changes[-1] < .001 + + # Independent Python dynamics use only the archived source HDF5/STL and + # configured wave/PTO. The last arm is a post-screen physical diagnostic. + python = _python(case, .05, SOURCE_INPUT) + body = python.bodies["ellipsoid"] + pto = python.ptos["PTO1"] + physical = { + "heave": body.position[:, 2], + "speed": body.velocity[:, 2], + "pto_stroke": pto.stroke, + "pto_speed": pto.velocity, + "pto_force": pto.force, + "pto_power": pto.absorbed_power, + } + assert _peak(python.time, time) < 1e-9 + assert _peak(python.wave_elevation, base["wave"]) < 1e-10 + for tag, _ in STEPS: + assert _peak(records[tag]["wave"], python.wave_elevation) < 1e-10 + final = records["0p003125"] + gaps = {field: _peak(final[field], values) + for field, values in physical.items()} + limits = dict(zip(CHANNELS, PHYSICAL_LIMITS[case], strict=True)) + limits["pto_stroke"] = limits["heave"] + limits["pto_speed"] = limits["speed"] + for field, limit in limits.items(): + assert gaps[field] < limit, (case, field, gaps[field], limit) + python_work = float(trapezoid(physical["pto_power"], time)) + final_work_error = abs(work[-1] / python_work - 1) + assert final_work_error < .001 + print(f"{case} post-screen 0.003125 s Python gaps: {gaps}; " + f"signed-work fraction: {final_work_error}") diff --git a/tests/test_ellipsoid_ode45_restoring_timing_portable.py b/tests/test_ellipsoid_ode45_restoring_timing_portable.py new file mode 100644 index 00000000..32f6d25e --- /dev/null +++ b/tests/test_ellipsoid_ode45_restoring_timing_portable.py @@ -0,0 +1,204 @@ +"""Portable 150 s screen for the copied ellipsoid source force-timing arms. + +The unchanged published 0.05 s source remains the separate original-input +gate. These edited-source arms are diagnostics; their inherited physical +motion, PTO, and work limits still fail. +""" + +from pathlib import Path + +import numpy as np +import pytest +from scipy.integrate import trapezoid + +from test_ellipsoid_ode45_convergence import PYTHON_LIMITS, _python +from wecsim.hydrodynamics.bodyClass import BodyClass +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro + + +ROOT = Path(__file__).resolve().parents[1] +FIXTURE = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_timing" +SOURCE_INPUT = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_convergence" +ORIGINAL = (ROOT / "tests/fixtures/matlab/ellipsoid_ode45_nonlinear_clock" + / "source_maxstep_common_grid.npz") +MODES = ("baseline", "continuous_restoring", "continuous_nonlinear") +CHANNELS = ("heave", "speed", "pto_stroke", "pto_speed", "pto_force", + "pto_power") +FIELDS = ("time", "heave", "speed", "restoring", "excitation", + "acceleration", "pto_stroke", "pto_speed", "pto_force", + "pto_power", "wave", "settings") +SAMPLE_PERIODS = { + "baseline": (.05, .05, .05), + "continuous_restoring": (0., .05, .05), + "continuous_nonlinear": (0., 0., 0.), +} + + +def _peak_gap(actual, expected): + assert actual.shape == expected.shape + assert np.isfinite(actual).all() and np.isfinite(expected).all() + return float(np.max(np.abs(actual - expected))) + + +@pytest.mark.parametrize("case", ("Regular", "RegularCIC")) +def test_full_duration_copied_source_force_timing_diagnostic(case): + name = "regular_cic" if case == "RegularCIC" else "regular" + with np.load(FIXTURE / f"{name}.npz", allow_pickle=False) as archive: + expected_keys = {"case"} | { + f"{mode}_{field}" for mode in MODES for field in FIELDS + } + assert set(archive.files) == expected_keys + assert str(archive["case"]) == case + records = { + mode: {field: archive[f"{mode}_{field}"] for field in FIELDS} + for mode in MODES + } + + times = np.arange(3001) * .05 + for mode, record in records.items(): + np.testing.assert_allclose(record["settings"], + [.05, .05, .05, .05, .05, 150, + *SAMPLE_PERIODS[mode]], + rtol=0, atol=1e-12) + assert record["settings"].shape == (9,) + for field in FIELDS[:-1]: + assert record[field].shape == (3001,), (case, mode, field) + assert np.isfinite(record[field]).all(), (case, mode, field) + np.testing.assert_allclose(record["time"], times, rtol=0, atol=1e-9) + assert _peak_gap(record["pto_speed"], record["speed"]) < 1e-9 + assert _peak_gap(record["pto_stroke"], record["heave"] + 2) < 1e-10 + assert _peak_gap(record["pto_force"], + -1_200_000 * record["pto_speed"]) < 1e-6 + assert _peak_gap(record["pto_power"], + -record["pto_force"] * record["pto_speed"]) < 1e-6 + assert _peak_gap(record["wave"], records["baseline"]["wave"]) < 1e-10 + + # The fresh copied-model baseline reproduces the unchanged published run, + # independently archived before this timing intervention. + with np.load(ORIGINAL, allow_pickle=False) as original: + baseline = records["baseline"] + for field, limit in (("time", 1e-10), ("heave", 1e-7), + ("speed", 1e-7), ("restoring", .1), + ("excitation", .1), ("pto_force", .1), + ("wave", 1e-10)): + assert _peak_gap( + baseline[field], original[f"{name}_published_{field}" + ]) < limit, (case, field) + np.testing.assert_allclose(baseline["settings"][:5], + original[f"{name}_published_settings"], + rtol=0, atol=1e-12) + + mesh = HeaveMeshHydro.from_stl( + SOURCE_INPUT / "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, + ) + bem = BodyClass(str(SOURCE_INPUT / "hydroData/ellipsoid.h5")) + bem.bodyNumber = bem.bodyTotal = 1 + bem.readH5file() + bem.mass = mesh.mass + bem.inertia = [1.375264e6, 1.375264e6, 1.341721e6] + bem.hydroStiffness = np.zeros((6, 6)) + bem.viscDrag = {"Drag": np.zeros((6, 6)), "cd": np.zeros(6), + "characteristicArea": np.zeros(6)} + bem.linearDamping = np.zeros((6, 6)) + omega = 2 * np.pi / 6 + cic = case == "RegularCIC" + convolution_time = np.arange(1201) * .05 if cic else np.array([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 times]) + linear_excitation = 2 * ramp * ( + re * np.cos(omega * times) - im * np.sin(omega * times) + ) + incident_elevation = mesh.amplitude * ramp * np.cos(omega * times) + for mode, record in records.items(): + mesh_forces = np.array([ + mesh.forces(t, z + 2, v) + for t, z, v in zip(record["time"], record["heave"], + record["speed"], strict=True) + ]) + applied = -record["restoring"] + previous_error = _peak_gap(applied[1:], mesh_forces[:-1, 0]) + current_error = _peak_gap(applied[1:], mesh_forces[1:, 0]) + if mode == "baseline": + assert previous_error < 1e-6 and current_error > 500 + else: + assert current_error < 1e-6 and previous_error > 500 + + # The logged excitation and wave use the current output state even + # in the baseline with 0.05 s sampled nonlinear blocks. This rules + # out a one-output-sample lag in the saved signals, but says nothing + # about force evaluation at ode45 minor steps between saved samples. + assert _peak_gap(record["wave"], incident_elevation) < 1e-10 + assert _peak_gap(record["wave"][1:], incident_elevation[:-1]) > .05 + current_excitation = linear_excitation + mesh_forces[:, 1] + assert _peak_gap(record["excitation"], current_excitation) < 1e-5 + if mode == "continuous_nonlinear": + previous = np.maximum(np.arange(len(times)) - 1, 0) + previous_pose_fk = np.array([ + mesh.forces(t, record["heave"][i] + 2, + record["speed"][i])[1] + for t, i in zip(times, previous, strict=True) + ]) + assert _peak_gap( + record["excitation"][1:], + (linear_excitation + previous_pose_fk)[1:], + ) > 10_000 + assert _peak_gap( + record["excitation"][1:], + (linear_excitation + mesh_forces[previous, 1])[1:], + ) > 10_000 + + # Python advances its own dynamics using the archived source-generated + # HDF5/STL and published wave/PTO input; no source motion or force is fed in. + python = _python(case, .05, SOURCE_INPUT) + python_body = python.bodies["ellipsoid"] + python_pto = python.ptos["PTO1"] + physical = { + "heave": python_body.position[:, 2], + "speed": python_body.velocity[:, 2], + "pto_stroke": python_pto.stroke, + "pto_speed": python_pto.velocity, + "pto_force": python_pto.force, + "pto_power": python_pto.absorbed_power, + } + assert _peak_gap(python.time, times) < 1e-9 + assert _peak_gap(python.wave_elevation, baseline["wave"]) < 1e-10 + python_work = float(trapezoid(physical["pto_power"], python.time)) + assert python_work > 1e3 + + gaps = { + mode: {field: _peak_gap(record[field], physical[field]) + for field in CHANNELS} + for mode, record in records.items() + } + work_gaps = { + mode: abs(float(trapezoid(record["pto_power"], times)) + / python_work - 1) + for mode, record in records.items() + } + pointwise_limits = PYTHON_LIMITS[case] + for mode in MODES[1:]: + # These are regression screens for the observed source improvement, + # not prospective acceptance limits. Both edited source routes still + # fail every inherited full-duration physical gate. + for field in CHANNELS: + assert gaps[mode][field] < .5 * gaps["baseline"][field], ( + case, mode, field, gaps[mode][field], gaps["baseline"][field] + ) + assert work_gaps[mode] < .5 * work_gaps["baseline"] + for field, limit in zip( + ("heave", "speed", "pto_force", "pto_power"), + pointwise_limits, strict=True, + ): + assert gaps[mode][field] > limit, (case, mode, field) + assert gaps[mode]["pto_stroke"] > pointwise_limits[0] + assert gaps[mode]["pto_speed"] > pointwise_limits[1] + assert work_gaps[mode] > .001 + print(f"{case} edited-source peak gaps: {gaps}; work fractions: {work_gaps}") diff --git a/tools/check_ellipsoid_ode45_restoring_step.py b/tools/check_ellipsoid_ode45_restoring_step.py new file mode 100644 index 00000000..7ecb40a6 --- /dev/null +++ b/tools/check_ellipsoid_ode45_restoring_step.py @@ -0,0 +1,178 @@ +"""Screen copied-model continuous nonlinear ellipsoid source over ode45 steps. + +The source input, 0.05 s output, nonlinear, and CIC clocks are held fixed. +Only the copied model's ode45 MaxStep changes from 0.05 to 0.025 to 0.0125 s. +This derived screen does not qualify the unchanged published input. +""" + +import argparse +import json +from pathlib import Path + +import numpy as np +from scipy.integrate import trapezoid + +from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro + + +SOURCE_LIMITS = { + "Regular": (.001, .0013, 1600, 1300), + "RegularCIC": (.0011, .0015, 2000, 1600), +} +PHYSICAL_LIMITS = { + "Regular": (.003, .004, 5000, 4000), + "RegularCIC": (.0035, .0045, 6000, 5000), +} +CHANNELS = ("heave", "speed", "pto_force", "pto_power") + + +def _read(folder, name, columns): + values = np.loadtxt(folder / name, delimiter=",", ndmin=2) + assert values.shape == (3001, columns), (folder, name, values.shape) + assert np.isfinite(values).all(), (folder, name) + return values + + +def _max_error(left, right): + assert left.shape == right.shape + return float(np.max(np.abs(left - right))) + + +def _source(folder, step, mesh): + body = _read(folder, "continuous_nonlinear_body.csv", 55) + pto = _read(folder, "continuous_nonlinear_pto.csv", 25) + wave = _read(folder, "continuous_nonlinear_wave.csv", 2) + settings = np.loadtxt( + folder / "continuous_nonlinear_settings.csv", delimiter=",", ndmin=2) + np.testing.assert_allclose( + settings, [[.05, .05, .05, .05, step, 150, 0, 0, 0]], + rtol=0, atol=1e-12, + ) + time = np.arange(3001) * .05 + for array in (body, pto, wave): + assert _max_error(array[:, 0], time) < 1e-9 + heave, speed, restoring = body[:, 3], body[:, 9], -body[:, 39] + expected = np.array([ + mesh.forces(t, z + 2, v)[0] + for t, z, v in zip(time, heave, speed, strict=True) + ]) + assert _max_error(restoring, expected) < 1e-6, ( + "Copied-model restoring force no longer follows current state", folder + ) + assert _max_error(pto[:, 9], speed) < 1e-9 + assert _max_error(pto[:, 15], -1_200_000 * pto[:, 9]) < 1e-6 + return { + "time": time, + "heave": heave, + "speed": speed, + "pto_force": pto[:, 15], + "pto_power": -pto[:, 15] * pto[:, 9], + "wave": wave[:, 1], + } + + +def _python(case, app): + wec = WEC("ellipsoid continuous nonlinear source-step screen") + 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) + output = wec.run( + wave, dt=.05, end_time=150, ramp_time=50, rho=1025, + radiation_memory=60 if case == "RegularCIC" else None, + ) + return { + "time": output.time, + "heave": output.bodies["ellipsoid"].position[:, 2], + "speed": output.bodies["ellipsoid"].velocity[:, 2], + "pto_force": output.ptos["PTO1"].force, + "pto_power": output.ptos["PTO1"].absorbed_power, + "wave": output.wave_elevation, + } + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("case", choices=("Regular", "RegularCIC")) + parser.add_argument("base", type=Path) + parser.add_argument("step025", type=Path) + parser.add_argument("step0125", type=Path) + parser.add_argument("--applications", type=Path, default=Path("applications")) + args = parser.parse_args() + app = args.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, + ) + labels = ("maxstep_0p05", "maxstep_0p025", "maxstep_0p0125") + source = { + label: _source(folder, step, mesh) + for label, folder, step in zip( + labels, (args.base, args.step025, args.step0125), + (.05, .025, .0125), strict=True + ) + } + reference = _python(args.case, app) + first, middle, final = (source[label] for label in labels) + for history in (first, middle, final): + assert _max_error(history["wave"], first["wave"]) < 1e-10 + assert _max_error(history["wave"], reference["wave"]) < 1e-10 + summary = {"case": args.case, "source_step": {}, "physical": {}} + source_ok = True + physical_ok = True + for name, source_limit, physical_limit in zip( + CHANNELS, SOURCE_LIMITS[args.case], PHYSICAL_LIMITS[args.case], + strict=True, + ): + coarse = _max_error(first[name], middle[name]) + fine = _max_error(middle[name], final[name]) + physical = _max_error(final[name], reference[name]) + summary["source_step"][name] = { + "0p05_to_0p025": coarse, + "0p025_to_0p0125": fine, + "final_limit": source_limit, + "contracts": fine < .7 * coarse, + "passes_final_limit": fine < source_limit, + } + summary["physical"][name] = { + "final_error": physical, + "limit": physical_limit, + "passes": physical < physical_limit, + } + source_ok &= fine < source_limit and fine < .7 * coarse + physical_ok &= physical < physical_limit + source_work = [float(trapezoid(history["pto_power"], history["time"])) + for history in (first, middle, final)] + python_work = float(trapezoid(reference["pto_power"], reference["time"])) + work_change = [abs(later - earlier) / later + for earlier, later in zip(source_work, source_work[1:])] + work_error = abs(source_work[-1] / python_work - 1) + summary["source_step"]["signed_work"] = { + "work_J": source_work, + "relative_changes": work_change, + "passes_final": work_change[-1] < .001, + "contracts": work_change[-1] < .7 * work_change[0], + } + summary["physical"]["signed_work"] = { + "final_error_fraction": work_error, + "limit": .001, + "passes": work_error < .001, + } + source_ok &= work_change[-1] < .001 and work_change[-1] < .7 * work_change[0] + physical_ok &= work_error < .001 + summary["passes_source_step_screen"] = bool(source_ok) + summary["passes_inherited_physical_limits"] = bool(physical_ok) + print(json.dumps(summary, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/tools/check_ellipsoid_ode45_restoring_timing.py b/tools/check_ellipsoid_ode45_restoring_timing.py new file mode 100644 index 00000000..75508e08 --- /dev/null +++ b/tools/check_ellipsoid_ode45_restoring_timing.py @@ -0,0 +1,191 @@ +"""Check the causal force timing in a disposable ellipsoid source-model run. + +This experiment changes only sampled subsystems in a copied MATLAB model. +It does not change or qualify the published source input. +""" + +import argparse +import json +from pathlib import Path + +import numpy as np +from scipy.integrate import trapezoid + +from wecsim import RegularCICWave, RegularWave, WEC, WorldPoint +from wecsim.hydrodynamics.nonlinearHydro import HeaveMeshHydro + + +def _read(folder, mode, suffix, columns, rows): + array = np.loadtxt(folder / f"{mode}_{suffix}.csv", delimiter=",", ndmin=2) + assert array.shape == (rows, columns), (mode, suffix, array.shape) + assert np.isfinite(array).all(), (mode, suffix) + return array + + +def _max_error(a, b): + assert a.shape == b.shape + return float(np.max(np.abs(a - b))) + + +def _python_physical(wave_case, app): + wec = WEC("ellipsoid force-timing source experiment") + 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 wave_case == "RegularCIC" + else RegularWave(4, 6)) + return wec.run( + wave, dt=.05, end_time=150, ramp_time=50, rho=1025, + radiation_memory=60 if wave_case == "RegularCIC" else None, + ) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("wave", choices=("Regular", "RegularCIC")) + parser.add_argument("folder", type=Path) + parser.add_argument("--duration", choices=("smoke", "full"), default="smoke") + parser.add_argument("--applications", type=Path, default=Path("applications")) + parser.add_argument("--modes", nargs="+", default=( + "baseline", "continuous_restoring", "continuous_nonlinear")) + args = parser.parse_args() + stop = 2 if args.duration == "smoke" else 150 + rows = round(stop / .05) + 1 + geometry = args.applications / "Nonlinear_Hydro/geometry/elipsoid.stl" + mesh = HeaveMeshHydro.from_stl( + geometry, 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, + ) + result = {} + for mode in args.modes: + body = _read(args.folder, mode, "body", 55, rows) + pto = _read(args.folder, mode, "pto", 25, rows) + wave = _read(args.folder, mode, "wave", 2, rows) + settings = _read(args.folder, mode, "settings", 9, 1)[0] + expected_periods = { + "baseline": (.05, .05, .05), + "continuous_restoring": (0, .05, .05), + "continuous_nonlinear": (0, 0, 0), + }[mode] + np.testing.assert_allclose( + settings, [.05, .05, .05, .05, .05, stop, *expected_periods], + rtol=0, atol=1e-12, + ) + times = np.arange(rows) * .05 + for name, source in (("body", body), ("pto", pto), ("wave", wave)): + assert _max_error(source[:, 0], times) < 1e-9, name + heave = body[:, 3] + speed = body[:, 9] + source_restoring = -body[:, 39] + current = np.array([ + mesh.forces(t, z + 2, v)[0] + for t, z, v in zip(times, heave, speed, strict=True) + ]) + previous = np.array([ + mesh.forces(t, z + 2, v)[0] + for t, z, v in zip(times[:-1], heave[:-1], speed[:-1], strict=True) + ]) + result[mode] = { + "compiled_sample_periods_s": tuple(float(x) for x in settings[-3:]), + "restoring_vs_current_N": _max_error(source_restoring[1:], current[1:]), + "restoring_vs_previous_N": _max_error(source_restoring[1:], previous), + "max_heave_m": float(np.max(np.abs(heave))), + "max_speed_mps": float(np.max(np.abs(speed))), + "wave": wave[:, 1], + "heave": heave, + "speed": speed, + "pto_force": pto[:, 15], + "pto_power": -pto[:, 15] * pto[:, 9], + } + baseline = result["baseline"] + assert baseline["restoring_vs_previous_N"] < 1e-6, ( + "The copied baseline no longer has the archived preceding-state law" + ) + assert baseline["restoring_vs_current_N"] > 5, ( + "The 2 s baseline is too short to distinguish the force clocks" + ) + for mode in args.modes: + if mode == "baseline": + continue + assert result[mode]["restoring_vs_current_N"] < 1e-6, ( + "The edited restoring force no longer follows the current state", + mode, + ) + assert result[mode]["restoring_vs_previous_N"] > 5, ( + "The edited arm is too short to distinguish the force clocks", + mode, + ) + archive = (Path(__file__).resolve().parents[1] + / "tests/fixtures/matlab/ellipsoid_ode45_nonlinear_clock" + / "source_maxstep_common_grid.npz") + prefix = "regular_cic" if args.wave == "RegularCIC" else "regular" + with np.load(archive, allow_pickle=False) as original: + for key, limit in (("heave", 1e-7), ("speed", 1e-7), + ("pto_force", .1), ("wave", 1e-10)): + error = _max_error( + baseline[key], original[f"{prefix}_published_{key}"][:rows]) + assert error < limit, ("baseline published trace", key, error) + summary = {} + if args.duration == "full": + python = _python_physical(args.wave, + args.applications / "Nonlinear_Hydro") + python_body = python.bodies["ellipsoid"] + python_pto = python.ptos["PTO1"] + physical = { + "heave": python_body.position[:, 2], + "speed": python_body.velocity[:, 2], + "pto_force": python_pto.force, + "pto_power": python_pto.absorbed_power, + "wave": python.wave_elevation, + } + assert _max_error(physical["wave"], baseline["wave"]) < 1e-10 + for mode, data in result.items(): + summary[mode] = {key: value for key, value in data.items() + if not isinstance(value, np.ndarray)} + if mode != "baseline": + summary[mode].update({ + "vs_baseline_heave_m": _max_error(data["heave"], baseline["heave"]), + "vs_baseline_speed_mps": _max_error(data["speed"], baseline["speed"]), + "vs_baseline_pto_force_N": _max_error( + data["pto_force"], baseline["pto_force"]), + "vs_baseline_pto_power_W": _max_error( + data["pto_power"], baseline["pto_power"]), + "vs_baseline_wave_m": _max_error(data["wave"], baseline["wave"]), + }) + assert summary[mode]["vs_baseline_wave_m"] < 1e-10, ( + "The incident wave changed in a force-timing arm", mode + ) + if args.duration == "full": + summary[mode].update({ + f"vs_python_{key}": _max_error(data[key], physical[key]) + for key in ("heave", "speed", "pto_force", "pto_power") + }) + source_work = float(trapezoid(data["pto_power"], times)) + python_work = float(trapezoid(physical["pto_power"], times)) + summary[mode]["source_signed_work_J"] = source_work + summary[mode]["vs_python_signed_work_fraction"] = abs( + source_work / python_work - 1) + motion_limits = ((.0035, .0045, 6000, 5000) + if args.wave == "RegularCIC" + else (.003, .004, 5000, 4000)) + summary[mode]["passes_inherited_motion_and_work_limits"] = ( + all(summary[mode][f"vs_python_{key}"] < limit + for key, limit in zip( + ("heave", "speed", "pto_force", "pto_power"), + motion_limits, strict=True)) + and summary[mode]["vs_python_signed_work_fraction"] < .001 + ) + print(json.dumps({"wave": args.wave, "duration_s": stop, + "results": summary}, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py b/tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py new file mode 100644 index 00000000..18a7ad5b --- /dev/null +++ b/tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py @@ -0,0 +1,128 @@ +"""Pack a full 150 s ellipsoid continuous-nonlinear mass-control artifact. + +Usage: + python tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control.py \ + CASE CSV_DIR OUTPUT.npz [--expect ARCHIVED.npz] [--source-hydro H5] + +The MATLAB workflow copies the application input and model. This script reads +its output, checks the matched clocks and no-active-heave-feedback identity, +and stores the source traces needed by the portable diagnostic gate. +""" + +import argparse +from pathlib import Path + +import h5py +import numpy as np + + +ROWS = 3001 +TIMES = np.arange(ROWS) * .05 +ROOT = Path(__file__).resolve().parents[1] +PINNED_HYDRO = (ROOT / "tests/fixtures/matlab/ellipsoid_ode45_convergence" + / "hydroData/ellipsoid.h5") + + +def _read(folder: Path, name: str, shape: tuple[int, int]) -> np.ndarray: + values = np.loadtxt(folder / f"{name}.csv", delimiter=",", ndmin=2) + assert values.shape == shape, (name, values.shape) + assert np.isfinite(values).all(), name + return values + + +def _datasets(group, prefix=""): + result = {} + for name, item in group.items(): + path = f"{prefix}/{name}" + if isinstance(item, h5py.Dataset): + result[path] = item[()] + else: + result.update(_datasets(item, path)) + return result + + +def _check_hydro(path: Path) -> None: + with h5py.File(path) as fresh, h5py.File(PINNED_HYDRO) as archived: + fresh_data, archived_data = _datasets(fresh), _datasets(archived) + assert len(fresh_data) == 38 and fresh_data.keys() == archived_data.keys() + for name in fresh_data: + np.testing.assert_equal(fresh_data[name], archived_data[name], + err_msg=f"source HDF5 dataset {name} changed") + print(f"Matched 38 pinned HDF5 datasets in {path}") + + +def pack(case: str, folder: Path, output: Path, + expected: Path | None = None, source_hydro: Path | None = None) -> None: + if case not in ("Regular", "RegularCIC"): + raise ValueError(f"Unsupported ellipsoid wave case: {case}") + if source_hydro is not None: + _check_hydro(source_hydro) + + body = _read(folder, "body", (ROWS, 55)) + pto = _read(folder, "pto", (ROWS, 25)) + wave = _read(folder, "wave", (ROWS, 2)) + applied = _read(folder, "applied_added", (ROWS, 7)) + settings = _read(folder, "settings", (1, 9))[0] + mass = _read(folder, "mass", (1, 12))[0] + np.testing.assert_allclose(settings, [.05]*5 + [150, 0, 0, 0], + rtol=0, atol=1e-12) + for name, values in (("body", body), ("pto", pto), ("wave", wave), + ("applied", applied)): + np.testing.assert_allclose(values[:, 0], TIMES, rtol=0, atol=1e-9, + err_msg=f"{name} changed its output grid") + np.testing.assert_allclose(mass[0], mass[2] / mass[1], rtol=0, atol=1e-13) + np.testing.assert_allclose(mass[4], mass[3] + mass[2], rtol=0, atol=1e-7) + assert abs(mass[5]) < 1e-7 and abs(mass[6] - .05) < 1e-12 + inactive_speed = body[:, [7, 8, 10, 11, 12]] + inactive_acceleration = body[:, [49, 50, 52, 53, 54]] + assert np.max(np.abs(inactive_speed)) < 1e-9 + assert np.max(np.abs(inactive_acceleration)) < 1e-9 + assert np.max(np.abs(inactive_acceleration @ mass[7:])) < 1e-6 + assert np.max(np.abs(applied[:, 3])) < 1e-6 + np.testing.assert_allclose(pto[:, 3], body[:, 3] + 2, rtol=0, atol=1e-9) + np.testing.assert_allclose(pto[:, 9], body[:, 9], rtol=0, atol=1e-9) + np.testing.assert_allclose(pto[:, 15], -1_200_000 * pto[:, 9], + rtol=0, atol=1e-6) + np.testing.assert_allclose(pto[:, 21], pto[:, 15] * pto[:, 9], + rtol=0, atol=1e-6) + + packed = { + "case": np.array(case), "time": body[:, 0], + "heave": body[:, 3], "speed": body[:, 9], + "total": body[:, 15], "excitation": body[:, 21], + "radiation": body[:, 27], "reported_added_mass": body[:, 33], + "restoring": body[:, 39], "viscous": body[:, 45], + "acceleration": body[:, 51], + "applied_heave_added": applied[:, 3], + "inactive_speed": inactive_speed, + "inactive_acceleration": inactive_acceleration, + "pto_stroke": pto[:, 3], "pto_speed": pto[:, 9], + "pto_force": pto[:, 15], "pto_power": -pto[:, 21], + "wave": wave[:, 1], "settings": settings, "mass": mass, + } + output.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(output, **packed) + print(f"Packed {case} full-duration mass-control traces to {output}") + if expected is not None: + with np.load(expected, allow_pickle=False) as archived: + assert set(archived.files) == set(packed) + for name, value in packed.items(): + if name == "case": + assert str(archived[name]) == case + else: + np.testing.assert_allclose( + value, archived[name], rtol=1e-9, atol=1e-9, + err_msg=f"{case} {name} differs from {expected}", + ) + print(f"Matched archived mass-control traces in {expected}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("case", choices=("Regular", "RegularCIC")) + parser.add_argument("csv_dir", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--expect", type=Path) + parser.add_argument("--source-hydro", type=Path) + args = parser.parse_args() + pack(args.case, args.csv_dir, args.output, args.expect, args.source_hydro) diff --git a/tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py b/tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py new file mode 100644 index 00000000..8bcbd71c --- /dev/null +++ b/tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py @@ -0,0 +1,124 @@ +"""Pack post-screen 150 s no-active-delay MaxStep controls by wave case. + +Usage: + python tools/pack_ellipsoid_ode45_continuous_nonlinear_mass_control_step.py \ + CASE ARTIFACT_ROOT OUTPUT.npz [--expect ARCHIVED.npz] + +ARTIFACT_ROOT contains CASE_0p0125 and CASE_0p00625 directories downloaded +from the four named artifacts in GitHub Actions run 40131806849. +""" + +import argparse +from pathlib import Path + +import numpy as np + +from pack_ellipsoid_ode45_continuous_nonlinear_mass_control import _check_hydro + + +ROOT = Path(__file__).resolve().parents[1] +BASE = (ROOT / "tests/fixtures/matlab" + / "ellipsoid_ode45_continuous_nonlinear_mass_control") +STEPS = (("0p0125", .0125), ("0p00625", .00625)) +ROWS = 3001 +TIME = np.arange(ROWS) * .05 + + +def _read(folder: Path, kind: str, shape: tuple[int, int]): + values = np.loadtxt(folder / f"{kind}.csv", delimiter=",", ndmin=2) + assert values.shape == shape and np.isfinite(values).all(), (folder, kind) + return values + + +def pack(case: str, artifact_root: Path, output: Path, + expected: Path | None = None) -> None: + if case not in ("Regular", "RegularCIC"): + raise ValueError(f"Unsupported ellipsoid wave case: {case}") + basename = "regular_cic" if case == "RegularCIC" else "regular" + with np.load(BASE / f"{basename}.npz", allow_pickle=False) as base: + base_mass = base["mass"] + base_wave = base["wave"] + np.testing.assert_allclose(base["time"], TIME, rtol=0, atol=1e-9) + np.testing.assert_allclose(base["settings"], + [.05, .05, .05, .05, .05, 150, 0, 0, 0], + rtol=0, atol=1e-12) + packed: dict[str, np.ndarray] = {"case": np.array(case), "time": TIME} + for tag, step in STEPS: + artifact = artifact_root / f"{case}_{tag}" + _check_hydro(artifact / "applications/Nonlinear_Hydro" + / "hydroData/ellipsoid.h5") + folder = (artifact + / "matlab-ellipsoid-ode45-continuous-nonlinear-mass-control-step" + / case / f"maxstep_{tag}") + body = _read(folder, "body", (ROWS, 55)) + pto = _read(folder, "pto", (ROWS, 25)) + wave = _read(folder, "wave", (ROWS, 2)) + applied = _read(folder, "applied_added", (ROWS, 7)) + settings = _read(folder, "settings", (1, 9))[0] + mass = _read(folder, "mass", (1, 12))[0] + np.testing.assert_allclose(settings, + [.05, .05, .05, .05, step, 150, 0, 0, 0], + rtol=0, atol=1e-12) + for name, values in (("body", body), ("pto", pto), ("wave", wave), + ("applied", applied)): + np.testing.assert_allclose(values[:, 0], TIME, rtol=0, atol=1e-9, + err_msg=f"{case} {tag} {name} time") + np.testing.assert_allclose(wave[:, 1], base_wave, rtol=0, atol=1e-10) + np.testing.assert_allclose(np.delete(mass, 6), + np.delete(base_mass, 6), rtol=0, atol=1e-7) + assert abs(mass[6] - step) < 1e-12 + inactive_speed = body[:, [7, 8, 10, 11, 12]] + inactive_acceleration = body[:, [49, 50, 52, 53, 54]] + assert np.max(np.abs(inactive_speed)) < 1e-9 + assert np.max(np.abs(inactive_acceleration)) < 1e-9 + assert np.max(np.abs(inactive_acceleration @ mass[7:])) < 1e-6 + assert np.max(np.abs(applied[:, 3])) < 1e-6 + np.testing.assert_allclose(pto[:, 3], body[:, 3] + 2, + rtol=0, atol=1e-9) + np.testing.assert_allclose(pto[:, 9], body[:, 9], rtol=0, atol=1e-9) + np.testing.assert_allclose(pto[:, 15], -1_200_000 * pto[:, 9], + rtol=0, atol=1e-6) + np.testing.assert_allclose(pto[:, 21], pto[:, 15] * pto[:, 9], + rtol=0, atol=1e-6) + prefix = f"maxstep_{tag}_" + packed.update({ + prefix + "heave": body[:, 3], + prefix + "speed": body[:, 9], + prefix + "total": body[:, 15], + prefix + "excitation": body[:, 21], + prefix + "radiation": body[:, 27], + prefix + "reported_added_mass": body[:, 33], + prefix + "restoring": body[:, 39], + prefix + "viscous": body[:, 45], + prefix + "acceleration": body[:, 51], + prefix + "applied_heave_added": applied[:, 3], + prefix + "inactive_speed": inactive_speed, + prefix + "inactive_acceleration": inactive_acceleration, + prefix + "pto_stroke": pto[:, 3], + prefix + "pto_speed": pto[:, 9], + prefix + "pto_force": pto[:, 15], + prefix + "pto_power": -pto[:, 21], + prefix + "wave": wave[:, 1], + prefix + "settings": settings, + prefix + "mass": mass, + }) + output.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(output, **packed) + print(f"Packed {case} no-active-delay MaxStep traces to {output}") + if expected is not None: + with np.load(expected, allow_pickle=False) as archived: + assert set(archived.files) == set(packed), expected + for name, value in packed.items(): + np.testing.assert_equal(value, archived[name], + err_msg=f"{case} {name} changed") + print(f"Matched archived no-active-delay traces in {expected}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("case", choices=("Regular", "RegularCIC")) + parser.add_argument("artifact_root", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--expect", type=Path) + args = parser.parse_args() + pack(args.case, args.artifact_root, args.output, args.expect) diff --git a/tools/pack_ellipsoid_ode45_restoring_step.py b/tools/pack_ellipsoid_ode45_restoring_step.py new file mode 100644 index 00000000..d854dffe --- /dev/null +++ b/tools/pack_ellipsoid_ode45_restoring_step.py @@ -0,0 +1,114 @@ +"""Pack copied-model ellipsoid MaxStep artifacts into a portable fixture. + +Usage: + python tools/pack_ellipsoid_ode45_restoring_step.py CASE ARTIFACT_ROOT OUT.npz + [--expect ARCHIVED.npz] + +ARTIFACT_ROOT contains the eight named artifacts from run 40131187182. +The archived 0.05 s continuous-nonlinear arm supplies the common base trace. +""" + +import argparse +from pathlib import Path + +import numpy as np + +from pack_ellipsoid_ode45_restoring_timing import _check_source_hydro + + +ROOT = Path(__file__).resolve().parents[1] +BASE = ROOT / "tests/fixtures/matlab/ellipsoid_ode45_restoring_timing" +STEPS = (("0.025", "0p025", .025), + ("0.0125", "0p0125", .0125), + ("0.00625", "0p00625", .00625), + ("0.003125", "0p003125", .003125)) +ROWS = 3001 +TIME = np.arange(ROWS) * .05 + + +def _read(folder: Path, kind: str, shape: tuple[int, int]) -> np.ndarray: + values = np.loadtxt(folder / f"continuous_nonlinear_{kind}.csv", + delimiter=",", ndmin=2) + assert values.shape == shape and np.isfinite(values).all(), (folder, kind) + return values + + +def pack(case: str, artifact_root: Path, output: Path, + expected: Path | None = None) -> None: + if case not in ("Regular", "RegularCIC"): + raise ValueError(f"Unsupported ellipsoid wave case: {case}") + basename = "regular_cic" if case == "RegularCIC" else "regular" + with np.load(BASE / f"{basename}.npz", allow_pickle=False) as base: + assert str(base["case"]) == case + base_time = base["continuous_nonlinear_time"] + base_wave = base["continuous_nonlinear_wave"] + np.testing.assert_allclose( + base["continuous_nonlinear_settings"], + [.05, .05, .05, .05, .05, 150, 0, 0, 0], rtol=0, atol=1e-12, + ) + np.testing.assert_allclose(base_time, TIME, rtol=0, atol=1e-9) + packed: dict[str, np.ndarray] = { + "case": np.array(case), + "time": base_time, + } + for step_text, tag, step in STEPS: + artifact = (artifact_root / f"ellipsoid-ode45-restoring-{case}" + f"-maxstep-{step_text}") + _check_source_hydro( + artifact / "applications/Nonlinear_Hydro/hydroData/ellipsoid.h5" + ) + folder = (artifact / "matlab-ellipsoid-ode45-restoring-timing-step" + / case / f"maxstep_{tag}") + body = _read(folder, "body", (ROWS, 55)) + pto = _read(folder, "pto", (ROWS, 25)) + wave = _read(folder, "wave", (ROWS, 2)) + settings = _read(folder, "settings", (1, 9))[0] + np.testing.assert_allclose(settings, + [.05, .05, .05, .05, step, 150, 0, 0, 0], + rtol=0, atol=1e-12, + err_msg=f"{case} {step_text} settings") + for name, values in (("body", body), ("pto", pto), ("wave", wave)): + np.testing.assert_allclose(values[:, 0], TIME, rtol=0, atol=1e-9, + err_msg=f"{case} {step_text} {name} time") + np.testing.assert_allclose(wave[:, 1], base_wave, rtol=0, atol=1e-10) + np.testing.assert_allclose(pto[:, 3], body[:, 3] + 2, + rtol=0, atol=1e-10) + np.testing.assert_allclose(pto[:, 9], body[:, 9], rtol=0, atol=1e-9) + np.testing.assert_allclose(pto[:, 15], -1_200_000 * pto[:, 9], + rtol=0, atol=1e-6) + np.testing.assert_allclose(pto[:, 21], pto[:, 15] * pto[:, 9], + rtol=0, atol=1e-6) + prefix = f"maxstep_{tag}_" + packed.update({ + prefix + "heave": body[:, 3], + prefix + "speed": body[:, 9], + prefix + "restoring": body[:, 39], + prefix + "pto_stroke": pto[:, 3], + prefix + "pto_speed": pto[:, 9], + prefix + "pto_force": pto[:, 15], + prefix + "pto_power": -pto[:, 21], + prefix + "wave": wave[:, 1], + prefix + "settings": settings, + }) + output.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(output, **packed) + print(f"Packed {case} MaxStep traces to {output}") + if expected is not None: + with np.load(expected, allow_pickle=False) as archived: + assert set(archived.files) == set(packed), expected + for name, value in packed.items(): + np.testing.assert_equal( + value, archived[name], + err_msg=f"{case} {name} differs from {expected}", + ) + print(f"Matched archived MaxStep traces in {expected}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("case", choices=("Regular", "RegularCIC")) + parser.add_argument("artifact_root", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--expect", type=Path) + args = parser.parse_args() + pack(args.case, args.artifact_root, args.output, args.expect) diff --git a/tools/pack_ellipsoid_ode45_restoring_timing.py b/tools/pack_ellipsoid_ode45_restoring_timing.py new file mode 100644 index 00000000..c49f43f6 --- /dev/null +++ b/tools/pack_ellipsoid_ode45_restoring_timing.py @@ -0,0 +1,135 @@ +"""Pack one full-duration ellipsoid force-timing source artifact for Python CI. + +Usage: + python tools/pack_ellipsoid_ode45_restoring_timing.py CASE CSV_DIR OUTPUT.npz + [--expect ARCHIVED.npz] [--source-hydro FRESH.h5] + +CSV_DIR contains the three 150 s arms uploaded by the repository-owned +ellipsoid-ode45-restoring-timing workflow. The original MATLAB repositories +are never modified by this tool. +""" + +import argparse +from pathlib import Path + +import h5py +import numpy as np + + +MODES = ("baseline", "continuous_restoring", "continuous_nonlinear") +ROWS = 3001 +DT = .05 +EXPECTED_RATES = { + "baseline": (.05, .05, .05), + "continuous_restoring": (0., .05, .05), + "continuous_nonlinear": (0., 0., 0.), +} + + +def _read(folder: Path, mode: str, name: str, shape: tuple[int, int]): + values = np.loadtxt(folder / f"{mode}_{name}.csv", delimiter=",", ndmin=2) + assert values.shape == shape, (mode, name, values.shape) + assert np.isfinite(values).all(), (mode, name) + return values + + +def _datasets(group, prefix=""): + result = {} + for name, item in group.items(): + path = f"{prefix}/{name}" + if isinstance(item, h5py.Dataset): + result[path] = item[()] + else: + result.update(_datasets(item, path)) + return result + + +def _check_source_hydro(path: Path) -> None: + reference = (Path(__file__).resolve().parents[1] + / "tests/fixtures/matlab/ellipsoid_ode45_convergence" + / "hydroData/ellipsoid.h5") + with h5py.File(path) as source, h5py.File(reference) as archived: + source_data, archived_data = _datasets(source), _datasets(archived) + assert len(source_data) == 38 and source_data.keys() == archived_data.keys() + for name in source_data: + np.testing.assert_equal(source_data[name], archived_data[name], + err_msg=f"source HDF5 dataset {name} changed") + print(f"Matched 38 pinned HDF5 datasets in {path}") + + +def pack(case: str, folder: Path, output: Path, + expected: Path | None = None, + source_hydro: Path | None = None) -> None: + if case not in ("Regular", "RegularCIC"): + raise ValueError(f"Unsupported ellipsoid wave case: {case}") + if source_hydro is not None: + _check_source_hydro(source_hydro) + times = np.arange(ROWS) * DT + packed: dict[str, np.ndarray] = {"case": np.array(case)} + waves = [] + for mode in MODES: + body = _read(folder, mode, "body", (ROWS, 55)) + pto = _read(folder, mode, "pto", (ROWS, 25)) + wave = _read(folder, mode, "wave", (ROWS, 2)) + settings = _read(folder, mode, "settings", (1, 9))[0] + np.testing.assert_allclose( + settings, [.05, .05, .05, .05, .05, 150, *EXPECTED_RATES[mode]], + rtol=0, atol=1e-12, err_msg=f"{mode} changed a pinned clock", + ) + for name, source in (("body", body), ("pto", pto), ("wave", wave)): + np.testing.assert_allclose(source[:, 0], times, rtol=0, + atol=1e-9, err_msg=f"{mode} {name} time") + np.testing.assert_allclose(pto[:, 9], body[:, 9], rtol=0, atol=1e-9, + err_msg=f"{mode} inactive PTO speed") + np.testing.assert_allclose(pto[:, 15], -1_200_000 * pto[:, 9], + rtol=0, atol=1e-6, + err_msg=f"{mode} PTO damping law") + np.testing.assert_allclose(pto[:, 21], pto[:, 15] * pto[:, 9], + rtol=0, atol=1e-6, + err_msg=f"{mode} PTO source power law") + waves.append(wave[:, 1]) + prefix = f"{mode}_" + packed.update({ + prefix + "time": body[:, 0], + prefix + "heave": body[:, 3], + prefix + "speed": body[:, 9], + prefix + "restoring": body[:, 39], + prefix + "excitation": body[:, 21], + prefix + "acceleration": body[:, 51], + prefix + "pto_stroke": pto[:, 3], + prefix + "pto_speed": pto[:, 9], + prefix + "pto_force": pto[:, 15], + prefix + "pto_power": -pto[:, 21], + prefix + "wave": wave[:, 1], + prefix + "settings": settings, + }) + for wave in waves[1:]: + np.testing.assert_allclose(wave, waves[0], rtol=0, atol=1e-10, + err_msg="force timing changed the incident wave") + output.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(output, **packed) + print(f"Packed {case} source timing traces to {output}") + if expected is not None: + with np.load(expected, allow_pickle=False) as archived: + assert set(archived.files) == set(packed), expected + for name, value in packed.items(): + if name == "case": + assert str(archived[name]) == case + else: + np.testing.assert_allclose( + value, archived[name], rtol=1e-9, atol=1e-9, + err_msg=f"{case} {name} differs from {expected}", + ) + print(f"Matched archived timing traces in {expected}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("case", choices=("Regular", "RegularCIC")) + parser.add_argument("csv_dir", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--expect", type=Path) + parser.add_argument("--source-hydro", type=Path) + args = parser.parse_args() + pack(args.case, args.csv_dir, args.output, args.expect, + args.source_hydro)