Skip to content

Commit 7e80c40

Browse files
dav-ogdav-og
committed
Add eta consistency regression test
New test that validates spectrum-generated and eta-file-imported irregular waves produce identical results (max diff: ~1e-11 m). This ensures both code paths handle the time range consistently and helps catch regressions in either wave generation method.
1 parent 7863073 commit 7e80c40

4 files changed

Lines changed: 352 additions & 1 deletion

File tree

‎data/reference_data/sphere/hc_ref_sphere_irreg_waves_eta_consistency.txt‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
# Reference file for sphere_irreg_waves_eta_consistency test
2+
# This test is self-validating: it compares spectrum-generated and eta-imported
3+
# irregular waves internally. The comparison script validates that the maximum
4+
# difference between the two code paths is within tolerance (1e-6 m).
5+
#
6+
# This file is a placeholder required by the build system.
7+
# The actual validation is performed by the test executable itself.
8+

‎tests/regression/CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,7 @@ build_regression_test(sphere sphere_decay FALSE)
7474
build_regression_test(sphere sphere_reg_waves TRUE)
7575
build_regression_test(sphere sphere_irreg_waves FALSE)
7676
build_regression_test(sphere sphere_irreg_waves_eta FALSE)
77+
build_regression_test(sphere sphere_irreg_waves_eta_consistency FALSE)
7778
build_regression_test(f3of f3of_dt1 FALSE)
7879
build_regression_test(f3of f3of_dt2 FALSE)
7980
build_regression_test(f3of f3of_dt3 FALSE)
@@ -102,5 +103,5 @@ set_tests_properties(
102103
LABELS "regression;report;documentation;core"
103104
# WORKING_DIRECTORY ${TEST_REGRESSION_OUTPUT_DIR}
104105
ENVIRONMENT "PATH=${Python3_ROOT_DIR};$ENV{PATH};HYDROCHRONO_BUILD_DIR=${CMAKE_BINARY_DIR}"
105-
DEPENDS "sphere_decay_regression;sphere_reg_waves_regression;sphere_irreg_waves_regression;sphere_irreg_waves_eta_regression;f3of_dt1_regression;f3of_dt2_regression;f3of_dt3_regression;oswec_decay_regression;oswec_reg_waves_regression;rm3_decay_regression;rm3_reg_waves_regression"
106+
DEPENDS "sphere_decay_regression;sphere_reg_waves_regression;sphere_irreg_waves_regression;sphere_irreg_waves_eta_regression;sphere_irreg_waves_eta_consistency_regression;f3of_dt1_regression;f3of_dt2_regression;f3of_dt3_regression;oswec_decay_regression;oswec_reg_waves_regression;rm3_decay_regression;rm3_reg_waves_regression"
106107
)

