Merge pull request #49 from tridelat/ex_interp

Interpolate when doing excitation force convolution

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/helper.h

@@ -3,40 +3,48 @@
 
 #pragma once
 
-#include <string>
-#include <iostream>
+#include <Eigen/Dense>  // Need for the container function
 #include <fstream>
-#include <Eigen/Dense> // Need for the container function
-
+#include <iostream>
+#include <string>
 
 #ifndef M_PI
     #define M_PI 3.14159265358979323846
 #endif
+
+/**@brief Returns last index of vector element below value.
+ *
+ * @param value Input value
+ * @param ticks Array of ticks from which to find lower-bound index (assuming ascending order)
+ *
+ */
+size_t get_lower_index(double value, const std::vector<double>& ticks);
+
 /**@brief Base namespace for HydroChrono library
- * 
-*/
+ *
+ */
 namespace hydroc {
 
 /**@brief Set initial environment
- * 
+ *
  * Use command line argument or env variables to set s
  * ome static data at initialization.
- * 
+ *
  * Set the main data directory ..
- * 
+ *
  * @param argc number of argument (same as for main function)
- * @param argv arguments of main function 
+ * @param argv arguments of main function
  * @return 1 on error 0 else
-*/
+ */
 int SetInitialEnvironment(int argc, char* argv[]) noexcept;
 
 /**@brief Get base name of data directory
- * 
+ *
  * @return the string containing the path in standard format
-*/
+ */
 std::string getDataDir() noexcept;
 template <typename T>
-void WriteDataToFile(const std::vector<T>& data, const std::string& filename){
+void WriteDataToFile(const std::vector<T>& data, const std::string& filename) {
     std::ofstream outFile(filename);
     if (outFile.is_open()) {
         for (const auto& item : data) {
@@ -66,7 +74,8 @@ void WriteContainerToFile(const Container& container, const std::string& file_na
  * @param file_name file to write container to
  */
 template <>
-void inline WriteContainerToFile<std::vector<double>>(const std::vector<double>& container, const std::string& file_name) {
+void inline WriteContainerToFile<std::vector<double>>(const std::vector<double>& container,
+                                                      const std::string& file_name) {
     std::ofstream output_file(file_name);
 
     if (!output_file) {
@@ -103,6 +112,6 @@ void inline WriteContainerToFile<Eigen::VectorXd>(const Eigen::VectorXd& contain
     output_file.close();
 };
 
-} // end namespace hydroc
+}  // end namespace hydroc
 
 #endif

include/hydroc/wave_types.h

@@ -13,15 +13,14 @@
 // todo move this helper function somewhere else?
 Eigen::VectorXd PiersonMoskowitzSpectrumHz(Eigen::VectorXd& f, double Hs, double Tp);
 
-Eigen::VectorXd JONSWAPSpectrumHz(Eigen::VectorXd& f, double Hs, double Tp, double gamma=3.3);
+Eigen::VectorXd JONSWAPSpectrumHz(Eigen::VectorXd& f, double Hs, double Tp, double gamma = 3.3);
 
 std::vector<double> FreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
                                          const Eigen::VectorXd& spectral_densities,
-                                         const Eigen::VectorXd& time_index,
+                                         const Eigen::VectorXd& time_array,
                                          double water_depth,
                                          int seed = 1);
 
-
 enum class WaveMode {
     /// @brief No waves
     noWaveCIC = 0,
@@ -187,7 +186,7 @@ class RegularWave : public WaveBase {
 };
 
 //// class to instantiate WaveBase for irregular waves
-//class IrregularWave : public WaveBase {
+// class IrregularWave : public WaveBase {
 //  public:
 //    IrregularWave();
 //    IrregularWave(unsigned int num_b);
@@ -207,7 +206,8 @@ class RegularWave : public WaveBase {
 //    double ramp_duration;
 //    Eigen::VectorXd eta;
 //
-//    void AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data, HydroData::SimulationParameters& sim_data);
+//    void AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data, HydroData::SimulationParameters&
+//    sim_data);
 //
 //  private:
 //    unsigned int num_bodies;
@@ -222,8 +222,8 @@ class RegularWave : public WaveBase {
 //
 //    // Eigen::MatrixXd GetExcitationIRF(int b) const;
 //    Eigen::VectorXd ResampleTime(const Eigen::VectorXd& t_old, const double dt_new);
-//    Eigen::MatrixXd ResampleVals(const Eigen::VectorXd& t_old, Eigen::MatrixXd& vals_old, const Eigen::VectorXd& t_new);
-//    double ExcitationConvolution(int body, int dof, double time);
+//    Eigen::MatrixXd ResampleVals(const Eigen::VectorXd& t_old, Eigen::MatrixXd& vals_old, const Eigen::VectorXd&
+//    t_new); double ExcitationConvolution(int body, int dof, double time);
 //
 //    void CreateSpectrum();
 //    void IrregularWave::CreateFreeSurfaceElevation();
@@ -236,10 +236,10 @@ struct IrregularWaveParams {
     double simulation_duration_;
     double ramp_duration_ = 0.0;
     std::string eta_file_path_;
-    double wave_height_ = 0.0;
-    double wave_period_ = 0.0;
+    double wave_height_             = 0.0;
+    double wave_period_             = 0.0;
     double peak_enhancement_factor_ = 1.0;
-    int seed_ = 1;
+    int seed_                       = 1;
 };
 
 class IrregularWaves : public WaveBase {
@@ -300,12 +300,13 @@ class IrregularWaves : public WaveBase {
     Eigen::Vector3<double> GetWaveMeshVelocity();
 
     // user input // TODO add default values in case user doesn't initialize these?
-    //double wave_height_;
-    //double wave_period_;
-    //double simulation_duration_;
-    //double simulation_dt_;
-    //double ramp_duration_;
-    //Eigen::VectorXd eta_;  // public for mesh printing functions, TODO maybe make those friends, so this can be private?
+    // double wave_height_;
+    // double wave_period_;
+    // double simulation_duration_;
+    // double simulation_dt_;
+    // double ramp_duration_;
+    // Eigen::VectorXd eta_;  // public for mesh printing functions, TODO maybe make those friends, so this can be
+    // private?
 
     /**
      * @brief Initializes other member variables for timestep calculations later.
@@ -329,16 +330,18 @@ class IrregularWaves : public WaveBase {
     int seed_;
     std::vector<double> spectrum_;
     std::vector<double> time_data_;
-    std::vector<double> free_surface_elevation_;
+    std::vector<double> free_surface_elevation_sampled_;
+    std::vector<double> free_surface_time_sampled_;
     bool spectrumCreated_;
 
     const WaveMode mode_ = WaveMode::irregular;
-    //unsigned int num_bodies_;
-    //const WaveMode mode_ = WaveMode::irregular;
+    // unsigned int num_bodies_;
+    // const WaveMode mode_ = WaveMode::irregular;
     std::vector<HydroData::IrregularWaveInfo> wave_info_;
     HydroData::SimulationParameters sim_data_;
-    std::vector<Eigen::MatrixXd> ex_irf_resampled_;
-    std::vector<Eigen::VectorXd> ex_irf_time_resampled_;
+    std::vector<Eigen::MatrixXd> ex_irf_sampled_;
+    std::vector<Eigen::VectorXd> ex_irf_time_sampled_;
+    std::vector<Eigen::VectorXd> ex_irf_width_sampled_;
     Eigen::VectorXd spectrum_frequencies_;
     Eigen::VectorXd spectral_densities_;
     std::string mesh_file_name_;
@@ -350,34 +353,24 @@ class IrregularWaves : public WaveBase {
 
     Eigen::MatrixXd GetExcitationIRF(int b) const;
 
-    /**
-     * @brief resamples irf time vector from old time vector and new timestep.
-     *
-     * Creates a resized time vector that starts and ends at the same values as t_old and has timestep t_new.
-     *
-     * @param t_old reference to old time vector from h5 file, the first and last element are transfered to the first
-     * and last element of t_new (return val)
-     * @param dt_new the time step to use in resampled vector, typically the timestep of chrono simulation
+    /** @brief Resamples IRF time, widths, and values.
      *
-     * @return newly sized vector that looks like \
-     * (t_old[0], t_old[0] + dt_new, t_old[0] + 2*dt_new, ... , t_old[t_old.size()-1])
+     * @param dt Time step value to resample
      */
-    Eigen::VectorXd ResampleTime(const Eigen::VectorXd& t_old, const double dt_new);
+    void ResampleIRF(double dt);
 
-    /**
-     * @brief Creates a resized values vector to interpolate values for a new timestep.
-     *
-     * @param t_old reference to old time vector from h5 file
-     * @param vals_old the original values corresponding to the times in t_old from h5 file
-     * @param t_new resampled times (return value from ResampleTime()) to use in interpolation
-     *
-     * @return matrix for interpolated vals_old over t_new
+    /** @brief Calculates width (used for excitation convolution).
      */
-    Eigen::MatrixXd ResampleVals(const Eigen::VectorXd& t_old, Eigen::MatrixXd& vals_old, const Eigen::VectorXd& t_new);
+    void IrregularWaves::CalculateWidthIRF();
 
     /**
      * @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/helper.cpp

@@ -5,6 +5,22 @@
 #include <memory>
 #include <vector>
 
+size_t get_lower_index(double value, const std::vector<double>& ticks) {
+    auto it = std::upper_bound(ticks.begin(), ticks.end(), value);
+    // get nearest-below index
+    size_t idx = it - ticks.begin() - 1;
+    // remove one if equal to value
+    if (ticks[idx] == value) {
+        idx -= 1;
+    }
+    if (idx <= 0 || idx >= ticks.size() - 1) {
+        throw std::runtime_error("Could not find index for value " + std::to_string(value) + " in array with bounds (" +
+                                 std::to_string(ticks.front()) + ", " + std::to_string(ticks.back()) + ").");
+    }
+    // return index
+    return idx;
+}
+
 using std::filesystem::path;
 
 static path DATADIR{};