Merge pull request #75 from Project-SEA-Stack/refactor/simd-fix-parallel-convolution-savgol

Refactor: SIMD/vectorization fixes, parallel convolution integral (~5× speedup), Savitzky–Golay smoothing, and profiling options

Author: dav-og

app/run_hydrochrono.cpp

@@ -83,6 +83,7 @@ struct CLIArgs {
     bool h5_verbose = false;            // NEW: HDF5 verbose diagnostics
     std::string h5_tag;                 // NEW: Optional tag appended to filename
     bool fail_fast = false;             // NEW: Stop on first failure during sweep
+    bool profile = false;               // NEW: Enable runtime profiling summary
 };
 
 static CLIArgs ParseArguments(int argc, char* argv[]) {
@@ -136,6 +137,8 @@ static CLIArgs ParseArguments(int argc, char* argv[]) {
             }
         } else if (arg == "--fail-fast") {
             args.fail_fast = true;
+        } else if (arg == "--profile") {
+            args.profile = true;
         } else if (arg.substr(0, 1) != "-") {
             // This is a positional argument (input directory)
             if (args.input_directory.empty()) {
@@ -286,6 +289,10 @@ int main(int argc, char* argv[]) {
         runner_args.push_back("--trace");
     }
 
+    if (args.profile) {
+        runner_args.push_back("--profile");
+    }
+    
     if (!args.model_file.empty()) {
         runner_args.push_back("--model_file");
         runner_args.push_back(args.model_file);

build.ps1

@@ -603,11 +603,20 @@ function Get-BuildArguments {
     $irrlicht   = ($Config.IrrlichtDir -replace '\\','/')
     $irrlichtDllWin = Join-Path $Config.IrrlichtDir "bin\Win64-VisualStudio\Irrlicht.dll"
     $irrlichtDll = ($irrlichtDllWin -replace '\\','/')
+    $irrlichtFlag = if ($Config.IrrlichtDir -and (Test-Path $Config.IrrlichtDir)) { "ON" } else { "OFF" }
 
     # Deterministic prefix path so our roots win over ambient environment
     $prefixParts = @()
     foreach ($p in @($Config.ChronoDir, $Config.Hdf5Dir, $Config.EigenDir, $Config.IrrlichtDir, $Config.PythonRoot)) {
-        if ($p -and (Test-Path $p)) { $prefixParts += ($p -replace '\\','/') }
+        if ($p) {
+            $isValid = $false
+            try { $null = [System.IO.Path]::GetFullPath($p); $isValid = $true } catch { $isValid = $false }
+            if ($isValid) {
+                try {
+                    if (Test-Path -LiteralPath $p) { $prefixParts += ($p -replace '\\','/') }
+                } catch { }
+            }
+        }
     }
     $prefixPath = ($prefixParts -join ';')
     $testsFlag = if ($NoTests) { "OFF" } else { "ON" }
@@ -623,6 +632,7 @@ function Get-BuildArguments {
         "-DEIGEN3_INCLUDE_DIR=`"$eigenDir`"",
         "-DIrrlicht_ROOT=`"$irrlicht`"",
         "-DIRRLICHT_DLL_PATH=`"$irrlichtDll`"",
+        "-DHYDROCHRONO_ENABLE_IRRLICHT=`"$irrlichtFlag`"",
         "-DCMAKE_PREFIX_PATH=`"$prefixPath`"",
         "-DCMAKE_FIND_PACKAGE_PREFER_CONFIG=ON",
         "-DCMAKE_FIND_USE_PACKAGE_REGISTRY=OFF",

include/hydroc/hydro_forces.h

@@ -149,6 +149,16 @@ class ForceFunc6d {
 
 class ChLoadAddedMass;
 
+// Lightweight hydrodynamics profiling stats
+struct HydroProfileStats {
+    double hydrostatics_seconds = 0.0;
+    double radiation_seconds    = 0.0;
+    double waves_seconds        = 0.0;
+    int hydrostatics_calls      = 0;
+    int radiation_calls         = 0;
+    int waves_calls             = 0;
+};
+
 // TODO: Rename TestHydro for clarity, perhaps to HydroForces?
 // TODO: Split TestHydro class from its helper classes for clearer code structure.
 class TestHydro {
@@ -220,6 +230,42 @@ class TestHydro {
      */
     double GetRIRFval(int row, int col, int st);
 
+    // Convolution mode selection
+    enum class RadiationConvolutionMode {
+        Baseline,
+        TaperedDirect
+    };
+
+    /**
+     * @brief Set the radiation convolution mode. Default is Baseline.
+     */
+    void SetRadiationConvolutionMode(RadiationConvolutionMode mode) { convolution_mode_ = mode; }
+
+        struct TaperedDirectOptions {
+            // smoothing: "sg" (Savitzky–Golay) or "moving_average"
+            std::string smoothing = "sg";
+            int window_length = 5;                 // odd, >= 3
+            
+            // RIRF truncation
+            double rirf_end_time = -1.0;           // end RIRF at this time (seconds), -1.0 = use full length
+            
+            // Simple taper control - sensible defaults for improved stability
+            double taper_start_percent = 0.8;      // start taper at 80% (taper last 20%)
+            double taper_end_percent = 1.0;        // end taper at 100% of total time series  
+            double taper_final_amplitude = 0.0;    // final amplitude as fraction of original (0.0 = zero, 1.0 = no change)
+            bool export_plot_csv = false;          // dump before/after CSV summaries (false by default)
+        };
+
+    /**
+     * @brief Set options for TaperedDirect preprocessing.
+     */
+    void SetTaperedDirectOptions(const TaperedDirectOptions& opts) { tapered_opts_ = opts; }
+
+    /**
+     * @brief Set the directory where diagnostics (e.g., CSVs) should be written.
+     */
+    void SetDiagnosticsOutputDirectory(const std::string& dir) { diagnostics_output_dir_ = dir; }
+
     /**
      * @brief Calculates or retrieves the total force on a specific body in a particular degree of freedom.
      *
@@ -233,6 +279,9 @@ class TestHydro {
      */
     double CoordinateFuncForBody(int b, int i);
 
+    // Hydrodynamics profiling accessors
+    HydroProfileStats GetProfileStats() const { return profile_stats_; }
+
   private:
     // Class properties related to the body and hydrodynamics
     std::vector<std::shared_ptr<ChBody>> bodies_;
@@ -281,6 +330,18 @@ class TestHydro {
      * @param index The DOF index, ranging from [0,1,...,5].
      */
     double SetVelHistory(double val, int step, int b_num, int index);
+
+    // Hydrodynamics profiling data (accumulated over run)
+    HydroProfileStats profile_stats_;
+
+    // Convolution kernel preprocessing (optional)
+    RadiationConvolutionMode convolution_mode_ = RadiationConvolutionMode::Baseline;
+    bool rirf_processed_ready_ = false;
+    std::vector<Eigen::Tensor<double, 3>> rirf_processed_; // per body [dof x col x step]
+    TaperedDirectOptions tapered_opts_;
+    std::string diagnostics_output_dir_;
+
+    void EnsureProcessedRIRF();
 };
 
 #endif