remove dead comparison helpers and document TestHydro as HydroSystem adapter

Remove dead internal helpers EvaluateHydroSystem() and CompareWithHydroSystem() from TestHydro.

Update header and implementation comments in hydro_forces to reflect the new architecture: TestHydro is a thin adapter over HydroSystem/ChronoHydroCoupler, and the Create*Component helpers are the single sources of truth for component configuration.

Legacy ComputeForce* methods and compare_mode_ are retained for API compatibility and clearly marked as legacy.

Author: dav-og

include/hydroc/hydro_forces.h

@@ -7,16 +7,29 @@
  * @brief Header file of TestHydro main class and helper classes
  * ComponentFunc and ForceFunc6d.
  *
- * OVERVIEW:
- * This header defines the public interface for hydrodynamic force computation
- * in multibody marine systems. It provides Chrono-compatible force callbacks
- * that integrate hydrostatics, radiation damping, and wave excitation.
+ * ARCHITECTURE:
+ * TestHydro is a thin adapter over HydroSystem + ChronoHydroCoupler.
+ * Force computation is handled internally by HydroSystem, which owns:
+ *   - HydrostaticsComponent (restoring forces + buoyancy)
+ *   - RadiationComponent (RIRF convolution, owns velocity history)
+ *   - ExcitationComponent (wave forcing)
+ *
+ * The sign convention is: total = hydrostatics - radiation + waves.
+ * (RadiationComponent applies the negative sign internally since damping
+ * opposes motion.)
  *
  * MAIN RESPONSIBILITIES:
- * - TestHydro: Main orchestrator class for all hydrodynamic force components
+ * - TestHydro: Façade over HydroSystem; provides Chrono force callbacks
  * - ForceFunc6d: Wraps 6-DOF force/torque callbacks for Chrono bodies
  * - ComponentFunc: Per-DOF force function for Chrono's ChForce system
  *
+ * COMPONENT CONSTRUCTION:
+ * Force components are constructed via factory methods:
+ *   - CreateHydrostaticsComponent()
+ *   - CreateRadiationComponent()
+ *   - CreateExcitationComponent()
+ * These are the single source of truth for component configuration.
+ *
  * INTERACTIONS:
  * - Used by setup_hydro_from_yaml to create and configure TestHydro instances
  * - Called by Chrono during simulation via ChForce callbacks
@@ -27,21 +40,9 @@
  * - Bodies are 1-indexed in CoordinateFuncForBody (legacy, from ForceFunc6d)
  * - All bodies share same H5 file (multibody data in single file)
  * - 6 DOF per body (surge, sway, heave, roll, pitch, yaw)
- * - Forces computed once per time step and cached
- *
- * KNOWN LIMITATIONS:
- * - Monolithic design: all physics components mixed in TestHydro
- * - Tight coupling to Chrono types (hard to test in isolation)
- * - Body indexing inconsistency (1-indexed vs 0-indexed)
- * - No per-body enable/disable of force components
- *
- * FUTURE REFACTORING:
- * This class will be replaced by HydroSystem + ChronoHydroCoupler.
- * TestHydro will become a thin legacy adapter for API stability.
+ * - Forces computed once per time step and cached via prev_time check
  *********************************************************************/
 
-// TODO: clean up include statements
-
 // Include hydro_types.h FIRST to ensure BodyForces and GeneralizedForce are available
 // before any other includes that might conflict (e.g., hydro_config_types.h)
 #include <hydroc/hydro_types.h>
