Use MPI groups to parallelize data assimilation

source/DataAssimilator.hh

@@ -66,7 +66,9 @@ public:
   /**
    * Constructor.
    */
-  DataAssimilator(boost::property_tree::ptree const &database);
+  DataAssimilator(MPI_Comm const &global_communicator,
+                  MPI_Comm const &local_communicator, int my_color,
+                  boost::property_tree::ptree const &database);
 
   /**
    * This is the main public interface for the class and is called to perform
@@ -83,8 +85,7 @@ public:
    * u is a perturbation to the observed solution
    * H is the observation matrix
    */
-  void update_ensemble(MPI_Comm const &communicator,
-                       std::vector<dealii::LA::distributed::BlockVector<double>>
+  void update_ensemble(std::vector<dealii::LA::distributed::BlockVector<double>>
                            &augmented_state_ensemble,
                        std::vector<double> const &expt_data,
                        dealii::SparseMatrix<double> const &R);
@@ -159,25 +160,34 @@ private:
       std::vector<dealii::LA::distributed::BlockVector<double>> const
           &vec_ensemble) const;
 
+  // TODO
+  MPI_Comm _global_communicator;
+
+  MPI_Comm _local_communicator;
+
+  int _global_rank = -1;
+
+  int _color = -1;
+
   /**
    * The number of ensemble members in the simulation.
    */
-  unsigned int _num_ensemble_members;
+  unsigned int _num_ensemble_members = 0;
 
   /**
    * The length of the data vector for each simulation ensemble member.
    */
-  unsigned int _sim_size;
+  unsigned int _sim_size = 0;
 
   /**
    * The length of the parameter vector for each simulation ensemble member.
    */
-  unsigned int _parameter_size;
+  unsigned int _parameter_size = 0;
 
   /**
    * The length of the data vector the experimental observations.
    */
-  unsigned int _expt_size;
+  unsigned int _expt_size = 0;
 
   /**
    * The sparsity pattern for the localized covariance matrix.

source/ensemble_management.cc

@@ -1,35 +1,33 @@
-/* Copyright (c) 2021, the adamantine authors.
+/* Copyright (c) 2021-2023, the adamantine authors.
  *
  * This file is subject to the Modified BSD License and may not be distributed
  * without copyright and license information. Please refer to the file LICENSE
  * for the text and further information on this license.
  */
 
 #include <ensemble_management.hh>
+#include <utils.hh>
 
 namespace adamantine
 {
-std::vector<double> fill_and_sync_random_vector(unsigned int length,
-                                                double mean, double stddev)
+std::vector<double> get_normal_random_vector(unsigned int length,
+                                             unsigned int n_rejected_draws,
+                                             double mean, double stddev)
 {
-  std::vector<double> output_vector(length);
-
-  unsigned int rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  ASSERT(stddev >= 0., "Internal Error");
 
-  if (rank == 0)
+  std::mt19937 pseudorandom_number_generator;
+  std::normal_distribution<> normal_dist_generator(mean, stddev);
+  for (unsigned int i = 0; i < n_rejected_draws; ++i)
   {
-    std::mt19937 pseudorandom_number_generator;
-    std::normal_distribution<> normal_dist_generator(0.0, 1.0);
-
-    for (unsigned int member = 0; member < length; ++member)
-    {
-      output_vector[member] =
-          mean + stddev * normal_dist_generator(pseudorandom_number_generator);
-    }
+    normal_dist_generator(pseudorandom_number_generator);
   }
 
-  output_vector =
-      dealii::Utilities::MPI::broadcast(MPI_COMM_WORLD, output_vector, 0);
+  std::vector<double> output_vector(length);
+  for (unsigned int i = 0; i < length; ++i)
+  {
+    output_vector[i] = normal_dist_generator(pseudorandom_number_generator);
+  }
 
   return output_vector;
 }

source/ensemble_management.hh

@@ -1,4 +1,4 @@
-/* Copyright (c) 2021, the adamantine authors.
+/* Copyright (c) 2021-2023, the adamantine authors.
  *
  * This file is subject to the Modified BSD License and may not be distributed
  * without copyright and license information. Please refer to the file LICENSE
@@ -14,8 +14,14 @@
 #include <vector>
 namespace adamantine
 {
-std::vector<double> fill_and_sync_random_vector(unsigned int length,
-                                                double mean, double stddev);
+/**
+ * Return a vector of size @p length, with random values drawn following a
+ * normal distribution of average @p mean and standard deviation @stddev. The
+ * first @p n_rejected_draws are rejected.
+ */
+std::vector<double> get_normal_random_vector(unsigned int length,
+                                             unsigned int n_rejected_draws,
+                                             double mean, double stddev);
 } // namespace adamantine
 
 #endif

