|
| 1 | +"""Pair WEC-Sim's pressure-integration spectral mean-drift option.""" |
| 2 | + |
| 3 | +from pathlib import Path |
| 4 | + |
| 5 | +import h5py |
| 6 | +import numpy as np |
| 7 | +from scipy.interpolate import CubicSpline |
| 8 | + |
| 9 | +from wecsim import PMWave, WEC |
| 10 | +from wecsim.environment.irregularWave import pm_equal_energy_components |
| 11 | + |
| 12 | + |
| 13 | +FIXTURE = Path(__file__).parent / "fixtures/matlab/sphere_pressure_drift" |
| 14 | +HYDRO = Path(__file__).parent / "fixtures/matlab/sphere_mean_drift/sphere.h5" |
| 15 | + |
| 16 | + |
| 17 | +def test_pressure_integration_drift_against_pinned_matlab(): |
| 18 | + with np.load(FIXTURE / "source.npz") as source: |
| 19 | + body, wave, recorded, mass = (source[key] for key in |
| 20 | + ("body", "wave", "components", "mass")) |
| 21 | + assert (body.shape, wave.shape, recorded.shape, mass.shape) == ( |
| 22 | + (3001, 19), (3001, 2), (500, 4), (4,), |
| 23 | + ) |
| 24 | + assert all(np.isfinite(values).all() for values in |
| 25 | + (body, wave, recorded, mass)) |
| 26 | + |
| 27 | + components = pm_equal_energy_components( |
| 28 | + HYDRO, significant_height=.5, peak_period=8, |
| 29 | + directions=[0], spreading=[1], seed=1, phase_generator="matlab", |
| 30 | + ) |
| 31 | + for actual, expected in ( |
| 32 | + (components.omega, recorded[:, 0]), |
| 33 | + (components.spectral_amplitude, recorded[:, 1]), |
| 34 | + (components.d_omega, recorded[:, 2]), |
| 35 | + (components.phase[:, 0], recorded[:, 3]), |
| 36 | + ): |
| 37 | + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-12) |
| 38 | + |
| 39 | + # Compute the selected mean-drift term directly from the HDF5 rather |
| 40 | + # than asking the production irregular-wave builder for its result. |
| 41 | + with h5py.File(HYDRO) as h5: |
| 42 | + frequency = np.asarray(h5["simulation_parameters/w"]).ravel() |
| 43 | + coefficient = np.asarray(h5[ |
| 44 | + "body1/hydro_coeffs/mean_drift/pressure_integration/val" |
| 45 | + ])[:, 0, :] |
| 46 | + expected_drift = ( |
| 47 | + CubicSpline(frequency, coefficient, axis=1)(components.omega) |
| 48 | + @ (components.spectral_amplitude * components.d_omega) |
| 49 | + ) * 1000 * 9.81 |
| 50 | + assert np.max(abs(expected_drift)) > 1 |
| 51 | + |
| 52 | + wec = WEC("Sphere pressure-integration mean drift") |
| 53 | + sphere = wec.body("sphere", HYDRO, inertia=mass[1:], |
| 54 | + mean_drift="pressure_integration") |
| 55 | + for axis in ("surge", "heave", "pitch"): |
| 56 | + wec.coordinate(axis, sphere.move(axis)) |
| 57 | + result = wec.run( |
| 58 | + PMWave(.5, 8, seed=1, phase_generator="matlab"), |
| 59 | + dt=.01, end_time=30, ramp_time=10, radiation_memory=10, |
| 60 | + ) |
| 61 | + output = dict(result.raw.extra_outputs) |
| 62 | + ramp = np.where(result.time < 10, |
| 63 | + (1 - np.cos(np.pi * result.time / 10)) / 2, 1) |
| 64 | + np.testing.assert_allclose(result.time, body[:, 0], rtol=0, atol=1e-9) |
| 65 | + np.testing.assert_allclose(result.wave_elevation, wave[:, 1], |
| 66 | + rtol=0, atol=1e-10) |
| 67 | + np.testing.assert_allclose(output["body1_mean_drift_force"], |
| 68 | + ramp[:, None] * expected_drift, |
| 69 | + rtol=0, atol=1e-6) |
| 70 | + assert np.max(abs(output["body1_excitation_force"] |
| 71 | + - body[:, 13:19])) < 1e-5 |
| 72 | + for axis, position_limit, speed_limit in ( |
| 73 | + (0, .003, .004), |
| 74 | + (2, .003, .004), |
| 75 | + (4, 1e-4, 5e-5), |
| 76 | + ): |
| 77 | + position = result.bodies["sphere"].position[:, axis] |
| 78 | + speed = result.bodies["sphere"].velocity[:, axis] |
| 79 | + assert np.max(abs(position - body[:, 1 + axis])) < position_limit |
| 80 | + assert np.max(abs(speed - body[:, 7 + axis])) < speed_limit |
0 commit comments