Skip to content

Commit 014e33b

Browse files
dav-ogdav-og
committed
Add RadiationStateSpaceModel math core with visual verification
- Implement state-space radiation model with exponential decay modes - Use exact exponential integrator for unconditional stability - Add 7 unit tests covering validation, steady-state, transient, multi-mode - Add diagnostic tool + Python plotter for visual verification (tests/manual/) - Reorganize test output: unit/, manual/, regression/ subdirectories - Model is Chrono-free, ready for integration in future stages
1 parent 2c88d25 commit 014e33b

7 files changed

Lines changed: 1482 additions & 0 deletions

File tree

‎src/hydro/radiation/radiation_ss_model.cpp‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
/*********************************************************************
2+
* @file radiation_ss_model.cpp
3+
* @brief Implementation of RadiationStateSpaceModel.
4+
*********************************************************************/
5+
6+
#include "radiation_ss_model.h"
7+
8+
#include <cmath>
9+
#include <stdexcept>
10+
#include <string>
11+
12+
namespace hydrochrono::hydro {
13+
14+
RadiationStateSpaceModel::RadiationStateSpaceModel(int num_dofs, int num_modes)
15+
: num_dofs_(num_dofs), num_modes_(num_modes) {
16+
17+
if (num_dofs <= 0) {
18+
throw std::invalid_argument(
19+
"RadiationStateSpaceModel: num_dofs must be > 0 (got " +
20+
std::to_string(num_dofs) + ")");
21+
}
22+
if (num_modes <= 0) {
23+
throw std::invalid_argument(
24+
"RadiationStateSpaceModel: num_modes must be > 0 (got " +
25+
std::to_string(num_modes) + ")");
26+
}
27+
28+
// Initialize storage
29+
alphas_.setZero(num_modes_);
30+
H_.setZero(num_dofs_, num_modes_);
31+
z_.setZero(num_dofs_, num_modes_);
32+
}
33+
34+
void RadiationStateSpaceModel::SetModeParameters(
35+
int mode_index,
36+
double alpha,
37+
const Eigen::VectorXd& h_column) {
38+
39+
if (mode_index < 0 || mode_index >= num_modes_) {
40+
throw std::out_of_range(
41+
"RadiationStateSpaceModel::SetModeParameters: mode_index " +
42+
std::to_string(mode_index) + " out of range [0, " +
43+
std::to_string(num_modes_) + ")");
44+
}
45+
if (alpha <= 0.0) {
46+
throw std::invalid_argument(
47+
"RadiationStateSpaceModel::SetModeParameters: alpha must be > 0 (got " +
48+
std::to_string(alpha) + ")");
49+
}
50+
if (h_column.size() != num_dofs_) {
51+
throw std::invalid_argument(
52+
"RadiationStateSpaceModel::SetModeParameters: h_column size mismatch "
53+
"(expected " + std::to_string(num_dofs_) + ", got " +
54+
std::to_string(h_column.size()) + ")");
55+
}
56+
57+
alphas_(mode_index) = alpha;
58+
H_.col(mode_index) = h_column;
59+
}
60+
61+
void RadiationStateSpaceModel::Reset() {
62+
z_.setZero();
63+
}
64+
65+
void RadiationStateSpaceModel::Step(double dt, const Eigen::VectorXd& v) {
66+
if (dt <= 0.0) {
67+
throw std::invalid_argument(
68+
"RadiationStateSpaceModel::Step: dt must be > 0 (got " +
69+
std::to_string(dt) + ")");
70+
}
71+
if (v.size() != num_dofs_) {
72+
throw std::invalid_argument(
73+
"RadiationStateSpaceModel::Step: velocity vector size mismatch "
74+
"(expected " + std::to_string(num_dofs_) + ", got " +
75+
std::to_string(v.size()) + ")");
76+
}
77+
78+
// Exact exponential integration for each mode:
79+
// z_new = exp(-α dt) * z_old + [1 - exp(-α dt)] / α * v
80+
//
81+
// This is unconditionally stable for any α > 0 and dt > 0.
82+
// The coefficient [1 - exp(-α dt)] / α approaches dt as α → 0,
83+
// but we assume α > 0 (validated in SetModeParameters).
84+
85+
for (int m = 0; m < num_modes_; ++m) {
86+
const double alpha = alphas_(m);
87+
88+
// Skip uninitialized modes (alpha == 0 shouldn't happen after proper setup)
89+
if (alpha <= 0.0) {
90+
continue;
91+
}
92+
93+
const double exp_factor = std::exp(-alpha * dt);
94+
const double input_coeff = (1.0 - exp_factor) / alpha;
95+
96+
// Update state for mode m:
97+
// z_.col(m) = exp_factor * z_.col(m) + input_coeff * v
98+
z_.col(m) = exp_factor * z_.col(m) + input_coeff * v;
99+
}
100+
}
101+
102+
Eigen::VectorXd RadiationStateSpaceModel::GetForces() const {
103+
// Radiation force: f_rad = Σₘ Hₘ ⊙ zₘ
104+
//
105+
// With H_ and z_ both [num_dofs × num_modes]:
106+
// f_rad[i] = Σₘ H_[i,m] * z_[i,m]
107+
//
108+
// This is the row-wise sum of element-wise product.
109+
110+
// Element-wise product, then sum across columns (modes)
111+
return (H_.array() * z_.array()).rowwise().sum();
112+
}
113+
114+
} // namespace hydrochrono::hydro
115+

