Included an auxiliary function to check data

GiuseppeOrlando878776 · GiuseppeOrlando878776 · commit 2351bbd6b87d · 2025-10-31T16:02:18.000+01:00
diff --git a/1D_RIEMANN/All_topology_solver.cpp b/1D_RIEMANN/All_topology_solver.cpp
@@ -120,26 +120,22 @@ int main(int argc, char* argv[]) {
   app.add_option("--p1L", Riemann_param.p1L, "Initial pressure phase 1 at left")->capture_default_str()->group("Initial conditions");
   app.add_option("--T1L", Riemann_param.T1L, "Initial temperature phase 1 at left")->capture_default_str()->group("Initial conditions");
   app.add_option("--u1L", Riemann_param.u1L, "Initial horizontal velocity phase 1 at left")->capture_default_str()->group("Initial conditions");
-  app.add_option("--v1L", Riemann_param.v1L, "Initial vertical velocity phase 1 at left")->capture_default_str()->group("Initial conditions");
 
   app.add_option("--rho2L", Riemann_param.rho2L, "Initial density phase 2 at left")->capture_default_str()->group("Initial conditions");
   app.add_option("--p2L", Riemann_param.p2L, "Initial pressure phase 2 at left")->capture_default_str()->group("Initial conditions");
   app.add_option("--T2L", Riemann_param.T2L, "Initial temperature phase 2 at left")->capture_default_str()->group("Initial conditions");
   app.add_option("--u2L", Riemann_param.u2L, "Initial horizontal velocity phase 2 at left")->capture_default_str()->group("Initial conditions");
-  app.add_option("--v2L", Riemann_param.v2L, "Initial vertical velocity phase 2 at left")->capture_default_str()->group("Initial conditions");
 
   app.add_option("--alpha1R", Riemann_param.alpha1R, "Initial volume fraction at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--rho1R", Riemann_param.rho1R, "Initial density phase 1 at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--p1R", Riemann_param.p1R, "Initial pressure phase 1 at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--T1R", Riemann_param.T1R, "Initial temperature phase 1 at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--u1R", Riemann_param.u1R, "Initial horizotnal velocity phase 1 at right")->capture_default_str()->group("Initial conditions");
-  app.add_option("--v1R", Riemann_param.v1R, "Initial vertical velocity phase 1 at right")->capture_default_str()->group("Initial conditions");
 
   app.add_option("--rho2R", Riemann_param.rho2R, "Initial density phase 2 at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--p2R", Riemann_param.p2R, "Initial pressure phase 2 at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--T2R", Riemann_param.T2R, "Initial temperature phase 2 at right")->capture_default_str()->group("Initial conditions");
   app.add_option("--u2R", Riemann_param.u2R, "Initial horizontal velocity phase 2 at right")->capture_default_str()->group("Initial conditions");
-  app.add_option("--v2R", Riemann_param.v2R, "Initial vertical velocity phase 2 at right")->capture_default_str()->group("Initial conditions");
 
   /*--- Create the instance of the class to perform the simulation ---*/
   CLI11_PARSE(app, argc, argv);
diff --git a/1D_RIEMANN/test_configuration/All_topology_solver.hpp b/1D_RIEMANN/test_configuration/All_topology_solver.hpp
@@ -21,6 +21,9 @@ namespace fs = std::filesystem;
 #include "schemes/Rusanov_flux.hpp"
 #include "schemes/non_conservative_flux.hpp"
 
+/*--- Define preprocessor to check whether to control data or not ---*/
+#define VERBOSE
+
 // This is the class for the simulation of the all-topology model
 //
 template<std::size_t dim>
@@ -118,6 +121,8 @@ class All_Topology_Solver {
 
   Number get_max_lambda() const; /*--- Compute the estimate of the maximum eigenvalue ---*/
 
+  void check_data(); /*--- Auxiliary routine to check if (small) spurious negative values are present ---*/
+
   void update_auxiliary_fields(); /*--- Routine to update auxiliary fields for output and time step update ---*/
 };
 
@@ -424,6 +429,57 @@ All_Topology_Solver<dim>::get_max_lambda() const {
   return global_res;
 }
 
+// Auxiliary routine to check if spurious negative values arise
+//
+template<std::size_t dim>
+void All_Topology_Solver<dim>::check_data() {
+  samurai::for_each_cell(mesh,
+                         [&](const auto& cell)
+                            {
+                              // Start with the volume fraction
+                              if(conserved_variables[cell][Indices::ALPHA1_INDEX] < static_cast<Number>(0.0)) {
+                                std::cerr << "Negative volume fraction for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(conserved_variables[cell][Indices::ALPHA1_INDEX] > static_cast<Number>(1.0)) {
+                                std::cerr << "Exceeding volume fraction for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_INDEX])) {
+                                std::cerr << "NaN volume fraction for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+
+                              // Sanity check for m1
+                              if(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX] < static_cast<Number>(0.0)) {
+                                std::cerr << "Negative mass for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX])) {
+                                std::cerr << "NaN mass for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+
+                              // Sanity check for m2
+                              if(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX] < static_cast<Number>(0.0)) {
+                                std::cerr << "Negative mass for phase 2" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(std::isnan(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX])){
+                                std::cerr << "NaN mass for phase 2" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                            }
+                        );
+}
+
 // Update auxiliary fields after solution of the system
 //
 template<std::size_t dim>
@@ -589,10 +645,20 @@ void All_Topology_Solver<dim>::run(const std::size_t nfiles) {
     }
 
     // Apply the numerical scheme
-    conserved_variables_np1 = conserved_variables
-                            - dt*Rusanov_flux(conserved_variables)
-                            - dt*NonConservative_flux(conserved_variables);
+    try {
+      conserved_variables_np1 = conserved_variables
+                              - dt*Rusanov_flux(conserved_variables)
+                              - dt*NonConservative_flux(conserved_variables);
+    }
+    catch(const std::exception& e) {
+      std::cerr << e.what() << std::endl;
+      save("_diverged", conserved_variables);
+      exit(1);
+    }
     samurai::swap(conserved_variables, conserved_variables_np1);
+    #ifdef VERBOSE
+      check_data();
+    #endif
 
     // Save the results
     if(t >= static_cast<Number>(nsave + 1)*dt_save || t == Tf) {
diff --git a/1D_RIEMANN/test_configuration/containers.hpp b/1D_RIEMANN/test_configuration/containers.hpp
@@ -60,23 +60,19 @@ struct Riemann_Parameters {
   T p1L;
   T T1L;
   T u1L;
-  T v1L;
   T rho2L;
   T p2L;
   T T2L;
   T u2L;
-  T v2L;
 
   /*--- Right state ---*/
   T alpha1R;
   T rho1R;
   T p1R;
   T T1R;
   T u1R;
-  T v1R;
   T rho2R;
   T p2R;
   T T2R;
   T u2R;
-  T v2R;
 };
diff --git a/2D_RIEMANN/test_configuration/All_topology_solver.hpp b/2D_RIEMANN/test_configuration/All_topology_solver.hpp
@@ -24,6 +24,9 @@ namespace fs = std::filesystem;
 #include "schemes/Rusanov_flux.hpp"
 #include "schemes/non_conservative_flux.hpp"
 
+/*--- Define preprocessor to check whether to control data or not ---*/
+#define VERBOSE
+
 // This is the class for the simulation of the all-topology model
 //
 template<std::size_t dim>
@@ -118,6 +121,8 @@ class All_Topology_Solver {
 
   Number get_max_lambda() const; /*--- Compute the estimate of the maximum eigenvalue ---*/
 
+  void check_data(); /*--- Auxiliary routine to check if (small) spurious negative values are present ---*/
+
   void update_auxiliary_fields(); /*--- Routine to update auxiliary fields for output and time step update ---*/
 };
 
@@ -433,6 +438,57 @@ All_Topology_Solver<dim>::get_max_lambda() const {
   return global_res;
 }
 
+// Auxiliary routine to check if spurious negative values arise
+//
+template<std::size_t dim>
+void All_Topology_Solver<dim>::check_data() {
+  samurai::for_each_cell(mesh,
+                         [&](const auto& cell)
+                            {
+                              // Start with the volume fraction
+                              if(conserved_variables[cell][Indices::ALPHA1_INDEX] < static_cast<Number>(0.0)) {
+                                std::cerr << "Negative volume fraction for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(conserved_variables[cell][Indices::ALPHA1_INDEX] > static_cast<Number>(1.0)) {
+                                std::cerr << "Exceeding volume fraction for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_INDEX])) {
+                                std::cerr << "NaN volume fraction for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+
+                              // Sanity check for m1
+                              if(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX] < static_cast<Number>(0.0)) {
+                                std::cerr << "Negative mass for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX])) {
+                                std::cerr << "NaN mass for phase 1" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+
+                              // Sanity check for m2
+                              if(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX] < static_cast<Number>(0.0)) {
+                                std::cerr << "Negative mass for phase 2" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                              else if(std::isnan(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX])){
+                                std::cerr << "NaN mass for phase 2" << std::endl;
+                                save("_diverged", conserved_variables);
+                                exit(1);
+                              }
+                            }
+                        );
+}
+
 // Update auxiliary fields after solution of the system
 //
 template<std::size_t dim>
@@ -595,10 +651,20 @@ void All_Topology_Solver<dim>::run(const std::size_t nfiles) {
     }
 
     // Apply the numerical scheme
-    conserved_variables_np1 = conserved_variables
-                            - dt*Rusanov_flux(conserved_variables)
-                            - dt*NonConservative_flux(conserved_variables);
+    try {
+      conserved_variables_np1 = conserved_variables
+                              - dt*Rusanov_flux(conserved_variables)
+                              - dt*NonConservative_flux(conserved_variables);
+    }
+    catch(const std::exception& e) {
+      std::cerr << e.what() << std::endl;
+      save("_diverged", conserved_variables);
+      exit(1);
+    }
     samurai::swap(conserved_variables, conserved_variables_np1);
+    #ifdef VERBOSE
+      check_data();
+    #endif
 
     // Save the results
     if(t >= static_cast<Number>(nsave + 1)*dt_save || t == Tf) {
diff --git a/include/schemes/Rusanov_flux.hpp b/include/schemes/Rusanov_flux.hpp
@@ -8,6 +8,8 @@
 
 #include "flux_base.hpp"
 
+#define VERBOSE_FLUX
+
 namespace samurai {
   /**
     * Implementation of a Rusanov flux
@@ -56,6 +58,22 @@ namespace samurai {
     const auto m2_L     = qL(Indices::ALPHA2_RHO2_INDEX);
     const auto m2E2_L   = qL(Indices::ALPHA2_RHO2_E2_INDEX);
 
+    // Verify if it is admissible
+    #ifdef VERBOSE_FLUX
+      if(m1_L < static_cast<Number>(0.0)) {
+        throw std::runtime_error(std::string("Negative mass phase 1 left state: " + std::to_string(m1_L)));
+      }
+      if(m2_L < static_cast<Number>(0.0)) {
+        throw std::runtime_error(std::string("Negative mass phase 2 left state: " + std::to_string(m2_L)));
+      }
+      if(alpha1_L < static_cast<Number>(0.0)) {
+        throw std::runtime_error(std::string("Negative volume fraction phase 1 left state: " + std::to_string(alpha1_L)));
+      }
+      else if(alpha1_L > static_cast<Number>(1.0)) {
+        throw std::runtime_error(std::string("Exceeding volume fraction phase 1 left state: " + std::to_string(alpha1_L)));
+      }
+    #endif
+
     // Phase 1
     const auto inv_m1_L = static_cast<Number>(1.0)/m1_L; /*--- TODO: Add treatment for vanishing volume fraction ---*/
     const auto vel1_L_d = qL(Indices::ALPHA1_RHO1_U1_INDEX + curr_d)*inv_m1_L; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -68,6 +86,11 @@ namespace samurai {
     }
     const auto p1_L = this->EOS_phase1.pres_value_Rhoe(rho1_L, e1_L);
     const auto c1_L = this->EOS_phase1.c_value_RhoP(rho1_L, p1_L);
+    #ifdef VERBOSE_FLUX
+      if(std::isnan(c1_L)) {
+        throw std::runtime_error(std::string("NaN speed of sound phase 1 left state"));
+      }
+    #endif
 
     // Phase 2
     const auto inv_m2_L = static_cast<Number>(1.0)/m2_L; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -81,6 +104,11 @@ namespace samurai {
     }
     const auto p2_L = this->EOS_phase2.pres_value_Rhoe(rho2_L, e2_L);
     const auto c2_L = this->EOS_phase2.c_value_RhoP(rho2_L, p2_L);
+    #ifdef VERBOSE_FLUX
+      if(std::isnan(c2_L)) {
+        throw std::runtime_error(std::string("NaN speed of sound phase 2 left state"));
+      }
+    #endif
 
     /*--- Right state ---*/
     // Pre-fetch variables that will be used several times so as to exploit possible vectorization
@@ -91,6 +119,22 @@ namespace samurai {
     const auto m2_R     = qR(Indices::ALPHA2_RHO2_INDEX);
     const auto m2E2_R   = qR(Indices::ALPHA2_RHO2_E2_INDEX);
 
+    // Verify if it is admissible
+    #ifdef VERBOSE_FLUX
+      if(m1_R < static_cast<Number>(0.0)) {
+        throw std::runtime_error(std::string("Negative mass phase 1 right state: " + std::to_string(m1_R)));
+      }
+      if(m2_R < static_cast<Number>(0.0)) {
+        throw std::runtime_error(std::string("Negative mass phase 2 right state: " + std::to_string(m2_R)));
+      }
+      if(alpha1_R < static_cast<Number>(0.0)) {
+        throw std::runtime_error(std::string("Negative volume fraction phase 1 right state: " + std::to_string(alpha1_R)));
+      }
+      else if(alpha1_R > static_cast<Number>(1.0)) {
+        throw std::runtime_error(std::string("Exceeding volume fraction phase 1 right state: " + std::to_string(alpha1_R)));
+      }
+    #endif
+
     // Phase 1
     const auto inv_m1_R = static_cast<Number>(1.0)/m1_R; /*--- TODO: Add treatment for vanishing volume fraction ---*/
     const auto vel1_R_d = qR(Indices::ALPHA1_RHO1_U1_INDEX + curr_d)*inv_m1_R; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -103,6 +147,11 @@ namespace samurai {
     }
     const auto p1_R = this->EOS_phase1.pres_value_Rhoe(rho1_R, e1_R);
     const auto c1_R = this->EOS_phase1.c_value_RhoP(rho1_R, p1_R);
+    #ifdef VERBOSE_FLUX
+      if(std::isnan(c1_R)) {
+        throw std::runtime_error(std::string("NaN speed of sound phase 1 right state"));
+      }
+    #endif
 
     // Phase 2
     const auto inv_m2_R = static_cast<Number>(1.0)/m2_R; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -116,6 +165,11 @@ namespace samurai {
     }
     const auto p2_R = this->EOS_phase2.pres_value_Rhoe(rho2_R, e2_R);
     const auto c2_R = this->EOS_phase2.c_value_RhoP(rho2_R, p2_R);
+    #ifdef VERBOSE_FLUX
+      if(std::isnan(c2_R)) {
+        throw std::runtime_error(std::string("NaN speed of sound phase 2 right state"));
+      }
+    #endif
 
     /*--- Compute the flux ---*/
     const auto lambda = std::max(std::max(std::abs(vel1_L_d) + c1_L, std::abs(vel1_R_d) + c1_R),
diff --git a/include/schemes/non_conservative_flux.hpp b/include/schemes/non_conservative_flux.hpp
@@ -6,6 +6,8 @@
 //
 #pragma once
 
+#include "flux_base.hpp"
+
 namespace samurai {
   /**
     * Implementation of the non-conservative flux
