Add documentation for radiation damping convolution

Author: tridelat

doc/user/_theory/theory.rst

@@ -114,6 +114,7 @@ The :math:`K_{rad}(t)` function is derived by implementing the inverse continuou
 
 This transform allows the frequency domain coefficients, :math:`B(\omega)`, to be remapped into the time domain, thus producing the radiation impulse response function, :math:`K_{rad}(t)`. The :math:`B(\omega)` values can be sourced using open-source BEM software.
 
+In HydroChrono, the force is computed through trapezoidal integration at the time values given by the RIRF time array relative to the current simulation time step. Linear interpolation is done for the velocity history if a given time value is between two values of the time series of the stored velocity history.
 Wave excitation force, :math:`F_{exc}(t)`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

include/hydroc/hydro_forces.h

@@ -193,6 +193,13 @@ class TestHydro {
     /**
      * @brief Computes the Radiation Damping force with convolution history for a 6N dimensional system.
      *
+     * The discretization uses the time series of the the RIRF relative to the current step.
+     * Linear interpolation is done on the velocity history if time_sim-time_rirf is between two values of the time
+     * history. Trapezoidal integration is used to compute the force.
+     *
+     * Time history is automatically added in this function (so it should only be called once per time step), and
+     * history that is older than the maximum RIRF time value is automatically removed.
+     *
      * @return 6N dimensional force for 6 DOF and N bodies in system.
      */
     std::vector<double> ComputeForceRadiationDampingConv();

src/hydro_forces.cpp

@@ -306,14 +306,12 @@ std::vector<double> TestHydro::ComputeForceRadiationDampingConv() {
 
     assert(numRows * size > 0 && numCols > 0);
 
-    std::vector<double> tmp_s(numRows * size, 0.0);
-
-    // Helper function for tmp_s indexing
-    auto TmpSIndex = [&](int row, int step) { return (row * size) + step; };
-
     // time history
     auto t_sim = bodies_[0]->GetChTime();
     auto t_min = t_sim - rirf_time_vector.tail<1>()[0];
+    if (time_history_.size() > 0 && t_sim == time_history_.front()) {
+        throw std::runtime_error("Tried to compute the radiation damping convolution twice within the same time step!");
+    }
     time_history_.insert(time_history_.begin(), t_sim);
 
     // velocity history