‎tests/regression/sphere/compare_sphere_irreg_waves_eta_consistency.py‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
#!/usr/bin/env python3
2+
"""
3+
HydroChrono Sphere Irregular Waves Eta Consistency Test Comparison
4+
5+
This script validates that spectrum-generated and eta-imported irregular waves
6+
produce identical results using the standardized comparison template.
7+
8+
Usage:
9+
python compare_sphere_irreg_waves_eta_consistency.py <reference_file> <results_file>
10+
Note: reference_file is unused for this self-validating test
11+
"""
12+
13+
import sys
14+
import os
15+
from pathlib import Path
16+
17+
# Import the comparison template
18+
sys.path.append(os.path.join(os.path.dirname(__file__), '../utilities'))
19+
from compare_template import run_consistency_comparison
20+
21+
if __name__ == '__main__':
22+
"""
23+
Compare spectrum-generated vs eta-imported irregular waves for consistency
24+
"""
25+
if len(sys.argv) != 3:
26+
print("Usage: python compare_sphere_irreg_waves_eta_consistency.py <reference_file> <results_file>")
27+
print("Note: reference_file is unused for this self-validating test")
28+
sys.exit(1)
29+
30+
# Reference file is unused - this is a self-validating test
31+
fname_rst = sys.argv[2]
32+
33+
if not os.path.exists(fname_rst):
34+
print(f"ERROR: Results file not found: {fname_rst}")
35+
sys.exit(1)
36+
37+
# Show where the plot will be saved
38+
test_file_path = Path(fname_rst)
39+
plots_dir = test_file_path.parent / "plots"
40+
plots_dir.mkdir(parents=True, exist_ok=True)
41+
print(f"Plot will be saved to: {plots_dir}")
42+
43+
# Test-specific configuration
44+
test_name = "Sphere Irreg Waves Eta Consistency"
45+
safe_test_name = test_name.lower().replace(' ', '_').replace('-', '_')
46+
print(f"Plot filename: {plots_dir}/{safe_test_name}_comparison.png")
47+
y_label = "Heave (m)"
48+
tolerance = 1e-6
49+
50+
# Run the consistency comparison using the template
51+
max_diff, passed = run_consistency_comparison(
52+
fname_rst, test_name, y_label,
53+
tolerance=tolerance,
54+
ref_label="Spectrum",
55+
sim_label="Eta-Import"
56+
)
57+
58+
sys.exit(0 if passed else 1)

‎tests/regression/sphere/test_sphere_irreg_waves_eta_consistency.cpp‎

