From 1fdf365498066299972c4e0e93f47546f57a4a12 Mon Sep 17 00:00:00 2001 From: Carlos Aguero Date: Thu, 18 Jun 2026 22:26:01 +0200 Subject: [PATCH 1/3] gz_waves_provider_gerstner: Gerstner wave engine, source system, and GUI registrar Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Carlos Aguero --- gz_waves_provider_gerstner/CMakeLists.txt | 119 ++++++ .../gz/sim/waves/GerstnerWaveSimulation.hh | 136 +++++++ gz_waves_provider_gerstner/package.xml | 29 ++ .../src/GerstnerWaveSimulation.cc | 353 ++++++++++++++++++ .../src/GerstnerWavesGui.cc | 53 +++ .../src/GerstnerWavesSystem.cc | 102 +++++ .../test/gerstner_test.cc | 292 +++++++++++++++ 7 files changed, 1084 insertions(+) create mode 100644 gz_waves_provider_gerstner/CMakeLists.txt create mode 100644 gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh create mode 100644 gz_waves_provider_gerstner/package.xml create mode 100644 gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc create mode 100644 gz_waves_provider_gerstner/src/GerstnerWavesGui.cc create mode 100644 gz_waves_provider_gerstner/src/GerstnerWavesSystem.cc create mode 100644 gz_waves_provider_gerstner/test/gerstner_test.cc diff --git a/gz_waves_provider_gerstner/CMakeLists.txt b/gz_waves_provider_gerstner/CMakeLists.txt new file mode 100644 index 00000000..85a40d81 --- /dev/null +++ b/gz_waves_provider_gerstner/CMakeLists.txt @@ -0,0 +1,119 @@ +cmake_minimum_required(VERSION 3.16) +project(gz_waves_provider_gerstner) + +if(NOT CMAKE_CXX_STANDARD) + set(CMAKE_CXX_STANDARD 17) + set(CMAKE_CXX_STANDARD_REQUIRED ON) + set(CMAKE_CXX_EXTENSIONS OFF) +endif() + +if(CMAKE_COMPILER_IS_GNUCXX OR CMAKE_CXX_COMPILER_ID MATCHES "Clang") + add_compile_options(-Wall -Wextra -Wpedantic) +endif() + +find_package(ament_cmake REQUIRED) +# gz_waves's exported library links gz-sim PUBLIC (the Wavefield component), so +# gz-sim::core must be resolvable before find_package(gz_waves). The core also +# provides WavesSystemBase, the base for the gerstner system plugin below. +find_package(gz_sim_vendor REQUIRED) +find_package(gz-sim REQUIRED) +find_package(gz_waves REQUIRED) # IWaveField + WaveParameters + base +find_package(gz_plugin_vendor REQUIRED) +find_package(gz-plugin REQUIRED COMPONENTS register) +find_package(Eigen3 REQUIRED) + +# -------------------------------------------------------------------------- +# Gerstner wave-field engine (libgz-waves-provider-gerstner): a plain analytic +# IWaveField implementation. It is LINKED directly by the gerstner system +# plugin (server) and the water visual (GUI), not discovered by a runtime +# loader, so it carries no GZ_ADD_PLUGIN and is exported as a normal target. +# MakeGerstnerWaveField() is the factory consumers register under "gerstner". +# -------------------------------------------------------------------------- +add_library(gz-waves-provider-gerstner SHARED + src/GerstnerWaveSimulation.cc +) +target_include_directories(gz-waves-provider-gerstner PUBLIC + $ + $) +target_include_directories(gz-waves-provider-gerstner PRIVATE + ${gz_waves_INCLUDE_DIRS}) +target_link_libraries(gz-waves-provider-gerstner + PUBLIC Eigen3::Eigen) + +# -------------------------------------------------------------------------- +# Gerstner wave source system (gz-sim-waves-gerstner-system): a WavesSystemBase +# subclass that builds the Gerstner engine. This is the gz-plugin SDF loads by +# filename; it links the engine directly, so no runtime engine loader is used. +# -------------------------------------------------------------------------- +add_library(gz-sim-waves-gerstner-system SHARED + src/GerstnerWavesSystem.cc +) +target_include_directories(gz-sim-waves-gerstner-system PRIVATE + $ + ${gz_waves_INCLUDE_DIRS}) +target_link_libraries(gz-sim-waves-gerstner-system + PRIVATE + gz_waves::gz_waves + gz-waves-provider-gerstner + gz-sim::core + gz-plugin::register) + +# -------------------------------------------------------------------------- +# Gerstner GUI registrar (gz-sim-waves-gerstner-gui): loaded in the GUI process +# alongside WaterVisual. It registers the "gerstner" factory so the visual can +# rebuild its own thread-private engine from the replicated recipe via +# CreateWaveSimulation. This is what lets gz_waves_rendering stay +# provider-agnostic (core-only): the engine dependency lives here, the GUI-side +# mirror of gz-sim-waves-gerstner-system. +# -------------------------------------------------------------------------- +add_library(gz-sim-waves-gerstner-gui SHARED + src/GerstnerWavesGui.cc +) +target_include_directories(gz-sim-waves-gerstner-gui PRIVATE + $ + ${gz_waves_INCLUDE_DIRS}) +target_link_libraries(gz-sim-waves-gerstner-gui + PRIVATE + gz_waves::gz_waves + gz-waves-provider-gerstner + gz-sim::core + gz-plugin::register) + +install( + TARGETS gz-waves-provider-gerstner + EXPORT export_gz_waves_provider_gerstner + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib + RUNTIME DESTINATION bin +) +install( + TARGETS gz-sim-waves-gerstner-system gz-sim-waves-gerstner-gui + LIBRARY DESTINATION lib + RUNTIME DESTINATION bin +) +install(DIRECTORY include/ DESTINATION include) + +ament_export_targets(export_gz_waves_provider_gerstner HAS_LIBRARY_TARGET) +ament_export_include_directories(include) +ament_export_dependencies(Eigen3) + +if(BUILD_TESTING) + find_package(ament_lint_auto REQUIRED) + set(ament_cmake_copyright_FOUND TRUE) + # Gazebo (Allman) style, matching the upstream gz-sim source. + set(ament_cmake_cpplint_FOUND TRUE) + set(ament_cmake_uncrustify_FOUND TRUE) + ament_lint_auto_find_test_dependencies() + + # Links the engine directly and registers its "gerstner" factory in-binary + # (see kGerstnerRegistered in the test), so the CreateWaveSimulation tests + # resolve via the registry — no dlopen, no GzPluginHook. + find_package(ament_cmake_gtest REQUIRED) + ament_add_gtest(gerstner_test test/gerstner_test.cc) + target_include_directories(gerstner_test PRIVATE + ${gz_waves_INCLUDE_DIRS}) + target_link_libraries(gerstner_test + gz_waves::gz_waves gz-waves-provider-gerstner) +endif() + +ament_package() diff --git a/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh b/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh new file mode 100644 index 00000000..f0485d55 --- /dev/null +++ b/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh @@ -0,0 +1,136 @@ +/* + * Copyright (C) 2026 Honu Robotics + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + */ + +#ifndef GZ_SIM_WAVES_GERSTNERWAVESIMULATION_HH_ +#define GZ_SIM_WAVES_GERSTNERWAVESIMULATION_HH_ + +#include +#include + +#include + +#include "gz/sim/waves/WaveSimulation.hh" +#include "gz/sim/waves/Wavefield.hh" + +namespace gz::sim::waves +{ + +/// \brief Analytic sum-of-Gerstners wave model. Closed-form for each of up +/// to N component waves; deterministic; unbounded in space. +class GerstnerWaveSimulation final : public IWaveField +{ + /// \brief Default-construct an unconfigured field. Call `SetParameters` + /// before sampling. Used by the engine factory (`MakeGerstnerWaveField`). + public: GerstnerWaveSimulation(); + + /// \brief Construct and sample the spectrum from `_params` (convenience; + /// equivalent to default-construct followed by `SetParameters`). + /// \param[in] _params Wave parameters to build the field from. + public: explicit GerstnerWaveSimulation(const WaveParameters &_params); + + /// \brief Destructor. + public: ~GerstnerWaveSimulation() override = default; + + // Documentation inherited + public: void SetParameters(const WaveParameters &_params) override; + // Documentation inherited + public: double Elevation(double _x, double _y, double _t) const override; + // Documentation inherited + public: Eigen::Vector3d ParticleVelocity( + double _x, double _y, double _t) const override; + // Documentation inherited + public: Eigen::Vector3d Normal( + double _x, double _y, double _t) const override; + // Documentation inherited + public: double Jacobian(double _x, double _y, double _t) const override; + // Documentation inherited + public: void Update(double _simTime) override; + // Documentation inherited + public: const WaveField2D *Field() const override; + // Documentation inherited + public: std::string_view Kind() const override { return "gerstner"; } + + // ---- Backend-specific introspection accessors (exercised by the unit + // tests; not part of the IWaveField interface — WaterVisual consumes + // Field() instead) ---- + + /// \brief Per-component wave amplitudes [m]. + public: const std::vector &Amplitudes() const + { + return this->amplitudes_; + } + /// \brief Per-component wavenumbers |k| [rad/m]. + public: const std::vector &Wavenumbers() const + { + return this->wavenumbers_; + } + /// \brief Per-component angular frequencies ω [rad/s]. + public: const std::vector &AngularFrequencies() const + { + return this->angularFrequencies_; + } + /// \brief Per-component Gerstner steepness in [0, 1]. + public: const std::vector &Steepnesses() const + { + return this->steepnesses_; + } + /// \brief Per-component unit propagation directions. + public: const std::vector &Directions() const + { + return this->directions_; + } + + /// \brief Per-component wave amplitudes [m]. + private: std::vector amplitudes_; + /// \brief Per-component wavenumbers |k| [rad/m]. + private: std::vector wavenumbers_; + /// \brief Per-component angular frequencies ω [rad/s]. + private: std::vector angularFrequencies_; + /// \brief Per-component Gerstner steepness in [0, 1]. + private: std::vector steepnesses_; + /// \brief Per-component unit propagation directions. + private: std::vector directions_; + /// \brief Startup-ramp time constant τ [s]. + private: double tau_{2.0}; + /// \brief Common phase offset φ [rad]. + private: double phase_{0.0}; + + // Render grid: the analytic field sampled onto an N×N tile each Update, + // exposed via Field() as the backend-agnostic rendering contract. Buffers + // are column-major; Update overwrites them in place, so field_'s data() + // pointers (bound in SetParameters) stay valid. + + /// \brief Render-grid resolution per axis (N). + private: std::size_t fieldN_{128}; + /// \brief Render-grid tile extent [m]. + private: double fieldTile_{0.0}; + /// \brief Sim time the render grid was last sampled at. + private: double currentTime_{-1.0}; + /// \brief Column-major vertical-displacement buffer for `field_.dz`. + private: std::vector dzBuf_; + /// \brief Column-major x-chop buffer for `field_.dx`. + private: std::vector dxBuf_; + /// \brief Column-major y-chop buffer for `field_.dy`. + private: std::vector dyBuf_; + /// \brief Column-major folding/foam buffer for `field_.foam`. + private: std::vector foamBuf_; + /// \brief The rendering contract returned by Field(). + private: WaveField2D field_; +}; + +/// \brief Factory: a default-constructed Gerstner wave-field engine (apply +/// `SetParameters` before use). Registered under the "gerstner" token so +/// `CreateWaveSimulation` can rebuild the engine from a serialized `Wavefield` +/// component on the GUI side. +std::shared_ptr MakeGerstnerWaveField(); + +} // namespace gz::sim::waves + +#endif // GZ_SIM_WAVES_GERSTNERWAVESIMULATION_HH_ diff --git a/gz_waves_provider_gerstner/package.xml b/gz_waves_provider_gerstner/package.xml new file mode 100644 index 00000000..2333a198 --- /dev/null +++ b/gz_waves_provider_gerstner/package.xml @@ -0,0 +1,29 @@ + + + + gz_waves_provider_gerstner + 0.0.0 + + Analytic Gerstner (sum-of-sines) wave-field engine for gz_waves. A plain + IWaveField implementation registered in the in-process engine registry and + linked by the gz-sim-waves-gerstner-system plugin — and the simplest + template for writing your own provider package. + + Carlos Agüero + Apache-2.0 + + ament_cmake + + gz_waves + gz_sim_vendor + gz_math_vendor + gz_plugin_vendor + + ament_cmake_gtest + ament_lint_auto + ament_lint_common + + + ament_cmake + + diff --git a/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc new file mode 100644 index 00000000..a60e946a --- /dev/null +++ b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc @@ -0,0 +1,353 @@ +/* + * Copyright (C) 2026 Honu Robotics + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + */ + +#include "gz/sim/waves/GerstnerWaveSimulation.hh" + +#include +#include +#include + + +#include "gz/sim/waves/Wavefield.hh" + +namespace gz::sim::waves +{ + +namespace +{ + +////////////////////////////////////////////////// +double DispersionOmega(double _k, double _g) +{ + return std::sqrt(_g * _k); +} + +////////////////////////////////////////////////// +double DispersionWavenumber(double _omega, double _g) +{ + return _omega * _omega / _g; +} + +////////////////////////////////////////////////// +double PiersonMoskowitz(double _omega, double _omegaP, double _g) +{ + constexpr double kAlpha = 0.0081; + const double ratio = _omegaP / _omega; + return kAlpha * _g * _g / std::pow(_omega, 5.0) * + std::exp(-1.25 * std::pow(ratio, 4.0)); +} + +} // namespace + +////////////////////////////////////////////////// +GerstnerWaveSimulation::GerstnerWaveSimulation() = default; + +////////////////////////////////////////////////// +GerstnerWaveSimulation::GerstnerWaveSimulation(const WaveParameters &_p) +{ + this->SetParameters(_p); +} + +////////////////////////////////////////////////// +void GerstnerWaveSimulation::SetParameters(const WaveParameters &_params) +{ + // Resolve (if set) into period/gain before configuring. + const WaveParameters p = WithSeaState(_params); + this->tau_ = p.tau; + this->phase_ = p.phase; + const double omegaMean = 2.0 * M_PI / p.period; + const std::size_t n = p.number; + + amplitudes_.resize(n); + wavenumbers_.resize(n); + angularFrequencies_.resize(n); + steepnesses_.resize(n); + directions_.resize(n); + + // Pre-compute spectral-bin widths once. The 3-component PMS sampler is the + // standard layout from asv_wave_sim / Tessendorf; if N != 3 we fall back + // to uniform spacing around the mean. + std::vector dOmega(n); + if (p.model == "PMS") + { + if (n == 3) + { + dOmega[0] = omegaMean * (1.0 - 1.0 / p.scale); + dOmega[1] = omegaMean * (p.scale - 1.0 / p.scale) * 0.5; + dOmega[2] = omegaMean * (p.scale - 1.0); + } + else + { + const double w = + omegaMean * (p.scale - 1.0 / p.scale) / static_cast(n); + std::fill(dOmega.begin(), dOmega.end(), w); + } + } + + for (std::size_t i = 0; i < n; ++i) + { + const int nIdx = static_cast(i) - static_cast(n / 2); + const double scaleFactor = std::pow(p.scale, nIdx); + + double a = 0.0; + double k = 0.0; + double omega = 0.0; + + if (p.model == "PMS") + { + omega = omegaMean * scaleFactor; + const double pms = PiersonMoskowitz(omega, omegaMean, p.gravity); + a = p.gain * std::sqrt(2.0 * pms * dOmega[i]); + k = DispersionWavenumber(omega, p.gravity); + } + else if (p.model == "CWR") + { + a = scaleFactor * p.amplitude; + k = DispersionWavenumber(omegaMean, p.gravity) / scaleFactor; + omega = DispersionOmega(k, p.gravity); + } + else + { + std::cerr << "[GerstnerWaveSimulation] unknown spectrum model '" + << p.model << "'; expected 'PMS' or 'CWR'." << '\n'; + amplitudes_.clear(); + wavenumbers_.clear(); + angularFrequencies_.clear(); + steepnesses_.clear(); + directions_.clear(); + return; + } + + double q = 0.0; + if (std::abs(a) > 1e-9) + { + q = std::min(1.0, p.steepness / (a * k * static_cast(n))); + } + + amplitudes_[i] = a; + wavenumbers_[i] = k; + angularFrequencies_[i] = omega; + steepnesses_[i] = q; + + const double theta = nIdx * p.angle + p.direction; + directions_[i] = Eigen::Vector2d(std::cos(theta), std::sin(theta)); + } + + // Render grid: sample the analytic field onto a tile sized to the longest + // component wavelength. Resolution reuses the gridSize knob; the buffers are + // filled by Update(). + this->fieldN_ = p.gridSize > 0 ? p.gridSize : 128; + double kMin = 0.0; + for (const double k : this->wavenumbers_) + if (k > 0.0 && (kMin == 0.0 || k < kMin)) kMin = k; + this->fieldTile_ = (kMin > 0.0) ? (2.0 * M_PI / kMin) : 100.0; + + // Quantize every component's wave vector to the tile lattice so each one + // completes a whole number of cycles over fieldTile_ (not just the longest). + // This makes the analytic field exactly periodic over the tile, so the + // rendered surface tiles with no edge seams. The fundamental is + // k0 = 2*pi/fieldTile_ (== kMin); snapping (kx, ky) to the nearest integer + // multiples of k0 shifts each component by at most half a lattice step. + // (|k| >= kMin guarantees the snapped vector is never zero.) omega and the + // steepness clamp are recomputed from the snapped wavenumber. + if (kMin > 0.0) + { + const double k0 = 2.0 * M_PI / this->fieldTile_; + for (std::size_t i = 0; i < n; ++i) + { + const Eigen::Vector2d kVec = this->wavenumbers_[i] * this->directions_[i]; + const Eigen::Vector2d kSnap(std::round(kVec.x() / k0) * k0, + std::round(kVec.y() / k0) * k0); + const double kMag = kSnap.norm(); + if (kMag <= 0.0) + continue; + this->wavenumbers_[i] = kMag; + this->directions_[i] = kSnap / kMag; + this->angularFrequencies_[i] = DispersionOmega(kMag, p.gravity); + if (std::abs(this->amplitudes_[i]) > 1e-9) + this->steepnesses_[i] = std::min(1.0, + p.steepness / (this->amplitudes_[i] * kMag * static_cast(n))); + } + } + const std::size_t cells = this->fieldN_ * this->fieldN_; + this->dzBuf_.assign(cells, 0.0); + this->dxBuf_.assign(cells, 0.0); + this->dyBuf_.assign(cells, 0.0); + this->foamBuf_.assign(cells, 0.0); + this->field_.n = this->fieldN_; + this->field_.tile = this->fieldTile_; + this->field_.dz = this->dzBuf_.data(); + this->field_.dx = this->dxBuf_.data(); + this->field_.dy = this->dyBuf_.data(); + this->field_.foam = this->foamBuf_.data(); + this->currentTime_ = -1.0; // force the first Update to sample + this->Update(0.0); +} + +////////////////////////////////////////////////// +double GerstnerWaveSimulation::Elevation(double _x, double _y, double _t) const +{ + double eta = 0.0; + const std::size_t n = amplitudes_.size(); + for (std::size_t i = 0; i < n; ++i) + { + const auto &d = directions_[i]; + const double theta = + wavenumbers_[i] * (d.x() * _x + d.y() * _y) - + angularFrequencies_[i] * _t + phase_; + eta += amplitudes_[i] * std::cos(theta); + } + return eta * StartupRamp(_t, this->tau_); +} + +////////////////////////////////////////////////// +Eigen::Vector3d GerstnerWaveSimulation::ParticleVelocity( + double _x, double _y, double _t) const +{ + double vx = 0.0, vy = 0.0, vz = 0.0; + const double r = StartupRamp(_t, this->tau_); + const std::size_t n = amplitudes_.size(); + for (std::size_t i = 0; i < n; ++i) + { + const auto &d = directions_[i]; + const double aw = amplitudes_[i] * angularFrequencies_[i]; + const double theta = + wavenumbers_[i] * (d.x() * _x + d.y() * _y) - + angularFrequencies_[i] * _t + phase_; + const double s = std::sin(theta); + const double c = std::cos(theta); + vx += aw * d.x() * s; + vy += aw * d.y() * s; + vz += aw * c; + } + return Eigen::Vector3d(vx * r, vy * r, vz * r); +} + +////////////////////////////////////////////////// +Eigen::Vector3d GerstnerWaveSimulation::Normal( + double _x, double _y, double _t) const +{ + double dhdx = 0.0, dhdy = 0.0; + const double r = StartupRamp(_t, this->tau_); + const std::size_t n = amplitudes_.size(); + for (std::size_t i = 0; i < n; ++i) + { + const auto &d = directions_[i]; + const double ak = amplitudes_[i] * wavenumbers_[i]; + const double theta = + wavenumbers_[i] * (d.x() * _x + d.y() * _y) - + angularFrequencies_[i] * _t + phase_; + const double s = std::sin(theta); + dhdx += -ak * d.x() * s; + dhdy += -ak * d.y() * s; + } + dhdx *= r; + dhdy *= r; + Eigen::Vector3d nvec(-dhdx, -dhdy, 1.0); + nvec.normalize(); + return nvec; +} + +////////////////////////////////////////////////// +double GerstnerWaveSimulation::Jacobian(double _x, double _y, double _t) const +{ + double dsxdx = 0.0, dsydy = 0.0, dsxdy = 0.0; + const double r = StartupRamp(_t, this->tau_); + const std::size_t n = amplitudes_.size(); + for (std::size_t i = 0; i < n; ++i) + { + const auto &d = directions_[i]; + const double qak = steepnesses_[i] * amplitudes_[i] * wavenumbers_[i]; + const double theta = + wavenumbers_[i] * (d.x() * _x + d.y() * _y) - + angularFrequencies_[i] * _t + phase_; + const double c = std::cos(theta); + dsxdx += -qak * d.x() * d.x() * c; + dsydy += -qak * d.y() * d.y() * c; + dsxdy += -qak * d.x() * d.y() * c; + } + dsxdx *= r; + dsydy *= r; + dsxdy *= r; + return (1.0 + dsxdx) * (1.0 + dsydy) - dsxdy * dsxdy; +} + +////////////////////////////////////////////////// +void GerstnerWaveSimulation::Update(double _simTime) +{ + // Analytic backend: the point queries (Elevation/Normal/...) stay + // closed-form. Update only (re)samples the render grid at the new time, + // matching the "Update advances, Field reads" model the FFT backend uses. + if (_simTime == this->currentTime_) + return; + this->currentTime_ = _simTime; + + const int N = static_cast(this->fieldN_); + if (N <= 0 || this->dzBuf_.empty()) + return; + const double T = this->fieldTile_; + const double r = StartupRamp(_simTime, this->tau_); + const std::size_t nc = this->amplitudes_.size(); + for (int j = 0; j < N; ++j) + { + const double y = static_cast(j) * T / static_cast(N); + for (int i = 0; i < N; ++i) + { + const double x = static_cast(i) * T / static_cast(N); + double dz = 0.0, dx = 0.0, dy = 0.0; + double dsxdx = 0.0, dsydy = 0.0, dsxdy = 0.0; + for (std::size_t c = 0; c < nc; ++c) + { + const auto &d = this->directions_[c]; + const double theta = + this->wavenumbers_[c] * (d.x() * x + d.y() * y) - + this->angularFrequencies_[c] * _simTime + this->phase_; + const double s = std::sin(theta); + const double co = std::cos(theta); + dz += this->amplitudes_[c] * co; + // Gerstner horizontal displacement (chop): -q·a·dir·sin(θ). + const double qa = this->steepnesses_[c] * this->amplitudes_[c]; + dx -= qa * d.x() * s; + dy -= qa * d.y() * s; + // Folding/foam terms reuse co — equivalent to Jacobian(x, y, t) but + // without a second pass over the components. + const double qak = qa * this->wavenumbers_[c]; + dsxdx -= qak * d.x() * d.x() * co; + dsydy -= qak * d.y() * d.y() * co; + dsxdy -= qak * d.x() * d.y() * co; + } + const std::size_t idx = static_cast(i) + + static_cast(j) * + static_cast(N); + this->dzBuf_[idx] = dz * r; + this->dxBuf_[idx] = dx * r; + this->dyBuf_[idx] = dy * r; + const double jxx = dsxdx * r, jyy = dsydy * r, jxy = dsxdy * r; + this->foamBuf_[idx] = (1.0 + jxx) * (1.0 + jyy) - jxy * jxy; + } + } +} + +////////////////////////////////////////////////// +const WaveField2D *GerstnerWaveSimulation::Field() const +{ + return &this->field_; +} + +////////////////////////////////////////////////// +/// \brief Factory used to register the Gerstner engine under the "gerstner" +/// token (see RegisterWaveEngineFactory). Returns a default-constructed engine; +/// the caller applies SetParameters. +std::shared_ptr MakeGerstnerWaveField() +{ + return std::make_shared(); +} + +} // namespace gz::sim::waves diff --git a/gz_waves_provider_gerstner/src/GerstnerWavesGui.cc b/gz_waves_provider_gerstner/src/GerstnerWavesGui.cc new file mode 100644 index 00000000..04ab24ab --- /dev/null +++ b/gz_waves_provider_gerstner/src/GerstnerWavesGui.cc @@ -0,0 +1,53 @@ +/* + * Copyright (C) 2026 Honu Robotics + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + */ + +#include + +#include "gz/sim/System.hh" +#include "gz/sim/waves/GerstnerWaveSimulation.hh" // MakeGerstnerWaveField +#include "gz/sim/waves/WaveSimulation.hh" // RegisterWaveEngineFactory + +namespace gz::sim::systems +{ +/// \brief GUI-side registrar for the Gerstner wave engine. Loaded in the GUI +/// process alongside WaterVisual so the in-process engine registry knows the +/// "gerstner" token, letting WaterVisual rebuild its own (thread-private) engine +/// from the replicated recipe via `CreateWaveSimulation`. The server builds its +/// engine directly (`GerstnerWaves::MakeEngine`) and does not need this. Keeping +/// registration here — rather than in `gz_waves_rendering` — lets the rendering +/// package stay provider-agnostic (it depends only on the gz_waves core). The +/// GUI-side mirror of the server's `gz-sim-waves-gerstner-system`. +/// +/// The factory is registered by the file-scope initializer below, which runs +/// the moment the GUI dlopens this library — so the plugin deliberately carries +/// no ISystem interface and parses no SDF, avoiding gz-sim's empty-plugin SDF +/// re-parse warning. The empty System body exists only so gz-sim has a plugin +/// to load (which is what triggers the dlopen). +class GerstnerWavesGui : public System +{ +}; +} // namespace gz::sim::systems + +namespace +{ +/// \brief Register the "gerstner" engine factory when this library is loaded. +const bool kGerstnerGuiRegistered = [] +{ + gz::sim::waves::RegisterWaveEngineFactory( + "gerstner", &gz::sim::waves::MakeGerstnerWaveField); + return true; +}(); +} // namespace + +////////////////////////////////////////////////// +GZ_ADD_PLUGIN(gz::sim::systems::GerstnerWavesGui, gz::sim::System) + +GZ_ADD_PLUGIN_ALIAS(gz::sim::systems::GerstnerWavesGui, + "gz::sim::systems::GerstnerWavesGui") diff --git a/gz_waves_provider_gerstner/src/GerstnerWavesSystem.cc b/gz_waves_provider_gerstner/src/GerstnerWavesSystem.cc new file mode 100644 index 00000000..291247f3 --- /dev/null +++ b/gz_waves_provider_gerstner/src/GerstnerWavesSystem.cc @@ -0,0 +1,102 @@ +/* + * Copyright (C) 2026 Honu Robotics + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + */ + +#include +#include + +#include + +#include "gz/sim/systems/WavesSystemBase.hh" +#include "gz/sim/waves/GerstnerWaveSimulation.hh" // MakeGerstnerWaveField +#include "gz/sim/waves/Wavefield.hh" // WaveParameters + +namespace gz::sim::systems +{ +/// \brief Wave source system backed by the analytic Gerstner engine (a +/// sum-of-Gerstners trochoidal wave field). Loaded by SDF as +/// `gz-sim-waves-gerstner-system`. The ECM/source plumbing lives in +/// WavesSystemBase; this names the engine and owns the Gerstner parameter +/// surface below. (It ignores the FFT-only tags: tile_size, grid_size, seed, +/// choppiness, spectrum, spreading, dispersion.) +/// +/// ## SDF parameters +/// +/// `` is at plugin level; the rest live in a `` block. All +/// are also runtime-settable via the `/world//wave/set_parameters` +/// service (`gz.msgs.Param`). +/// +/// | Tag | Type | Default | Meaning | +/// |---|---|---|---| +/// | `` | double [Hz] | 30.0 | Throttle for the backend `Update()` (plugin level) | +/// | `` | string | PMS | Sampling model: `PMS` (wind-derived components) or `CWR` (constant amplitude) | +/// | `` | uint | 3 | Number of component waves | +/// | `` | double [s] | 5.0 | Mean wave period | +/// | `` | double [m] | 0.0 | Mean amplitude (CWR model only) | +/// | `` | double [rad] | 0.0 | Mean propagation direction from +X | +/// | `` | double [rad] | 0.4 | Angular spread between components | +/// | `` | double | 1.1 | Amplitude/length ratio between mean and extreme components | +/// | `` | double [0,1] | 0.0 | Crest sharpness (0 = round sine, 1 = pinched) | +/// | `` | double [rad] | 0.0 | Common phase offset | +/// | `` | double [s] | 2.0 | Startup-ramp time constant, (1 - exp(-t/tau)) | +/// | `` | double | 1.0 | Amplitude multiplier | +/// | `` | int [0-9] | -1 (off) | WMO sea-state code; when >= 0 it OVERRIDES ``/`` — see "Sea state" below | +/// +/// ### Sea state vs. manual height +/// +/// `` is the one-knob control: a WMO code 0-9 (table 3700) that sets +/// the sea by deriving its significant wave height and peak period. Precedence: +/// - When `>= 0` it **overrides `` and ``** — setting those +/// alongside it has no effect, so use one approach or the other. +/// - `0` is calm/glassy (flat) water; `-1` (the default, or simply omitting the +/// tag) means "manual": drive the sea with `` + `` instead. +/// - Every other knob (``, ``, ``, ``, +/// ``, ...) is independent of `` and applies in both modes. +/// +/// \verbatim +/// +/// 30 +/// +/// PMS +/// 3 +/// 3.2 +/// 2.356 +/// 0.4 +/// 1.1 +/// 0.5 +/// 1.0 +/// 2.0 +/// +/// +/// \endverbatim +class GerstnerWaves : public WavesSystemBase +{ + // Documentation inherited + protected: std::string EngineToken() const override { return "gerstner"; } + + // Documentation inherited + protected: std::shared_ptr MakeEngine( + const waves::WaveParameters &_params) const override + { + auto engine = waves::MakeGerstnerWaveField(); + engine->SetParameters(_params); + return engine; + } +}; +} // namespace gz::sim::systems + +GZ_ADD_PLUGIN(gz::sim::systems::GerstnerWaves, + gz::sim::System, + gz::sim::systems::GerstnerWaves::ISystemConfigure, + gz::sim::systems::GerstnerWaves::ISystemPreUpdate, + gz::sim::systems::GerstnerWaves::ISystemReset) + +GZ_ADD_PLUGIN_ALIAS(gz::sim::systems::GerstnerWaves, + "gz::sim::systems::GerstnerWaves") diff --git a/gz_waves_provider_gerstner/test/gerstner_test.cc b/gz_waves_provider_gerstner/test/gerstner_test.cc new file mode 100644 index 00000000..7045354c --- /dev/null +++ b/gz_waves_provider_gerstner/test/gerstner_test.cc @@ -0,0 +1,292 @@ +/* + * Copyright (C) 2026 Honu Robotics + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + */ + +#include +#include +#include + +#include + +#include "gz/sim/waves/Eval.hh" +#include "gz/sim/waves/GerstnerWaveSimulation.hh" +#include "gz/sim/waves/WaveSimulation.hh" +#include "gz/sim/waves/Wavefield.hh" + +namespace gsw = gz::sim::waves; + +namespace +{ +// The gerstner engine is no longer a dlopen'd gz-plugin; this binary links it +// directly, so register its factory for the CreateWaveSimulation tests. +const bool kGerstnerRegistered = [] { + gsw::RegisterWaveEngineFactory("gerstner", &gsw::MakeGerstnerWaveField); + return true; +}(); + +/// \brief Build a PMS-model Gerstner WavefieldData with a configured engine. +/// \return The populated WavefieldData (algorithm "gerstner"). +gsw::WavefieldData MakePmsField() +{ + gsw::WavefieldData data; + data.algorithm = "gerstner"; + data.params.model = "PMS"; + data.params.number = 3; + data.params.period = 5.0; + data.params.gain = 0.3; + data.params.direction = 0.0; + data.params.angle = 0.4; + data.params.scale = 1.1; + data.params.steepness = 0.0; + data.params.tau = 2.0; + data.simulation = gsw::CreateWaveSimulation(data.algorithm, data.params); + return data; +} +} // namespace + +////////////////////////////////////////////////// +TEST(Factory, GerstnerInstantiates) +{ + auto wf = MakePmsField(); + EXPECT_NE(wf.simulation, nullptr); + EXPECT_EQ(wf.simulation->Kind(), "gerstner"); +} + +////////////////////////////////////////////////// +TEST(Factory, UnknownAlgorithmReturnsNullptr) +{ + gsw::WaveParameters p; + auto sim = gsw::CreateWaveSimulation("bogus", p); + EXPECT_EQ(sim, nullptr); +} + +// NOTE: the FFT serialization/decoupling test lives in fft_test.cc, which +// links the fft engine and registers its factory. Engines are plain libraries +// now (no GzPluginHook), so a binary links + registers whichever engine it +// exercises; this one covers gerstner. + +////////////////////////////////////////////////// +TEST(Gerstner, AccessorsArePopulated) +{ + auto wf = MakePmsField(); + const auto *g = + dynamic_cast(wf.simulation.get()); + ASSERT_NE(g, nullptr); + EXPECT_EQ(g->Amplitudes().size(), 3u); + EXPECT_EQ(g->Wavenumbers().size(), 3u); + EXPECT_EQ(g->AngularFrequencies().size(), 3u); + EXPECT_EQ(g->Steepnesses().size(), 3u); + EXPECT_EQ(g->Directions().size(), 3u); +} + +////////////////////////////////////////////////// +TEST(Gerstner, DispersionRelationHolds) +{ + auto wf = MakePmsField(); + const auto *g = + dynamic_cast(wf.simulation.get()); + ASSERT_NE(g, nullptr); + for (std::size_t i = 0; i < g->Amplitudes().size(); ++i) + { + const double omega = g->AngularFrequencies()[i]; + const double k = g->Wavenumbers()[i]; + EXPECT_NEAR(omega * omega, gsw::WaveParameters{}.gravity * k, + 1e-6 * omega * omega); + } +} + +////////////////////////////////////////////////// +TEST(SurfaceElevation, ZeroAmplitudeProducesZeroElevation) +{ + gsw::WavefieldData data; + data.algorithm = "gerstner"; + data.params.model = "CWR"; + data.params.number = 3; + data.params.period = 5.0; + data.params.amplitude = 0.0; + data.simulation = gsw::CreateWaveSimulation(data.algorithm, data.params); + + EXPECT_NEAR(gsw::SurfaceElevation(data, 0.0, 0.0, 100.0), 0.0, 1e-12); + EXPECT_NEAR(gsw::SurfaceElevation(data, 10.0, -3.0, 50.0), 0.0, 1e-12); +} + +////////////////////////////////////////////////// +TEST(SurfaceElevation, IsBounded) +{ + auto wf = MakePmsField(); + const auto *g = + dynamic_cast(wf.simulation.get()); + double bound = 0.0; + for (double a : g->Amplitudes()) + bound += std::abs(a); + + for (double t = 5.0; t < 30.0; t += 1.7) + { + for (double x = -50.0; x <= 50.0; x += 7.3) + { + for (double y = -50.0; y <= 50.0; y += 6.1) + { + const double eta = gsw::SurfaceElevation(wf, x, y, t); + EXPECT_LE(std::abs(eta), bound + 1e-9); + } + } + } +} + +////////////////////////////////////////////////// +TEST(SurfaceElevation, RampUpFromZero) +{ + auto wf = MakePmsField(); + // The startup ramp (1 - exp(-t/tau)) multiplies the whole field: exactly zero + // at t = 0, growing with time. Probe it via the spatial RMS rather than a + // single point — the RMS is phase-independent (a single point can sit on a + // node at one instant and a crest at another), so RMS(t) tracks the ramp + // directly. With tau = 2 s, ramp(0.5)/ramp(100) ~ 0.22, well under 0.5. + EXPECT_NEAR(gsw::SurfaceElevation(wf, 1.0, 2.0, 0.0), 0.0, 1e-12); + const auto spatialRms = [&](double t) + { + double sumSq = 0.0; + int count = 0; + for (double x = 0.0; x < 150.0; x += 7.0) + for (double y = 0.0; y < 150.0; y += 7.0) + { + const double e = gsw::SurfaceElevation(wf, x, y, t); + sumSq += e * e; + ++count; + } + return std::sqrt(sumSq / count); + }; + EXPECT_LT(spatialRms(0.5), spatialRms(100.0) * 0.5); +} + +////////////////////////////////////////////////// +TEST(ParticleVelocity, FiniteOnPmsField) +{ + auto wf = MakePmsField(); + for (double t = 5.0; t < 20.0; t += 1.0) + { + const auto v = gsw::ParticleVelocity(wf, 0.0, 0.0, t); + EXPECT_TRUE(std::isfinite(v.x())); + EXPECT_TRUE(std::isfinite(v.y())); + EXPECT_TRUE(std::isfinite(v.z())); + EXPECT_LT(v.norm(), 10.0); + } +} + +////////////////////////////////////////////////// +TEST(Normal, IsUnitVector) +{ + auto wf = MakePmsField(); + for (double t = 5.0; t < 20.0; t += 1.0) + { + const auto n = gsw::Normal(wf, 1.0, 2.0, t); + EXPECT_NEAR(n.norm(), 1.0, 1e-9); + EXPECT_GT(n.z(), 0.0); + } +} + +////////////////////////////////////////////////// +TEST(Jacobian, OneAtZeroSteepness) +{ + auto wf = MakePmsField(); + EXPECT_NEAR(gsw::Jacobian(wf, 0.0, 0.0, 10.0), 1.0, 1e-12); + EXPECT_NEAR(gsw::Jacobian(wf, 5.0, -3.0, 7.0), 1.0, 1e-12); +} + +////////////////////////////////////////////////// +TEST(FoamMask, ZeroAtZeroSteepness) +{ + auto wf = MakePmsField(); + EXPECT_EQ(gsw::FoamMask(wf, 0.0, 0.0, 10.0), 0.0); +} + +// The component serializes only the recipe; operator>> restores the params and +// leaves `simulation` null. A consumer rebuilds its own engine from the recipe. +////////////////////////////////////////////////// +TEST(Serialization, RoundTripPreservesRecipe) +{ + auto wf = MakePmsField(); + std::stringstream ss; + ss << wf; + gsw::WavefieldData rebuilt; + ss >> rebuilt; + + EXPECT_EQ(rebuilt.simulation, nullptr) + << "operator>> must not build a (process-global) engine"; + EXPECT_EQ(rebuilt.algorithm, "gerstner"); + + // Rebuild from the recipe; same params → same analytic field. + rebuilt.simulation = + gsw::CreateWaveSimulation(rebuilt.algorithm, rebuilt.params); + ASSERT_NE(rebuilt.simulation, nullptr); + EXPECT_NEAR( + gsw::SurfaceElevation(wf, 10.0, 5.0, 7.0), + gsw::SurfaceElevation(rebuilt, 10.0, 5.0, 7.0), + 1e-9); +} + +// The analytic backend samples itself onto a grid for the unified rendering +// contract. The sampled values must match the closed-form queries at the same +// world points, so the renderer draws the same field the buoyancy sees. +////////////////////////////////////////////////// +TEST(Gerstner, FieldSamplesAnalyticGrid) +{ + auto wf = MakePmsField(); + auto *g = dynamic_cast(wf.simulation.get()); + ASSERT_NE(g, nullptr); + g->Update(3.0); + const auto *f = g->Field(); + ASSERT_NE(f, nullptr); + EXPECT_GT(f->n, 0u); + EXPECT_GT(f->tile, 0.0); + ASSERT_NE(f->dz, nullptr); + ASSERT_NE(f->dx, nullptr); + ASSERT_NE(f->dy, nullptr); + ASSERT_NE(f->foam, nullptr); + + const int N = static_cast(f->n); + const double T = f->tile; + double maxAbs = 0.0; + for (int j = 0; j < N; j += N / 4) + { + for (int i = 0; i < N; i += N / 4) + { + const double x = static_cast(i) * T / N; + const double y = static_cast(j) * T / N; + EXPECT_NEAR(f->dz[i + j * N], g->Elevation(x, y, 3.0), 1e-9); + EXPECT_NEAR(f->foam[i + j * N], g->Jacobian(x, y, 3.0), 1e-9); + EXPECT_TRUE(std::isfinite(f->dx[i + j * N])); + EXPECT_TRUE(std::isfinite(f->dy[i + j * N])); + const double a = std::abs(f->dz[i + j * N]); + if (a > maxAbs) maxAbs = a; + } + } + EXPECT_GT(maxAbs, 0.0) << "sampled field is flat — Update didn't fill it"; +} + +// scales the Gerstner field: a rougher sea -> bigger waves. +////////////////////////////////////////////////// +TEST(Gerstner, SeaStateScalesWaveHeight) +{ + auto sigma = [](int code) + { + gsw::WaveParameters p; + p.model = "PMS"; + p.number = 3; + p.seaState = code; + gsw::GerstnerWaveSimulation sim; + sim.SetParameters(p); + double s2 = 0.0; + for (double a : sim.Amplitudes()) + s2 += a * a; + return std::sqrt(0.5 * s2); // surface RMS = sqrt(sum a_i^2 / 2) + }; + EXPECT_GT(sigma(2), 0.0); + EXPECT_GT(sigma(6), sigma(2)) << "rougher sea state should produce bigger waves"; +} From 165a67b811a54c00c201bcc79b7f9736e86880e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20Ag=C3=BCero?= Date: Mon, 22 Jun 2026 18:14:14 +0200 Subject: [PATCH 2/3] gz_waves_provider_gerstner: gz::math vectors, dep wiring, and style cleanup MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Drop Eigen in favor of gz::math vectors (matching the core interface), wire gz-math into the engine target, switch members to gz-style (this->member, no trailing underscore), and factor the per-component phase into a shared Phase() helper. Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Carlos Agüero --- gz_waves_provider_gerstner/CMakeLists.txt | 7 +- .../gz/sim/waves/GerstnerWaveSimulation.hh | 70 +++--- .../src/GerstnerWaveSimulation.cc | 220 +++++++++--------- .../test/gerstner_test.cc | 12 +- 4 files changed, 161 insertions(+), 148 deletions(-) diff --git a/gz_waves_provider_gerstner/CMakeLists.txt b/gz_waves_provider_gerstner/CMakeLists.txt index 85a40d81..35d232bc 100644 --- a/gz_waves_provider_gerstner/CMakeLists.txt +++ b/gz_waves_provider_gerstner/CMakeLists.txt @@ -20,7 +20,8 @@ find_package(gz-sim REQUIRED) find_package(gz_waves REQUIRED) # IWaveField + WaveParameters + base find_package(gz_plugin_vendor REQUIRED) find_package(gz-plugin REQUIRED COMPONENTS register) -find_package(Eigen3 REQUIRED) +find_package(gz_math_vendor REQUIRED) +find_package(gz-math REQUIRED) # gz::math vectors + GZ_PI # -------------------------------------------------------------------------- # Gerstner wave-field engine (libgz-waves-provider-gerstner): a plain analytic @@ -38,7 +39,7 @@ target_include_directories(gz-waves-provider-gerstner PUBLIC target_include_directories(gz-waves-provider-gerstner PRIVATE ${gz_waves_INCLUDE_DIRS}) target_link_libraries(gz-waves-provider-gerstner - PUBLIC Eigen3::Eigen) + PUBLIC gz-math::gz-math) # -------------------------------------------------------------------------- # Gerstner wave source system (gz-sim-waves-gerstner-system): a WavesSystemBase @@ -95,7 +96,7 @@ install(DIRECTORY include/ DESTINATION include) ament_export_targets(export_gz_waves_provider_gerstner HAS_LIBRARY_TARGET) ament_export_include_directories(include) -ament_export_dependencies(Eigen3) +ament_export_dependencies(gz_math_vendor) if(BUILD_TESTING) find_package(ament_lint_auto REQUIRED) diff --git a/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh b/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh index f0485d55..6008c1a8 100644 --- a/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh +++ b/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh @@ -14,7 +14,8 @@ #include #include -#include +#include +#include #include "gz/sim/waves/WaveSimulation.hh" #include "gz/sim/waves/Wavefield.hh" @@ -43,10 +44,10 @@ class GerstnerWaveSimulation final : public IWaveField // Documentation inherited public: double Elevation(double _x, double _y, double _t) const override; // Documentation inherited - public: Eigen::Vector3d ParticleVelocity( + public: gz::math::Vector3d ParticleVelocity( double _x, double _y, double _t) const override; // Documentation inherited - public: Eigen::Vector3d Normal( + public: gz::math::Vector3d Normal( double _x, double _y, double _t) const override; // Documentation inherited public: double Jacobian(double _x, double _y, double _t) const override; @@ -64,65 +65,76 @@ class GerstnerWaveSimulation final : public IWaveField /// \brief Per-component wave amplitudes [m]. public: const std::vector &Amplitudes() const { - return this->amplitudes_; + return this->amplitudes; } /// \brief Per-component wavenumbers |k| [rad/m]. public: const std::vector &Wavenumbers() const { - return this->wavenumbers_; + return this->wavenumbers; } /// \brief Per-component angular frequencies ω [rad/s]. public: const std::vector &AngularFrequencies() const { - return this->angularFrequencies_; + return this->angularFrequencies; } /// \brief Per-component Gerstner steepness in [0, 1]. public: const std::vector &Steepnesses() const { - return this->steepnesses_; + return this->steepnesses; } /// \brief Per-component unit propagation directions. - public: const std::vector &Directions() const + public: const std::vector &Directions() const { - return this->directions_; + return this->directions; } + /// \brief Phase θ_i(x, y, t) = k_i·(d_i·(x, y)) − ω_i·t + φ of component `_i` + /// at world point (x, y) and time t. Shared by every per-component sum + /// (Elevation / ParticleVelocity / Normal / Jacobian / the render grid) so + /// the phase is computed one way only. + /// \param[in] _i Component index (< number of components). + /// \param[in] _x World-frame x coordinate [m]. + /// \param[in] _y World-frame y coordinate [m]. + /// \param[in] _t Simulation time [s]. + /// \return The component's phase [rad]. + private: double Phase(std::size_t _i, double _x, double _y, double _t) const; + /// \brief Per-component wave amplitudes [m]. - private: std::vector amplitudes_; + private: std::vector amplitudes; /// \brief Per-component wavenumbers |k| [rad/m]. - private: std::vector wavenumbers_; + private: std::vector wavenumbers; /// \brief Per-component angular frequencies ω [rad/s]. - private: std::vector angularFrequencies_; + private: std::vector angularFrequencies; /// \brief Per-component Gerstner steepness in [0, 1]. - private: std::vector steepnesses_; + private: std::vector steepnesses; /// \brief Per-component unit propagation directions. - private: std::vector directions_; + private: std::vector directions; /// \brief Startup-ramp time constant τ [s]. - private: double tau_{2.0}; + private: double tau{2.0}; /// \brief Common phase offset φ [rad]. - private: double phase_{0.0}; + private: double phase{0.0}; // Render grid: the analytic field sampled onto an N×N tile each Update, // exposed via Field() as the backend-agnostic rendering contract. Buffers - // are column-major; Update overwrites them in place, so field_'s data() + // are column-major; Update overwrites them in place, so field's data() // pointers (bound in SetParameters) stay valid. /// \brief Render-grid resolution per axis (N). - private: std::size_t fieldN_{128}; + private: std::size_t fieldN{128}; /// \brief Render-grid tile extent [m]. - private: double fieldTile_{0.0}; + private: double fieldTile{0.0}; /// \brief Sim time the render grid was last sampled at. - private: double currentTime_{-1.0}; - /// \brief Column-major vertical-displacement buffer for `field_.dz`. - private: std::vector dzBuf_; - /// \brief Column-major x-chop buffer for `field_.dx`. - private: std::vector dxBuf_; - /// \brief Column-major y-chop buffer for `field_.dy`. - private: std::vector dyBuf_; - /// \brief Column-major folding/foam buffer for `field_.foam`. - private: std::vector foamBuf_; + private: double currentTime{-1.0}; + /// \brief Column-major vertical-displacement buffer for `field.dz`. + private: std::vector dzBuf; + /// \brief Column-major x-chop buffer for `field.dx`. + private: std::vector dxBuf; + /// \brief Column-major y-chop buffer for `field.dy`. + private: std::vector dyBuf; + /// \brief Column-major folding/foam buffer for `field.foam`. + private: std::vector foamBuf; /// \brief The rendering contract returned by Field(). - private: WaveField2D field_; + private: WaveField2D field; }; /// \brief Factory: a default-constructed Gerstner wave-field engine (apply diff --git a/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc index a60e946a..683576a1 100644 --- a/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc +++ b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc @@ -14,6 +14,7 @@ #include #include +#include // GZ_PI #include "gz/sim/waves/Wavefield.hh" @@ -60,16 +61,16 @@ void GerstnerWaveSimulation::SetParameters(const WaveParameters &_params) { // Resolve (if set) into period/gain before configuring. const WaveParameters p = WithSeaState(_params); - this->tau_ = p.tau; - this->phase_ = p.phase; - const double omegaMean = 2.0 * M_PI / p.period; + this->tau = p.tau; + this->phase = p.phase; + const double omegaMean = 2.0 * GZ_PI / p.period; const std::size_t n = p.number; - amplitudes_.resize(n); - wavenumbers_.resize(n); - angularFrequencies_.resize(n); - steepnesses_.resize(n); - directions_.resize(n); + this->amplitudes.resize(n); + this->wavenumbers.resize(n); + this->angularFrequencies.resize(n); + this->steepnesses.resize(n); + this->directions.resize(n); // Pre-compute spectral-bin widths once. The 3-component PMS sampler is the // standard layout from asv_wave_sim / Tessendorf; if N != 3 we fall back @@ -117,11 +118,11 @@ void GerstnerWaveSimulation::SetParameters(const WaveParameters &_params) { std::cerr << "[GerstnerWaveSimulation] unknown spectrum model '" << p.model << "'; expected 'PMS' or 'CWR'." << '\n'; - amplitudes_.clear(); - wavenumbers_.clear(); - angularFrequencies_.clear(); - steepnesses_.clear(); - directions_.clear(); + this->amplitudes.clear(); + this->wavenumbers.clear(); + this->angularFrequencies.clear(); + this->steepnesses.clear(); + this->directions.clear(); return; } @@ -131,127 +132,129 @@ void GerstnerWaveSimulation::SetParameters(const WaveParameters &_params) q = std::min(1.0, p.steepness / (a * k * static_cast(n))); } - amplitudes_[i] = a; - wavenumbers_[i] = k; - angularFrequencies_[i] = omega; - steepnesses_[i] = q; + this->amplitudes[i] = a; + this->wavenumbers[i] = k; + this->angularFrequencies[i] = omega; + this->steepnesses[i] = q; const double theta = nIdx * p.angle + p.direction; - directions_[i] = Eigen::Vector2d(std::cos(theta), std::sin(theta)); + this->directions[i] = gz::math::Vector2d(std::cos(theta), std::sin(theta)); } // Render grid: sample the analytic field onto a tile sized to the longest // component wavelength. Resolution reuses the gridSize knob; the buffers are // filled by Update(). - this->fieldN_ = p.gridSize > 0 ? p.gridSize : 128; + this->fieldN = p.gridSize > 0 ? p.gridSize : 128; double kMin = 0.0; - for (const double k : this->wavenumbers_) + for (const double k : this->wavenumbers) if (k > 0.0 && (kMin == 0.0 || k < kMin)) kMin = k; - this->fieldTile_ = (kMin > 0.0) ? (2.0 * M_PI / kMin) : 100.0; + this->fieldTile = (kMin > 0.0) ? (2.0 * GZ_PI / kMin) : 100.0; // Quantize every component's wave vector to the tile lattice so each one - // completes a whole number of cycles over fieldTile_ (not just the longest). + // completes a whole number of cycles over this->fieldTile (not just the longest). // This makes the analytic field exactly periodic over the tile, so the // rendered surface tiles with no edge seams. The fundamental is - // k0 = 2*pi/fieldTile_ (== kMin); snapping (kx, ky) to the nearest integer + // k0 = 2*pi/this->fieldTile (== kMin); snapping (kx, ky) to the nearest integer // multiples of k0 shifts each component by at most half a lattice step. // (|k| >= kMin guarantees the snapped vector is never zero.) omega and the // steepness clamp are recomputed from the snapped wavenumber. if (kMin > 0.0) { - const double k0 = 2.0 * M_PI / this->fieldTile_; + const double k0 = 2.0 * GZ_PI / this->fieldTile; for (std::size_t i = 0; i < n; ++i) { - const Eigen::Vector2d kVec = this->wavenumbers_[i] * this->directions_[i]; - const Eigen::Vector2d kSnap(std::round(kVec.x() / k0) * k0, - std::round(kVec.y() / k0) * k0); - const double kMag = kSnap.norm(); + const gz::math::Vector2d kVec = + this->wavenumbers[i] * this->directions[i]; + const gz::math::Vector2d kSnap(std::round(kVec.X() / k0) * k0, + std::round(kVec.Y() / k0) * k0); + const double kMag = kSnap.Length(); if (kMag <= 0.0) continue; - this->wavenumbers_[i] = kMag; - this->directions_[i] = kSnap / kMag; - this->angularFrequencies_[i] = DispersionOmega(kMag, p.gravity); - if (std::abs(this->amplitudes_[i]) > 1e-9) - this->steepnesses_[i] = std::min(1.0, - p.steepness / (this->amplitudes_[i] * kMag * static_cast(n))); + this->wavenumbers[i] = kMag; + this->directions[i] = kSnap / kMag; + this->angularFrequencies[i] = DispersionOmega(kMag, p.gravity); + if (std::abs(this->amplitudes[i]) > 1e-9) + this->steepnesses[i] = std::min(1.0, + p.steepness / (this->amplitudes[i] * kMag * static_cast(n))); } } - const std::size_t cells = this->fieldN_ * this->fieldN_; - this->dzBuf_.assign(cells, 0.0); - this->dxBuf_.assign(cells, 0.0); - this->dyBuf_.assign(cells, 0.0); - this->foamBuf_.assign(cells, 0.0); - this->field_.n = this->fieldN_; - this->field_.tile = this->fieldTile_; - this->field_.dz = this->dzBuf_.data(); - this->field_.dx = this->dxBuf_.data(); - this->field_.dy = this->dyBuf_.data(); - this->field_.foam = this->foamBuf_.data(); - this->currentTime_ = -1.0; // force the first Update to sample + const std::size_t cells = this->fieldN * this->fieldN; + this->dzBuf.assign(cells, 0.0); + this->dxBuf.assign(cells, 0.0); + this->dyBuf.assign(cells, 0.0); + this->foamBuf.assign(cells, 0.0); + this->field.n = this->fieldN; + this->field.tile = this->fieldTile; + this->field.dz = this->dzBuf.data(); + this->field.dx = this->dxBuf.data(); + this->field.dy = this->dyBuf.data(); + this->field.foam = this->foamBuf.data(); + this->currentTime = -1.0; // force the first Update to sample this->Update(0.0); } +////////////////////////////////////////////////// +double GerstnerWaveSimulation::Phase( + std::size_t _i, double _x, double _y, double _t) const +{ + const auto &d = this->directions[_i]; + return this->wavenumbers[_i] * (d.X() * _x + d.Y() * _y) - + this->angularFrequencies[_i] * _t + this->phase; +} + ////////////////////////////////////////////////// double GerstnerWaveSimulation::Elevation(double _x, double _y, double _t) const { double eta = 0.0; - const std::size_t n = amplitudes_.size(); + const std::size_t n = this->amplitudes.size(); for (std::size_t i = 0; i < n; ++i) { - const auto &d = directions_[i]; - const double theta = - wavenumbers_[i] * (d.x() * _x + d.y() * _y) - - angularFrequencies_[i] * _t + phase_; - eta += amplitudes_[i] * std::cos(theta); + eta += this->amplitudes[i] * std::cos(this->Phase(i, _x, _y, _t)); } - return eta * StartupRamp(_t, this->tau_); + return eta * StartupRamp(_t, this->tau); } ////////////////////////////////////////////////// -Eigen::Vector3d GerstnerWaveSimulation::ParticleVelocity( +gz::math::Vector3d GerstnerWaveSimulation::ParticleVelocity( double _x, double _y, double _t) const { double vx = 0.0, vy = 0.0, vz = 0.0; - const double r = StartupRamp(_t, this->tau_); - const std::size_t n = amplitudes_.size(); + const double r = StartupRamp(_t, this->tau); + const std::size_t n = this->amplitudes.size(); for (std::size_t i = 0; i < n; ++i) { - const auto &d = directions_[i]; - const double aw = amplitudes_[i] * angularFrequencies_[i]; - const double theta = - wavenumbers_[i] * (d.x() * _x + d.y() * _y) - - angularFrequencies_[i] * _t + phase_; + const auto &d = this->directions[i]; + const double aw = this->amplitudes[i] * this->angularFrequencies[i]; + const double theta = this->Phase(i, _x, _y, _t); const double s = std::sin(theta); const double c = std::cos(theta); - vx += aw * d.x() * s; - vy += aw * d.y() * s; + vx += aw * d.X() * s; + vy += aw * d.Y() * s; vz += aw * c; } - return Eigen::Vector3d(vx * r, vy * r, vz * r); + return gz::math::Vector3d(vx * r, vy * r, vz * r); } ////////////////////////////////////////////////// -Eigen::Vector3d GerstnerWaveSimulation::Normal( +gz::math::Vector3d GerstnerWaveSimulation::Normal( double _x, double _y, double _t) const { double dhdx = 0.0, dhdy = 0.0; - const double r = StartupRamp(_t, this->tau_); - const std::size_t n = amplitudes_.size(); + const double r = StartupRamp(_t, this->tau); + const std::size_t n = this->amplitudes.size(); for (std::size_t i = 0; i < n; ++i) { - const auto &d = directions_[i]; - const double ak = amplitudes_[i] * wavenumbers_[i]; - const double theta = - wavenumbers_[i] * (d.x() * _x + d.y() * _y) - - angularFrequencies_[i] * _t + phase_; + const auto &d = this->directions[i]; + const double ak = this->amplitudes[i] * this->wavenumbers[i]; + const double theta = this->Phase(i, _x, _y, _t); const double s = std::sin(theta); - dhdx += -ak * d.x() * s; - dhdy += -ak * d.y() * s; + dhdx += -ak * d.X() * s; + dhdy += -ak * d.Y() * s; } dhdx *= r; dhdy *= r; - Eigen::Vector3d nvec(-dhdx, -dhdy, 1.0); - nvec.normalize(); + gz::math::Vector3d nvec(-dhdx, -dhdy, 1.0); + nvec.Normalize(); return nvec; } @@ -259,19 +262,18 @@ Eigen::Vector3d GerstnerWaveSimulation::Normal( double GerstnerWaveSimulation::Jacobian(double _x, double _y, double _t) const { double dsxdx = 0.0, dsydy = 0.0, dsxdy = 0.0; - const double r = StartupRamp(_t, this->tau_); - const std::size_t n = amplitudes_.size(); + const double r = StartupRamp(_t, this->tau); + const std::size_t n = this->amplitudes.size(); for (std::size_t i = 0; i < n; ++i) { - const auto &d = directions_[i]; - const double qak = steepnesses_[i] * amplitudes_[i] * wavenumbers_[i]; - const double theta = - wavenumbers_[i] * (d.x() * _x + d.y() * _y) - - angularFrequencies_[i] * _t + phase_; + const auto &d = this->directions[i]; + const double qak = + this->steepnesses[i] * this->amplitudes[i] * this->wavenumbers[i]; + const double theta = this->Phase(i, _x, _y, _t); const double c = std::cos(theta); - dsxdx += -qak * d.x() * d.x() * c; - dsydy += -qak * d.y() * d.y() * c; - dsxdy += -qak * d.x() * d.y() * c; + dsxdx += -qak * d.X() * d.X() * c; + dsydy += -qak * d.Y() * d.Y() * c; + dsxdy += -qak * d.X() * d.Y() * c; } dsxdx *= r; dsydy *= r; @@ -285,16 +287,16 @@ void GerstnerWaveSimulation::Update(double _simTime) // Analytic backend: the point queries (Elevation/Normal/...) stay // closed-form. Update only (re)samples the render grid at the new time, // matching the "Update advances, Field reads" model the FFT backend uses. - if (_simTime == this->currentTime_) + if (_simTime == this->currentTime) return; - this->currentTime_ = _simTime; + this->currentTime = _simTime; - const int N = static_cast(this->fieldN_); - if (N <= 0 || this->dzBuf_.empty()) + const int N = static_cast(this->fieldN); + if (N <= 0 || this->dzBuf.empty()) return; - const double T = this->fieldTile_; - const double r = StartupRamp(_simTime, this->tau_); - const std::size_t nc = this->amplitudes_.size(); + const double T = this->fieldTile; + const double r = StartupRamp(_simTime, this->tau); + const std::size_t nc = this->amplitudes.size(); for (int j = 0; j < N; ++j) { const double y = static_cast(j) * T / static_cast(N); @@ -305,32 +307,30 @@ void GerstnerWaveSimulation::Update(double _simTime) double dsxdx = 0.0, dsydy = 0.0, dsxdy = 0.0; for (std::size_t c = 0; c < nc; ++c) { - const auto &d = this->directions_[c]; - const double theta = - this->wavenumbers_[c] * (d.x() * x + d.y() * y) - - this->angularFrequencies_[c] * _simTime + this->phase_; + const auto &d = this->directions[c]; + const double theta = this->Phase(c, x, y, _simTime); const double s = std::sin(theta); const double co = std::cos(theta); - dz += this->amplitudes_[c] * co; + dz += this->amplitudes[c] * co; // Gerstner horizontal displacement (chop): -q·a·dir·sin(θ). - const double qa = this->steepnesses_[c] * this->amplitudes_[c]; - dx -= qa * d.x() * s; - dy -= qa * d.y() * s; + const double qa = this->steepnesses[c] * this->amplitudes[c]; + dx -= qa * d.X() * s; + dy -= qa * d.Y() * s; // Folding/foam terms reuse co — equivalent to Jacobian(x, y, t) but // without a second pass over the components. - const double qak = qa * this->wavenumbers_[c]; - dsxdx -= qak * d.x() * d.x() * co; - dsydy -= qak * d.y() * d.y() * co; - dsxdy -= qak * d.x() * d.y() * co; + const double qak = qa * this->wavenumbers[c]; + dsxdx -= qak * d.X() * d.X() * co; + dsydy -= qak * d.Y() * d.Y() * co; + dsxdy -= qak * d.X() * d.Y() * co; } const std::size_t idx = static_cast(i) + static_cast(j) * static_cast(N); - this->dzBuf_[idx] = dz * r; - this->dxBuf_[idx] = dx * r; - this->dyBuf_[idx] = dy * r; + this->dzBuf[idx] = dz * r; + this->dxBuf[idx] = dx * r; + this->dyBuf[idx] = dy * r; const double jxx = dsxdx * r, jyy = dsydy * r, jxy = dsxdy * r; - this->foamBuf_[idx] = (1.0 + jxx) * (1.0 + jyy) - jxy * jxy; + this->foamBuf[idx] = (1.0 + jxx) * (1.0 + jyy) - jxy * jxy; } } } @@ -338,7 +338,7 @@ void GerstnerWaveSimulation::Update(double _simTime) ////////////////////////////////////////////////// const WaveField2D *GerstnerWaveSimulation::Field() const { - return &this->field_; + return &this->field; } ////////////////////////////////////////////////// diff --git a/gz_waves_provider_gerstner/test/gerstner_test.cc b/gz_waves_provider_gerstner/test/gerstner_test.cc index 7045354c..2191bfca 100644 --- a/gz_waves_provider_gerstner/test/gerstner_test.cc +++ b/gz_waves_provider_gerstner/test/gerstner_test.cc @@ -172,10 +172,10 @@ TEST(ParticleVelocity, FiniteOnPmsField) for (double t = 5.0; t < 20.0; t += 1.0) { const auto v = gsw::ParticleVelocity(wf, 0.0, 0.0, t); - EXPECT_TRUE(std::isfinite(v.x())); - EXPECT_TRUE(std::isfinite(v.y())); - EXPECT_TRUE(std::isfinite(v.z())); - EXPECT_LT(v.norm(), 10.0); + EXPECT_TRUE(std::isfinite(v.X())); + EXPECT_TRUE(std::isfinite(v.Y())); + EXPECT_TRUE(std::isfinite(v.Z())); + EXPECT_LT(v.Length(), 10.0); } } @@ -186,8 +186,8 @@ TEST(Normal, IsUnitVector) for (double t = 5.0; t < 20.0; t += 1.0) { const auto n = gsw::Normal(wf, 1.0, 2.0, t); - EXPECT_NEAR(n.norm(), 1.0, 1e-9); - EXPECT_GT(n.z(), 0.0); + EXPECT_NEAR(n.Length(), 1.0, 1e-9); + EXPECT_GT(n.Z(), 0.0); } } From bb957dc34d6339c974fc22a6a9af679f65783d05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20Ag=C3=BCero?= Date: Mon, 22 Jun 2026 19:16:26 +0200 Subject: [PATCH 3/3] gz_waves_provider_gerstner: gzerr logging and a non-positive period guard MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Switch the engine's diagnostics from std::cerr to gzerr (wiring gz-common into the engine target), and reject a non-positive in SetParameters — which would otherwise make omega=2*pi/period infinite and the PMS spectrum blow up — leaving the field flat instead. Adds a test. Co-Authored-By: Claude Opus 4.8 (1M context) Signed-off-by: Carlos Agüero --- gz_waves_provider_gerstner/CMakeLists.txt | 5 ++++- gz_waves_provider_gerstner/package.xml | 1 + .../src/GerstnerWaveSimulation.cc | 22 ++++++++++++++++--- .../test/gerstner_test.cc | 21 ++++++++++++++++++ 4 files changed, 45 insertions(+), 4 deletions(-) diff --git a/gz_waves_provider_gerstner/CMakeLists.txt b/gz_waves_provider_gerstner/CMakeLists.txt index 35d232bc..573a9dca 100644 --- a/gz_waves_provider_gerstner/CMakeLists.txt +++ b/gz_waves_provider_gerstner/CMakeLists.txt @@ -22,6 +22,8 @@ find_package(gz_plugin_vendor REQUIRED) find_package(gz-plugin REQUIRED COMPONENTS register) find_package(gz_math_vendor REQUIRED) find_package(gz-math REQUIRED) # gz::math vectors + GZ_PI +find_package(gz_common_vendor REQUIRED) +find_package(gz-common REQUIRED) # gzerr/gzwarn console # -------------------------------------------------------------------------- # Gerstner wave-field engine (libgz-waves-provider-gerstner): a plain analytic @@ -39,7 +41,8 @@ target_include_directories(gz-waves-provider-gerstner PUBLIC target_include_directories(gz-waves-provider-gerstner PRIVATE ${gz_waves_INCLUDE_DIRS}) target_link_libraries(gz-waves-provider-gerstner - PUBLIC gz-math::gz-math) + PUBLIC gz-math::gz-math + PRIVATE gz-common::gz-common) # -------------------------------------------------------------------------- # Gerstner wave source system (gz-sim-waves-gerstner-system): a WavesSystemBase diff --git a/gz_waves_provider_gerstner/package.xml b/gz_waves_provider_gerstner/package.xml index 2333a198..ed79c2d5 100644 --- a/gz_waves_provider_gerstner/package.xml +++ b/gz_waves_provider_gerstner/package.xml @@ -17,6 +17,7 @@ gz_waves gz_sim_vendor gz_math_vendor + gz_common_vendor gz_plugin_vendor ament_cmake_gtest diff --git a/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc index 683576a1..87161586 100644 --- a/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc +++ b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc @@ -12,8 +12,8 @@ #include #include -#include +#include #include // GZ_PI #include "gz/sim/waves/Wavefield.hh" @@ -61,6 +61,22 @@ void GerstnerWaveSimulation::SetParameters(const WaveParameters &_params) { // Resolve (if set) into period/gain before configuring. const WaveParameters p = WithSeaState(_params); + + // A non-positive period makes omegaMean (2*pi/period) infinite/garbage and + // the Pierson-Moskowitz spectrum (∝ 1/omega^5) blow up. Reject it and leave + // the field flat rather than emit NaNs into elevation/normals downstream. + if (!(p.period > 0.0)) + { + gzerr << "[GerstnerWaveSimulation] non-positive (" << p.period + << "); leaving the wave field flat." << '\n'; + this->amplitudes.clear(); + this->wavenumbers.clear(); + this->angularFrequencies.clear(); + this->steepnesses.clear(); + this->directions.clear(); + return; + } + this->tau = p.tau; this->phase = p.phase; const double omegaMean = 2.0 * GZ_PI / p.period; @@ -116,8 +132,8 @@ void GerstnerWaveSimulation::SetParameters(const WaveParameters &_params) } else { - std::cerr << "[GerstnerWaveSimulation] unknown spectrum model '" - << p.model << "'; expected 'PMS' or 'CWR'." << '\n'; + gzerr << "[GerstnerWaveSimulation] unknown spectrum model '" + << p.model << "'; expected 'PMS' or 'CWR'." << '\n'; this->amplitudes.clear(); this->wavenumbers.clear(); this->angularFrequencies.clear(); diff --git a/gz_waves_provider_gerstner/test/gerstner_test.cc b/gz_waves_provider_gerstner/test/gerstner_test.cc index 2191bfca..906b34e0 100644 --- a/gz_waves_provider_gerstner/test/gerstner_test.cc +++ b/gz_waves_provider_gerstner/test/gerstner_test.cc @@ -290,3 +290,24 @@ TEST(Gerstner, SeaStateScalesWaveHeight) EXPECT_GT(sigma(2), 0.0); EXPECT_GT(sigma(6), sigma(2)) << "rougher sea state should produce bigger waves"; } + +// A non-positive would make omega = 2*pi/period infinite and the PMS +// spectrum (∝ 1/omega^5) blow up. SetParameters must reject it and leave the +// field flat rather than emit NaNs into elevation/normals. +////////////////////////////////////////////////// +TEST(Gerstner, NonPositivePeriodLeavesFieldFlat) +{ + for (double badPeriod : {0.0, -1.0}) + { + gsw::WaveParameters p; + p.model = "PMS"; + p.number = 3; + p.period = badPeriod; + gsw::GerstnerWaveSimulation sim; + sim.SetParameters(p); // must not crash + EXPECT_EQ(sim.Amplitudes().size(), 0u) << "period=" << badPeriod; + const double eta = sim.Elevation(1.0, 2.0, 3.0); + EXPECT_TRUE(std::isfinite(eta)) << "period=" << badPeriod; + EXPECT_DOUBLE_EQ(eta, 0.0) << "period=" << badPeriod; + } +}