fix irregular waves & export of data to .h5 file

Author: dav-og

include/hydroc/hydro_forces.h

@@ -218,6 +218,8 @@ class TestHydro {
      * @return 6N dimensional force for 6 DOF and N bodies in system (already Eigen type).
      */
     Eigen::VectorXd ComputeForceWaves();
+    // Expose the wave object (read-only) so callers can query inputs if needed
+    std::shared_ptr<WaveBase> GetWave() const { return user_waves_; }
 
     /**
      * @brief Fetches the RIRF value from the h5 file based on the provided indices.

include/hydroc/simulation_exporter.h

@@ -133,6 +133,16 @@ class SimulationExporter {
      */
     void Finalize();
 
+    /**
+     * @brief Write irregular wave input arrays under /inputs/simulation/waves/irregular.
+     *
+     * All arrays are 1D and written with attributes documenting units.
+     */
+    void WriteIrregularInputs(const std::vector<double>& frequencies_hz,
+                              const std::vector<double>& spectral_densities,
+                              const std::vector<double>& free_surface_time,
+                              const std::vector<double>& free_surface_eta);
+
 
 
     // Runtime metadata (written into /meta/run at Finalize)

include/hydroc/wave_types.h

@@ -297,9 +297,14 @@ class IrregularWaves : public WaveBase {
     void Initialize() override {}
 
     void CreateSpectrum();
+    // Returns spectral density values S(f) (one per frequency bin)
     std::vector<double> GetSpectrum();
+    // Returns precomputed free-surface elevation eta(t) samples
     std::vector<double> GetFreeSurfaceElevation();
-    std::vector<double> GetEtaTimeData();
+    // Returns time samples corresponding to eta(t)
+    std::vector<double> GetFreeSurfaceTime() const;
+    // Returns the frequency bins (Hz)
+    std::vector<double> GetFrequenciesHz() const;
 
     Eigen::VectorXd GetForceAtTime(double t) override;
 

src/hydro_forces.cpp

@@ -236,6 +236,9 @@ TestHydro::TestHydro(std::vector<std::shared_ptr<ChBody>> user_bodies,
     // Set up hydro inputs
     user_waves_ = waves;
     AddWaves(user_waves_);
+
+    // If irregular waves are active, publish spectrum and free-surface inputs into HDF5 if an exporter is present.
+    // The exporter is managed by the runner; here we only expose getters via the waves object.
 }
 
 void TestHydro::AddWaves(std::shared_ptr<WaveBase> waves) {

src/hydrochrono_runner/run_hydrochrono_from_yaml.cpp

@@ -5,6 +5,7 @@
 #include "../hydro_yaml_parser.h"
 #include <hydroc/hydro_forces.h>
 #include <hydroc/simulation_exporter.h>
+#include <hydroc/wave_types.h>
 
 #include <chrono_parsers/yaml/ChParserMbsYAML.h>
 #include <chrono/physics/ChSystem.h>
@@ -450,8 +451,10 @@ int RunHydroChronoFromYAML(int argc, char* argv[]) {
 
                     // Setup hydrodynamic forces
                     hydroc::debug::LogDebug("Initializing TestHydro...");
-                    // Provide a neutral horizon (Chrono governs runtime via YAML/UI)
-                    test_hydro = SetupHydroFromYAML(hydro_data, bodies, loop_dt, /*time_end_hint*/ 0.0, 0.0);
+                    // Provide simulation horizon from YAML end_time (important for irregular waves spectrum)
+                    double sim_duration_hint = 0.0;
+                    TryFindYamlDouble(sim_file, "end_time", sim_duration_hint);
+                    test_hydro = SetupHydroFromYAML(hydro_data, bodies, loop_dt, sim_duration_hint, 0.0);
                     hydroc::debug::LogDebug("Hydrodynamic forces initialized successfully");
 
                     // Inform location for diagnostics CSVs: write to hydro file directory
@@ -657,9 +660,23 @@ int RunHydroChronoFromYAML(int argc, char* argv[]) {
                 exporter = std::make_unique<hydroc::SimulationExporter>(exp_opts);
 
                 // Write static info and model before stepping
-                exporter->WriteSimulationInfo(system.get(), std::string(""), std::filesystem::path(model_file).filename().generic_string(), loop_dt, /*duration_seconds*/ 0.0);
+                double duration_hint = 0.0; TryFindYamlDouble(sim_file, "end_time", duration_hint);
+                exporter->WriteSimulationInfo(system.get(), std::string(""), std::filesystem::path(model_file).filename().generic_string(), loop_dt, duration_hint);
                 exporter->WriteModel(system.get());
                 exporter->BeginResults(system.get(), /*expected_steps*/ 0);
+
+                // If irregular waves are configured, persist spectrum and eta(t) inputs to HDF5
+                if (test_hydro) {
+                    auto wave_ptr = test_hydro->GetWave();
+                    if (wave_ptr && wave_ptr->GetWaveMode() == WaveMode::irregular) {
+                        auto irreg = std::static_pointer_cast<IrregularWaves>(wave_ptr);
+                        std::vector<double> f = irreg->GetFrequenciesHz();
+                        std::vector<double> S = irreg->GetSpectrum();
+                        std::vector<double> tvec = irreg->GetFreeSurfaceTime();
+                        std::vector<double> eta = irreg->GetFreeSurfaceElevation();
+                        exporter->WriteIrregularInputs(f, S, tvec, eta);
+                    }
+                }
             } catch (const std::exception& e) {
                 hydroc::cli::LogWarning(std::string("HDF5 exporter disabled: ") + e.what());
                 exporter.reset();

src/setup_hydro_from_yaml.cpp

@@ -56,7 +56,8 @@ std::shared_ptr<WaveBase> CreateWaveFromSettings(const WaveSettings& wave_settin
         params.ramp_duration_ = ramp_duration;
         params.wave_height_ = wave_settings.height;
         params.wave_period_ = wave_settings.period;
-        params.seed_ = 1;  // Default seed, could be made configurable
+        // Use YAML-provided seed if available; fall back to a default deterministic seed
+        params.seed_ = (wave_settings.seed > 0 ? wave_settings.seed : 1);
 
         auto irregular_wave = std::make_shared<IrregularWaves>(params);
 

src/simulation_exporter.cpp

@@ -5,6 +5,7 @@
 
 #include <hydroc/simulation_exporter.h>
 #include <hydroc/h5_writer.h>
+#include <hydroc/wave_types.h>
 #include <hydroc/version.h>
 #include <hydroc/logging.h>
 
@@ -81,6 +82,8 @@ struct SimulationExporter::Impl {
     H5Writer::Group g_results_model_tsdas;
     H5Writer::Group g_results_model_rsdas;
     H5Writer::Group g_results_model_joints;
+    // irregular waves datasets (inputs)
+    H5Writer::Group g_inputs_waves_irregular;
     // meta subtree
     H5Writer::Group g_meta;
 
@@ -184,6 +187,7 @@ struct SimulationExporter::Impl {
         g_inputs_sim_time = writer.RequireGroup("/inputs/simulation/time");
         g_inputs_sim_env = writer.RequireGroup("/inputs/simulation/environment");
         g_inputs_sim_waves = writer.RequireGroup("/inputs/simulation/waves");
+        g_inputs_waves_irregular = writer.RequireGroup("/inputs/simulation/waves/irregular");
         // results
         g_results = writer.RequireGroup("/results");
         g_results_model = writer.RequireGroup("/results/model");
@@ -347,6 +351,7 @@ void SimulationExporter::WriteSimulationInfo(chrono::ChSystem* system,
         g_waves.WriteAttribute("Hs", impl_->options.scenario_Hs);
         g_waves.WriteAttribute("Tp", impl_->options.scenario_Tp);
         if (impl_->options.scenario_seed >= 0) g_waves.WriteAttribute("seed", static_cast<double>(impl_->options.scenario_seed));
+        // Group already reserved in ctor; datasets will be written later if available
     }
 }
 
@@ -357,6 +362,36 @@ void SimulationExporter::WriteInitialConditions(chrono::ChSystem* ,
                                                 double ) {
     // Deprecated in schema v0.3: handled in WriteSimulationInfo under inputs/simulation/*
 }
+void SimulationExporter::WriteIrregularInputs(const std::vector<double>& frequencies_hz,
+                                             const std::vector<double>& spectral_densities,
+                                             const std::vector<double>& free_surface_time,
+                                             const std::vector<double>& free_surface_eta) {
+    auto g = impl_->g_inputs_waves_irregular;
+    // Write datasets if provided (skip empty vectors)
+    if (!frequencies_hz.empty()) {
+        std::array<hsize_t,1> d1 = {static_cast<hsize_t>(frequencies_hz.size())};
+        g.WriteDataset("frequencies_hz", frequencies_hz, d1);
+        g.WriteAttribute("frequencies_hz.units", std::string("Hz"));
+    }
+    if (!spectral_densities.empty()) {
+        std::array<hsize_t,1> d1 = {static_cast<hsize_t>(spectral_densities.size())};
+        g.WriteDataset("spectral_densities", spectral_densities, d1);
+        g.WriteAttribute("spectral_densities.units", std::string("m^2/Hz"));
+        g.WriteAttribute("spectral_densities.convention", std::string("JONSWAP (if gamma>1), else PM"));
+    }
+    if (!free_surface_time.empty()) {
+        std::array<hsize_t,1> d1 = {static_cast<hsize_t>(free_surface_time.size())};
+        g.WriteDataset("free_surface_time", free_surface_time, d1);
+        g.WriteAttribute("free_surface_time.units", std::string("s"));
+    }
+    if (!free_surface_eta.empty()) {
+        std::array<hsize_t,1> d1 = {static_cast<hsize_t>(free_surface_eta.size())};
+        g.WriteDataset("free_surface_eta", free_surface_eta, d1);
+        g.WriteAttribute("free_surface_eta.units", std::string("m"));
+        g.WriteAttribute("free_surface_eta.location", std::string("x=0,y=0,z=0 (assumed)"));
+    }
+}
+
 
 void SimulationExporter::WriteModel(chrono::ChSystem* system) {
     if (system == nullptr) {

src/wave_types.cpp

@@ -470,8 +470,11 @@ std::vector<double> IrregularWaves::GetFreeSurfaceElevation() {
     return free_surface_elevation_sampled_;
 }
 
-std::vector<double> IrregularWaves::GetEtaTimeData() {
-    return time_data_;
+std::vector<double> IrregularWaves::GetFreeSurfaceTime() const { return free_surface_time_sampled_; }
+std::vector<double> IrregularWaves::GetFrequenciesHz() const {
+    std::vector<double> out(spectrum_frequencies_.size());
+    for (int i = 0; i < spectrum_frequencies_.size(); ++i) out[i] = spectrum_frequencies_[i];
+    return out;
 }
 
 void IrregularWaves::ReadEtaFromFile() {
@@ -669,22 +672,7 @@ void IrregularWaves::CreateSpectrum() {
     auto omegas  = 2 * M_PI * spectrum_frequencies_;
     wavenumbers_ = ComputeWaveNumbers(omegas, water_depth_, g_);
 
-    // Open a file stream for writing
-    std::ofstream outputFile("spectral_densities.txt");
-    outputFile.precision(9);
-
-    // Check if the file stream is open
-    if (outputFile.is_open()) {
-        // Write the spectral densities and their corresponding frequencies to the file
-        for (size_t i = 0; i < spectral_densities_.size(); ++i) {
-            outputFile << spectrum_frequencies_[i] << " : " << spectral_densities_[i] << std::endl;
-        }
-
-        // Close the file stream
-        outputFile.close();
-    } else {
-        hydroc::cli::LogError("Unable to open file for writing.");
-    }
+    // Removed text file dump; this data will be written into HDF5 by the exporter when available.
 }
 
 // TODO put spectrum functions in a new namespace (when we have more options?)
@@ -780,21 +768,7 @@ void IrregularWaves::CreateFreeSurfaceElevation() {
         }
     }
 
-    // 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()) {
-        // 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;
-        }
-        // Close the file stream
-        eta_output.close();
-    } else {
-        hydroc::cli::LogError("Unable to open file for writing.");
-    }
+    // Removed text file dump; this data will be written into HDF5 by the exporter when available.
 
     hydroc::debug::LogDebug("Finished precalculating free surface elevation.");
 }
@@ -836,9 +810,9 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
             // get free surface elevation
             double eta_val;
             if (t_tau == t1) {
-                eta_val = free_surface_time_sampled_[idx];
+                eta_val = free_surface_elevation_sampled_[idx];
             } else if (t_tau == t2) {
-                eta_val = free_surface_time_sampled_[idx + 1];
+                eta_val = free_surface_elevation_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];