‎src/hydro/radiation/radiation_ss_model.h‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
/*********************************************************************
2+
* @file radiation_ss_model.h
3+
* @brief State-space model for radiation damping (exponential kernel approximation).
4+
*
5+
* This class implements the mathematical core of state-space radiation:
6+
* - Approximates the radiation impulse response as a sum of exponentials
7+
* - Integrates internal states using an exact exponential scheme
8+
* - Returns radiation forces given velocity inputs
9+
*
10+
* THEORY:
11+
* The radiation kernel K(τ) is approximated as:
12+
* K(τ) ≈ Σₘ Hₘ exp(−αₘ τ)
13+
*
14+
* This leads to ODEs for internal states zₘ:
15+
* żₘ(t) = v(t) − αₘ zₘ(t)
16+
*
17+
* And radiation forces:
18+
* f_rad(t) = Σₘ Hₘ ⊙ zₘ(t)
19+
*
20+
* where ⊙ denotes element-wise multiplication.
21+
*
22+
* USAGE:
23+
* 1. Construct with number of DOFs and modes
24+
* 2. Set parameters for each mode via SetModeParameters()
25+
* 3. Call Reset() to zero states
26+
* 4. In simulation loop: Step(dt, v), then GetForces()
27+
*
28+
* NOTE: This class is Chrono-free and operates on pure Eigen vectors.
29+
*********************************************************************/
30+
31+
#ifndef HYDRO_RADIATION_SS_MODEL_H
32+
#define HYDRO_RADIATION_SS_MODEL_H
33+
34+
#include <Eigen/Dense>
35+
36+
namespace hydrochrono::hydro {
37+
38+
/**
39+
* @brief State-space model for radiation damping.
40+
*
41+
* Models radiation forces using a sum of exponential decay modes.
42+
* Each mode m has:
43+
* - A decay rate αₘ > 0
44+
* - A gain vector Hₘ of size num_dofs
45+
* - An internal state vector zₘ of size num_dofs
46+
*
47+
* The model uses an exact exponential integrator for unconditional stability.
48+
*/
49+
class RadiationStateSpaceModel {
50+
public:
51+
/**
52+
* @brief Construct a state-space model.
53+
*
54+
* @param num_dofs Number of degrees of freedom (e.g., 6 for single body, 12 for two bodies)
55+
* @param num_modes Number of exponential modes in the approximation
56+
*
57+
* @throws std::invalid_argument if num_dofs <= 0 or num_modes <= 0
58+
*/
59+
RadiationStateSpaceModel(int num_dofs, int num_modes);
60+
61+
/**
62+
* @brief Set parameters for a single mode.
63+
*
64+
* Must be called for each mode (0 to num_modes-1) before stepping.
65+
*
66+
* @param mode_index Index of the mode (0-based)
67+
* @param alpha Decay rate (must be > 0 for stability)
68+
* @param h_column Gain vector mapping this mode's state to forces (size num_dofs)
69+
*
70+
* @throws std::out_of_range if mode_index is invalid
71+
* @throws std::invalid_argument if alpha <= 0 or h_column size mismatches
72+
*/
73+
void SetModeParameters(int mode_index, double alpha, const Eigen::VectorXd& h_column);
74+
75+
/**
76+
* @brief Reset all internal states to zero.
77+
*
78+
* Call this at the start of a simulation or to reinitialize.
79+
*/
80+
void Reset();
81+
82+
/**
83+
* @brief Advance the model by one time step.
84+
*
85+
* Uses exact exponential integration:
86+
* zₘ(t+dt) = exp(−αₘ dt) zₘ(t) + [1 − exp(−αₘ dt)] / αₘ · v
87+
*
88+
* This scheme is unconditionally stable for any dt > 0 and αₘ > 0.
89+
*
90+
* @param dt Time step (seconds, must be > 0)
91+
* @param v Velocity vector for all DOFs (size num_dofs)
92+
*
93+
* @throws std::invalid_argument if dt <= 0 or v size mismatches
94+
*/
95+
void Step(double dt, const Eigen::VectorXd& v);
96+
97+
/**
98+
* @brief Get the current radiation force vector.
99+
*
100+
* Computes: f_rad = Σₘ Hₘ ⊙ zₘ
101+
* where ⊙ is element-wise multiplication.
102+
*
103+
* @return Radiation force vector (size num_dofs)
104+
*/
105+
Eigen::VectorXd GetForces() const;
106+
107+
// Accessors for testing and diagnostics
108+
int num_dofs() const { return num_dofs_; }
109+
int num_modes() const { return num_modes_; }
110+
const Eigen::VectorXd& alphas() const { return alphas_; }
111+
const Eigen::MatrixXd& H() const { return H_; }
112+
const Eigen::MatrixXd& z() const { return z_; }
113+
114+
private:
115+
int num_dofs_; ///< Number of degrees of freedom
116+
int num_modes_; ///< Number of exponential modes
117+
118+
Eigen::VectorXd alphas_; ///< Decay rates [num_modes]
119+
Eigen::MatrixXd H_; ///< Gain matrix [num_dofs × num_modes]
120+
Eigen::MatrixXd z_; ///< Internal states [num_dofs × num_modes]
121+
};
122+
123+
} // namespace hydrochrono::hydro
124+
125+
#endif // HYDRO_RADIATION_SS_MODEL_H
126+