@@ -252,29 +253,35 @@ class TestHydro {
     void AddWaves(std::shared_ptr<WaveBase> waves);
 
     /**
-     * @brief Computes the Hydrostatic stiffness force plus buoyancy force for a 6N dimensional system.
+     * @brief Legacy API: Computes hydrostatic forces for a 6N dimensional system.
+     *
+     * NOTE: Retained for backward compatibility. The main force path now uses HydroSystem
+     * (via CoordinateFuncForBody). This method is not called during normal simulation.
      *
      * @return 6N dimensional force for 6 DOF and N bodies in system.
      */
     std::vector<double> ComputeForceHydrostatics();
 
     /**
-     * @brief Computes the Radiation Damping force with convolution history for a 6N dimensional system.
+     * @brief Legacy API: Computes radiation damping forces for a 6N dimensional system.
      *
-     * The discretization uses the time series of the the RIRF relative to the current step.
-     * Linear interpolation is done on the velocity history if time_sim-time_rirf is between two values of the time
-     * history. Trapezoidal integration is used to compute the force.
+     * NOTE: Retained for backward compatibility. The main force path now uses HydroSystem
+     * (via CoordinateFuncForBody). This method is not called during normal simulation.
      *
-     * Time history is automatically added in this function (so it should only be called once per time step), and
-     * history that is older than the maximum RIRF time value is automatically removed.
+     * Returns positive damping magnitudes (legacy convention); the main path applies
+     * the correct sign internally via RadiationComponent.
      *
      * @return 6N dimensional force for 6 DOF and N bodies in system.
      */
     std::vector<double> ComputeForceRadiationDampingConv();
 
     /**
-     * @brief Computes the 6N dimensional force from any waves applied to the system.
-     * @return 6N dimensional force for 6 DOF and N bodies in system (already Eigen type).
+     * @brief Legacy API: Computes wave excitation forces for a 6N dimensional system.
+     *
+     * NOTE: Retained for backward compatibility. The main force path now uses HydroSystem
+     * (via CoordinateFuncForBody). This method is not called during normal simulation.
+     *
+     * @return 6N dimensional force for 6 DOF and N bodies in system.
      */
     Eigen::VectorXd ComputeForceWaves();
     // Expose the wave object (read-only) so callers can query inputs if needed
@@ -352,8 +359,8 @@ class TestHydro {
     // Hydrodynamics profiling accessors
     HydroProfileStats GetProfileStats() const { return profile_stats_; }
 
-    // Compare mode: when enabled, compares legacy force path vs HydroSystem path
-    // and logs any discrepancies. Default is off.
+    // Compare mode: legacy debugging feature, retained for API compatibility.
+    // No longer affects behavior since the main path now uses HydroSystem exclusively.
     void SetCompareMode(bool enable) { compare_mode_ = enable; }
     bool GetCompareMode() const { return compare_mode_; }
 
@@ -416,7 +423,7 @@ class TestHydro {
     std::unique_ptr<hydrochrono::hydro::HydroSystem> hydro_system_;
     std::unique_ptr<hydrochrono::hydro::ChronoHydroCoupler> chrono_coupler_;
 
-    // Compare mode flag: when true, compares legacy vs HydroSystem forces
+    // Legacy compare mode flag: retained for API compatibility; no longer used.
     bool compare_mode_ = false;
 
     // Convolution kernel preprocessing (optional)
@@ -451,15 +458,9 @@ class TestHydro {
     // Note: Each instance owns its own velocity history (they are independent).
     std::unique_ptr<hydrochrono::hydro::RadiationComponent> CreateRadiationComponent() const;
 
-    // Internal helpers for HydroSystem + ChronoHydroCoupler path
-    // Constructs hydro_system_ and chrono_coupler_ once; subsequent calls are no-ops.
+    // Internal helper: constructs hydro_system_ and chrono_coupler_ once.
+    // Subsequent calls are no-ops. Called automatically by CoordinateFuncForBody().
     void EnsureHydroSystemAndCoupler();
-    // Evaluate forces via HydroSystem (used by comparison harness)
-    hydrochrono::hydro::BodyForces EvaluateHydroSystem(double time);
-
-    // Comparison harness: compares legacy forces vs HydroSystem forces and logs discrepancies.
-    // Only called when compare_mode_ is true.
-    void CompareWithHydroSystem(double time, int total_dofs);
 };
 
 #endif

src/hydro_forces.cpp

@@ -4,14 +4,34 @@
  * @brief Implementation of TestHydro main class and helper classes
  * ComponentFunc and ForceFunc6d.
  *
- * OVERVIEW:
- * This file implements the core hydrodynamic force computation for multibody
- * marine systems. It connects Chrono physics bodies to hydrodynamic components
- * (hydrostatics, radiation damping, wave excitation) and applies the
- * resulting forces during simulation.
+ * ARCHITECTURE:
+ * TestHydro is a thin adapter over HydroSystem + ChronoHydroCoupler.
+ *
+ * Force computation flow:
+ *   1. Chrono calls CoordinateFuncForBody(body, dof) via ChForce callbacks.
+ *   2. CoordinateFuncForBody detects new timesteps and calls chrono_coupler_->Evaluate().
+ *   3. ChronoHydroCoupler extracts system state from Chrono bodies and invokes HydroSystem.
+ *   4. HydroSystem evaluates all force components and returns BodyForces.
+ *   5. Forces are flattened into the legacy total_force_ array for Chrono consumption.
+ *
+ * Sign convention: total = hydrostatics - radiation + waves
+ *   (RadiationComponent applies the negative sign internally since damping opposes motion.)
+ *
+ * COMPONENT CONSTRUCTION:
+ * Force components are built via factory methods (single source of truth):
+ *   - CreateHydrostaticsComponent()
+ *   - CreateRadiationComponent()
+ *   - CreateExcitationComponent()
+ *
+ * LEGACY API:
+ * The following methods are retained for backward compatibility but are NOT used
+ * by the main force path (CoordinateFuncForBody):
+ *   - ComputeForceHydrostatics()
+ *   - ComputeForceRadiationDampingConv()
+ *   - ComputeForceWaves()
  *
  * MAIN RESPONSIBILITIES:
- * - TestHydro: Orchestrates all hydrodynamic force components for multiple bodies
+ * - TestHydro: Façade over HydroSystem; provides Chrono force callbacks
  * - ForceFunc6d: Wraps Chrono ChForce/ChTorque callbacks for each body
  * - ComponentFunc: Provides per-DOF force values to Chrono's force system
  *
@@ -26,24 +46,11 @@
  * - Bodies are 1-indexed in ForceFunc6d interface (legacy)
  * - 6 DOF per body (surge, sway, heave, roll, pitch, yaw)
  * - Forces computed once per time step, cached via prev_time check
- * - RadiationComponent is the single source of truth for radiation history/forces
- *
- * KNOWN LIMITATIONS:
- * - Monolithic design: all force components mixed in one class
- * - Tight coupling to Chrono types (hard to test without Chrono)
- * - Body indexing inconsistency (1-indexed in some places, 0-indexed in others)
- * - No per-body enable/disable of radiation or excitation (system-wide only)
  *
  * DEBUG INSTRUMENTATION:
  * To enable debug output, compile with -DHYDROCHRONO_DEBUG.
- * This will print force components and timing info.
- * Baseline outputs to record:
- * - Total force vector per time step (total_force_)
- * - Individual components (hydrostatic, radiation, waves)
- * - RIRF kernel access patterns
  *********************************************************************/
 
-// TODO minimize include statements, move all to header file hydro_forces.h?
 #include "hydroc/hydro_forces.h"
 #include <hydroc/chloadaddedmass.h>
 #include <hydroc/h5fileinfo.h>
@@ -547,59 +554,6 @@ void TestHydro::EnsureHydroSystemAndCoupler() {
         hydro_system_shared, bodies_);
 }
 
-hydrochrono::hydro::BodyForces TestHydro::EvaluateHydroSystem(double time) {
-    EnsureHydroSystemAndCoupler();
-    return chrono_coupler_->Evaluate(time);
-}
-
-void TestHydro::CompareWithHydroSystem(double time, int total_dofs) {
-    // DEPRECATED: This method is no longer called from CoordinateFuncForBody().
-    // The main force path now routes through HydroSystem, so there's nothing to compare.
-    // This method is kept for potential external use but will likely be removed in a future cleanup.
-
-    // Evaluate forces via the persistent HydroSystem path
-    hydrochrono::hydro::BodyForces new_forces = EvaluateHydroSystem(time);
-
-    // Flatten new_forces into legacy ordering: [body0_dof0, ..., body0_dof5, body1_dof0, ...]
-    std::vector<double> new_total(total_dofs, 0.0);
-    for (int body_idx = 0; body_idx < num_bodies_; ++body_idx) {
-        const int offset = kDofPerBody * body_idx;
-        for (int dof = 0; dof < kDofPerBody; ++dof) {
-            new_total[offset + dof] = new_forces[body_idx][dof];
-        }
-    }
-
-    // Comparison tolerances
-    const double abs_tol = 1e-6;  // Absolute tolerance (N or N·m)
-    const double rel_tol = 1e-4;  // Relative tolerance (0.01%)
-
-    // Compare entry-by-entry
-    int mismatch_count = 0;
-    for (int idx = 0; idx < total_dofs; ++idx) {
-        const double legacy_val = total_force_[idx];
-        const double new_val = new_total[idx];
-        const double diff = std::abs(legacy_val - new_val);
-        const double max_abs = std::max(std::abs(legacy_val), std::abs(new_val));
-        const double rel_diff = (max_abs > abs_tol) ? diff / max_abs : 0.0;
-
-        if (diff > abs_tol && rel_diff > rel_tol) {
-            const int body_idx = idx / kDofPerBody;
-            const int dof = idx % kDofPerBody;
-            std::cout << "[COMPARE] Mismatch at t=" << time
-                      << " body=" << body_idx << " dof=" << dof
-                      << ": legacy=" << legacy_val
-                      << " new=" << new_val
-                      << " diff=" << diff
-                      << " rel=" << rel_diff << std::endl;
-            ++mismatch_count;
-        }
-    }
-
-    if (mismatch_count == 0) {
-        // Optionally log success at intervals (every 100 steps or so)
-        // For now, silent when all forces match
-    }
-}
 
 // Legacy radiation convolution computation.
 // RadiationComponent is the single source of truth for radiation history and forces.
@@ -776,9 +730,6 @@ double TestHydro::CoordinateFuncForBody(int b, int dof_index) {
         }
     }
 
-    // Note: compare_mode_ is no longer used here because the main path IS HydroSystem now.
-    // The legacy per-component path is preserved for external callers but not used by this method.
-
 #ifdef HYDROCHRONO_DEBUG
     std::cout << "[HYDRO_DEBUG] CoordinateFuncForBody: body=" << b << ", dof=" << dof_index 
               << ", time=" << prev_time << ", force=" << total_force_[body_num_offset + dof_index] << std::endl;