Precalculate width for radiation convolution

Author: tridelat

include/hydroc/hydro_forces.h

@@ -244,8 +244,8 @@ class TestHydro {
     // Additional properties related to equilibrium and hydrodynamics
     std::vector<double> equilibrium_;
     std::vector<double> cb_minus_cg_;
-    double rirf_timestep_;
     Eigen::VectorXd rirf_time_vector;  // Assumed consistent for each body
+    Eigen::VectorXd rirf_width_vector;
 
     // Properties for velocity history management and time tracking
     std::vector<std::vector<std::vector<double>>> velocity_history_;

src/h5fileinfo.cpp

@@ -7,12 +7,11 @@
 // TODO: this include statement list looks good
 #include <H5Cpp.h>
 #include <hydroc/h5fileinfo.h>
-#include <filesystem> // std::filesystem::absolute
+#include <filesystem>  // std::filesystem::absolute
 
 using namespace chrono;  // TODO narrow this using namespace to specify what we use from chrono or put chrono:: in front
                          // of it all?
 
-
 H5FileInfo::H5FileInfo(std::string file, int num_bod) {
     h5_file_name_ = file;
     num_bodies_   = num_bod;
@@ -41,7 +40,7 @@ HydroData H5FileInfo::ReadH5Data() {
     for (int i = 0; i < num_bodies_; i++) {
         // body data
         data_to_init.body_data_[i].body_name = "body" + std::to_string(i + 1);
-        std::string bodyName                = data_to_init.body_data_[i].body_name;  // shortcut for reading later
+        std::string bodyName                 = data_to_init.body_data_[i].body_name;  // shortcut for reading later
         data_to_init.body_data_[i].body_num  = i;
 
         InitScalar(userH5File, bodyName + "/properties/disp_vol", data_to_init.body_data_[i].disp_vol);
@@ -54,7 +53,8 @@ HydroData H5FileInfo::ReadH5Data() {
 
         Init1D(userH5File, bodyName + "/properties/cg", data_to_init.body_data_[i].cg);
         Init1D(userH5File, bodyName + "/properties/cb", data_to_init.body_data_[i].cb);
-        Init2D(userH5File, bodyName + "/hydro_coeffs/linear_restoring_stiffness", data_to_init.body_data_[i].lin_matrix);
+        Init2D(userH5File, bodyName + "/hydro_coeffs/linear_restoring_stiffness",
+               data_to_init.body_data_[i].lin_matrix);
         Init2D(userH5File, bodyName + "/hydro_coeffs/added_mass/inf_freq", data_to_init.body_data_[i].inf_added_mass);
         data_to_init.body_data_[i].inf_added_mass *= rho;
         Init3D(userH5File, bodyName + "/hydro_coeffs/radiation_damping/impulse_response_fun/K",
@@ -242,5 +242,17 @@ int HydroData::GetRIRFDims(int i) const {
 }
 
 Eigen::VectorXd HydroData::GetRIRFTimeVector() const {
-    return body_data_[0].rirf_time_vector;
+    double tol = 1e-10;
+    // check if all time vectors are the same within tolerance
+    auto& rirf_time_vector = body_data_[0].rirf_time_vector;
+    for (size_t ii = 1; ii < body_data_.size(); ii++) {
+        for (size_t jj = 0; jj < body_data_[ii].rirf_time_vector.size(); jj++) {
+            if (abs(body_data_[ii].rirf_time_vector[jj] - rirf_time_vector[jj]) > tol) {
+                throw std::runtime_error(
+                    "RIRF time vectors have to be exactly the same for all bodies. Difference found in body " +
+                    std::to_string(jj) + " at time index " + std::to_string(jj) + ".");
+            }
+        }
+    }
+    return rirf_time_vector;
 }

src/hydro_forces.cpp

@@ -170,7 +170,17 @@ TestHydro::TestHydro(std::vector<std::shared_ptr<ChBody>> user_bodies,
 
     // Set up time vector
     rirf_time_vector = file_info_.GetRIRFTimeVector();
-    rirf_timestep_   = rirf_time_vector[1] - rirf_time_vector[0];
+    // width array
+    rirf_width_vector.resize(rirf_time_vector.size());
+    for (int ii = 0; ii < rirf_width_vector.size(); ii++) {
+        rirf_width_vector[ii] = 0.0;
+        if (ii < rirf_time_vector.size() - 1) {
+            rirf_width_vector[ii] += 0.5 * abs(rirf_time_vector[ii + 1] - rirf_time_vector[ii]);
+        }
+        if (ii > 0) {
+            rirf_width_vector[ii] += 0.5 * abs(rirf_time_vector[ii] - rirf_time_vector[ii - 1]);
+        }
+    }
 
     // Total degrees of freedom
     int total_dofs = kDofPerBody * num_bodies_;
@@ -303,7 +313,7 @@ std::vector<double> TestHydro::ComputeForceRadiationDampingConv() {
 
     // time history
     auto t_sim = bodies_[0]->GetChTime();
-    auto t_min = t_sim - rirf_timestep_ * size;
+    auto t_min = t_sim - rirf_time_vector.tail<1>()[0];
     time_history_.insert(time_history_.begin(), t_sim);
 
     // velocity history
@@ -333,7 +343,7 @@ std::vector<double> TestHydro::ComputeForceRadiationDampingConv() {
 
         // iterate over RIRF steps
         for (int step = 0; step < size; step++) {
-            auto t_rirf = t_sim - rirf_timestep_ * step;
+            auto t_rirf = t_sim - rirf_time_vector[step];
             while (time_history_[idx_history + 1] > t_rirf && idx_history < time_history_.size() - 1) {
                 idx_history += 1;
             }
@@ -368,19 +378,11 @@ std::vector<double> TestHydro::ComputeForceRadiationDampingConv() {
 
                     // iterate over rows
                     for (int row = 0; row < numRows; row++) {
-                        tmp_s[TmpSIndex(row, step)] += GetRIRFval(row, col, step) * vel_weighted;
+                        force_radiation_damping_[row] +=
+                            GetRIRFval(row, col, step) * vel_weighted * rirf_width_vector[step];
                     }
                 }
             }
-
-            if (step > 0) {
-                // cummulate values for current step on all rows
-                for (int row = 0; row < numRows; row++) {
-                    // Integrate tmp_s
-                    force_radiation_damping_[row] += (tmp_s[TmpSIndex(row, step - 1)] + tmp_s[TmpSIndex(row, step)]) /
-                                                     2.0 * (rirf_time_vector[step] - rirf_time_vector[step - 1]);
-                }
-            }
         }
     }
     return force_radiation_damping_;