Decouple irregular wave model from simulation parameters; add frequency-domain excitation transfer function

Remove simulation_dt_ and simulation_duration_ from IrregularWaveParams, eliminating the dependency that prevented adaptive time integration.

The excitation convolution is replaced by a pre-computed frequency-domain transfer function that evaluates in O(N_freq) per DOF per timestep instead of O(N_irf × N_freq). The eta-file-import path retains direct interpolation-based convolution. nfrequencies_ is now a required parameter (no fallback to simulation_duration_). A new ComputeElevationTimeSeries() utility replaces the old pre-computed arrays. Theory documentation updated with the derivation.

Author: dav-og

demos/sphere/demo_sphere_irreg_waves.cpp

@@ -103,8 +103,6 @@ int main(int argc, char* argv[]) {
     bodies.push_back(sphereBody);
 
     IrregularWaveParams wave_inputs;
-    wave_inputs.simulation_dt_       = timestep;
-    wave_inputs.simulation_duration_ = simulationDuration;
     wave_inputs.ramp_duration_       = 60.0;
     wave_inputs.wave_height_         = 2.0;
     wave_inputs.wave_period_         = 12.0;

demos/sphere/demo_sphere_irreg_waves_eta_import.cpp

@@ -106,8 +106,6 @@ int main(int argc, char* argv[]) {
     std::cout << "Defining irregular wave input parameters..." << std::endl;
     IrregularWaveParams params;
     std::cout << "bodies.size() = " << bodies.size() << std::endl;
-    params.simulation_dt_       = timestep;
-    params.simulation_duration_ = simulationDuration;
     params.ramp_duration_       = 0.0;
     params.eta_file_path_ = (DATADIR / "demos" / "sphere" / "eta" / "eta.txt").lexically_normal().generic_string();
     params.frequency_min_ = 0.001;

docs/_main_pages/theory.md

@@ -102,15 +102,50 @@ In HydroChrono, the force is computed through trapezoidal integration at the tim
 
 ### Wave excitation force, \\( F_{exc}(t) \\)
 
-The method to compute the wave excitation force involves convolution between the excitation Impulse Response Function (IRF) \\( K_{exc}(t) \\) and the wave elevation time sequence \\( \eta(t) \\):
+The wave excitation force is defined by the convolution of the excitation Impulse Response Function (IRF) \\( K_{exc}(\tau) \\) with the free-surface elevation \\( \eta(t) \\):
 
 $$
-F_{exc}(t) = \int_{-\infty}^{+\infty} K_{exc}(\tau) \eta(x, y, t-\tau) d\tau
+F_{exc}(t) = \int_{-\infty}^{+\infty} K_{exc}(\tau) \, \eta(x, y, t-\tau) \, d\tau
 $$
 
 By amalgamating these forces into the equation of motion, one can effectively model the behavior of a multibody oceanic system influenced by hydrodynamic forces.
 
-In HydroChrono, the force is computed through trapezoidal integration by discretizing at the time values given by the excitation IRF time array relative to the current simulation time step. Linear interpolation is done for the free surface elevation if a given time value is between two values of the time series of the precomputed free surface elevation.
+#### Numerical evaluation
+
+The continuous convolution integral is discretized using the quadrature points \\( \tau_j \\) and widths \\( \Delta\tau_j \\) provided by the excitation IRF time array from the BEM data:
+
+$$
+F_{exc,d}(t) = \sum_j K_{exc,d}(\tau_j) \, \eta(t - \tau_j) \, \Delta\tau_j
+$$
+
+where \\( d \\) denotes a particular degree of freedom.
+
+HydroChrono supports two evaluation strategies, depending on how the sea state is defined.
+
+**Frequency-domain transfer function (spectrum-based waves).** When the irregular sea state is represented as a linear superposition of \\( N_f \\) wave components,
+
+$$
+\eta(t) = \sum_{i=1}^{N_f} A_i \cos(\varphi_i - \omega_i \, t)
+$$
+
+the double summation over IRF points and frequency components can be factored by swapping the summation order and applying the angle-sum identity. This yields:
+
+$$
+F_{exc,d}(t) = \sum_{i=1}^{N_f} A_i \left[ C_{d,i} \cos(\varphi_i - \omega_i \, t) \;-\; S_{d,i} \sin(\varphi_i - \omega_i \, t) \right]
+$$
+
+where the *excitation transfer function* coefficients are pre-computed once from the IRF and the spectrum:
+
+$$
+C_{d,i} = \sum_j K_{exc,d}(\tau_j) \, \Delta\tau_j \, \cos(\omega_i \, \tau_j), \qquad
+S_{d,i} = \sum_j K_{exc,d}(\tau_j) \, \Delta\tau_j \, \sin(\omega_i \, \tau_j)
+$$
+
+This reduces the per-timestep cost from \\( \mathcal{O}(N_{irf} \times N_f) \\) to \\( \mathcal{O}(N_f) \\) — the same cost as a single free-surface elevation evaluation — while producing results identical to the time-domain convolution to machine precision. Crucially, because the wave elevation is evaluated analytically, this formulation has no dependency on the simulation time step or duration, enabling the use of adaptive time integrators.
+
+During the ramp-up period (when a cosine or linear ramp modifies \\( \eta \\)), the frequency-domain factoring does not apply. HydroChrono falls back to a batched matrix evaluation that computes \\( \eta \\) at all IRF quadrature points simultaneously via pre-computed trigonometric matrices, keeping the ramp startup efficient.
+
+**Interpolation-based evaluation (imported elevation time series).** When the free-surface elevation is provided as an external time series (e.g., from a wave tank measurement or a coupled model), the convolution is evaluated directly by linearly interpolating the stored \\( \eta \\) data at each quadrature point \\( t - \tau_j \\). This path retains the traditional \\( \mathcal{O}(N_{irf}) \\) per-timestep cost.
 
 ## Links to related open-source software:
 

include/hydroc/waves/irregular_wave.h

@@ -13,19 +13,17 @@
 #include <vector>
 
 struct IrregularWaveParams {
-    double simulation_dt_ = 0.0;
-    double simulation_duration_ = 0.0;
-    double ramp_duration_ = 0.0;
-    std::string eta_file_path_;
     double wave_height_             = 0.0;
     double wave_period_             = 0.0;
     double frequency_min_           = 0.001;
     double frequency_max_           = 1.0;
-    double nfrequencies_            = 0;
+    int    nfrequencies_            = 0;
     double peak_enhancement_factor_ = 1.0;
-    bool is_normalized_             = false;
-    int seed_                       = 1;
-    bool wave_stretching_           = true;
+    bool   is_normalized_           = false;
+    int    seed_                    = 1;
+    bool   wave_stretching_         = true;
+    double ramp_duration_           = 0.0;
+    std::string eta_file_path_;
 };
 
 class IrregularWaves : public WaveBase {
@@ -39,13 +37,18 @@ class IrregularWaves : public WaveBase {
     std::vector<double> GetFreeSurfaceTime() const;
     std::vector<double> GetFrequenciesHz() const;
 
+    std::pair<std::vector<double>, std::vector<double>>
+    ComputeElevationTimeSeries(double t_start, double t_end, double dt) const;
+
     Eigen::VectorXd GetForceAtTime(double t) override;
     WaveMode GetWaveMode() override { return mode_; }
 
     Eigen::VectorXd SetSpectrumFrequencies(double start, double end, int num_steps);
 
     void AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data, HydroData::SimulationParameters& sim_data);
 
+    double GetIRFTimeMax() const { return irf_time_max_; }
+
     double GetElevation(const Eigen::Vector3d& position, double time) override;
     Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time, double elevation) override;
     Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time, double elevation) override;
@@ -62,10 +65,8 @@ class IrregularWaves : public WaveBase {
   private:
     IrregularWaveParams params_;
     std::vector<double> spectrum_;
-    std::vector<double> time_data_;
-    std::vector<double> free_surface_elevation_sampled_;
-    std::vector<double> free_surface_time_sampled_;
     bool spectrumCreated_ = false;
+    bool use_eta_from_file_ = false;
 
     const WaveMode mode_ = WaveMode::irregular;
     std::vector<HydroData::IrregularWaveInfo> wave_info_;
@@ -78,31 +79,37 @@ class IrregularWaves : public WaveBase {
     Eigen::VectorXd wavenumbers_;
     Eigen::VectorXd wave_phases_;
 
-    // -------------------------------------------------------------------------
-    // Pre-computed arrays for fast elevation calculation.
-    // These are computed once from the spectrum and reused for every GetElevation() call.
-    // This optimization is critical for real-time visualization where GetElevation()
-    // may be called thousands of times per frame (once per water surface grid point).
-    // -------------------------------------------------------------------------
-    
-    /// Wave amplitude for each frequency component [m].
-    /// Computed as: A_i = sqrt(2 * S(f_i) * df_i), where S is the spectral density.
+    // Stored eta data — only populated when importing from file.
+    std::vector<double> time_data_;
+    std::vector<double> free_surface_elevation_sampled_;
+    std::vector<double> free_surface_time_sampled_;
+
+    // Pre-computed arrays for fast elevation calculation (SIMD-optimized).
+    // Computed once from the spectrum and reused for every GetElevation() call.
     Eigen::VectorXd amplitudes_;
-    
-    /// Angular frequency for each component [rad/s].
-    /// Computed as: omega_i = 2 * pi * f_i
     Eigen::VectorXd angular_freqs_;
 
+    // Pre-computed excitation transfer function coefficients.
+    // Converts the O(N_irf * N_freq) time-domain convolution into an
+    // O(N_freq) frequency-domain evaluation per DOF per timestep:
+    //   F_ex(dof, t) = sum_i A_i * [C(dof,i)*cos(theta_i) - S(dof,i)*sin(theta_i)]
+    // where theta_i = phi_i - omega_i * t.
+    std::vector<Eigen::MatrixXd> ex_transfer_cos_;  // [body](6, N_freq)
+    std::vector<Eigen::MatrixXd> ex_transfer_sin_;  // [body](6, N_freq)
+
+    // Pre-computed trig matrices for batch eta evaluation during ramp period.
+    // irf_cos_wt_[b](j, i) = cos(omega_i * tau_j), similarly for sin.
+    // Allows computing eta at all IRF points via matrix-vector product.
+    std::vector<Eigen::MatrixXd> irf_cos_wt_;  // [body](N_irf, N_freq)
+    std::vector<Eigen::MatrixXd> irf_sin_wt_;  // [body](N_irf, N_freq)
+    double irf_time_max_ = 0.0;
+
     void InitializeIRFVectors();
-    
-    /// Pre-compute amplitudes_ and angular_freqs_ arrays from the spectrum.
-    /// Called automatically after CreateSpectrum(). Must be called before GetElevation().
     void PrecomputeAmplitudes();
+    void PrecomputeExcitationTransfer();
     void ReadEtaFromFile();
-    void CreateFreeSurfaceElevation();
 
     Eigen::MatrixXd GetExcitationIRF(int b) const;
-    void ResampleIRF(double dt);
     void CalculateWidthIRF();
     double ExcitationConvolution(int body, int dof, double time);
 };

src/hydro/config/setup_from_yaml.cpp

@@ -60,11 +60,10 @@ std::shared_ptr<WaveBase> CreateWaveFromSettings(const WaveSettings& wave_settin
 
     } else if (type == "irregular") {
         IrregularWaveParams params;
-        params.simulation_dt_ = timestep;
-        params.simulation_duration_ = sim_duration;
         params.ramp_duration_ = ramp_duration;
         params.wave_height_ = wave_settings.height;
         params.wave_period_ = wave_settings.period;
+        params.nfrequencies_ = 1000;
         params.seed_ = (wave_settings.seed > 0 ? wave_settings.seed : 1);
 
         auto irregular_wave = std::make_shared<IrregularWaves>(params);

src/hydro/runner/run_from_yaml.cpp

@@ -748,8 +748,7 @@ int RunHydroChronoFromYAML(int argc, char* argv[]) {
                         auto irreg = std::static_pointer_cast<IrregularWaves>(wave_ptr);
                         std::vector<double> f = irreg->GetFrequenciesHz();
                         std::vector<double> S = irreg->GetSpectrum();
-                        std::vector<double> tvec = irreg->GetFreeSurfaceTime();
-                        std::vector<double> eta = irreg->GetFreeSurfaceElevation();
+                        auto [tvec, eta] = irreg->ComputeElevationTimeSeries(0.0, duration_hint, loop_dt);
                         exporter->WriteIrregularInputs(f, S, tvec, eta);
                     }
                 }