Merge pull request #353 from Rombur/fix_random

Rework ensembles

application/adamantine.hh

@@ -1430,35 +1430,59 @@ run_ensemble(MPI_Comm const &global_communicator,
   const double initial_temperature_stddev =
       ensemble_database.get("initial_temperature_stddev", 0.0);
 
+  unsigned int n_rejected_draws = first_local_member;
   std::vector<double> initial_temperature =
       adamantine::get_normal_random_vector(
-          local_ensemble_size, first_local_member, initial_temperature_mean,
-          initial_temperature_stddev);
+          local_ensemble_size, n_rejected_draws, initial_temperature_mean,
+          initial_temperature_stddev, verbose_output);
 
   // PropertyTreeInput ensemble.new_material_temperature_stddev
   const double new_material_temperature_stddev =
       ensemble_database.get("new_material_temperature_stddev", 0.0);
 
+  // Update the number of rejected draws to make sure that the variables are not
+  // correlated.
+  n_rejected_draws += global_ensemble_size;
   std::vector<double> new_material_temperature =
       adamantine::get_normal_random_vector(
-          local_ensemble_size, first_local_member,
-          new_material_temperature_mean, new_material_temperature_stddev);
+          local_ensemble_size, n_rejected_draws, new_material_temperature_mean,
+          new_material_temperature_stddev, verbose_output);
 
   // PropertyTreeInput ensemble.beam_0_max_power_stddev
   const double beam_0_max_power_stddev =
       ensemble_database.get("beam_0_max_power_stddev", 0.0);
 
+  n_rejected_draws += global_ensemble_size;
   std::vector<double> beam_0_max_power = adamantine::get_normal_random_vector(
-      local_ensemble_size, first_local_member, beam_0_max_power_mean,
-      beam_0_max_power_stddev);
+      local_ensemble_size, n_rejected_draws, beam_0_max_power_mean,
+      beam_0_max_power_stddev, verbose_output);
 
   // PropertyTreeInput ensemble.beam_0_absorption_stddev
   const double beam_0_absorption_stddev =
       ensemble_database.get("beam_0_absorption_stddev", 0.0);
 
+  n_rejected_draws += global_ensemble_size;
   std::vector<double> beam_0_absorption = adamantine::get_normal_random_vector(
-      local_ensemble_size, first_local_member, beam_0_absorption_mean,
-      beam_0_absorption_stddev);
+      local_ensemble_size, n_rejected_draws, beam_0_absorption_mean,
+      beam_0_absorption_stddev, verbose_output);
+
+  // Write the ensemble variables in files to simplify data analysis.
+  if (dealii::Utilities::MPI::this_mpi_process(local_communicator) == 0)
+  {
+    for (unsigned int member = 0; member < local_ensemble_size; ++member)
+    {
+
+      std::string member_data_filename =
+          post_processor_database.get<std::string>("filename_prefix") + "_m" +
+          std::to_string(first_local_member + member) + "_data.txt";
+      std::ofstream file(member_data_filename);
+      file << "Initial temperature: " << initial_temperature[member] << "\n";
+      file << "New material temperature: " << new_material_temperature[member]
+           << "\n";
+      file << "Beam_0 max power: " << beam_0_max_power[member] << "\n";
+      file << "Beam_0 absorption: " << beam_0_absorption[member] << "\n";
+    }
+  }
 
   // Create a new property tree database for each ensemble member
   std::vector<boost::property_tree::ptree> database_ensemble(

source/ensemble_management.cc

@@ -9,7 +9,8 @@ namespace adamantine
 {
 std::vector<double> get_normal_random_vector(unsigned int length,
                                              unsigned int n_rejected_draws,
-                                             double mean, double stddev)
+                                             double mean, double stddev,
+                                             bool verbose)
 {
   ASSERT(stddev >= 0., "Internal Error");
 
@@ -23,7 +24,20 @@ std::vector<double> get_normal_random_vector(unsigned int length,
   std::vector<double> output_vector(length);
   for (unsigned int i = 0; i < length; ++i)
   {
-    output_vector[i] = normal_dist_generator(pseudorandom_number_generator);
+    // We reject negative values because physical quantities we care about are
+    // all positive and we cannot guarantee that the normal distribution will
+    // always be positive.
+    do
+    {
+      output_vector[i] = normal_dist_generator(pseudorandom_number_generator);
+
+      if (verbose && output_vector[i] < 0.)
+      {
+        std::cout << "Random value rejected because it was negative: "
+                  << output_vector[i] << std::endl;
+      }
+
+    } while (output_vector[i] < 0.);
   }
 
   return output_vector;

source/ensemble_management.hh

@@ -17,11 +17,14 @@ namespace adamantine
  * first @p n_rejected_draws are rejected. @p n_rejected_draws is used to
  * ensure that wether we use one MPI rank or ten, we use the same random
  * numbers. If we don't reject the first few draws, all the processors will use
- * the same "random" numbers since they have the same seed.
+ * the same "random" numbers since they have the same seed. We also reject
+ * negative values even because our physical quantities must be positive. If @p
+ * verbose is true, we output a message when a negative value was rejected.
  */
 std::vector<double> get_normal_random_vector(unsigned int length,
                                              unsigned int n_rejected_draws,
-                                             double mean, double stddev);
+                                             double mean, double stddev,
+                                             bool verbose);
 } // namespace adamantine
 
 #endif

tests/test_ensemble_management.cc

@@ -20,11 +20,11 @@ BOOST_AUTO_TEST_CASE(fill_and_sync_random_vector, *utf::tolerance(10.))
 {
 
   // Create the random vector
-  double mean = -1.2;
+  double mean = 12.2;
   double stddev = 0.25;
   unsigned int ensemble_size = 5000;
-  std::vector<double> vec =
-      adamantine::get_normal_random_vector(ensemble_size, 0, mean, stddev);
+  std::vector<double> vec = adamantine::get_normal_random_vector(
+      ensemble_size, 0, mean, stddev, false);
 
   // Check vector size
   BOOST_TEST(vec.size() == ensemble_size);