Port PR #77 - wave stretching fixes to refactored wave modules (#80)

* Port Wheeler stretching fixes from PR #77 (by tridelat)

PR #77 was made before a major refactor - I've ported the changes over into the relevant new locations.

* Improve robustness in wave kinematics utilities

- Treat non-positive/very large/infinite water depth as “effectively infinite” in is_in_deep_water() (consistent with ComputeWaveNumber()), and clarify the 89.4 kh threshold is a numerical switch (tanh(kh)≈1).
- Avoid Wheeler stretching when depth is non-physical/unbounded (≤0, inf, or very large), and guard against h + eta near zero to prevent division-by-zero.
- Add checks and small‑kh limiting forms in water velocity/acceleration to avoid sinh(kh) division‑by‑zero / NaNs as k -> 0.

Author: dav-og

include/hydroc/waves/regular_wave.h

@@ -26,6 +26,7 @@ class RegularWave : public WaveBase {
     double regular_wave_amplitude_;
     double regular_wave_omega_;
     double regular_wave_phase_ = 0.0;
+    bool wave_stretching_      = true;
 
     void AddH5Data(std::vector<HydroData::RegularWaveInfo>& reg_h5_data, HydroData::SimulationParameters& sim_data);
 

src/hydro/waves/irregular_wave.cpp

@@ -110,9 +110,7 @@ void IrregularWaves::AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_
 Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, double time) {
     auto position_stretched = position;
     if (params_.wave_stretching_) {
-        auto eta = GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, wave_phases_, wavenumbers_);
-        auto z_pos = position.z() - mwl_;
-        position_stretched[2] = water_depth_ * (z_pos - eta) / (water_depth_ + eta);
+        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
     }
 
     return GetWaterVelocityIrregular(position_stretched,
@@ -129,9 +127,7 @@ Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, dou
 Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, double time) {
     auto position_stretched = position;
     if (params_.wave_stretching_) {
-        auto eta = GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, wave_phases_, wavenumbers_);
-        auto z_pos = position.z() - mwl_;
-        position_stretched[2] = water_depth_ * (z_pos - eta) / (water_depth_ + eta);
+        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
     }
 
     return GetWaterAccelerationIrregular(position_stretched,

src/hydro/waves/regular_wave.cpp

@@ -42,7 +42,11 @@ void RegularWave::AddH5Data(std::vector<HydroData::RegularWaveInfo>& reg_h5_data
 }
 
 Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double time) {
-    return GetWaterVelocity(position,
+    auto position_stretched = position;
+    if (wave_stretching_) {
+        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
+    }
+    return GetWaterVelocity(position_stretched,
                             time,
                             regular_wave_omega_,
                             regular_wave_amplitude_,
@@ -53,7 +57,11 @@ Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double
 }
 
 Eigen::Vector3d RegularWave::GetAcceleration(const Eigen::Vector3d& position, double time) {
-    return GetWaterAcceleration(position,
+    auto position_stretched = position;
+    if (wave_stretching_) {
+        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
+    }
+    return GetWaterAcceleration(position_stretched,
                                 time,
                                 regular_wave_omega_,
                                 regular_wave_amplitude_,

src/hydro/waves/wave_utilities.cpp

@@ -15,9 +15,43 @@
 #include <fstream>
 #include <iomanip>
 #include <iostream>
+#include <limits>
 
 #include <hydroc/logging.h>
 
+bool is_in_deep_water(double wavenumber, double water_depth) {
+    // Keep behavior consistent with ComputeWaveNumber(), which treats non-positive depth, very large
+    // depth, or infinite depth as "effectively infinite".
+    constexpr double EFFECTIVE_INFINITE_DEPTH = 1000.0;
+    if (water_depth <= 0.0 || water_depth > EFFECTIVE_INFINITE_DEPTH || std::isinf(water_depth)) {
+        return true;
+    }
+
+    // Numerical threshold for using deep-water asymptotic form (tanh(kh) ≈ 1).
+    return (std::abs(wavenumber) * water_depth > 89.4);
+}
+
+Eigen::Vector3d GetWheelerStretchedPosition(const Eigen::Vector3d& position, double eta, double water_depth, double mwl) {
+    // Wheeler stretching assumes finite depth. If depth is being used as a sentinel for "infinite",
+    // fall back to no stretching to avoid undefined behavior.
+    if (water_depth <= 0.0 || std::isinf(water_depth)) {
+        return position;
+    }
+
+    // Position relative to mean water level
+    double z_pos = position.z() - mwl;
+    // Wheeler stretching
+    double denom = water_depth + eta;
+    const double scale = std::max(1.0, std::abs(water_depth));
+    const double eps   = std::max(1e-12 * scale, 100.0 * std::numeric_limits<double>::epsilon() * scale);
+    if (std::abs(denom) < eps) {
+        // Pathological case (e.g., eta ≈ -water_depth): avoid division-by-zero blow-up.
+        return position;
+    }
+    double z_stretched = water_depth * (z_pos - eta) / denom;
+    return Eigen::Vector3d(position.x(), position.y(), z_stretched + mwl);
+}
+
 double GetEta(const Eigen::Vector3d& position,
               double time,
               double omega,
@@ -76,16 +110,31 @@ Eigen::Vector3d GetWaterVelocity(const Eigen::Vector3d& position,
     double z_pos = position.z() - mwl;
 
     Eigen::Vector3d water_velocity(0.0, 0.0, 0.0);
-    if (2 * M_PI / wavenumber > water_depth || wavenumber * water_depth > 500.0) {
+    const double k_mag = std::abs(wavenumber);
+    if (!std::isfinite(k_mag) || k_mag == 0.0 || !std::isfinite(omega) || omega == 0.0 || !std::isfinite(amplitude) || amplitude == 0.0) {
+        return water_velocity;
+    }
+
+    const double phase_arg = wavenumber * x_pos - omega * time + phase;
+
+    if (is_in_deep_water(wavenumber, water_depth)) {
         water_velocity[0] =
-            omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_pos - omega * time + phase);
+            omega * amplitude * std::exp(k_mag * z_pos) * cos(phase_arg);
         water_velocity[2] =
-            omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_pos - omega * time + phase);
+            omega * amplitude * std::exp(k_mag * z_pos) * sin(phase_arg);
     } else {
-        water_velocity[0] = omega * amplitude * std::cosh(wavenumber * (z_pos + water_depth)) /
-                            std::sinh(wavenumber * water_depth) * cos(wavenumber * x_pos - omega * time + phase);
-        water_velocity[2] = omega * amplitude * std::sinh(wavenumber * (z_pos + water_depth)) /
-                            std::sinh(wavenumber * water_depth) * sin(wavenumber * x_pos - omega * time + phase);
+        const double kh = k_mag * water_depth;
+        if (std::abs(kh) < 1e-8) {
+            // Small-kh safe form (avoids sinh(kh)→0 leading to inf*0 -> NaN).
+            // Uses the limiting behavior: cosh(k(z+h))/sinh(kh) ~ 1/(kh) and sinh(k(z+h))/sinh(kh) ~ (z+h)/h.
+            const double omega_over_k = omega / k_mag;
+            water_velocity[0]         = omega_over_k * amplitude / water_depth * cos(phase_arg);
+            water_velocity[2]         = omega * amplitude * ((z_pos + water_depth) / water_depth) * sin(phase_arg);
+        } else {
+            const double denom = std::sinh(kh);
+            water_velocity[0] = omega * amplitude * std::cosh(k_mag * (z_pos + water_depth)) / denom * cos(phase_arg);
+            water_velocity[2] = omega * amplitude * std::sinh(k_mag * (z_pos + water_depth)) / denom * sin(phase_arg);
+        }
     }
     return water_velocity;
 }
@@ -102,16 +151,30 @@ Eigen::Vector3d GetWaterAcceleration(const Eigen::Vector3d& position,
     double z_pos = position.z() - mwl;
 
     Eigen::Vector3d water_acceleration(0.0, 0.0, 0.0);
-    if (2 * M_PI / wavenumber > water_depth || wavenumber * water_depth > 500.0) {
+    const double k_mag = std::abs(wavenumber);
+    if (!std::isfinite(k_mag) || k_mag == 0.0 || !std::isfinite(omega) || omega == 0.0 || !std::isfinite(amplitude) || amplitude == 0.0) {
+        return water_acceleration;
+    }
+
+    const double phase_arg = wavenumber * x_pos - omega * time + phase;
+
+    if (is_in_deep_water(wavenumber, water_depth)) {
         water_acceleration[0] =
-            omega * omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_pos - omega * time + phase);
+            omega * omega * amplitude * std::exp(k_mag * z_pos) * sin(phase_arg);
         water_acceleration[2] =
-            -omega * omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_pos - omega * time + phase);
+            -omega * omega * amplitude * std::exp(k_mag * z_pos) * cos(phase_arg);
     } else {
-        water_acceleration[0] = omega * omega * amplitude * std::cosh(wavenumber * (z_pos + water_depth)) /
-                                std::sinh(wavenumber * water_depth) * sin(wavenumber * x_pos - omega * time + phase);
-        water_acceleration[2] = -omega * omega * amplitude * std::sinh(wavenumber * (z_pos + water_depth)) /
-                                std::sinh(wavenumber * water_depth) * cos(wavenumber * x_pos - omega * time + phase);
+        const double kh = k_mag * water_depth;
+        if (std::abs(kh) < 1e-8) {
+            // Small-kh safe form (avoids sinh(kh)→0 leading to inf*0 -> NaN).
+            const double omega_over_k = omega / k_mag;
+            water_acceleration[0]     = omega * omega_over_k * amplitude / water_depth * sin(phase_arg);
+            water_acceleration[2]     = -omega * omega * amplitude * ((z_pos + water_depth) / water_depth) * cos(phase_arg);
+        } else {
+            const double denom = std::sinh(kh);
+            water_acceleration[0] = omega * omega * amplitude * std::cosh(k_mag * (z_pos + water_depth)) / denom * sin(phase_arg);
+            water_acceleration[2] = -omega * omega * amplitude * std::sinh(k_mag * (z_pos + water_depth)) / denom * cos(phase_arg);
+        }
     }
     return water_acceleration;
 }

src/hydro/waves/wave_utilities.h

@@ -6,11 +6,19 @@
 #ifndef HYDRO_WAVES_WAVE_UTILITIES_H
 #define HYDRO_WAVES_WAVE_UTILITIES_H
 
+#include <cstddef>
 #include <Eigen/Dense>
 #include <array>
 #include <string>
 #include <vector>
 
+bool is_in_deep_water(double wavenumber, double water_depth);
+
+Eigen::Vector3d GetWheelerStretchedPosition(const Eigen::Vector3d& position,
+                                            double                eta,
+                                            double                water_depth,
+                                            double                mwl);
+
 double GetFreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
                                const Eigen::VectorXd& spectral_densities,
                                const Eigen::VectorXd& spectral_widths,