Skip to content

Commit 2705a67

Browse files
tridelattridelat
committed
Add water_depth .h5 of sphere
1 parent 30302f4 commit 2705a67

3 files changed

Lines changed: 87 additions & 49 deletions

File tree

‎demos/sphere/hydroData/sphere.h5‎

2.01 KB
Binary file not shown.

‎include/hydroc/wave_types.h‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -19,12 +19,23 @@ Eigen::VectorXd JONSWAPSpectrumHz(Eigen::VectorXd& f,
1919
double gamma = 3.3,
2020
bool is_normalized = false);
2121

22-
std::vector<double> FreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
23-
const Eigen::VectorXd& spectral_densities,
24-
const Eigen::VectorXd& time_array,
25-
const Eigen::VectorXd& wave_phases,
26-
double position,
27-
double water_depth);
22+
double GetFreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
23+
const Eigen::VectorXd& spectral_densities,
24+
const Eigen::VectorXd& spectral_widths,
25+
const Eigen::VectorXd& wave_phases,
26+
const Eigen::VectorXd& wavenumbers,
27+
const Eigen::Vector3d& position,
28+
double time_value,
29+
double water_depth);
30+
31+
std::vector<double> GetFreeSurfaceElevationTimeSeries(const Eigen::VectorXd& freqs_hz,
32+
const Eigen::VectorXd& spectral_densities,
33+
const Eigen::VectorXd& spectral_widths,
34+
const Eigen::VectorXd& wave_phases,
35+
const Eigen::VectorXd& wavenumbers,
36+
const Eigen::Vector3d& position,
37+
const Eigen::VectorXd& time_array,
38+
double water_depth);
2839

