Nan check

include/aether.h

@@ -14,6 +14,9 @@
 /// This is for generating random numbers:
 #include <random>
 
+//This is to look into bits for flagging nans and infs
+#include <cstdint>
+
 // Types
 // Precision compile-time aliasing
 #ifdef AETHER_USE_PRECISION_DOUBLE

include/inputs.h

@@ -56,6 +56,10 @@ class Inputs {
   bool get_do_lat_dependent_radius();
   bool get_do_J2();
 
+  bool get_check_for_nans();
+  bool get_nan_test();
+  std::string get_nan_test_variable();
+
   bool get_is_cubesphere();
 
   bool get_NO_cooling();

include/ions.h

@@ -101,7 +101,8 @@ class Ions {
 							     Grid grid);
   void calc_ion_temperature(Neutrals neutrals, Grid grid, Times time);
   void calc_electron_temperature(Neutrals neutrals, Grid grid);
-
+  bool check_for_nonfinites();
+  void nan_test(std::string variable);
   bool restart_file(std::string dir, bool DoRead);
 };
 #endif  // INCLUDE_IONS_H_

include/neutrals.h

@@ -372,6 +372,16 @@ class Neutrals {
    **/
   int get_species_id(std::string name);
 
+  /*****************************************************************************
+      \brief  Checks for nans and +/- infinities in density, temp, and velocity
+    **/
+  bool check_for_nonfinites();
+
+  /**********************************************************************
+      \brief Checks for nans in the specified variable
+    **/
+  void nan_test(std::string variable);
+
   /**********************************************************************
      \brief Read / Write restart files for the neutral variables
      \param dir directory to write restart files

include/tools.h

@@ -4,72 +4,114 @@
 #ifndef INCLUDE_TOOLS_H_
 #define INCLUDE_TOOLS_H_
 
-// ----------------------------------------------------------------------------
-// 
-// ----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
+// Display an armadillo vector
+// ----------------------------------------------------------------------
 
 void display_vector(arma_vec vec);
 
-// -----------------------------------------------------------------------------
-// synchronize a variable across all processors
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
+// synchronize a (boolean) variable across all processors
+// ----------------------------------------------------------------------
 
 bool sync_across_all_procs(bool value);
 
-// ----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 // Calculate the average value across all processors
-// ----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 
 precision_t sync_mean_across_all_procs(precision_t value);
 
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 // Generate a vector of normally distributed random doubles
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 
 std::vector<double> get_normal_random_vect(double mean,
 					   double std,
 					   int64_t nValues,
 					   int seed);
 
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 // Generate a vector of uniformly distributed random unsigned ints
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 
 std::vector<unsigned int> get_random_unsigned_vect(int64_t nValues,
 						   int seed);
 
+// ----------------------------------------------------------------------
+// Make a vector of arma cubes:
+// ----------------------------------------------------------------------
+
 std::vector<arma_cube> make_cube_vector(int64_t nLons,
-				    int64_t nLats,
-				    int64_t nAlts,
-				    int64_t nComps);
+					int64_t nLats,
+					int64_t nAlts,
+					int64_t nComps);
+
+// ----------------------------------------------------------------------
+// Take the dot product between two armadilo cubes
+// ----------------------------------------------------------------------
 
 arma_cube dot_product(std::vector<arma_cube> vec1,
 		      std::vector<arma_cube> vec2);
 
+// ----------------------------------------------------------------------
+// Take the cross product between two arma cubes
+// ----------------------------------------------------------------------
+
 std::vector<arma_cube> cross_product(std::vector<arma_cube> vec1,
 				     std::vector<arma_cube> vec2);
 
+// ----------------------------------------------------------------------
+// Convert an armadillo vector to a c++ vector
+// ----------------------------------------------------------------------
+
 std::vector<precision_t> make_vector_from_fvec(arma_vec in_fvec);
+
+// ----------------------------------------------------------------------
+// Convert a c++ vector to an armadillo vector
+// ----------------------------------------------------------------------
+
 arma_vec make_fvec_from_vector(std::vector<precision_t> in_vector);
 
+// ----------------------------------------------------------------------
+// Convert an integer to a 0-padded string
+// ----------------------------------------------------------------------
+
 std::string tostr(int64_t num_to_convert, int64_t zero_padding_len);
 
+// ----------------------------------------------------------------------
+// Convert a string to a number
+// ----------------------------------------------------------------------
+
 precision_t str_to_num(std::string input);
 
+// ----------------------------------------------------------------------
+// Read a json from a file
+// ----------------------------------------------------------------------
+
 json read_json(std::string json_file);
+
+// ----------------------------------------------------------------------
+// Write a json to a file
+// ----------------------------------------------------------------------
+
 bool write_json(std::string json_file, json json_output);
 
+// ----------------------------------------------------------------------
+// Compare two numbers and fail if the difference is too large
+// ----------------------------------------------------------------------
+
 bool compare(precision_t value1, precision_t value2);
 
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 // calculate mean of vector
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 
 precision_t mean(std::vector<precision_t> values);
 
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 // calculate standard deviation of vector
-// -----------------------------------------------------------------------------
+// ----------------------------------------------------------------------
 
 precision_t standard_deviation(std::vector<precision_t> values);
 
@@ -80,7 +122,8 @@ precision_t standard_deviation(std::vector<precision_t> values);
 std::vector<precision_t> get_min_mean_max(const arma_cube &value);
 
 //-------------------------------------------------------------
-// Find the name of given species in neutrals and ions. Throw exception if not found
+// Find the name of given species in neutrals and ions. Throw
+// exception if not found
 //-------------------------------------------------------------
 
 const arma_cube& find_species_density(const std::string &name,
@@ -95,4 +138,49 @@ std::vector<precision_t> get_min_mean_max_density(const std::string &name,
                                                   Neutrals &neutrals,
                                                   Ions &ions);
 
+//-----------------------------------------------------------------------
+// Checks if armacube(s) has all finite values, if not, adds them to
+// errors in report class
+//-----------------------------------------------------------------------
+
+bool all_finite(arma_cube cube, std::string name);
+bool all_finite(std::vector<arma_cube> cubes, std::string name);
+
+
+//-----------------------------------------------------------------------
+// Takes an index of an armacube and converts it to latitude,
+// longitude, and altitude
+//-----------------------------------------------------------------------
+
+std::vector<int> index_to_ijk(arma_cube cube, int index);
+
+//-----------------------------------------------------------------------
+// Inserts 3 nans and 3 inf into the specified armacube(s)
+//-----------------------------------------------------------------------
+
+std::vector<int> insert_indefinites(arma_cube &cube);
+
+//-----------------------------------------------------------------------
+// Returns a vector of ints as a string with spaces in between
+//-----------------------------------------------------------------------
+std::string print_nan_vector(std::vector<int> input, arma_cube cube);
+
+//-----------------------------------------------------------------------
+// Returns whether a given arma cube has all finite values
+//-----------------------------------------------------------------------
+
+bool is_finite(arma_cube &cube);
+
+//-----------------------------------------------------------------------
+//Returns whether the double/float value is a nan or infinity
+//-----------------------------------------------------------------------
+
+bool is_nan_inf(double value);
+
+//-----------------------------------------------------------------------
+// Returns vector of indefinite values
+//-----------------------------------------------------------------------
+
+std::vector<int> indef_vector(arma_cube cube);
+
 #endif  // INCLUDE_TOOLS_H_

share/run/UA/inputs/defaults.json

@@ -6,7 +6,14 @@
 	"dt" : 60.0,
 	"TimingPercent" : 1.0,
 	"iTimingDepth" : 5,
-	"iProc" : 0},
+	"iProc" : 0,
+        "iFunctionVerbose" : {
+            "Grid::create_altitudes": 0},	
+	"check_for_nans" : false,
+	"nan_test" : {
+            "insert" : false,
+            "variable" : "temperature_scgc"}
+    },
 
     "InitialConditions" : {
 	"type" : "Planet"},    

share/run/aether.json

@@ -1,14 +1,15 @@
 
 {
-
     "Ensembles" : {
-		"nMembers" : 1},
+	"nMembers" : 1},
 
-	"Debug" : {
+    "Debug" : {
 	"iVerbose" : 0,
         "iFunctionVerbose" : {
             "Grid::create_altitudes": 0},	
-	"dt" : 10.0},
+	"dt" : 10.0,
+	"check_for_nans" : false
+    },
 
     "EndTime" : [2011, 3, 20, 0, 10, 0],
 

src/advance.cpp

@@ -96,6 +96,7 @@ int advance(Planets &planet,
     time.restart_file(input.get_restartout_dir(), DoWrite);
   }
 
+
   iErr = output(neutrals, ions, gGrid, time, planet);
 
   logfile.write_logfile(indices, neutrals, ions, gGrid, time);

src/inputs.cpp

@@ -751,6 +751,38 @@ json Inputs::get_perturb_values() {
   return values;
 }
 
+
+// -----------------------------------------------------------------------
+// Flag to check neutral and ions for nans and infinites
+// -----------------------------------------------------------------------
+
+
+bool Inputs::get_check_for_nans() {
+  return settings.at("Debug").at("check_for_nans");;
+}
+
+
+
+
+// -----------------------------------------------------------------------
+// Checks to see if nan_test is needed
+// -----------------------------------------------------------------------
+
+
+bool Inputs::get_nan_test() {
+  return settings.at("Debug").at("nan_test").at("insert");
+}
+
+
+// -----------------------------------------------------------------------
+// Returns which variable is being tested for nans
+// -----------------------------------------------------------------------
+
+
+std::string Inputs::get_nan_test_variable() {
+  return settings.at("Debug").at("nan_test").at("variable");
+}
+
 // -----------------------------------------------------------------------
 // Flag to have a latitude dependent radius, and by extension gravity
 // -----------------------------------------------------------------------

src/ions.cpp

@@ -153,6 +153,49 @@ int Ions::read_planet_file(Planets planet) {
 }
 
 
+//----------------------------------------------------------------------
+// Reports location of nans inserted into specified variable
+//----------------------------------------------------------------------
+void Ions::nan_test(std::string variable){
+  std::vector<int> locations;
+  std::string message = ("For Ions " + variable + " ");
+  if (variable == "temperature_scgc") {
+    locations = insert_indefinites(temperature_scgc);
+    message += print_nan_vector(locations, temperature_scgc);
+  }
+  if (variable == "density_scgc") {
+    locations = insert_indefinites(density_scgc);
+    message += print_nan_vector(locations, density_scgc);
+  }
+  if (variable == "velocity_vcgc") {
+    locations = insert_indefinites(velocity_vcgc[0]);
+    message +=
+      "at the x loc " + print_nan_vector(locations, velocity_vcgc[0]);
+    locations = insert_indefinites(velocity_vcgc[1]);
+    message +=
+      "at the y loc " + print_nan_vector(locations, velocity_vcgc[1]);
+    locations = insert_indefinites(velocity_vcgc[2]);
+    message +=
+      "at the z loc " + print_nan_vector(locations, velocity_vcgc[2]);
+  }
+  std::cout << message;
+}
+
+//----------------------------------------------------------------------
+// Checks for nans and +/- infinities in density, temp, and velocity
+//----------------------------------------------------------------------
+
+bool Ions::check_for_nonfinites() {
+  bool non_finites_exist = false;
+  if (!all_finite(density_scgc, "density_scgc") ||
+      !all_finite(temperature_scgc, "temperature_scgc") ||
+      !all_finite(velocity_vcgc, "velocity_vcgc"))
+    non_finites_exist = true;
+  if (non_finites_exist)
+    throw std::string("Check for nonfinites failed!!!\n");
+  return non_finites_exist;
+}
+
 // -----------------------------------------------------------------------------
 // Calculate the electron density from the sum of all ion species
 // -----------------------------------------------------------------------------