Skip to content

Commit 62f99ec

Browse files
committed
Use ChLoadHydrodynamics for added mass; retain legacy as toggleable option
Switch default added-mass path to Chrono's built-in ChLoadHydrodynamics. The original ChLoadAddedMass (ChLoadCustomMultiple-based) implementation is kept in the repo and can be re-activated by defining HYDROCHRONO_USE_LEGACY_ADDED_MASS in hydro_system.h. Also properly documents the regular wave excitation phase indexing fix (body_offset was missing for body 2+ in GetForceAtTime), which affected multibody regular wave results such as the RM3 plate heave response.
1 parent 0748d12 commit 62f99ec

10 files changed

Lines changed: 327 additions & 82 deletions

File tree

CHANGELOG.md

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,21 +21,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
2121

2222
### Changed
2323

24+
- **Switched added-mass implementation to Chrono's built-in `ChLoadHydrodynamics`** — replaces HydroChrono's legacy `ChLoadAddedMass` (`ChLoadCustomMultiple`-based) as the default. The legacy implementation is retained in the codebase and can be re-activated by uncommenting `#define HYDROCHRONO_USE_LEGACY_ADDED_MASS` in `hydro_system.h`.
2425
- Default solver changed from GMRES to SPARSE_QR for most regression tests (faster, deterministic)
2526
- Irregular wave surface evaluation performance improved (parallelism and caching)
27+
- Regression test report generator now correctly parses CTest logs for PASS/FAIL status instead of assuming PASS when plots exist
2628

2729
### Fixed
2830

2931
- **YAML runner: `LoadSolverData` never called** — YAML structure mismatch prevented solver data from being loaded
3032
- **YAML runner: mesh file paths broken**`m_script_directory` was empty, breaking relative `model_file:` paths
33+
- **Regression report image paths** — report generator now computes correct relative paths to plot images instead of using hardcoded paths
3134
- OSWEC solver switched from SPARSE_QR to GMRES to prevent divergence (see Known Issues)
3235
- RM3 decay test fixes and cleanup
3336
- Sphere irregular wave test default arguments corrected
34-
- Regular wave bug fixes
37+
- **Regular wave excitation phase indexing for multi-body systems**`RegularWave::GetForceAtTime()` used `excitation_force_phase_[rowEx]` instead of `excitation_force_phase_[body_offset + rowEx]`, causing body 2+ to use body 1's excitation phase. This affected the RM3 plate (and any second+ body) heave response in regular waves. Single-body models (e.g., sphere) were unaffected.
3538

3639
### Known Issues
3740

38-
- **SPARSE_QR solver causes OSWEC simulations to diverge.** The constrained multi-body OSWEC model (revolute joints) diverges under the SPARSE_QR solver. OSWEC demos and tests use GMRES as a workaround. Other models (sphere, RM3, F3OF) are unaffected. Root cause not yet diagnosed.
41+
- **OSWEC tests use GMRES solver.** OSWEC demos and tests use GMRES rather than SPARSE_QR. SPARSE_QR may work for OSWEC but has not yet been validated.
3942
- **PSOR solver cannot handle added-mass assembly.** The added-mass determinism unit test solver sweep reports that PSOR errors out because it cannot handle stiffness/damping matrices. All other swept solvers (SPARSE_QR, SPARSE_LU, MINRES, GMRES, BICGSTAB, BARZILAIBORWEIN, APGD) pass assembly.
4043

4144
## [0.4.0] — 2025

CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -303,6 +303,7 @@ set(HC_SOURCES
303303
src/hydro/radiation/radiation_rirf_convolution.cpp
304304
src/hydro/core/hydro_forces.cpp
305305
# Chrono coupling layer (bridges Chrono physics with Chrono-free core)
306+
src/hydro/chrono/added_mass.cpp
306307
src/hydro/chrono/chrono_state_utils.cpp
307308
src/hydro/chrono/chrono_hydro_coupler.cpp
308309
src/hydro/force_components/hydrostatics_component.cpp
@@ -325,6 +326,7 @@ set(HC_HEADERS
325326
include/hydroc/core/hydro_forces.h
326327
include/hydroc/core/hydro_types.h
327328
include/hydroc/core/system_state.h
329+
include/hydroc/coupling/added_mass.h
328330
include/hydroc/coupling/chrono_coupler.h
329331
include/hydroc/io/h5_reader.h
330332
include/hydroc/io/h5_writer.h

docs/_main_pages/developer_docs/unit_tests.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ ctest -C Release -L unit --test-dir build -V
2929

3030
## Test: added-mass assembly (`test_added_mass_determinism`)
3131

32-
Verifies that Chrono's `ChLoadHydrodynamics` correctly assembles a cross-coupled added-mass matrix for a multi-body system.
32+
Verifies that Chrono's `ChLoadHydrodynamics` (the default added-mass path) correctly assembles a cross-coupled added-mass matrix for a multi-body system.
3333

3434
### What it does
3535

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
/*********************************************************************
2+
* @file added_mass.h
3+
* @brief ChLoadAddedMass: Chrono load for infinite-frequency added mass.
4+
*
5+
* STATUS: LEGACY IMPLEMENTATION -- retained for reference and as a
6+
* switchable alternative. The default added-mass path now uses
7+
* Chrono's built-in ChLoadHydrodynamics.
8+
*
9+
* To activate this implementation, define HYDROCHRONO_USE_LEGACY_ADDED_MASS
10+
* in hydro_system.h (see toggle comment near the top of that file).
11+
*
12+
* MAIN TYPES:
13+
* - ChLoadAddedMass: ChLoadCustomMultiple subclass that applies
14+
* added mass from HDF5 hydrodynamic data to Chrono bodies.
15+
*
16+
* ROLE: Chrono coupling layer. Provides a ChLoadCustomMultiple-based
17+
* approach to applying added mass forces to the Chrono simulation system.
18+
*********************************************************************/
19+
20+
#ifndef HYDROC_COUPLING_ADDED_MASS_H
21+
#define HYDROC_COUPLING_ADDED_MASS_H
22+
23+
#include <chrono/core/ChMatrix.h>
24+
#include <chrono/physics/ChBody.h>
25+
#include <hydroc/io/h5_reader.h>
26+
#include <vector>
27+
28+
#include <chrono/physics/ChLoad.h>
29+
#include <chrono/physics/ChSystem.h>
30+
31+
// =============================================================================
32+
class ChLoadAddedMass : public chrono::ChLoadCustomMultiple {
33+
public:
34+
/**
35+
* @brief Initializes body to have load applied to and added mass matrix from h5 file initialized object.
36+
*
37+
* Vector of bodies with hydro forces (and added mass) need to be listed in same order as they are added to the
38+
* system. Also need to be separate from any bodies without hydro forces (or added mass) and hydro bodies need to be
39+
* added to system before any bodies without hydro forces applied.
40+
*
41+
* @param body_info_struct HydroData::BodyInfo for each body with h5 information including added mass matrix
42+
* @param bodies vector of Project Chrono bodies to apply added mass to. Must be added to system in same order as in
43+
* this matrix.
44+
* @param system pointer to system containing the bodies, used for getting system mass matrix size at any time.
45+
*/
46+
ChLoadAddedMass(const std::vector<HydroData::BodyInfo>& body_info_struct,
47+
std::vector<std::shared_ptr<chrono::ChLoadable>>& bodies,
48+
chrono::ChSystem* system);
49+
50+
/**
51+
* @brief "Virtual" copy constructor (covariant return type). Required from chrono inheritance.
52+
*/
53+
virtual ChLoadAddedMass* Clone() const override { return new ChLoadAddedMass(*this); }
54+
55+
/**
56+
* @brief Compute Q, the generalized load.
57+
*
58+
* From Chrono documentation:
59+
* In this case, it computes the quadratic (centrifugal, gyroscopic) terms.
60+
* Signs are negative as Q assumed at right hand side, so Q= -Fgyro -Fcentrifugal
61+
* Called automatically at each Update().
62+
* The M*a term is not added: to this end one could use LoadIntLoadResidual_Mv afterward.
63+
*/
64+
virtual void ComputeQ(chrono::ChState* state_x, ///< state position to evaluate Q
65+
chrono::ChStateDelta* state_w ///< state speed to evaluate Q
66+
) override {}
67+
68+
/**
69+
* @brief This is the function that sets the infinite added mass matrix every timestep.
70+
*
71+
* From Chrono docs:
72+
* Compute the K=-dQ/dx, R=-dQ/dv, M=-dQ/da Jacobians.
73+
* Implementation in a derived class should load the Jacobian matrices K, R, M in the structure 'm_jacobians'.
74+
* Note the sign that is flipped because we assume equations are written with Q moved to left-hand side.
75+
*
76+
* @param state_x state position to evaluate jacobians
77+
* @param state_w state speed to evaluate jacobians
78+
*/
79+
virtual void ComputeJacobian(chrono::ChState* state_x, chrono::ChStateDelta* state_w) override;
80+
81+
/**
82+
* @brief Computes LoadIntLoadResidual_Mv for vector w, const c, and vector R. Also carried over from chrono
83+
* inheritance.
84+
*
85+
* Note R here is vector, and is not R gyroscopic damping matrix from ComputeJacobian.
86+
* Just for efficiency, override the default LoadIntLoadResidual_Mv, because we can do this in a simplified way.
87+
*
88+
* @param R result: the R residual, R += c*M*w
89+
* @param w the w vector
90+
* @param c a scaling factor
91+
*/
92+
virtual void LoadIntLoadResidual_Mv(chrono::ChVectorDynamic<>& R,
93+
const chrono::ChVectorDynamic<>& w,
94+
const double c) override;
95+
96+
private:
97+
chrono::ChSystem* system;
98+
chrono::ChMatrixDynamic<double> infinite_added_mass; ///< added mass at infinite frequency in global coordinates
99+
chrono::ChMatrixDynamic<double>
100+
infinite_added_mass_system; ///< added mass at infinite frequency in global coordinates (system matrix)
101+
virtual bool IsStiff() override { return true; } // this to force the use of the inertial M, R and K matrices
102+
};
103+
104+
#endif // HYDROC_COUPLING_ADDED_MASS_H
105+

include/hydroc/hydro_system.h

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,16 @@
1717
#ifndef HYDRO_SYSTEM_H
1818
#define HYDRO_SYSTEM_H
1919

20+
// ─────────────────────────────────────────────────────────────────────────────
21+
// Added Mass Implementation Toggle
22+
// ─────────────────────────────────────────────────────────────────────────────
23+
// By default, HydroChrono uses Chrono's built-in ChLoadHydrodynamics for
24+
// infinite-frequency added mass. To switch to HydroChrono's legacy
25+
// ChLoadAddedMass (ChLoadCustomMultiple-based) implementation, uncomment:
26+
//
27+
// #define HYDROCHRONO_USE_LEGACY_ADDED_MASS
28+
//
29+
2030
// Include hydro_types.h FIRST to ensure BodyForces and GeneralizedForce are available
2131
// before any other includes that might conflict (e.g., config/hydro_config.h)
2232
#include <hydroc/core/hydro_types.h>
@@ -37,7 +47,9 @@
3747
// ChForce.h includes ChFunction internally, so we don't need a separate include.
3848
#include <chrono/physics/ChBody.h>
3949
#include <chrono/physics/ChForce.h>
50+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
4051
#include <chrono/physics/ChLoadContainer.h>
52+
#endif
4153

4254
// Hydroc library includes
4355
#include <hydroc/io/h5_reader.h>
@@ -185,8 +197,10 @@ class ForceFunc6d {
185197
HydroSystem* all_hydro_forces_ = nullptr; ///< Pointer to HydroSystem for calculations.
186198
};
187199

200+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
188201
// Forward declaration of ChLoadAddedMass (defined in added_mass.h, included in .cpp)
189202
class ChLoadAddedMass;
203+
#endif
190204

191205
// Lightweight hydrodynamics profiling stats
192206
struct HydroProfileStats {
@@ -436,9 +450,11 @@ class HydroSystem {
436450
hydrochrono::hydro::SystemState cached_state_;
437451
double cached_state_time_ = std::numeric_limits<double>::quiet_NaN();
438452

439-
// Added mass related properties
453+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
454+
// Added mass related properties (legacy ChLoadAddedMass path)
440455
std::shared_ptr<chrono::ChLoadContainer> my_loadcontainer;
441456
std::shared_ptr<ChLoadAddedMass> my_loadbodyinertia;
457+
#endif
442458

443459
/**
444460
* @brief Returns the cached SystemState for the given time.

src/hydro/chrono/added_mass.cpp

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
/*********************************************************************
2+
* @file added_mass.cpp
3+
* @brief Implementation of ChLoadAddedMass: Chrono load for added mass.
4+
*
5+
* STATUS: LEGACY IMPLEMENTATION -- retained for reference and as a
6+
* switchable alternative. The default added-mass path now uses
7+
* Chrono's built-in ChLoadHydrodynamics.
8+
*
9+
* To activate this implementation, define HYDROCHRONO_USE_LEGACY_ADDED_MASS
10+
* in hydro_system.h (see toggle comment near the top of that file).
11+
*
12+
* ROLE: Chrono coupling layer. Provides a ChLoadCustomMultiple subclass
13+
* that applies infinite-frequency added mass from HDF5 data to bodies.
14+
*********************************************************************/
15+
16+
#include <hydroc/coupling/added_mass.h>
17+
18+
#include <utility>
19+
20+
#include "chrono/physics/ChBody.h"
21+
22+
using namespace chrono;
23+
24+
ChLoadAddedMass::ChLoadAddedMass(const std::vector<HydroData::BodyInfo>& user_h5_body_data,
25+
std::vector<std::shared_ptr<ChLoadable>>& bodies,
26+
ChSystem* system)
27+
: ChLoadCustomMultiple(bodies), system(system) {
28+
auto nBodies = bodies.size();
29+
30+
infinite_added_mass.setZero(6 * nBodies, 6 * nBodies);
31+
for (int i = 0; i < nBodies; i++) {
32+
infinite_added_mass.block(i * 6, 0, 6, nBodies * 6) = user_h5_body_data[i].inf_added_mass;
33+
}
34+
35+
// initialize added mass matrix for whole system
36+
infinite_added_mass_system = infinite_added_mass;
37+
}
38+
39+
void ChLoadAddedMass::ComputeJacobian(ChState* state_x, ///< state position to evaluate jacobians
40+
ChStateDelta* state_w ///< state speed to evaluate jacobians
41+
) {
42+
// The following ensures that the added mass matrix matches the size of the ChSystem mass matrix. It is necessary
43+
// for systems that have both hydro and non-hydro bodies when adding a system-wide load.
44+
// @todo if possible, remove hack by using initialiazer function called AFTER the ChSystem is assembled.
45+
46+
// get mass matrix size
47+
auto mmrows = system->GetNumCoordsVelLevel(); // number of rows in system mass matrix
48+
// check if ChSystem mass matrix size different from added mass matrix size
49+
if (mmrows != infinite_added_mass_system.rows() && mmrows > 0) {
50+
// initialize/update system matrix;
51+
infinite_added_mass_system.setZero(mmrows, mmrows);
52+
auto amrows = infinite_added_mass.rows();
53+
infinite_added_mass_system.block(0, 0, amrows, amrows) = infinite_added_mass;
54+
}
55+
// set mass matrix here
56+
m_jacobians->M = infinite_added_mass_system;
57+
58+
// R gyroscopic damping matrix terms (6Nx6N)
59+
// 0 for added mass
60+
m_jacobians->R.setZero();
61+
62+
// K inertial stiffness matrix terms (6Nx6N)
63+
// 0 for added mass
64+
m_jacobians->K.setZero();
65+
}
66+
67+
void ChLoadAddedMass::LoadIntLoadResidual_Mv(ChVectorDynamic<>& R, const ChVectorDynamic<>& w, const double c) {
68+
if (!this->m_jacobians) return;
69+
70+
// R+=c*M*a
71+
// segment gives the chunk of vector starting at the first argument, and going for as many elements as the second
72+
// argument... in this case, segment gets the 3vector starting at the 0th DOF's offset (ie 0)
73+
// R.segment(loadable->GetSubBlockOffset(0), 3) += c * (this->mass * (a_x + chrono::Vcross(a_w,
74+
// this->c_m))).eigen();
75+
// in this case, segment gets the 3vector starting at the 0th DOF's + 3 offset (ie 3)
76+
// R.segment(loadable->GetSubBlockOffset(0) + 3, 3) += c * (this->mass * chrono::Vcross(this->c_m, a_x) + this->I *
77+
// a_w).eigen();
78+
// since R is a vector, we can probably just do R += C*M*a with no need to separate w into a_x and a_w above
79+
R += c * m_jacobians->M * w;
80+
}

src/hydro/hydro_system.cpp

Lines changed: 38 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,11 @@
1010
*********************************************************************/
1111

1212
#include "hydroc/hydro_system.h"
13+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
14+
#include <hydroc/coupling/added_mass.h>
15+
#else
16+
#include <chrono/physics/ChLoadHydrodynamics.h>
17+
#endif
1318
#include <hydroc/io/h5_reader.h>
1419
#include <hydroc/waves/wave_base.h>
1520
#include <hydroc/waves/regular_wave.h>
@@ -25,13 +30,14 @@
2530
#include "radiation/radiation_rirf_processing.h"
2631

2732
// Chrono includes needed for implementation (not exposed in public header)
33+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
2834
#include <chrono/physics/ChLoad.h>
2935
#include <chrono/physics/ChLoadable.h>
36+
#endif
3037
#include <chrono/physics/ChBodyEasy.h>
3138
#include <chrono/physics/ChLoadsBody.h>
3239
#include <chrono/physics/ChSystemNSC.h>
3340
#include <chrono/physics/ChSystemSMC.h>
34-
#include <chrono/physics/ChLoadHydrodynamics.h>
3541
#include <chrono/solver/ChIterativeSolverLS.h>
3642
#include <chrono/timestepper/ChTimestepper.h>
3743
#include <chrono/timestepper/ChTimestepperHHT.h>
@@ -59,9 +65,16 @@
5965
// Local using declarations for Chrono types (implementation only - not leaked to users)
6066
// These provide convenience within this translation unit while keeping the public header clean.
6167
using chrono::ChBody;
68+
using chrono::ChBodyEasyMesh;
6269
using chrono::ChForce;
63-
using chrono::ChBodyAddedMassBlocks;
64-
using chrono::ChLoadHydrodynamics;
70+
using chrono::ChFunction;
71+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
72+
using chrono::ChLoadable;
73+
using chrono::ChLoadContainer;
74+
#endif
75+
using chrono::ChSolver;
76+
using chrono::ChSystem;
77+
using chrono::ChVector3d;
6578

6679
// kDofPerBody is already defined in hydro_system.h as constexpr
6780
const int kDofLinOrRot = 3;
@@ -277,19 +290,37 @@ HydroSystem::HydroSystem(std::vector<std::shared_ptr<ChBody>> user_bodies,
277290
force_per_body_.emplace_back(bodies_[b], this);
278291
}
279292

280-
// Handle added mass info (applied via Chrono's ChLoadHydrodynamics).
293+
// Handle added mass info (infinite-frequency added mass applied to Chrono system).
294+
// Toggle between implementations via HYDROCHRONO_USE_LEGACY_ADDED_MASS in hydro_system.h.
295+
#ifdef HYDROCHRONO_USE_LEGACY_ADDED_MASS
296+
// Legacy path: HydroChrono's ChLoadAddedMass (ChLoadCustomMultiple-based).
297+
my_loadcontainer = std::make_shared<ChLoadContainer>();
298+
299+
std::vector<std::shared_ptr<ChLoadable>> loadables(bodies_.size());
300+
for (int i = 0; i < static_cast<int>(bodies_.size()); ++i) {
301+
loadables[i] = bodies_[i];
302+
}
303+
304+
my_loadbodyinertia =
305+
std::make_shared<ChLoadAddedMass>(file_info_.GetBodyInfos(), loadables, bodies_[0]->GetSystem());
306+
307+
bodies_[0]->GetSystem()->Add(my_loadcontainer);
308+
my_loadcontainer->Add(my_loadbodyinertia);
309+
#else
310+
// Default path: Chrono's built-in ChLoadHydrodynamics.
281311
// ChBodyAddedMassBlocks is a std::vector<ChBodyAddedMassBlock>; each entry
282-
// pairs a body with its 6×6N row of the infinite-frequency added mass matrix.
312+
// pairs a body with its 6x6N row of the infinite-frequency added mass matrix.
283313
const auto& body_info = file_info_.GetBodyInfos();
284-
ChBodyAddedMassBlocks body_blocks;
314+
chrono::ChBodyAddedMassBlocks body_blocks;
285315
for (size_t i = 0; i < num_bodies_; i++) {
286316
body_blocks.push_back({bodies_[i], body_info[i].inf_added_mass});
287317
}
288318
if (num_bodies_ > 0) {
289-
auto hydro_load = chrono_types::make_shared<ChLoadHydrodynamics>(body_blocks);
319+
auto hydro_load = chrono_types::make_shared<chrono::ChLoadHydrodynamics>(body_blocks);
290320
hydro_load->SetVerbose(false);
291321
bodies_[0]->GetSystem()->Add(hydro_load);
292322
}
323+
#endif
293324

294325
// Set up hydro inputs
295326
user_waves_ = waves;

tests/regression/rm3/test_rm3_decay.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ int main(int argc, char* argv[]) {
5050

5151
double timestep = 0.01;
5252
system.SetTimestepperType(ChTimestepper::Type::HHT);
53-
system.SetSolverType(ChSolver::Type::SPARSE_QR);
53+
system.SetSolverType(ChSolver::Type::GMRES);
5454
double simulationDuration = 40.0;
5555

5656
// Output timeseries

tests/regression/rm3/test_rm3_reg_waves.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ int main(int argc, char* argv[]) {
4949
system.SetGravitationalAcceleration(ChVector3d(0.0, 0.0, -9.81));
5050
double timestep = 0.01;
5151
system.SetTimestepperType(ChTimestepper::Type::HHT);
52-
system.SetSolverType(ChSolver::Type::SPARSE_QR);
52+
system.SetSolverType(ChSolver::Type::GMRES);
5353
double simulationDuration = 40.0;
5454

5555
// Output timeseries

0 commit comments

Comments
 (0)