2940
enum class WaveMode {
3041
/// @brief No waves
@@ -345,6 +356,8 @@ class IrregularWaves : public WaveBase {
345356
std::vector<Eigen::VectorXd> ex_irf_width_sampled_;
346357
Eigen::VectorXd spectrum_frequencies_;
347358
Eigen::VectorXd spectral_densities_;
359+
Eigen::VectorXd spectral_widths_;
360+
Eigen::VectorXd wavenumbers_;
348361
Eigen::VectorXd wave_phases_;
349362
std::string mesh_file_name_;
350363

‎src/wave_types.cpp‎

Lines changed: 68 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -86,12 +86,12 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
8686
return excitationPhase;
8787
}
8888

89-
std::vector<double> ComputeWaveNumbers(const std::vector<double>& omegas,
90-
double water_depth,
91-
double g = 9.81,
92-
double tolerance = 1e-6,
93-
int max_iterations = 100) {
94-
std::vector<double> wave_numbers(omegas.size());
89+
Eigen::VectorXd ComputeWaveNumbers(const Eigen::VectorXd& omegas,
90+
double water_depth,
91+
double g,
92+
double tolerance = 1e-6,
93+
int max_iterations = 100) {
94+
Eigen::VectorXd wavenumbers(omegas.size());
9595

9696
for (size_t i = 0; i < omegas.size(); ++i) {
9797
double omega = omegas[i];
@@ -112,33 +112,43 @@ std::vector<double> ComputeWaveNumbers(const std::vector<double>& omegas,
112112
iterations++;
113113
}
114114

115-
wave_numbers[i] = k;
115+
wavenumbers[i] = k;
116116
}
117117

118-
return wave_numbers;
118+
return wavenumbers;
119119
}
120120

121-
std::vector<double> FreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
122-
const Eigen::VectorXd& spectral_densities,
123-
const Eigen::VectorXd& wave_phases,
124-
const Eigen::VectorXd& time_index,
125-
double position,
126-
double water_depth) {
127-
int nf = freqs_hz.size();
128-
std::vector<double> omegas(nf);
121+
double GetFreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
122+
const Eigen::VectorXd& spectral_densities,
123+
const Eigen::VectorXd& spectral_widths,
124+
const Eigen::VectorXd& wave_phases,
125+
const Eigen::VectorXd& wavenumbers,
126+
const Eigen::Vector3d& position,
127+
double time_value,
128+
double water_depth) {
129+
// x position assuming wave direction along global X axis
130+
double x_pos = position.x();
131+
132+
double eta = 0.0;
129133
for (size_t i = 0; i < freqs_hz.size(); ++i) {
130-
omegas[i] = 2 * M_PI * freqs_hz[i];
134+
eta += std::sqrt(2 * spectral_densities[i] * spectral_widths[i]) *
135+
std::cos(wavenumbers[i] * x_pos - 2 * M_PI * freqs_hz[i] * time_value + wave_phases[i]);
131136
}
132-
double delta_f = freqs_hz(Eigen::last) / nf;
133-
134-
std::vector<double> wave_numbers = ComputeWaveNumbers(omegas, water_depth);
137+
return eta;
138+
}
135139

140+
std::vector<double> GetFreeSurfaceElevationTimeSeries(const Eigen::VectorXd& freqs_hz,
141+
const Eigen::VectorXd& spectral_densities,
142+
const Eigen::VectorXd& spectral_widths,
143+
const Eigen::VectorXd& wave_phases,
144+
const Eigen::VectorXd& wavenumbers,
145+
const Eigen::Vector3d& position,
146+
const Eigen::VectorXd& time_index,
147+
double water_depth) {
136148
std::vector<double> eta(time_index.size(), 0.0);
137-
for (size_t i = 0; i < nf; ++i) {
138-
for (size_t j = 0; j < time_index.size(); ++j) {
139-
eta[j] += std::sqrt(2 * spectral_densities[i] * delta_f) *
140-
std::cos(wave_numbers[i] * position - omegas[i] * time_index[j] + wave_phases[i]);
141-
}
149+
for (size_t j = 0; j < time_index.size(); ++j) {
150+
eta[j] += GetFreeSurfaceElevation(freqs_hz, spectral_densities, spectral_widths, wave_phases, wavenumbers,
151+
position, time_index[j], water_depth);
142152
}
143153

144154
return eta;
@@ -285,6 +295,7 @@ Eigen::MatrixXd IrregularWaves::GetExcitationIRF(int b) const {
285295
void IrregularWaves::AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data,
286296
HydroData::SimulationParameters& sim_data) {
287297
wave_info_ = irreg_h5_data;
298+
sim_data_ = sim_data;
288299

289300
InitializeIRFVectors();
290301
}
@@ -345,21 +356,25 @@ void IrregularWaves::ResampleIRF(double dt) {
345356
}
346357
}
347358

359+
Eigen::VectorXd GetWidthArray(const Eigen::VectorXd& input_array) {
360+
auto width_array = Eigen::VectorXd(input_array.size());
361+
for (int ii = 0; ii < width_array.size(); ii++) {
362+
width_array[ii] = 0.0;
363+
if (ii < input_array.size() - 1) {
364+
width_array[ii] += 0.5 * abs(input_array[ii + 1] - input_array[ii]);
365+
}
366+
if (ii > 0) {
367+
width_array[ii] += 0.5 * abs(input_array[ii] - input_array[ii - 1]);
368+
}
369+
}
370+
return width_array;
371+
}
372+
348373
void IrregularWaves::CalculateWidthIRF() {
349374
for (unsigned int b = 0; b < params_.num_bodies_; b++) {
350375
auto& time_array = ex_irf_time_sampled_[b];
351376
auto& width_array = ex_irf_width_sampled_[b];
352-
353-
width_array.resize(time_array.size());
354-
for (int ii = 0; ii < width_array.size(); ii++) {
355-
width_array[ii] = 0.0;
356-
if (ii < time_array.size() - 1) {
357-
width_array[ii] += 0.5 * abs(time_array[ii + 1] - time_array[ii]);
358-
}
359-
if (ii > 0) {
360-
width_array[ii] += 0.5 * abs(time_array[ii] - time_array[ii - 1]);
361-
}
362-
}
377+
width_array = GetWidthArray(time_array);
363378
}
364379
}
365380

@@ -389,17 +404,24 @@ void IrregularWaves::CreateSpectrum() {
389404
}
390405
spectrum_frequencies_ = Eigen::VectorXd::LinSpaced(nf, params_.frequency_min_, params_.frequency_max_);
391406

392-
// define random phases
407+
// Calculate the Pierson-Moskowitz Spectrum
408+
spectral_densities_ = JONSWAPSpectrumHz(spectrum_frequencies_, params_.wave_height_, params_.wave_period_,
409+
params_.peak_enhancement_factor_, params_.is_normalized_);
410+
411+
// precompute spectral widths
412+
spectral_widths_ = GetWidthArray(spectrum_frequencies_);
413+
414+
// precompute random phases
393415
wave_phases_ = Eigen::VectorXd(nf);
394416
std::mt19937 rng(params_.seed_);
395417
std::uniform_real_distribution<double> dist(0.0, 2 * M_PI);
396418
for (size_t i = 0; i < nf; ++i) {
397419
wave_phases_[i] = dist(rng);
398420
}
399421

400-
// Calculate the Pierson-Moskowitz Spectrum
401-
spectral_densities_ = JONSWAPSpectrumHz(spectrum_frequencies_, params_.wave_height_, params_.wave_period_,
402-
params_.peak_enhancement_factor_, params_.is_normalized_);
422+
// precompute wavenumbers
423+
auto omegas = 2 * M_PI * spectrum_frequencies_;
424+
wavenumbers_ = ComputeWaveNumbers(omegas, sim_data_.water_depth, sim_data_.g);
403425

404426
// Open a file stream for writing
405427
std::ofstream outputFile("spectral_densities.txt");
@@ -494,9 +516,12 @@ void IrregularWaves::CreateFreeSurfaceElevation() {
494516
<< std::endl;
495517

496518
// Calculate the free surface elevation
497-
double position = 0.0;
498-
free_surface_elevation_sampled_ = FreeSurfaceElevation(spectrum_frequencies_, spectral_densities_, wave_phases_,
499-
time_array, position, sim_data_.water_depth);
519+
// position assumed at (0.0, 0.0, 0.0)
520+
auto position = Eigen::Vector3d(0.0, 0.0, 0.0);
521+
// get timeseries
522+
free_surface_elevation_sampled_ =
523+
GetFreeSurfaceElevationTimeSeries(spectrum_frequencies_, spectral_densities_, spectral_widths_, wave_phases_,
524+
wavenumbers_, position, time_array, sim_data_.water_depth);
500525

501526
// Apply ramp if ramp_duration is greater than 0
502527
if (params_.ramp_duration_ > 0.0) {

0 commit comments

Comments
 (0)