diff --git a/gz_waves_provider_gerstner/CMakeLists.txt b/gz_waves_provider_gerstner/CMakeLists.txt new file mode 100644 index 00000000..573a9dca --- /dev/null +++ b/gz_waves_provider_gerstner/CMakeLists.txt @@ -0,0 +1,123 @@ +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(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 +# 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 gz-math::gz-math + PRIVATE gz-common::gz-common) + +# -------------------------------------------------------------------------- +# 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(gz_math_vendor) + +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..6008c1a8 --- /dev/null +++ b/gz_waves_provider_gerstner/include/gz/sim/waves/GerstnerWaveSimulation.hh @@ -0,0 +1,148 @@ +/* + * 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 + +#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: gz::math::Vector3d ParticleVelocity( + double _x, double _y, double _t) const override; + // Documentation inherited + 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; + // 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 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; + /// \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..ed79c2d5 --- /dev/null +++ b/gz_waves_provider_gerstner/package.xml @@ -0,0 +1,30 @@ + + + + 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_common_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..87161586 --- /dev/null +++ b/gz_waves_provider_gerstner/src/GerstnerWaveSimulation.cc @@ -0,0 +1,369 @@ +/* + * 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_PI + +#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); + + // 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; + const std::size_t n = p.number; + + 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 + // 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 + { + gzerr << "[GerstnerWaveSimulation] unknown spectrum model '" + << p.model << "'; expected 'PMS' or 'CWR'." << '\n'; + this->amplitudes.clear(); + this->wavenumbers.clear(); + this->angularFrequencies.clear(); + this->steepnesses.clear(); + this->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))); + } + + this->amplitudes[i] = a; + this->wavenumbers[i] = k; + this->angularFrequencies[i] = omega; + this->steepnesses[i] = q; + + const double theta = nIdx * p.angle + p.direction; + 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; + 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 * 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 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/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 * GZ_PI / this->fieldTile; + for (std::size_t i = 0; i < n; ++i) + { + 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))); + } + } + 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 = this->amplitudes.size(); + for (std::size_t i = 0; i < n; ++i) + { + eta += this->amplitudes[i] * std::cos(this->Phase(i, _x, _y, _t)); + } + return eta * StartupRamp(_t, this->tau); +} + +////////////////////////////////////////////////// +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 = this->amplitudes.size(); + for (std::size_t i = 0; i < n; ++i) + { + 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; + vz += aw * c; + } + return gz::math::Vector3d(vx * r, vy * r, vz * r); +} + +////////////////////////////////////////////////// +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 = this->amplitudes.size(); + for (std::size_t i = 0; i < n; ++i) + { + 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 *= r; + dhdy *= r; + gz::math::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 = this->amplitudes.size(); + for (std::size_t i = 0; i < n; ++i) + { + 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 *= 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->Phase(c, x, y, _simTime); + 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..906b34e0 --- /dev/null +++ b/gz_waves_provider_gerstner/test/gerstner_test.cc @@ -0,0 +1,313 @@ +/* + * 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.Length(), 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.Length(), 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"; +} + +// 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; + } +}