Modify test_rm3_decay to perform validation against reference file and plot validation results using GnuPlot

Author: rserban

CMakeLists.txt

@@ -431,6 +431,11 @@ target_include_directories(HydroChronoGUI
         "$<INSTALL_INTERFACE:include>"
 )
 
+if(MSVC)
+  target_compile_options(HydroChronoGUI PRIVATE $<$<COMPILE_LANGUAGE:CXX>:/wd4996>)  # deprecated function or class member
+  target_compile_options(HydroChronoGUI PRIVATE $<$<COMPILE_LANGUAGE:CXX>:/wd4458>)  # declaration hides class member
+endif()
+
 set_target_properties(HydroChronoGUI PROPERTIES POSITION_INDEPENDENT_CODE ON)
 target_link_libraries(HydroChronoGUI
     PUBLIC

data/reference_data/rm3/hc_ref_rm3_decay.txt

@@ -1,4 +1,4 @@
-Time (s)   Float Heave (m) Plate Heave (m)
+Time(s)        Float_Heave(m) Plate_Heave(m)
 0.01           -0.62000515    -21.29000007
 0.02           -0.62002233    -21.29000032
 0.03           -0.62005382    -21.29000078

tests/regression/CMakeLists.txt

@@ -22,6 +22,7 @@ function(build_regression_test TEST_GROUP TEST_NAME MULTIPLE)
 
    target_compile_definitions(${TEST_EXE} 
 	   PRIVATE
+       "REFERENCE_FILE_NAME=\"${FILE_REF}\""
        "RESULTS_DIR_NAME=\"${TEST_GROUP}\""
        "RESULTS_FILE_NAME=\"results_${TEST_NAME}\""
    )

tests/regression/rm3/test_rm3_decay.cpp

@@ -6,14 +6,17 @@
 #include <chrono/physics/ChBodyEasy.h>
 #include <chrono/physics/ChSystemNSC.h>
 
-#include "chrono_postprocess/ChGnuPlot.h"
+#include <chrono/utils/ChUtils.h>
+
+#include "../test_utils.h"
 
 #include <chrono>   // std::chrono::high_resolution_clock::now
 #include <iomanip>  // std::setprecision
 #include <vector>   // std::vector<double>
 
 // Use the namespaces of Chrono
 using namespace chrono;
+using namespace chrono::utils;
 
 int main(int argc, char* argv[]) {
     std::cout << "Chrono version: " << CHRONO_VERSION << "\n\n";
@@ -53,14 +56,8 @@ int main(int argc, char* argv[]) {
 
     // Create user interface
     std::shared_ptr<hydroc::gui::UI> pui = hydroc::gui::CreateUI(visualizationOn);
-
     hydroc::gui::UI& ui = *pui.get();
 
-    // some io/viz options
-    std::vector<double> time_vector;
-    std::vector<double> float_heave_position;
-    std::vector<double> plate_heave_position;
-
     // set up body from a mesh
     std::cout << "Attempting to open mesh file: " << body1_meshfame << std::endl;
     std::shared_ptr<ChBody> float_body1 = chrono_types::make_shared<ChBodyEasyMesh>(  //
@@ -127,10 +124,10 @@ int main(int argc, char* argv[]) {
 
     HydroForces hydroForces(bodies, h5fname, default_dont_add_waves);
 
-    //// Debug printing added mass matrix and system mass matrix
-    // ChSparseMatrix M;
-    // system.GetMassMatrix(&M);
-    // std::cout << M << std::endl;
+    // Result arrays
+    std::vector<double> time_vector;
+    std::vector<double> float_heave_position;
+    std::vector<double> plate_heave_position;
 
     // for profiling
     auto start = std::chrono::high_resolution_clock::now();
@@ -190,14 +187,39 @@ int main(int argc, char* argv[]) {
         }
     }
 
+    // Read reference data
+    ChValidation::Headers col_headers;
+    ChValidation::Data ref_data = ChValidation::ReadDataFile(REFERENCE_FILE_NAME, col_headers);
+    ChAssertAlways(ref_data.size() == 3);
+
+    // Simulation data
+    ChValidation::Data res_data = ChValidation::CreateData({time_vector, float_heave_position, plate_heave_position});
+
+    // Perform validation
+    auto norm_type   = ChValidation::NormType::RMS;
+    double tolerance = 1e-8;
+    ChValidation::DataVector error_norms;
+    bool passed = ChValidation::Test(res_data, ref_data, ChValidation::NormType::RMS, tolerance, error_norms);
+
+    std::cout << "\nValidation";
+    std::cout << "\n  Reference file:   " << REFERENCE_FILE_NAME;
+    std::cout << "\n  Data series:      ";
+    for (const auto& c : col_headers)
+        std::cout << c << "  ";
+    std::cout << "\n  Ref. data points: " << ref_data[0].size();
+    std::cout << "\n  Sim. data points: " << res_data[0].size();
+    std::cout << "\n  Validation norm:  " << ChValidation::GetNormTypeAsString(norm_type);
+    std::cout << "\n  Tolerance:        " << tolerance;
+    std::cout << "\n  " << (passed ? "Passed" : "Failed");
+    std::cout << "            [ ";
+    for (const auto& nrm : error_norms)
+        std::cout << nrm << " ";
+    std::cout << "]" << std::endl;
+
+    // Plot simulation and reference results
     if (plotOn) {
-        postprocess::ChGnuPlot gplot(out_dir + "/rm3_decay.gpl");
-        gplot.SetGrid();
-        gplot.SetLabelX("time (s)");
-        gplot.SetLabelY("heave (m)");
-        gplot.SetTitle("RM3 decay");
-        gplot.Plot(time_vector, plate_heave_position, "", " with lines lt rgb '#FF5500' lw 2");
+        PlotValidation(out_dir + "/rm3_decay.gpl", "RM3 decay", col_headers, ref_data, res_data, simulationDuration);
     }
 
-    return 0;
-}
+    return !passed;
+}

tests/regression/test_utils.h

@@ -0,0 +1,45 @@
+#ifndef TEST_UTILS_H
+#define TEST_UTILS_H
+
+#include <chrono/utils/ChValidation.h>
+#include <chrono_postprocess/ChGnuPlot.h>
+
+void PlotValidation(const std::string& out_script,
+                    const std::string& title,
+                    const chrono::utils::ChValidation::Headers& headers,
+                    chrono::utils::ChValidation::Data& ref_data,
+                    chrono::utils::ChValidation::Data& res_data,
+                    double simulationDuration) {
+    int n_plots = (int)res_data.size() - 1;
+
+    chrono::postprocess::ChGnuPlot gplot(out_script);
+    gplot.SetCanvasSize(1200, 500 * n_plots);
+    gplot.SetTitle("{/:Bold " + title + "}");
+    gplot.SetCommand("set multiplot layout " + std::to_string(n_plots) + ", 1");
+
+    for (int i = 1; i <= n_plots; i++) {
+        gplot.SetRangeX(0, simulationDuration);
+        gplot.SetCommand("set xlabel '" + headers[0] + "' noenhanced");
+
+        gplot.SetCommand("set ylabel '" + headers[i] + "' noenhanced");        
+        gplot.Plot(ref_data[0], ref_data[i], "Reference", "with lines lt rgb '#0055FF' lw 4");
+        gplot.Plot(res_data[0], res_data[i], "Simulation", "with lines lt rgb '#FF3300' lw 4 dashtype 2");
+
+        auto n_err                      = std::min(res_data[0].size(), ref_data[0].size());
+        chrono::ChVectorDynamic<> t_err = ref_data[0].head(n_err);
+        chrono::ChVectorDynamic<> err   = ref_data[i].head(n_err) - res_data[i].head(n_err);
+        if (err.norm() < 1e-4) {
+            gplot.SetRangeY2(-0.01, +0.01);
+            gplot.SetCommand("set y2tics -0.01 , 0.01, 0.01");
+        }
+        gplot.SetLabelY2("error");
+        gplot.Plot(t_err, err, "Error (ref - sim)", "axes x1y2 with lines lt rgb '#44CC44' lw 4");
+
+        gplot.SetLegend("box lt -1 lw 1");
+        gplot.SetGrid();
+
+        gplot.FlushPlots();
+    }
+}
+
+#endif