Add doc for excitation convolution + check bounds

Author: tridelat

doc/user/_theory/theory.rst

@@ -117,7 +117,7 @@ This transform allows the frequency domain coefficients, :math:`B(\omega)`, to b
 Wave excitation force, :math:`F_{exc}(t)`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-The following method to compute the wave excitation force involves convolution between the excitation impulse response function :math:`K_{exc}(t)` and the wave elevation time sequence :math:`\eta(t)`:
+The following method to compute the wave excitation force involves convolution between the excitation Impulse Response Function (IRF) :math:`K_{exc}(t)` and the wave elevation time sequence :math:`\eta(t)`:
 
 .. math::
     :label: f_ex
@@ -126,6 +126,8 @@ The following method to compute the wave excitation force involves convolution b
 
 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.
+
 
 .. rubric:: Footnotes
 

include/hydroc/wave_types.h

@@ -366,6 +366,11 @@ class IrregularWaves : public WaveBase {
     /**
      * @brief Calculates the component of force from Convolution integral for specified body, dof, time.
      *
+     * The discretization uses the time series of the of the IRF relative to the current time step.
+     * Linear interpolation is done for the free surface elevation if time_sim-time_irf is between two
+     * values of the time series of the precomputed free surface elevation.
+     * Trapezoidal integration is used to compute the force.
+     *
      * @param body which body currently calculating for
      * @param dof which degree of freedom to calculate force value for
      * @param time the time to compute force for

src/wave_types.cpp

@@ -530,7 +530,7 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
     auto tmin    = free_surface_time_sampled_.front();
     auto tmax    = free_surface_time_sampled_.back();
     double t_tau = time - irf_time_array[0];
-    int idx   = 0;
+    int idx      = 0;
     if (t_tau <= tmin) {
         idx = 0;
     } else if (t_tau >= tmax) {
@@ -549,26 +549,38 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
                 idx -= 1;
             }
 
-            // linearly interpolate free surface elevation between bounds
             // free surface time values
             auto t1 = free_surface_time_sampled_[idx];
             auto t2 = free_surface_time_sampled_[idx + 1];
-            // double check that t_tau is between bounds
-            if (t_tau < t1 || t_tau > t2) {
+
+            // get free surface elevation
+            double eta_val;
+            if (t_tau == t1) {
+                eta_val = free_surface_time_sampled_[idx];
+            } else if (t_tau == t2) {
+                eta_val = free_surface_time_sampled_[idx + 1];
+            } else if (t_tau > t1 && t_tau < t2) {
+                // linearly interpolate free surface elevation between bounds
+                auto eta1 = free_surface_elevation_sampled_[idx];
+                auto eta2 = free_surface_elevation_sampled_[idx + 1];
+                // weights
+                auto w1 = (t2 - t_tau) / (t2 - t1);
+                auto w2 = 1.0 - w1;
+                // weighted value
+                auto eta_val = w1 * eta1 + w2 * eta2;
+            } else {
                 throw std::runtime_error("Excitation convolution: wrong tau value " + std::to_string(tau) +
                                          " not between " + std::to_string(t1) + " and " + std::to_string(t2) + ".");
             }
-            // free surface elevation values
-            auto eta1 = free_surface_elevation_sampled_[idx];
-            auto eta2 = free_surface_elevation_sampled_[idx + 1];
-            // weights
-            auto w1 = (t2 - t_tau) / (t2 - t1);
-            auto w2 = 1.0 - w1;
-            // weighted value
-            auto eta_val = w1 * eta1 + w2 * eta2;
 
             // add to excitation force
             f_ex += irf_val_mat(dof, j) * eta_val * irf_width_array[j];
+
+        } else {
+            throw std::runtime_error(
+                "Excitation convolution: trying to find free surface elevation at a time out of bounds from "
+                "the precomputed free surface elevation (" +
+                std::to_string(t_tau) + "not in [" + std::to_string(tmin) + ", " + std::to_string(tmax) + "]).");
         }
     }