tests/CMakeLists.txt

@@ -22,8 +22,6 @@ list(APPEND
      test_timer
      test_utils
      test_integration_3d_amr
-     test_integration_da
-     test_integration_da_augmented
      test_validate_input_database
      )
 
@@ -33,6 +31,8 @@ list(APPEND
      test_experimental_data
      test_integration_2d
      test_integration_3d
+     test_integration_da
+     test_integration_da_augmented
      test_material_deposition
      test_thermal_physics
      test_ensemble_management

tests/test_data_assimilator.cc

@@ -27,7 +27,7 @@ class DataAssimilatorTester
     boost::property_tree::ptree database;
 
     // First checking the dealii default values
-    DataAssimilator da0(database);
+    DataAssimilator da0(MPI_COMM_WORLD, MPI_COMM_WORLD, 0, database);
 
     double tol = 1.0e-12;
     BOOST_TEST(da0._solver_control.tolerance() - 1.0e-10 == 0.,
@@ -39,7 +39,7 @@ class DataAssimilatorTester
     database.put("solver.convergence_tolerance", 1.0e-6);
     database.put("solver.max_iterations", 25);
     database.put("solver.max_number_of_temp_vectors", 4);
-    DataAssimilator da1(database);
+    DataAssimilator da1(MPI_COMM_WORLD, MPI_COMM_WORLD, 0, database);
     BOOST_TEST(da1._solver_control.tolerance() - 1.0e-6 == 0.,
                tt::tolerance(tol));
     BOOST_TEST(da1._solver_control.max_steps() == 25u);
@@ -81,7 +81,7 @@ class DataAssimilatorTester
     expt_to_dof_mapping.second[1] = 3;
 
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._sim_size = sim_size;
     da._expt_size = expt_size;
     da._parameter_size = 0;
@@ -183,7 +183,8 @@ class DataAssimilatorTester
     expt_to_dof_mapping.second[2] = 3;
 
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(MPI_COMM_WORLD, MPI_COMM_WORLD, 0,
+                       solver_settings_database);
     da._sim_size = sim_size;
     da._expt_size = expt_size;
     da.update_dof_mapping<2>(expt_to_dof_mapping);
@@ -228,7 +229,7 @@ class DataAssimilatorTester
     expt_to_dof_mapping.second[2] = 3;
 
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._sim_size = sim_size;
     da._expt_size = expt_size;
     da.update_dof_mapping<2>(expt_to_dof_mapping);
@@ -296,7 +297,7 @@ class DataAssimilatorTester
     expt_to_dof_mapping.second[2] = 3;
 
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._sim_size = sim_size;
     da._expt_size = expt_size;
     da.update_dof_mapping<2>(expt_to_dof_mapping);
@@ -328,7 +329,7 @@ class DataAssimilatorTester
 
     // Effectively a dense matrix
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._localization_cutoff_distance = 100.0;
     da._localization_cutoff_function = LocalizationCutoff::step_function;
     da.update_covariance_sparsity_pattern<2>(dof_handler, 0);
@@ -389,7 +390,7 @@ class DataAssimilatorTester
     solver_settings_database.put("localization_cutoff_distance", 100.0);
     solver_settings_database.put("localization_cutoff_function",
                                  "step_function");
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._sim_size = sim_vec.size();
     da._num_ensemble_members = vec_ensemble.size();
 
@@ -563,7 +564,8 @@ class DataAssimilatorTester
     if (R_is_diagonal)
     {
       boost::property_tree::ptree solver_settings_database;
-      DataAssimilator da(solver_settings_database);
+      DataAssimilator da(MPI_COMM_WORLD, MPI_COMM_WORLD, 0,
+                         solver_settings_database);
 
       dealii::SparsityPattern pattern(3, 3, 1);
       pattern.add(0, 0);
@@ -596,7 +598,8 @@ class DataAssimilatorTester
     else
     {
       boost::property_tree::ptree solver_settings_database;
-      DataAssimilator da(solver_settings_database);
+      DataAssimilator da(MPI_COMM_WORLD, MPI_COMM_WORLD, 0,
+                         solver_settings_database);
 
       dealii::SparsityPattern pattern(3, 3, 3);
       pattern.add(0, 0);
@@ -669,7 +672,7 @@ class DataAssimilatorTester
     expt_to_dof_mapping.second[1] = 3;
 
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._sim_size = sim_size;
     da._parameter_size = 0;
     da._expt_size = expt_size;
@@ -726,7 +729,7 @@ class DataAssimilatorTester
     sim_at_expt_pt_2_before.push_back(augmented_state_ensemble[2].block(0)[3]);
 
     // Update the simulation data
-    da.update_ensemble(communicator, augmented_state_ensemble, expt_vec, R);
+    da.update_ensemble(augmented_state_ensemble, expt_vec, R);
 
     // Save the data at the observation points after assimilation
     std::vector<double> sim_at_expt_pt_1_after(3);
@@ -787,7 +790,7 @@ class DataAssimilatorTester
     expt_to_dof_mapping.second[1] = 3;
 
     boost::property_tree::ptree solver_settings_database;
-    DataAssimilator da(solver_settings_database);
+    DataAssimilator da(communicator, communicator, 0, solver_settings_database);
     da._sim_size = sim_size;
     da._parameter_size = parameter_size;
     da._expt_size = expt_size;
@@ -855,7 +858,7 @@ class DataAssimilatorTester
     sim_at_expt_pt_2_before.push_back(augmented_state_ensemble[2].block(0)[3]);
 
     // Update the simulation data
-    da.update_ensemble(communicator, augmented_state_ensemble, expt_vec, R);
+    da.update_ensemble(augmented_state_ensemble, expt_vec, R);
 
     // Save the data at the observation points after assimilation
     std::vector<double> sim_at_expt_pt_1_after(3);

tests/test_ensemble_management.cc

@@ -27,7 +27,7 @@ BOOST_AUTO_TEST_CASE(fill_and_sync_random_vector, *utf::tolerance(10.))
   double stddev = 0.25;
   unsigned int ensemble_size = 5000;
   std::vector<double> vec =
-      adamantine::fill_and_sync_random_vector(ensemble_size, mean, stddev);
+      adamantine::get_normal_random_vector(ensemble_size, 0, mean, stddev);
 
   // Check vector size
   BOOST_TEST(vec.size() == ensemble_size);

tests/test_integration_2d.cc

@@ -56,18 +56,18 @@ BOOST_AUTO_TEST_CASE(integration_2D_ensemble, *utf::tolerance(0.1))
   boost::property_tree::ptree database;
   boost::property_tree::info_parser::read_info(filename, database);
 
-  auto result = run_ensemble<2, dealii::MemorySpace::Host>(communicator,
-                                                           database, timers);
+  auto result_ensemble = run_ensemble<2, dealii::MemorySpace::Host>(
+      communicator, database, timers);
 
-  for (unsigned int member = 0; member < 3; ++member)
+  for (auto &result_member : result_ensemble)
   {
     std::ifstream gold_file("integration_2d_gold.txt");
-    for (unsigned int i = 0; i < result[member].block(0).locally_owned_size();
+    for (unsigned int i = 0; i < result_member.block(0).locally_owned_size();
          ++i)
     {
       double gold_value = -1.;
       gold_file >> gold_value;
-      BOOST_TEST(result[member].block(0).local_element(i) == gold_value);
+      BOOST_TEST(result_member.block(0).local_element(i) == gold_value);
     }
   }
 }