Skip to content

Commit 654e942

Browse files
tridelattridelat
committed
Simplify IRF resampling for excitation force
1 parent cf20bb1 commit 654e942

2 files changed

Lines changed: 69 additions & 79 deletions

File tree

‎include/hydroc/wave_types.h‎

Lines changed: 8 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -339,8 +339,9 @@ class IrregularWaves : public WaveBase {
339339
// const WaveMode mode_ = WaveMode::irregular;
340340
std::vector<HydroData::IrregularWaveInfo> wave_info_;
341341
HydroData::SimulationParameters sim_data_;
342-
std::vector<Eigen::MatrixXd> ex_irf_resampled_;
343-
std::vector<Eigen::VectorXd> ex_irf_time_resampled_;
342+
std::vector<Eigen::MatrixXd> ex_irf_sampled_;
343+
std::vector<Eigen::VectorXd> ex_irf_time_sampled_;
344+
std::vector<Eigen::VectorXd> ex_irf_width_sampled_;
344345
Eigen::VectorXd spectrum_frequencies_;
345346
Eigen::VectorXd spectral_densities_;
346347
std::string mesh_file_name_;
@@ -352,30 +353,15 @@ class IrregularWaves : public WaveBase {
352353

353354
Eigen::MatrixXd GetExcitationIRF(int b) const;
354355

355-
/**
356-
* @brief resamples irf time vector from old time vector and new timestep.
357-
*
358-
* Creates a resized time vector that starts and ends at the same values as t_old and has timestep t_new.
359-
*
360-
* @param t_old reference to old time vector from h5 file, the first and last element are transfered to the first
361-
* and last element of t_new (return val)
362-
* @param dt_new the time step to use in resampled vector, typically the timestep of chrono simulation
356+
/** @brief Resamples IRF time, widths, and values.
363357
*
364-
* @return newly sized vector that looks like \
365-
* (t_old[0], t_old[0] + dt_new, t_old[0] + 2*dt_new, ... , t_old[t_old.size()-1])
358+
* @param dt Time step value to resample
366359
*/
367-
Eigen::VectorXd ResampleTime(const Eigen::VectorXd& t_old, const double dt_new);
360+
void ResampleIRF(double dt);
368361

369-
/**
370-
* @brief Creates a resized values vector to interpolate values for a new timestep.
371-
*
372-
* @param t_old reference to old time vector from h5 file
373-
* @param vals_old the original values corresponding to the times in t_old from h5 file
374-
* @param t_new resampled times (return value from ResampleTime()) to use in interpolation
375-
*
376-
* @return matrix for interpolated vals_old over t_new
362+
/** @brief Calculates width (used for excitation convolution).
377363
*/
378-
Eigen::MatrixXd ResampleVals(const Eigen::VectorXd& t_old, Eigen::MatrixXd& vals_old, const Eigen::VectorXd& t_new);
364+
void IrregularWaves::CalculateWidthIRF();
379365

380366
/**
381367
* @brief Calculates the component of force from Convolution integral for specified body, dof, time.

‎src/wave_types.cpp‎

Lines changed: 61 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -482,8 +482,9 @@ IrregularWaves::IrregularWaves(const IrregularWaveParams& params)
482482
ramp_duration_(params.ramp_duration_),
483483
seed_(params.seed_) {
484484
std::cout << "Creating IrregularWaves object..." << std::endl;
485-
ex_irf_resampled_.resize(num_bodies_);
486-
ex_irf_time_resampled_.resize(num_bodies_);
485+
ex_irf_sampled_.resize(num_bodies_);
486+
ex_irf_time_sampled_.resize(num_bodies_);
487+
ex_irf_width_sampled_.resize(num_bodies_);
487488
std::cout << "eta_file_path (1)" << eta_file_path_ << std::endl;
488489
if (!eta_file_path_.empty()) {
489490
std::cout << "Reading eta.txt..." << std::endl;
@@ -521,14 +522,14 @@ void IrregularWaves::InitializeIRFVectors() {
521522
std::vector<Eigen::VectorXd> ex_irf_time_old(num_bodies_);
522523

523524
for (unsigned int b = 0; b < num_bodies_; b++) {
524-
ex_irf_old[b] = GetExcitationIRF(b);
525-
ex_irf_time_old[b] = wave_info_[b].excitation_irf_time;
525+
ex_irf_sampled_[b] = GetExcitationIRF(b);
526+
ex_irf_time_sampled_[b] = wave_info_[b].excitation_irf_time;
527+
CalculateWidthIRF();
526528
}
527529

528530
// Resample excitation IRF time series
529-
for (unsigned int b = 0; b < num_bodies_; b++) {
530-
ex_irf_time_resampled_[b] = ResampleTime(ex_irf_time_old[b], simulation_dt_);
531-
ex_irf_resampled_[b] = ResampleVals(ex_irf_time_old[b], ex_irf_old[b], ex_irf_time_resampled_[b]);
531+
if (simulation_dt_ > 0.0) {
532+
ResampleIRF(simulation_dt_);
532533
}
533534
}
534535

@@ -616,44 +617,56 @@ Eigen::VectorXd IrregularWaves::GetForceAtTime(double t) {
616617
return f;
617618
}
618619

619-
Eigen::MatrixXd IrregularWaves::ResampleVals(const Eigen::VectorXd& t_old,
620-
Eigen::MatrixXd& vals_old,
621-
const Eigen::VectorXd& t_new) {
622-
assert(vals_old.rows() == 6);
623-
624-
Eigen::MatrixXd vals_new(6, t_new.size());
625-
// we need to scale t to be [0,1] for spline use
626-
// TODO: change this to accomodate variable dt instead of const dt
627-
Eigen::VectorXd t_old_scaled = Eigen::VectorXd::LinSpaced(t_old.size(), 0, 1);
628-
Eigen::VectorXd t_new_scaled = Eigen::VectorXd::LinSpaced(t_new.size(), 0, 1);
629-
630-
Eigen::Spline<double, 6> spline =
631-
Eigen::SplineFitting<Eigen::Spline<double, 6>>::Interpolate(vals_old, 3, t_old_scaled);
632-
for (int i = 0; i < t_new.rows(); i++) {
633-
vals_new.col(i) = spline(t_new_scaled[i]);
620+
void IrregularWaves::ResampleIRF(double dt) {
621+
for (unsigned int b = 0; b < num_bodies_; b++) {
622+
// 0) access arrays to modify at body index
623+
auto& time_array = ex_irf_time_sampled_[b];
624+
auto& width_array = ex_irf_width_sampled_[b];
625+
auto& val_array = ex_irf_sampled_[b];
626+
// copy time array
627+
auto time_array_old = time_array;
628+
629+
// 1) resample time
630+
auto t0 = time_array_old[0];
631+
auto t1 = time_array_old[time_array_old.size() - 1];
632+
time_array = Eigen::VectorXd::LinSpaced(static_cast<int>(ceil((t1 - t0) / dt)), t0, t1);
633+
634+
// 2) resample width
635+
CalculateWidthIRF();
636+
637+
// 3) resample values
638+
assert(val_array.rows() == 6);
639+
Eigen::MatrixXd vals_new(6, time_array.size());
640+
// we need to scale t to be [0,1] for spline use
641+
Eigen::VectorXd t_old_scaled = Eigen::VectorXd::LinSpaced(time_array_old.size(), 0, 1);
642+
Eigen::VectorXd t_new_scaled = Eigen::VectorXd::LinSpaced(time_array.size(), 0, 1);
643+
// use spline to get new values
644+
Eigen::Spline<double, 6> spline =
645+
Eigen::SplineFitting<Eigen::Spline<double, 6>>::Interpolate(val_array, 3, t_old_scaled);
646+
for (int i = 0; i < time_array.rows(); i++) {
647+
vals_new.col(i) = spline(t_new_scaled[i]);
648+
}
649+
val_array = vals_new;
634650
}
635-
636-
// std::cout << "Old time vector:\n" << t_old << std::endl;
637-
// std::cout << "Old values:\n" << vals_old << std::endl;
638-
// std::cout << "New time vector:\n" << t_new << std::endl;
639-
// std::cout << "New values:\n" << vals_new << std::endl; // assuming new_vals is the output of the function
640-
641-
return vals_new;
642651
}
643652

644-
Eigen::VectorXd IrregularWaves::ResampleTime(const Eigen::VectorXd& t_old, const double dt_new) {
645-
double dt_old = t_old[1] - t_old[0];
646-
int size_new = static_cast<int>(ceil(t_old.size() * dt_old / dt_new));
647-
double t_initial = t_old[0];
648-
double t_final = t_old[t_old.size() - 1];
649-
// t_0 and t_final should be the same in old/new versions, use new time step, which means a different number of
650-
// timesteps
651-
Eigen::VectorXd t_new = Eigen::VectorXd::LinSpaced(size_new, t_initial, t_final);
652-
653-
int N_old = t_old.size() - 1;
654-
int N_new = t_new.size() - 1;
655-
656-
return t_new;
653+
void IrregularWaves::CalculateWidthIRF() {
654+
for (unsigned int b = 0; b < num_bodies_; b++) {
655+
auto& time_array = ex_irf_time_sampled_[b];
656+
auto& width_array = ex_irf_width_sampled_[b];
657+
658+
// 2) calculate width
659+
width_array.resize(time_array.size());
660+
for (int ii = 0; ii < width_array.size(); ii++) {
661+
width_array[ii] = 0.0;
662+
if (ii < time_array.size() - 1) {
663+
width_array[ii] += 0.5 * abs(time_array[ii + 1] - time_array[ii]);
664+
}
665+
if (ii > 0) {
666+
width_array[ii] += 0.5 * abs(time_array[ii] - time_array[ii - 1]);
667+
}
668+
}
669+
}
657670
}
658671

659672
Eigen::VectorXd IrregularWaves::SetSpectrumFrequencies(double start, double end, int num_points) {
@@ -774,9 +787,10 @@ void IrregularWaves::CreateFreeSurfaceElevation() {
774787
}
775788

776789
double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
777-
double f_ex = 0.0;
778-
auto& irf_time_array = ex_irf_time_resampled_[body];
779-
auto& irf_val_mat = ex_irf_resampled_[body];
790+
double f_ex = 0.0;
791+
auto& irf_time_array = ex_irf_time_sampled_[body];
792+
auto& irf_val_mat = ex_irf_sampled_[body];
793+
auto& irf_width_array = ex_irf_width_sampled_[body];
780794

781795
// asumptions: irf_time_array in ascending order, free_surface_time_sampled_ in ascending order
782796
// get initial index
@@ -820,18 +834,8 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
820834
// weighted value
821835
auto eta_val = w1 * eta1 + w2 * eta2;
822836

823-
// calculate width
824-
double width = 0;
825-
if (j > 0) {
826-
width += 0.5 * (irf_time_array[j] - irf_time_array[j - 1]);
827-
}
828-
if (j < irf_time_array.size() - 1) {
829-
width += 0.5 * (irf_time_array[j + 1] - irf_time_array[j]);
830-
}
831-
832-
double ex_irf_val = irf_val_mat(dof, j);
833-
834-
f_ex += ex_irf_val * eta_val * width;
837+
// add to excitation force
838+
f_ex += irf_val_mat(dof, j) * eta_val * irf_width_array[j];
835839
}
836840
}
837841

0 commit comments

Comments
 (0)