Lines changed: 284 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,284 @@
1+
/**
2+
* @file test_sphere_irreg_waves_eta_consistency.cpp
3+
* @brief Validates that spectrum-generated and eta-file-imported irregular waves produce identical results.
4+
*
5+
* This test ensures the two code paths for irregular wave simulation are consistent:
6+
* 1. Generate waves from spectrum parameters (wave_height, wave_period)
7+
* 2. Import waves from an eta file
8+
*
9+
* The test runs both simulations with the same parameters and compares the heave response.
10+
* Any significant difference indicates a bug in one of the code paths.
11+
*/
12+
13+
#include <chrono>
14+
#include <cmath>
15+
#include <filesystem>
16+
#include <fstream>
17+
#include <iomanip>
18+
#include <iostream>
19+
#include <vector>
20+
21+
#include <chrono/core/ChRealtimeStep.h>
22+
#include <chrono/physics/ChBodyEasy.h>
23+
#include <chrono/physics/ChSystemNSC.h>
24+
25+
#include <hydroc/gui/guihelper.h>
26+
#include <hydroc/helper.h>
27+
#include <hydroc/hydro_forces.h>
28+
29+
using namespace chrono;
30+
31+
// Shared simulation parameters
32+
const double TIMESTEP = 0.015;
33+
const double DURATION = 60.0; // Short duration for faster testing
34+
const double WAVE_HEIGHT = 2.0;
35+
const double WAVE_PERIOD = 12.0;
36+
const double RAMP_DURATION = 0.0; // No ramp for exact comparison
37+
const int SEED = 42; // Fixed seed for reproducibility
38+
39+
int main(int argc, char* argv[]) {
40+
std::cout << "=== IRREGULAR WAVES ETA CONSISTENCY TEST ===" << std::endl;
41+
std::cout << "Validates spectrum-generated and eta-imported waves produce identical results.\n" << std::endl;
42+
43+
// Parse CLI arguments
44+
bool profilingOn = true;
45+
bool saveDataOn = true;
46+
bool visualizationOn = false;
47+
std::string data_dir;
48+
if (!hydroc::GetCLIArguments(argc, argv, "Sphere irregular waves eta consistency test",
49+
saveDataOn, profilingOn, visualizationOn, data_dir))
50+
return 1;
51+
if (!hydroc::SetInitialEnvironment(data_dir))
52+
return 1;
53+
54+
std::filesystem::path DATADIR(hydroc::getDataDir());
55+
auto meshfname = (DATADIR / "demos" / "sphere" / "geometry" / "oes_task10_sphere.obj")
56+
.lexically_normal().generic_string();
57+
auto h5fname = (DATADIR / "demos" / "sphere" / "hydroData" / "sphere.h5")
58+
.lexically_normal().generic_string();
59+
60+
// Output directory setup
61+
std::string out_dir = hydroc::getTestOutDir();
62+
std::filesystem::create_directories(out_dir + "/" + RESULTS_DIR_NAME);
63+
std::string eta_file = out_dir + "/" + RESULTS_DIR_NAME + "/temp_eta.txt";
64+
65+
std::vector<double> heave_spectrum;
66+
std::vector<double> heave_eta;
67+
68+
// ========== PHASE 1: Run with spectrum-generated waves ==========
69+
std::cout << "Phase 1: Running simulation with spectrum-generated waves..." << std::endl;
70+
{
71+
ChSystemNSC system;
72+
system.SetGravitationalAcceleration(ChVector3d(0.0, 0.0, -9.81));
73+
system.SetSolverType(ChSolver::Type::GMRES);
74+
system.GetSolver()->AsIterative()->SetMaxIterations(300);
75+
76+
// Ground
77+
auto ground = chrono_types::make_shared<ChBody>();
78+
system.AddBody(ground);
79+
ground->SetPos(ChVector3d(0, 0, -5));
80+
ground->SetFixed(true);
81+
ground->EnableCollision(false);
82+
83+
// Sphere body
84+
auto sphereBody = chrono_types::make_shared<ChBodyEasyMesh>(
85+
meshfname, 1000, false, true, false);
86+
system.Add(sphereBody);
87+
sphereBody->SetName("body1");
88+
sphereBody->SetPos(ChVector3d(0, 0, -2));
89+
sphereBody->SetMass(261.8e3);
90+
91+
// Prismatic joint (heave only)
92+
auto prismatic = chrono_types::make_shared<ChLinkLockPrismatic>();
93+
prismatic->Initialize(sphereBody, ground, false,
94+
ChFramed(ChVector3d(0, 0, -2)),
95+
ChFramed(ChVector3d(0, 0, -5)));
96+
system.AddLink(prismatic);
97+
98+
// Spring (zero stiffness/damping)
99+
auto spring = chrono_types::make_shared<ChLinkTSDA>();
100+
spring->Initialize(sphereBody, ground, false,
101+
ChVector3d(0, 0, -2), ChVector3d(0, 0, -5));
102+
spring->SetSpringCoefficient(0.0);
103+
spring->SetDampingCoefficient(0.0);
104+
system.AddLink(spring);
105+
106+
// Create spectrum-based waves
107+
IrregularWaveParams spectrum_params;
108+
spectrum_params.num_bodies_ = 1;
109+
spectrum_params.simulation_dt_ = TIMESTEP;
110+
spectrum_params.simulation_duration_ = DURATION;
111+
spectrum_params.ramp_duration_ = RAMP_DURATION;
112+
spectrum_params.wave_height_ = WAVE_HEIGHT;
113+
spectrum_params.wave_period_ = WAVE_PERIOD;
114+
spectrum_params.frequency_min_ = 0.001;
115+
spectrum_params.frequency_max_ = 1.0;
116+
spectrum_params.nfrequencies_ = 1000;
117+
spectrum_params.seed_ = SEED;
118+
119+
auto spectrum_waves = std::make_shared<IrregularWaves>(spectrum_params);
120+
121+
// Setup hydro forces - this initializes the waves via AddH5Data
122+
std::vector<std::shared_ptr<ChBody>> bodies = {sphereBody};
123+
HydroForces hydro_forces(bodies, h5fname);
124+
hydro_forces.AddWaves(spectrum_waves);
125+
126+
// NOW we can get the free surface data (after AddWaves initialized it)
127+
auto fse_time = spectrum_waves->GetFreeSurfaceTime();
128+
auto fse_elevation = spectrum_waves->GetFreeSurfaceElevation();
129+
130+
std::cout << " Generated " << fse_time.size() << " free surface elevation samples." << std::endl;
131+
if (!fse_time.empty()) {
132+
std::cout << " Time range: [" << fse_time.front() << ", " << fse_time.back() << "]" << std::endl;
133+
}
134+
135+
// Export eta data to file
136+
std::cout << " Exporting eta data to: " << eta_file << std::endl;
137+
{
138+
std::ofstream eta_out(eta_file);
139+
if (!eta_out.is_open()) {
140+
std::cerr << "ERROR: Could not create eta file: " << eta_file << std::endl;
141+
return 1;
142+
}
143+
eta_out << std::setprecision(10);
144+
for (size_t i = 0; i < fse_time.size(); ++i) {
145+
eta_out << fse_time[i] << ":" << fse_elevation[i] << "\n";
146+
}
147+
eta_out.close();
148+
}
149+
150+
// Run simulation
151+
while (system.GetChTime() <= DURATION) {
152+
system.DoStepDynamics(TIMESTEP);
153+
heave_spectrum.push_back(sphereBody->GetPos().z());
154+
}
155+
std::cout << " Completed. " << heave_spectrum.size() << " timesteps." << std::endl;
156+
}
157+
158+
// ========== PHASE 2: Run with eta-file-imported waves ==========
159+
std::cout << "\nPhase 2: Running simulation with eta-file-imported waves..." << std::endl;
160+
{
161+
ChSystemNSC system;
162+
system.SetGravitationalAcceleration(ChVector3d(0.0, 0.0, -9.81));
163+
system.SetSolverType(ChSolver::Type::GMRES);
164+
system.GetSolver()->AsIterative()->SetMaxIterations(300);
165+
166+
// Ground
167+
auto ground = chrono_types::make_shared<ChBody>();
168+
system.AddBody(ground);
169+
ground->SetPos(ChVector3d(0, 0, -5));
170+
ground->SetFixed(true);
171+
ground->EnableCollision(false);
172+
173+
// Sphere body
174+
auto sphereBody = chrono_types::make_shared<ChBodyEasyMesh>(
175+
meshfname, 1000, false, true, false);
176+
system.Add(sphereBody);
177+
sphereBody->SetName("body1");
178+
sphereBody->SetPos(ChVector3d(0, 0, -2));
179+
sphereBody->SetMass(261.8e3);
180+
181+
// Prismatic joint (heave only)
182+
auto prismatic = chrono_types::make_shared<ChLinkLockPrismatic>();
183+
prismatic->Initialize(sphereBody, ground, false,
184+
ChFramed(ChVector3d(0, 0, -2)),
185+
ChFramed(ChVector3d(0, 0, -5)));
186+
system.AddLink(prismatic);
187+
188+
// Spring (zero stiffness/damping)
189+
auto spring = chrono_types::make_shared<ChLinkTSDA>();
190+
spring->Initialize(sphereBody, ground, false,
191+
ChVector3d(0, 0, -2), ChVector3d(0, 0, -5));
192+
spring->SetSpringCoefficient(0.0);
193+
spring->SetDampingCoefficient(0.0);
194+
system.AddLink(spring);
195+
196+
// Create eta-import waves
197+
IrregularWaveParams eta_params;
198+
eta_params.num_bodies_ = 1;
199+
eta_params.simulation_dt_ = TIMESTEP;
200+
eta_params.simulation_duration_ = DURATION;
201+
eta_params.ramp_duration_ = RAMP_DURATION;
202+
eta_params.eta_file_path_ = eta_file;
203+
eta_params.frequency_min_ = 0.001;
204+
eta_params.frequency_max_ = 1.0;
205+
eta_params.nfrequencies_ = 1000;
206+
207+
auto eta_waves = std::make_shared<IrregularWaves>(eta_params);
208+
209+
// Setup hydro forces
210+
std::vector<std::shared_ptr<ChBody>> bodies = {sphereBody};
211+
HydroForces hydro_forces(bodies, h5fname);
212+
hydro_forces.AddWaves(eta_waves);
213+
214+
// Run simulation
215+
while (system.GetChTime() <= DURATION) {
216+
system.DoStepDynamics(TIMESTEP);
217+
heave_eta.push_back(sphereBody->GetPos().z());
218+
}
219+
std::cout << " Completed. " << heave_eta.size() << " timesteps." << std::endl;
220+
}
221+
222+
// ========== PHASE 3: Compare results ==========
223+
std::cout << "\nPhase 3: Comparing results..." << std::endl;
224+
225+
if (heave_spectrum.size() != heave_eta.size()) {
226+
std::cerr << "ERROR: Different number of timesteps! Spectrum: "
227+
<< heave_spectrum.size() << ", Eta: " << heave_eta.size() << std::endl;
228+
return 1;
229+
}
230+
231+
double max_diff = 0.0;
232+
double sum_diff_sq = 0.0;
233+
size_t max_diff_idx = 0;
234+
235+
for (size_t i = 0; i < heave_spectrum.size(); ++i) {
236+
double diff = std::abs(heave_spectrum[i] - heave_eta[i]);
237+
sum_diff_sq += diff * diff;
238+
if (diff > max_diff) {
239+
max_diff = diff;
240+
max_diff_idx = i;
241+
}
242+
}
243+
244+
double rms_diff = std::sqrt(sum_diff_sq / heave_spectrum.size());
245+
double time_at_max = max_diff_idx * TIMESTEP;
246+
247+
std::cout << " Max difference: " << max_diff << " m (at t=" << time_at_max << "s)" << std::endl;
248+
std::cout << " RMS difference: " << rms_diff << " m" << std::endl;
249+
250+
// ========== PHASE 4: Save results ==========
251+
if (saveDataOn) {
252+
std::string results_file = out_dir + "/" + RESULTS_DIR_NAME + "/" + RESULTS_FILE_NAME + ".txt";
253+
std::ofstream out(results_file);
254+
if (out.is_open()) {
255+
out << std::setprecision(10);
256+
out << "Time (s) Heave_Spectrum (m) Heave_Eta (m) Difference (m)\n";
257+
out << "-------- ------------------ ------------- --------------\n";
258+
for (size_t i = 0; i < heave_spectrum.size(); ++i) {
259+
double t = i * TIMESTEP;
260+
out << std::setw(10) << std::fixed << std::setprecision(3) << t
261+
<< std::setw(20) << std::fixed << std::setprecision(8) << heave_spectrum[i]
262+
<< std::setw(18) << std::fixed << std::setprecision(8) << heave_eta[i]
263+
<< std::setw(18) << std::scientific << std::setprecision(4)
264+
<< (heave_spectrum[i] - heave_eta[i]) << "\n";
265+
}
266+
out.close();
267+
std::cout << "\n Results saved to: " << results_file << std::endl;
268+
}
269+
}
270+
271+
// ========== PASS/FAIL determination ==========
272+
const double tolerance = 1e-6; // 1 micrometer tolerance
273+
bool passed = (max_diff < tolerance);
274+
275+
std::cout << "\n=== TEST " << (passed ? "PASSED" : "FAILED") << " ===" << std::endl;
276+
if (!passed) {
277+
std::cerr << "Spectrum and eta-imported waves produced different results!" << std::endl;
278+
std::cerr << "Max difference " << max_diff << " exceeds tolerance " << tolerance << std::endl;
279+
return 1;
280+
}
281+
282+
std::cout << "Spectrum-generated and eta-imported waves produce identical results." << std::endl;
283+
return 0;
284+
}

0 commit comments

Comments
 (0)