From 5187d01b824e58470425e7c37e14dea2e0cf5264 Mon Sep 17 00:00:00 2001 From: Shabara Date: Sat, 12 Jul 2025 00:09:11 -0600 Subject: [PATCH 1/4] Adds wave directionality to regular waves --- include/hydroc/h5fileinfo.h | 1 + include/hydroc/helper.h | 49 +++++++++++++++++++++++++++ include/hydroc/wave_types.h | 10 ++++-- src/h5fileinfo.cpp | 8 ++++- src/wave_types.cpp | 66 ++++++++++++++++++++++++++++++++++--- 5 files changed, 127 insertions(+), 7 deletions(-) diff --git a/include/hydroc/h5fileinfo.h b/include/hydroc/h5fileinfo.h index ea7bbaca..b1902f61 100644 --- a/include/hydroc/h5fileinfo.h +++ b/include/hydroc/h5fileinfo.h @@ -55,6 +55,7 @@ class HydroData { }; struct RegularWaveInfo { Eigen::VectorXd freq_list; + Eigen::VectorXd wave_directions_list; Eigen::Tensor excitation_mag_matrix; Eigen::Tensor excitation_phase_matrix; }; diff --git a/include/hydroc/helper.h b/include/hydroc/helper.h index 5f5adaf9..bc76c1ec 100644 --- a/include/hydroc/helper.h +++ b/include/hydroc/helper.h @@ -112,6 +112,55 @@ void inline WriteContainerToFile(const Eigen::VectorXd& contain output_file.close(); }; + +/// +/// Performs bilinear interpolation on a 3D matrix along the wave direction and frequency axes. +/// This is used to interpolate a scalar value (e.g., excitation magnitude or phase) +/// at a given wave direction and frequency index between grid points. +/// +/// +/// A callable that takes three integers (i, j, k) and returns the matrix value at that location. +/// +/// Function or lambda used to access elements of the 3D matrix. +/// Index of the degree of freedom (DOF) or row in the matrix. +/// Continuous wave direction index (used for interpolation along direction axis). +/// Continuous frequency index (used for interpolation along frequency axis). +/// Total number of discrete wave direction entries in the matrix. +/// Total number of discrete frequency entries in the matrix. +/// Interpolated scalar value at the specified DOF, direction, and frequency. +template +double BilinearInterp3D(const MatrixGetter& get_val, + int i, + double dir_index_des, + double freq_index_des, + int num_dirs, + int num_freqs) { + int freq_floor = static_cast(std::floor(freq_index_des)); + int freq_ceil = freq_floor + 1; + double freq_alpha = freq_index_des - freq_floor; + + int dir_floor = static_cast(std::floor(dir_index_des)); + int dir_ceil = dir_floor + 1; + double dir_alpha = dir_index_des - dir_floor; + + // Clamp to valid index ranges + freq_floor = std::max(0, std::min(freq_floor, num_freqs - 2)); + freq_ceil = freq_floor + 1; + + dir_floor = std::max(0, std::min(dir_floor, num_dirs - 2)); + dir_ceil = dir_floor + 1; + + // Access values + double v00 = get_val(i, dir_floor, freq_floor); + double v10 = get_val(i, dir_ceil, freq_floor); + double v01 = get_val(i, dir_floor, freq_ceil); + double v11 = get_val(i, dir_ceil, freq_ceil); + + return (1 - freq_alpha) * (1 - dir_alpha) * v00 + (1 - freq_alpha) * dir_alpha * v10 + + freq_alpha * (1 - dir_alpha) * v01 + freq_alpha * dir_alpha * v11; +} + + } // end namespace hydroc #endif \ No newline at end of file diff --git a/include/hydroc/wave_types.h b/include/hydroc/wave_types.h index ed7993c1..ca56d9b2 100644 --- a/include/hydroc/wave_types.h +++ b/include/hydroc/wave_types.h @@ -7,6 +7,7 @@ * @brief header file for Wavebase and classes inheriting from WaveBase. *********************************************************************/ #pragma once + #include #include @@ -165,6 +166,7 @@ class RegularWave : public WaveBase { double regular_wave_amplitude_; double regular_wave_omega_; double regular_wave_phase_ = 0.0; + double regular_wave_direction_ = 0.0; /** * @brief Initializes other member variables for timestep calculations later. @@ -200,6 +202,8 @@ class RegularWave : public WaveBase { */ double GetOmegaDelta() const; + double GetInterpolatedDirectionIndex(double dir_input) const; + /** * @brief gets interpolated excitation magnitude value. * @@ -212,7 +216,8 @@ class RegularWave : public WaveBase { * * @return excitation magnitudes for body b, row i, column j, frequency ix k */ - double GetExcitationMagInterp(int b, int i, int j, double freq_index_des) const; + //double GetExcitationMagInterp(int b, int i, int j, double freq_index_des) const; + double GetExcitationMagInterp(int b, int i, double dir_index_des, double freq_index_des) const; /** * @brief gets interpolated excitation phase value. @@ -226,7 +231,8 @@ class RegularWave : public WaveBase { * * @return excitation phases for body b, row i, column j, frequency ix k */ - double GetExcitationPhaseInterp(int b, int i, int j, double freq_index_des) const; + // double GetExcitationPhaseInterp(int b, int i, int j, double freq_index_des) const; + double GetExcitationPhaseInterp(int b, int i, double dir_index_des, double freq_index_des) const; }; //// class to instantiate WaveBase for irregular waves diff --git a/src/h5fileinfo.cpp b/src/h5fileinfo.cpp index 553fbeb2..561b53a6 100644 --- a/src/h5fileinfo.cpp +++ b/src/h5fileinfo.cpp @@ -5,9 +5,12 @@ * H5FileInfo. *********************************************************************/ // TODO: this include statement list looks good +#include "hydroc/h5fileinfo.h" +#include "chrono/utils/ChConstants.h" + #include -#include #include // std::filesystem::absolute +#include using namespace chrono; // TODO narrow this using namespace to specify what we use from chrono or put chrono:: in front // of it all? @@ -64,9 +67,12 @@ HydroData H5FileInfo::ReadH5Data() { // reg wave Init1D(userH5File, "simulation_parameters/w", data_to_init.reg_wave_data_[i].freq_list); + Init1D(userH5File, "simulation_parameters/wave_dir", data_to_init.reg_wave_data_[i].wave_directions_list); Init3D(userH5File, bodyName + "/hydro_coeffs/excitation/mag", data_to_init.reg_wave_data_[i].excitation_mag_matrix); + data_to_init.reg_wave_data_[i].wave_directions_list *= CH_DEG_TO_RAD; + // scale by rho * g data_to_init.reg_wave_data_[i].excitation_mag_matrix = data_to_init.reg_wave_data_[i].excitation_mag_matrix * diff --git a/src/wave_types.cpp b/src/wave_types.cpp index f9290cf7..b1195c8d 100644 --- a/src/wave_types.cpp +++ b/src/wave_types.cpp @@ -3,10 +3,13 @@ * * @brief implementation file for Wavebase and classes inheriting from WaveBase. *********************************************************************/ -#include -#include +#include "hydroc/helper.h" +#include "hydroc/wave_types.h" +#include "chrono/utils/ChConstants.h" + #include +using namespace chrono; double GetEta(const Eigen::Vector3d& position, double time, double omega, @@ -218,14 +221,22 @@ void RegularWave::AddH5Data(std::vector& reg_h5_data double wave_omega_delta = GetOmegaDelta(); double freq_index_des = (regular_wave_omega_ / wave_omega_delta) - 1; + + double dir_index_des = GetInterpolatedDirectionIndex(regular_wave_direction_); + + std::cout << "dir_index_des = " << dir_index_des << std::endl; + for (int b = 0; b < num_bodies_; b++) { for (int rowEx = 0; rowEx < 6; rowEx++) { int body_offset = 6 * b; // why are these always 0? vvv TODO check/change this - excitation_force_mag_[body_offset + rowEx] = GetExcitationMagInterp(b, rowEx, 0, freq_index_des); - excitation_force_phase_[body_offset + rowEx] = GetExcitationPhaseInterp(b, rowEx, 0, freq_index_des); + excitation_force_mag_[body_offset + rowEx] = GetExcitationMagInterp(b, rowEx, dir_index_des, freq_index_des); + excitation_force_phase_[body_offset + rowEx] = GetExcitationPhaseInterp(b, rowEx, dir_index_des, freq_index_des); } } + std::cout << "excitation_force_mag_ = " << excitation_force_mag_ << std::endl; + std::cout << "excitation_force_phase_ = " << excitation_force_phase_ << std::endl; + std::cout << "dir_index_des " << dir_index_des << "freq_index_des " << freq_index_des << std::endl; } Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double time) { @@ -262,6 +273,37 @@ double RegularWave::GetOmegaDelta() const { return omega_max / num_freqs; } +double RegularWave::GetInterpolatedDirectionIndex(double dir_input) const { + const auto& dirs = wave_info_[0].wave_directions_list; + const size_t N = dirs.size(); + + if (N < 2) { + throw std::runtime_error("Not enough wave directions to interpolate."); + } + + // Wrap input to [0, 2pi) + double wrapped_dir = fmod(fmod(dir_input, CH_PI) + CH_PI, CH_PI); + std::cout << "wrapped_dir = " << wrapped_dir << std::endl; + // Linear search for lower bracket + for (size_t i = 0; i < N - 1; ++i) { + if (wrapped_dir >= dirs[i] && wrapped_dir <= dirs[i + 1]) { + double delta = dirs[i + 1] - dirs[i]; + double frac = (wrapped_dir - dirs[i]) / delta; + return static_cast(i) + frac; + } + } + + // Wraparound interpolation between last and first (e.g., 2pi and 0) + if (wrapped_dir >= dirs[dirs.size() - 1]) { + double delta = (CH_PI)-dirs[dirs.size() - 1] + dirs[0]; + double frac = (wrapped_dir - dirs[dirs.size() - 1]) / delta; + return static_cast(N - 1) + frac; + } + + throw std::runtime_error("Could not find interpolation interval for direction input."); +} + +/* double RegularWave::GetExcitationMagInterp(int b, int i, int j, double freq_index_des) const { double freq_interp_val = freq_index_des - floor(freq_index_des); double excitationMagFloor = wave_info_[b].excitation_mag_matrix(i, j, (int)floor(freq_index_des)); @@ -280,6 +322,22 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in return excitationPhase; } +*/ + +double RegularWave::GetExcitationMagInterp(int b, int i, double dir_index_des, double freq_index_des) const { + const auto& mat = wave_info_[b].excitation_mag_matrix; + return hydroc::BilinearInterp3D( + [&](int i_, int j_, int k_) { return mat(i_, j_, k_); }, + i, dir_index_des, freq_index_des, mat.dimension(1), mat.dimension(2)); +} + +double RegularWave::GetExcitationPhaseInterp(int b, int i, double dir_index_des, double freq_index_des) const { + const auto& mat = wave_info_[b].excitation_phase_matrix; + return hydroc::BilinearInterp3D( + [&](int i_, int j_, int k_) { return mat(i_, j_, k_); }, + i, dir_index_des, freq_index_des, mat.dimension(1), mat.dimension(2)); +} + Eigen::VectorXd ComputeWaveNumbers(const Eigen::VectorXd& omegas, double water_depth, From c5ac24fcecb73005f9e7cae2808a283ef922823d Mon Sep 17 00:00:00 2001 From: Shabara Date: Sat, 12 Jul 2025 00:39:40 -0600 Subject: [PATCH 2/4] removes log lines --- include/hydroc/helper.h | 1 - src/h5fileinfo.cpp | 1 - src/wave_types.cpp | 8 -------- 3 files changed, 10 deletions(-) diff --git a/include/hydroc/helper.h b/include/hydroc/helper.h index bc76c1ec..107d81b5 100644 --- a/include/hydroc/helper.h +++ b/include/hydroc/helper.h @@ -160,7 +160,6 @@ double BilinearInterp3D(const MatrixGetter& get_val, freq_alpha * (1 - dir_alpha) * v01 + freq_alpha * dir_alpha * v11; } - } // end namespace hydroc #endif \ No newline at end of file diff --git a/src/h5fileinfo.cpp b/src/h5fileinfo.cpp index 561b53a6..82b62f83 100644 --- a/src/h5fileinfo.cpp +++ b/src/h5fileinfo.cpp @@ -72,7 +72,6 @@ HydroData H5FileInfo::ReadH5Data() { data_to_init.reg_wave_data_[i].excitation_mag_matrix); data_to_init.reg_wave_data_[i].wave_directions_list *= CH_DEG_TO_RAD; - // scale by rho * g data_to_init.reg_wave_data_[i].excitation_mag_matrix = data_to_init.reg_wave_data_[i].excitation_mag_matrix * diff --git a/src/wave_types.cpp b/src/wave_types.cpp index b1195c8d..abd9cb60 100644 --- a/src/wave_types.cpp +++ b/src/wave_types.cpp @@ -221,11 +221,7 @@ void RegularWave::AddH5Data(std::vector& reg_h5_data double wave_omega_delta = GetOmegaDelta(); double freq_index_des = (regular_wave_omega_ / wave_omega_delta) - 1; - double dir_index_des = GetInterpolatedDirectionIndex(regular_wave_direction_); - - std::cout << "dir_index_des = " << dir_index_des << std::endl; - for (int b = 0; b < num_bodies_; b++) { for (int rowEx = 0; rowEx < 6; rowEx++) { int body_offset = 6 * b; @@ -234,9 +230,6 @@ void RegularWave::AddH5Data(std::vector& reg_h5_data excitation_force_phase_[body_offset + rowEx] = GetExcitationPhaseInterp(b, rowEx, dir_index_des, freq_index_des); } } - std::cout << "excitation_force_mag_ = " << excitation_force_mag_ << std::endl; - std::cout << "excitation_force_phase_ = " << excitation_force_phase_ << std::endl; - std::cout << "dir_index_des " << dir_index_des << "freq_index_des " << freq_index_des << std::endl; } Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double time) { @@ -338,7 +331,6 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, double dir_index_des, i, dir_index_des, freq_index_des, mat.dimension(1), mat.dimension(2)); } - Eigen::VectorXd ComputeWaveNumbers(const Eigen::VectorXd& omegas, double water_depth, double g, From 3521fd1e0bdbc1b26f15ba07962ce77014a92c88 Mon Sep 17 00:00:00 2001 From: Shabara Date: Mon, 21 Jul 2025 15:43:55 -0600 Subject: [PATCH 3/4] Adds multiple wave headings for irregular waves --- include/hydroc/h5fileinfo.h | 5 +- include/hydroc/wave_types.h | 19 +- src/h5fileinfo.cpp | 10 +- src/wave_types.cpp | 436 +++++++++++++++++++++++++----------- 4 files changed, 328 insertions(+), 142 deletions(-) diff --git a/include/hydroc/h5fileinfo.h b/include/hydroc/h5fileinfo.h index b1902f61..a1f1b09d 100644 --- a/include/hydroc/h5fileinfo.h +++ b/include/hydroc/h5fileinfo.h @@ -55,7 +55,7 @@ class HydroData { }; struct RegularWaveInfo { Eigen::VectorXd freq_list; - Eigen::VectorXd wave_directions_list; + Eigen::VectorXd wave_direction_list; Eigen::Tensor excitation_mag_matrix; Eigen::Tensor excitation_phase_matrix; }; @@ -65,11 +65,12 @@ class HydroData { // Eigen::Tensor excitation_im_matrix; // Eigen::Vector3i im_dims; Eigen::VectorXd excitation_irf_time; - Eigen::MatrixXd excitation_irf_matrix; // TODO needs to be tensor? + Eigen::Tensor excitation_irf_matrix; // see std::optional documentation for how to use std::optional excitation_irf_resampled; // TODO needs to be tensor? std::optional excitation_irf_time_resampled; + Eigen::VectorXd wave_direction_list; }; private: diff --git a/include/hydroc/wave_types.h b/include/hydroc/wave_types.h index ca56d9b2..d2d29b81 100644 --- a/include/hydroc/wave_types.h +++ b/include/hydroc/wave_types.h @@ -79,6 +79,9 @@ class WaveBase { double g_ = 9.81; /// @brief Water depth double water_depth_ = 0.0; + +protected: + double GetInterpolatedDirectionIndex(const Eigen::VectorXd& dirs, double dir_input) const; }; /** @@ -167,6 +170,7 @@ class RegularWave : public WaveBase { double regular_wave_omega_; double regular_wave_phase_ = 0.0; double regular_wave_direction_ = 0.0; + double regular_wave_spread_ = 1; /** * @brief Initializes other member variables for timestep calculations later. @@ -290,11 +294,13 @@ struct IrregularWaveParams { 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; + std::vector wave_direction_ = {0.0}; + std::vector wave_spread_ = {1.0}; }; class IrregularWaves : public WaveBase { @@ -304,7 +310,7 @@ class IrregularWaves : public WaveBase { void CreateSpectrum(); std::vector GetSpectrum(); - std::vector GetFreeSurfaceElevation(); + Eigen::MatrixXd GetFreeSurfaceElevation(); std::vector GetEtaTimeData(); Eigen::VectorXd GetForceAtTime(double t) override; @@ -384,7 +390,8 @@ class IrregularWaves : public WaveBase { IrregularWaveParams params_; std::vector spectrum_; std::vector time_data_; - std::vector free_surface_elevation_sampled_; + // std::vector free_surface_elevation_sampled_; + Eigen::MatrixXd free_surface_elevation_sampled_; std::vector free_surface_time_sampled_; bool spectrumCreated_; @@ -392,21 +399,21 @@ class IrregularWaves : public WaveBase { // unsigned int num_bodies_; // const WaveMode mode_ = WaveMode::irregular; std::vector wave_info_; - std::vector ex_irf_sampled_; + std::vector> ex_irf_sampled_; std::vector ex_irf_time_sampled_; std::vector ex_irf_width_sampled_; Eigen::VectorXd spectrum_frequencies_; Eigen::VectorXd spectral_densities_; Eigen::VectorXd spectral_widths_; Eigen::VectorXd wavenumbers_; - Eigen::VectorXd wave_phases_; + Eigen::MatrixXd wave_phases_; std::string mesh_file_name_; void InitializeIRFVectors(); void ReadEtaFromFile(); void CreateFreeSurfaceElevation(); - Eigen::MatrixXd GetExcitationIRF(int b) const; + Eigen::Tensor GetExcitationIRF(int b) const; /** @brief Resamples IRF time, widths, and values. * diff --git a/src/h5fileinfo.cpp b/src/h5fileinfo.cpp index 82b62f83..dab4e5dd 100644 --- a/src/h5fileinfo.cpp +++ b/src/h5fileinfo.cpp @@ -11,6 +11,7 @@ #include #include // std::filesystem::absolute #include +#include using namespace chrono; // TODO narrow this using namespace to specify what we use from chrono or put chrono:: in front // of it all? @@ -67,11 +68,11 @@ HydroData H5FileInfo::ReadH5Data() { // reg wave Init1D(userH5File, "simulation_parameters/w", data_to_init.reg_wave_data_[i].freq_list); - Init1D(userH5File, "simulation_parameters/wave_dir", data_to_init.reg_wave_data_[i].wave_directions_list); + Init1D(userH5File, "simulation_parameters/wave_dir", data_to_init.reg_wave_data_[i].wave_direction_list); Init3D(userH5File, bodyName + "/hydro_coeffs/excitation/mag", data_to_init.reg_wave_data_[i].excitation_mag_matrix); - data_to_init.reg_wave_data_[i].wave_directions_list *= CH_DEG_TO_RAD; + data_to_init.reg_wave_data_[i].wave_direction_list *= CH_DEG_TO_RAD; // scale by rho * g data_to_init.reg_wave_data_[i].excitation_mag_matrix = data_to_init.reg_wave_data_[i].excitation_mag_matrix * @@ -89,8 +90,9 @@ HydroData H5FileInfo::ReadH5Data() { // TODO look up Eigen resize and map to make this temp conversion better Eigen::Tensor temp; Init3D(userH5File, bodyName + "/hydro_coeffs/excitation/impulse_response_fun/f", temp); - data_to_init.irreg_wave_data_[i].excitation_irf_matrix = SqueezeMid(temp); - data_to_init.irreg_wave_data_[i].excitation_irf_matrix *= rho * g; + auto scaled_temp = temp * temp.constant(rho * g); + data_to_init.irreg_wave_data_[i].excitation_irf_matrix = (temp * temp.constant(rho * g)).eval(); + data_to_init.irreg_wave_data_[i].wave_direction_list = data_to_init.reg_wave_data_[i].wave_direction_list; } userH5File.close(); diff --git a/src/wave_types.cpp b/src/wave_types.cpp index abd9cb60..e909cf37 100644 --- a/src/wave_types.cpp +++ b/src/wave_types.cpp @@ -8,6 +8,7 @@ #include "chrono/utils/ChConstants.h" #include +#include using namespace chrono; double GetEta(const Eigen::Vector3d& position, @@ -15,11 +16,14 @@ double GetEta(const Eigen::Vector3d& position, double omega, double amplitude, double phase, - double wavenumber) { + double wavenumber, + double wave_direction) { // assuming wave direction along global X axis auto x_pos = position.x(); + auto y_pos = position.y(); + auto x_rot = x_pos * cos(wave_direction) + y_pos * sin(wave_direction); - auto eta = amplitude * cos(wavenumber * x_pos - omega * time + phase); + auto eta = amplitude * cos(wavenumber * x_rot - omega * time + phase); return eta; }; @@ -29,30 +33,33 @@ double GetEtaIrregular(const Eigen::Vector3d& position, const Eigen::VectorXd& spectral_densities, const Eigen::VectorXd& spectral_widths, const Eigen::VectorXd& wave_phases, - const Eigen::VectorXd& wavenumbers) { + const Eigen::VectorXd& wavenumbers, + const double wave_direction, + const double wave_spread) { // x position assuming wave direction along global X axis double x_pos = position.x(); - double eta = 0.0; for (size_t i = 0; i < freqs_hz.size(); ++i) { - auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i]); + auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i] * wave_spread); auto omega = 2 * M_PI * freqs_hz[i]; - eta += GetEta(position, time, omega, amplitude, wave_phases[i], wavenumbers[i]); + eta += GetEta(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], wave_direction); } return eta; } -std::vector GetEtaIrregularTimeSeries(const Eigen::Vector3d& position, +Eigen::RowVectorXd GetEtaIrregularTimeSeries(const Eigen::Vector3d& position, const std::vector time_index, const Eigen::VectorXd& freqs_hz, const Eigen::VectorXd& spectral_densities, const Eigen::VectorXd& spectral_widths, const Eigen::VectorXd& wave_phases, - const Eigen::VectorXd& wavenumbers) { - std::vector eta(time_index.size(), 0.0); + const Eigen::VectorXd& wavenumbers, + const double wave_direction, + const double wave_spread) { + Eigen::RowVectorXd eta(time_index.size()); for (size_t j = 0; j < time_index.size(); ++j) { eta[j] = GetEtaIrregular(position, time_index[j], freqs_hz, spectral_densities, spectral_widths, wave_phases, - wavenumbers); + wavenumbers, wave_direction, wave_spread); } return eta; } @@ -64,9 +71,13 @@ Eigen::Vector3d GetWaterVelocity(const Eigen::Vector3d& position, double phase, double wavenumber, double water_depth, - double mwl) { + double mwl, + double wave_direction, + double wave_spread) { // assuming wave along global X axis position auto x_pos = position.x(); + auto y_pos = position.y(); + auto x_rot = x_pos * cos(wave_direction) + y_pos * sin(wave_direction); // position relative to mean water level auto z_pos = position.z() - mwl; @@ -75,15 +86,15 @@ Eigen::Vector3d GetWaterVelocity(const Eigen::Vector3d& position, if (2 * M_PI / wavenumber > water_depth || wavenumber * water_depth > 500.0) { // deep water water_velocity[0] = - omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_pos - omega * time + phase); + omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_rot - omega * time + phase); water_velocity[2] = - omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_pos - omega * time + phase); + omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_rot - omega * time + phase); } else { // shallow water water_velocity[0] = omega * amplitude * std::cosh(wavenumber * (z_pos + water_depth)) / - std::sinh(wavenumber * water_depth) * cos(wavenumber * x_pos - omega * time + phase); + std::sinh(wavenumber * water_depth) * cos(wavenumber * x_rot - omega * time + phase); water_velocity[2] = omega * amplitude * std::sinh(wavenumber * (z_pos + water_depth)) / - std::sinh(wavenumber * water_depth) * sin(wavenumber * x_pos - omega * time + phase); + std::sinh(wavenumber * water_depth) * sin(wavenumber * x_rot - omega * time + phase); } return water_velocity; @@ -96,9 +107,13 @@ Eigen::Vector3d GetWaterAcceleration(const Eigen::Vector3d& position, double phase, double wavenumber, double water_depth, - double mwl) { + double mwl, + double wave_direction, + double wave_spread) { // assuming wave along global X axis position auto x_pos = position.x(); + auto y_pos = position.y(); + auto x_rot = x_pos * cos(wave_direction) + y_pos * sin(wave_direction); // position relative to mean water level auto z_pos = position.z() - mwl; @@ -107,15 +122,15 @@ Eigen::Vector3d GetWaterAcceleration(const Eigen::Vector3d& position, if (2 * M_PI / wavenumber > water_depth || wavenumber * water_depth > 500.0) { // deep water water_acceleration[0] = - omega * omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_pos - omega * time + phase); + omega * omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_rot - omega * time + phase); water_acceleration[2] = - -omega * omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_pos - omega * time + phase); + -omega * omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_rot - omega * time + phase); } else { // shallow water water_acceleration[0] = omega * omega * amplitude * std::cosh(wavenumber * (z_pos + water_depth)) / - std::sinh(wavenumber * water_depth) * sin(wavenumber * x_pos - omega * time + phase); + std::sinh(wavenumber * water_depth) * sin(wavenumber * x_rot - omega * time + phase); water_acceleration[2] = -omega * omega * amplitude * std::sinh(wavenumber * (z_pos + water_depth)) / - std::sinh(wavenumber * water_depth) * cos(wavenumber * x_pos - omega * time + phase); + std::sinh(wavenumber * water_depth) * cos(wavenumber * x_rot - omega * time + phase); } return water_acceleration; } @@ -128,13 +143,15 @@ Eigen::Vector3d GetWaterVelocityIrregular(const Eigen::Vector3d& position, const Eigen::VectorXd& wave_phases, const Eigen::VectorXd& wavenumbers, double water_depth, - double mwl) { + double mwl, + double wave_direction, + double wave_spread) { auto water_velocity = Eigen::Vector3d(0.0, 0.0, 0.0); for (size_t i = 0; i < freqs_hz.size(); ++i) { auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i]); auto omega = 2 * M_PI * freqs_hz[i]; - water_velocity += - GetWaterVelocity(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], water_depth, mwl); + water_velocity += GetWaterVelocity(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], + water_depth, mwl, wave_direction, wave_spread); } return water_velocity; } @@ -147,13 +164,15 @@ Eigen::Vector3d GetWaterAccelerationIrregular(const Eigen::Vector3d& position, const Eigen::VectorXd& wave_phases, const Eigen::VectorXd& wavenumbers, double water_depth, - double mwl) { + double mwl, + double wave_direction, + double wave_spread) { auto water_acceleration = Eigen::Vector3d(0.0, 0.0, 0.0); for (size_t i = 0; i < freqs_hz.size(); ++i) { - auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i]); + auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i] * wave_spread); auto omega = 2 * M_PI * freqs_hz[i]; water_acceleration += - GetWaterAcceleration(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], water_depth, mwl); + GetWaterAcceleration(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], water_depth, mwl, wave_direction, wave_spread); } return water_acceleration; } @@ -234,16 +253,16 @@ void RegularWave::AddH5Data(std::vector& reg_h5_data Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double time) { return GetWaterVelocity(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_, - wavenumber_, water_depth_, mwl_); + wavenumber_, water_depth_, mwl_, regular_wave_direction_, regular_wave_spread_); }; Eigen::Vector3d RegularWave::GetAcceleration(const Eigen::Vector3d& position, double time) { return GetWaterAcceleration(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_, - wavenumber_, water_depth_, mwl_); + wavenumber_, water_depth_, mwl_, regular_wave_direction_, regular_wave_spread_); }; double RegularWave::GetElevation(const Eigen::Vector3d& position, double time) { - return GetEta(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_, wavenumber_); + return GetEta(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_, wavenumber_, regular_wave_direction_); }; Eigen::VectorXd RegularWave::GetForceAtTime(double t) { @@ -267,7 +286,7 @@ double RegularWave::GetOmegaDelta() const { } double RegularWave::GetInterpolatedDirectionIndex(double dir_input) const { - const auto& dirs = wave_info_[0].wave_directions_list; + const auto& dirs = wave_info_[0].wave_direction_list; const size_t N = dirs.size(); if (N < 2) { @@ -343,7 +362,7 @@ Eigen::VectorXd ComputeWaveNumbers(const Eigen::VectorXd& omegas, return wavenumbers; } -std::vector> CreateFreeSurface3DPts(const std::vector& eta, +std::vector> CreateFreeSurface3DPts(const Eigen::VectorXd& eta, const Eigen::VectorXd& t_vec) { std::vector> surface(t_vec.size() * 2); @@ -447,7 +466,7 @@ std::vector IrregularWaves::GetSpectrum() { return spectrum_; } -std::vector IrregularWaves::GetFreeSurfaceElevation() { +Eigen::MatrixXd IrregularWaves::GetFreeSurfaceElevation() { return free_surface_elevation_sampled_; } @@ -462,23 +481,106 @@ void IrregularWaves::ReadEtaFromFile() { throw std::runtime_error("Unable to open file at: " + params_.eta_file_path_ + "."); } - std::string line; - double time, eta; + std::vector times; + std::vector> elevations_rows; + std::string line; while (std::getline(file, line)) { std::stringstream ss(line); - char delimiter; - if (!(ss >> time >> delimiter >> eta) || delimiter != ':') { - throw std::runtime_error("Could not parse line: " + line + "."); + double val; + std::vector values; + + while (ss >> val) { + values.push_back(val); + // Skip non-numeric delimiters automatically + if (ss.peek() == ':' || ss.peek() == ',') ss.ignore(); + } + + if (values.size() < 2) { + throw std::runtime_error("Each line must have at least time and one elevation value."); } - time_data_.push_back(time); - free_surface_elevation_sampled_.push_back(eta); + + // first column is time + times.push_back(values[0]); + + // rest are elevations for different wave directions (columns) + values.erase(values.begin()); + elevations_rows.push_back(values); } - std::cout << "Finished reading eta file." << std::endl; + + // Check that all rows have the same number of directions (columns) + size_t num_directions = elevations_rows[0].size(); + for (const auto& row : elevations_rows) { + if (row.size() != num_directions) { + throw std::runtime_error("Inconsistent number of wave directions per row."); + } + } + + // Fill Eigen::MatrixXd: rows = time steps, cols = directions + free_surface_time_sampled_.resize(times.size()); + for (size_t i = 0; i < times.size(); ++i) { + free_surface_time_sampled_[i] = times[i]; + } + + free_surface_elevation_sampled_.resize(times.size(), num_directions); + for (size_t i = 0; i < times.size(); ++i) { + for (size_t j = 0; j < num_directions; ++j) { + free_surface_elevation_sampled_(i, j) = elevations_rows[i][j]; + } + } + + std::cout << "Finished reading eta file with " << times.size() << " time steps and " << num_directions + << " wave directions." << std::endl; } -Eigen::MatrixXd IrregularWaves::GetExcitationIRF(int b) const { - return wave_info_[b].excitation_irf_matrix; + +Eigen::Tensor IrregularWaves::GetExcitationIRF(int b) const { + if (wave_info_[b].excitation_irf_matrix.dimension(1) == 1) { + return wave_info_[b].excitation_irf_matrix; + } + Eigen::VectorXd dir_index_des(params_.wave_direction_.size()); // n_dir_target + + for (int i = 0; i < params_.wave_direction_.size(); ++i) { + dir_index_des[i] = GetInterpolatedDirectionIndex(wave_info_[0].wave_direction_list, params_.wave_direction_[i]); + } + std::cout << dir_index_des << std::endl; + const auto& irf = wave_info_[b].excitation_irf_matrix; + const int ndof = irf.dimension(0); // 6 + const int ntime = irf.dimension(2); + const int ndir_target = dir_index_des.size(); + Eigen::Tensor interpolated_irf(ndof, ndir_target, ntime); + + // Loop over each target direction + for (int i = 0; i < ndir_target; ++i) { + double idx = dir_index_des[i]; + int i0 = static_cast(std::floor(idx)); + int i1 = i0 + 1; + + double alpha = idx - i0; // interpolation weight + double w0 = 1.0 - alpha; // weight of the lower index + double w1 = alpha; // weight of the upper index + + // Clamp to valid bounds + i0 = std::max(0, std::min(i0, static_cast(irf.dimension(1)) - 1)); + i1 = std::max(0, std::min(i1, static_cast(irf.dimension(1)) - 1)); + + if (w1 <= 1e-6) { + // No interpolation needed, just copy IRF at i0 + for (int dof = 0; dof < ndof; ++dof) { + for (int t = 0; t < ntime; ++t) { + interpolated_irf(dof, i, t) = irf(dof, i0, t); + } + } + } else { + // Interpolate between irf(:, i0, :) and irf(:, i1, :) + for (int dof = 0; dof < ndof; ++dof) { + for (int t = 0; t < ntime; ++t) { + interpolated_irf(dof, i, t) = w0 * irf(dof, i0, t) + w1 * irf(dof, i1, t); + } + } + } + } + return interpolated_irf; } void IrregularWaves::AddH5Data(std::vector& irreg_h5_data, @@ -495,7 +597,7 @@ Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, dou auto position_stretched = position; if (params_.wave_stretching_) { auto eta = GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, - wave_phases_, wavenumbers_); + wave_phases_, wavenumbers_, params_.wave_direction_[0], params_.wave_spread_[0]); // position relative to mean water level auto z_pos = position.z() - mwl_; // Wheeler stretching @@ -503,7 +605,8 @@ Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, dou } return GetWaterVelocityIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_, - spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_); + spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_, + params_.wave_direction_[0], params_.wave_spread_[0]); }; Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, double time) { @@ -511,7 +614,7 @@ Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, auto position_stretched = position; if (params_.wave_stretching_) { auto eta = GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, - wave_phases_, wavenumbers_); + wave_phases_, wavenumbers_, params_.wave_direction_[0], params_.wave_spread_[0]); // position relative to mean water level auto z_pos = position.z() - mwl_; // Wheeler stretching @@ -519,12 +622,13 @@ Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, } return GetWaterAccelerationIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_, - spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_); + spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_, + params_.wave_direction_[0], params_.wave_spread_[0]); }; double IrregularWaves::GetElevation(const Eigen::Vector3d& position, double time) { return GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, wave_phases_, - wavenumbers_); + wavenumbers_, params_.wave_direction_[0], params_.wave_spread_[0]); }; Eigen::VectorXd IrregularWaves::GetForceAtTime(double t) { @@ -543,7 +647,6 @@ Eigen::VectorXd IrregularWaves::GetForceAtTime(double t) { f[b_offset + dof] = f_dof; } } - return f; } @@ -559,27 +662,43 @@ void IrregularWaves::ResampleIRF(double dt) { // 1) Resample time auto t0 = time_array_old[0]; auto t1 = time_array_old[time_array_old.size() - 1]; + int n_time_new = static_cast(ceil((t1 - t0) / dt)); time_array = Eigen::VectorXd::LinSpaced(static_cast(ceil((t1 - t0) / dt)), t0, t1); // 2) Resample width CalculateWidthIRF(); - // 3) Resample values - assert(val_array.rows() == 6); - Eigen::MatrixXd vals_new(6, time_array.size()); + // 3) Resample val_array [DOF, N_dir, N_time_old] -> [DOF, N_dir, N_time_new] + int ndof = val_array.dimension(0); + int ndir = val_array.dimension(1); + int nt_old = val_array.dimension(2); + int nt_new = n_time_new; + + // 4) Resample values + assert(ndof == 6); + Eigen::Tensor val_array_new(ndof, ndir, nt_new); // We need to scale t to be [0,1] for spline use Eigen::VectorXd t_old_scaled = Eigen::VectorXd::LinSpaced(time_array_old.size(), 0, 1); Eigen::VectorXd t_new_scaled = Eigen::VectorXd::LinSpaced(time_array.size(), 0, 1); - // Use spline to get new values - Eigen::Spline spline = - Eigen::SplineFitting>::Interpolate(val_array, 3, t_old_scaled); - for (int i = 0; i < time_array.rows(); i++) { - vals_new.col(i) = spline(t_new_scaled[i]); - } + for (int i = 0; i < ndof; ++i) { + for (int j = 0; j < ndir; ++j) { + // Extract the [i,j,:] time series into a VectorXd + Eigen::VectorXd f_old(nt_old); + for (int k = 0; k < nt_old; ++k) + f_old(k) = val_array(i, j, k); - val_array = vals_new; + // Fit spline + auto spline = + Eigen::SplineFitting>::Interpolate(f_old.transpose(), 1, t_old_scaled); + + // Evaluate spline at new time points + for (int k = 0; k < nt_new; ++k) + val_array_new(i, j, k) = spline(t_new_scaled[k])(0); + } + } + val_array = val_array_new; } } @@ -621,6 +740,7 @@ Eigen::VectorXd IrregularWaves::SetSpectrumFrequencies(double start, double end, void IrregularWaves::CreateSpectrum() { // Define the frequency vector int nf; + int ndir = 1; // params_.wave_direction_.size(); if (params_.nfrequencies_ == 0) { // automatically calculate number of frequencies necessary so that timeseries does not repeat itself double df = 1.0 / params_.simulation_duration_; @@ -639,11 +759,13 @@ void IrregularWaves::CreateSpectrum() { spectral_widths_ = GetWidthArray(spectrum_frequencies_); // precompute random phases - wave_phases_ = Eigen::VectorXd(nf); + wave_phases_ = Eigen::MatrixXd(ndir, nf); std::mt19937 rng(params_.seed_); std::uniform_real_distribution dist(0.0, 2 * M_PI); - for (size_t i = 0; i < nf; ++i) { - wave_phases_[i] = dist(rng); + for (size_t i = 0; i < ndir; ++i) { + for (size_t j = 0; j < nf; ++j) { + wave_phases_(i, j) = dist(rng); + } } // precompute wavenumbers @@ -743,40 +865,48 @@ void IrregularWaves::CreateFreeSurfaceElevation() { // Calculate the free surface elevation // position assumed at (0.0, 0.0, 0.0) auto position = Eigen::Vector3d(0.0, 0.0, 0.0); + const int num_directions = static_cast(params_.wave_direction_.size()); + const int num_times = static_cast(free_surface_time_sampled_.size()); + // get timeseries - free_surface_elevation_sampled_ = - GetEtaIrregularTimeSeries(position, free_surface_time_sampled_, spectrum_frequencies_, spectral_densities_, - spectral_widths_, wave_phases_, wavenumbers_); - - // Apply ramp if ramp_duration is greater than 0 - if (params_.ramp_duration_ > 0.0) { - for (size_t i = 0; i < free_surface_time_sampled_.size(); ++i) { - if (free_surface_time_sampled_[i] < params_.ramp_duration_) { - if (free_surface_time_sampled_[i] <= 0.0) { - free_surface_elevation_sampled_[i] *= 0.0; - } else { - free_surface_elevation_sampled_[i] *= free_surface_time_sampled_[i] / params_.ramp_duration_; + free_surface_elevation_sampled_.resize(num_directions, num_times); + + for (int i = 0; i < params_.wave_direction_.size(); i++) { // loop over wave directions + free_surface_elevation_sampled_.row(i) = GetEtaIrregularTimeSeries( + position, free_surface_time_sampled_, spectrum_frequencies_, spectral_densities_, spectral_widths_, + wave_phases_.row(i).transpose(), wavenumbers_, params_.wave_direction_[i], params_.wave_spread_[i]); + + // Apply ramp if ramp_duration is greater than 0 + if (params_.ramp_duration_ > 0.0) { + for (size_t j = 0; j < free_surface_time_sampled_.size(); ++j) { + if (free_surface_time_sampled_[j] < params_.ramp_duration_) { + if (free_surface_time_sampled_[j] <= 0.0) { + free_surface_elevation_sampled_.row(i)[j] *= 0.0; + } else { + free_surface_elevation_sampled_.row(i)[j] *= free_surface_time_sampled_[j] / params_.ramp_duration_; + } } } } - } - + } + // Open a file stream for writing std::ofstream eta_output("eta.txt"); eta_output.precision(9); // Check if the file stream is open if (eta_output.is_open()) { + Eigen::RowVectorXd sum_row = free_surface_elevation_sampled_.colwise().sum(); // Write the spectral densities and their corresponding frequencies to the file for (size_t i = 0; i < free_surface_elevation_sampled_.size(); ++i) { - eta_output << free_surface_time_sampled_[i] << " : " << free_surface_elevation_sampled_[i] << std::endl; + eta_output << free_surface_time_sampled_[i] << " : " << sum_row(i) << std::endl; } // Close the file stream eta_output.close(); } else { std::cerr << "Unable to open file for writing." << std::endl; } - + std::cout << "Finished precalculating free surface elevation." << std::endl; } @@ -786,67 +916,84 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) { auto& irf_val_mat = ex_irf_sampled_[body]; auto& irf_width_array = ex_irf_width_sampled_[body]; + + const int num_directions = params_.wave_direction_.size(); // Number of wave directions + + if (irf_val_mat.dimension(1) == 1 && num_directions > 1) { + throw std::runtime_error( + "Excitation IRF data contains only 1 wave direction, but multiple wave directions were requested (" + + std::to_string(params_.wave_direction_.size()) + + ").\n" + "To simulate multi-directional waves, the excitation IRF must include response data for each wave " + "direction.\n"); + } // asumptions: irf_time_array in ascending order, free_surface_time_sampled_ in ascending order // get initial index 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; - if (t_tau <= tmin) { - idx = 0; - } else if (t_tau >= tmax) { - idx = free_surface_time_sampled_.size() - 2; - } else { - idx = get_lower_index(t_tau, free_surface_time_sampled_); - } - // loop for all irf time values - for (size_t j = 0; j < irf_time_array.size(); ++j) { - double tau = irf_time_array[j]; - double t_tau = time - tau; - if (tmin <= t_tau && t_tau <= tmax) { - // find next index if current lower bound > t_tau - while (free_surface_time_sampled_[idx] > t_tau) { - idx -= 1; - } + for (size_t k = 0; k < num_directions; ++k) { + double direction_weight = params_.wave_spread_[k]; - // free surface time values - auto t1 = free_surface_time_sampled_[idx]; - auto t2 = free_surface_time_sampled_[idx + 1]; - - // 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 - 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) + "."); - } + double t_tau = time - irf_time_array[0]; + int idx = 0; + if (t_tau <= tmin) { + idx = 0; + } else if (t_tau >= tmax) { + idx = free_surface_time_sampled_.size() - 2; + } else { + idx = get_lower_index(t_tau, free_surface_time_sampled_); + } - // add to excitation force - f_ex += irf_val_mat(dof, j) * eta_val * irf_width_array[j]; + // loop for all irf time values + //#pragma omp parallel for // Parallelize over directions + for (size_t j = 0; j < irf_time_array.size(); ++j) { + double tau = irf_time_array[j]; + double t_tau = time - tau; + if (tmin <= t_tau && t_tau <= tmax) { + // find next index if current lower bound > t_tau + while (free_surface_time_sampled_[idx] > t_tau) { + idx -= 1; + } - } else { - // throw error if trying to compute convolution after the maximum precomputed free elevation time - 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) + - "]). Excitation force ignored at this time step."); + // free surface time values + auto t1 = free_surface_time_sampled_[idx]; + auto t2 = free_surface_time_sampled_[idx + 1]; + + // 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_.row(k)[idx]; + auto eta2 = free_surface_elevation_sampled_.row(k)[idx + 1]; + // weights + auto w1 = (t2 - t_tau) / (t2 - t1); + auto w2 = 1.0 - w1; + // weighted value + eta_val = w1 * eta1 + w2 * eta2; + } else { + throw std::runtime_error("Excitation convolution: wrong tau value " + std::to_string(t_tau) + + " not between " + std::to_string(t1) + " and " + std::to_string(t2) + "."); + } + + // add to excitation force + f_ex += irf_val_mat(static_cast(dof), static_cast(k), + static_cast(j)) * eta_val * irf_width_array[j]; + + } else { + // throw error if trying to compute convolution after the maximum precomputed free elevation time + 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) + + "]). Excitation force ignored at this time step."); + } } } - return f_ex; } @@ -855,8 +1002,9 @@ void IrregularWaves::SetUpWaveMesh(std::string filename) { int num_timesteps = static_cast(ceil(params_.simulation_duration_ / params_.simulation_dt_)); Eigen::VectorXd time_index = Eigen::VectorXd::LinSpaced(num_timesteps + 1, 0, num_timesteps * params_.simulation_dt_); - std::vector> free_surface_3d_pts = - CreateFreeSurface3DPts(free_surface_elevation_sampled_, time_index); + + Eigen::VectorXd sum_over_rows = free_surface_elevation_sampled_.colwise().sum(); + std::vector> free_surface_3d_pts = CreateFreeSurface3DPts(sum_over_rows, time_index); std::vector> free_surface_triangles = CreateFreeSurfaceTriangles(time_index.size()); WriteFreeSurfaceMeshObj(free_surface_3d_pts, free_surface_triangles, mesh_file_name_); @@ -869,3 +1017,31 @@ std::string IrregularWaves::GetMeshFile() { Eigen::Vector3 IrregularWaves::GetWaveMeshVelocity() { return Eigen::Vector3d(1.0, 0, 0); } + +double WaveBase::GetInterpolatedDirectionIndex(const Eigen::VectorXd& dirs, double dir_input) const { + const size_t N = dirs.size(); + + if (N < 2) { + throw std::runtime_error("Not enough wave directions to interpolate."); + } + + // Wrap input to [0, 2pi) + double wrapped_dir = fmod(fmod(dir_input, CH_PI) + CH_PI, CH_PI); + // Linear search for lower bracket + for (size_t i = 0; i < N - 1; ++i) { + if (wrapped_dir >= dirs[i] && wrapped_dir <= dirs[i + 1]) { + double delta = dirs[i + 1] - dirs[i]; + double frac = (wrapped_dir - dirs[i]) / delta; + return static_cast(i) + frac; + } + } + + // Wraparound interpolation between last and first (e.g., 2pi and 0) + if (wrapped_dir >= dirs[dirs.size() - 1]) { + double delta = (CH_PI)-dirs[dirs.size() - 1] + dirs[0]; + double frac = (wrapped_dir - dirs[dirs.size() - 1]) / delta; + return static_cast(N - 1) + frac; + } + + throw std::runtime_error("Could not find interpolation interval for direction input."); +} From 5062b4a1acc7b4be7452aa61a51c222402299e34 Mon Sep 17 00:00:00 2001 From: Shabara Date: Mon, 21 Jul 2025 17:13:48 -0600 Subject: [PATCH 4/4] Adds multiple wave headings for the stretching --- src/wave_types.cpp | 58 +++++++++++++++++++--------------------------- 1 file changed, 24 insertions(+), 34 deletions(-) diff --git a/src/wave_types.cpp b/src/wave_types.cpp index e909cf37..ea3e559e 100644 --- a/src/wave_types.cpp +++ b/src/wave_types.cpp @@ -144,14 +144,16 @@ Eigen::Vector3d GetWaterVelocityIrregular(const Eigen::Vector3d& position, const Eigen::VectorXd& wavenumbers, double water_depth, double mwl, - double wave_direction, - double wave_spread) { + const std::vector& wave_direction, + const std::vector& wave_spread) { auto water_velocity = Eigen::Vector3d(0.0, 0.0, 0.0); - for (size_t i = 0; i < freqs_hz.size(); ++i) { - auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i]); - auto omega = 2 * M_PI * freqs_hz[i]; - water_velocity += GetWaterVelocity(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], - water_depth, mwl, wave_direction, wave_spread); + for (size_t j = 0; j < freqs_hz.size(); ++j) { + for (size_t i = 0; i < freqs_hz.size(); ++i) { + auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i] * wave_spread[j]); + auto omega = 2 * M_PI * freqs_hz[i]; + water_velocity += GetWaterVelocity(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], + water_depth, mwl, wave_direction[j], wave_spread[j]); + } } return water_velocity; } @@ -315,27 +317,6 @@ double RegularWave::GetInterpolatedDirectionIndex(double dir_input) const { throw std::runtime_error("Could not find interpolation interval for direction input."); } -/* -double RegularWave::GetExcitationMagInterp(int b, int i, int j, double freq_index_des) const { - double freq_interp_val = freq_index_des - floor(freq_index_des); - double excitationMagFloor = wave_info_[b].excitation_mag_matrix(i, j, (int)floor(freq_index_des)); - double excitationMagCeil = wave_info_[b].excitation_mag_matrix(i, j, (int)floor(freq_index_des) + 1); - double excitationMag = (freq_interp_val * (excitationMagCeil - excitationMagFloor)) + excitationMagFloor; - - return excitationMag; -} - -double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_index_des) const { - double freq_interp_val = freq_index_des - floor(freq_index_des); // look into c++ modf TODO - double excitationPhaseFloor = wave_info_[b].excitation_phase_matrix( - i, j, (int)floor(freq_index_des)); // TODO check if freq_index_des is >0, if so just cast instead of floor - double excitationPhaseCeil = wave_info_[b].excitation_phase_matrix(i, j, (int)floor(freq_index_des) + 1); - double excitationPhase = (freq_interp_val * (excitationPhaseCeil - excitationPhaseFloor)) + excitationPhaseFloor; - - return excitationPhase; -} -*/ - double RegularWave::GetExcitationMagInterp(int b, int i, double dir_index_des, double freq_index_des) const { const auto& mat = wave_info_[b].excitation_mag_matrix; return hydroc::BilinearInterp3D( @@ -606,24 +587,32 @@ Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, dou return GetWaterVelocityIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_, - params_.wave_direction_[0], params_.wave_spread_[0]); + params_.wave_direction_, params_.wave_spread_); }; Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, double time) { // apply wave stretching (if enabled) auto position_stretched = position; + double eta = 0.0; if (params_.wave_stretching_) { - auto eta = GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, - wave_phases_, wavenumbers_, params_.wave_direction_[0], params_.wave_spread_[0]); + for (size_t i = 0; i < params_.wave_direction_.size(); i++) { + eta += GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, + wave_phases_, wavenumbers_, params_.wave_direction_[i], params_.wave_spread_[i]); + } // position relative to mean water level auto z_pos = position.z() - mwl_; // Wheeler stretching position_stretched[2] = water_depth_ * (z_pos - eta) / (water_depth_ + eta); } - return GetWaterAccelerationIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_, - spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_, - params_.wave_direction_[0], params_.wave_spread_[0]); + Eigen::Vector3d TotWaterAccelerationIrregular = Eigen::Vector3d::Zero(); + for (size_t i = 0; i < params_.wave_direction_.size(); i++) { + TotWaterAccelerationIrregular += GetWaterAccelerationIrregular( + position_stretched, time, spectrum_frequencies_, spectral_densities_, + spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_, + params_.wave_direction_[i], params_.wave_spread_[i]); + } + return TotWaterAccelerationIrregular; }; double IrregularWaves::GetElevation(const Eigen::Vector3d& position, double time) { @@ -916,6 +905,7 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) { auto& irf_val_mat = ex_irf_sampled_[body]; auto& irf_width_array = ex_irf_width_sampled_[body]; + std::cout << "irf_val_mat = " << irf_val_mat.dimension(1) << std::endl; const int num_directions = params_.wave_direction_.size(); // Number of wave directions