‎tests/manual/CMakeLists.txt‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
# Manual/Diagnostic Tests for HydroChrono
2+
# =========================================
3+
#
4+
# These are NOT automated tests - they are diagnostic tools for visual
5+
# verification and debugging. They output data files that can be plotted
6+
# using the accompanying Python scripts.
7+
#
8+
# Build: cmake --build build --target radiation_ss_visual_check
9+
# Run: ./radiation_ss_visual_check
10+
# Plot: python plot_radiation_ss.py
11+
#
12+
# Output structure:
13+
# build/bin/tests/manual/
14+
# radiation_ss_visual_check.exe
15+
# plot_radiation_ss.py
16+
# output/
17+
# radiation_ss_*.csv
18+
# plot_*.png
19+
20+
# Output to tests/manual/ subdirectory
21+
set(MANUAL_TEST_OUTPUT_DIR ${HYDROCHRONO_TEST_OUTPUT_DIR}/manual)
22+
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${MANUAL_TEST_OUTPUT_DIR})
23+
if(CMAKE_CONFIGURATION_TYPES)
24+
foreach(cfg ${CMAKE_CONFIGURATION_TYPES})
25+
string(TOUPPER "${cfg}" cfg_uc)
26+
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_${cfg_uc} ${MANUAL_TEST_OUTPUT_DIR})
27+
endforeach()
28+
endif()
29+
30+
# =============================================================================
31+
# State-Space Radiation Model Visual Verification
32+
# =============================================================================
33+
add_executable(radiation_ss_visual_check
34+
radiation_ss_visual_check.cpp
35+
../../src/hydro/radiation/radiation_ss_model.cpp
36+
)
37+
38+
target_include_directories(radiation_ss_visual_check
39+
PRIVATE
40+
../../src
41+
../../include
42+
)
43+
44+
target_link_libraries(radiation_ss_visual_check
45+
PRIVATE
46+
Eigen3::Eigen
47+
)
48+
49+
# Copy the Python plotting script to the build output directory
50+
configure_file(
51+
${CMAKE_CURRENT_SOURCE_DIR}/plot_radiation_ss.py
52+
${MANUAL_TEST_OUTPUT_DIR}/plot_radiation_ss.py
53+
COPYONLY
54+
)
55+

0 commit comments

Comments
 (0)