Add capability to store Monte Carlo sample data in parallel

ref #13879

framework/include/samplers/Sampler.h

@@ -137,6 +137,21 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
    */
   dof_id_type getLocalRowEnd();
 
+  /**
+   * Return true if the sample data is distributed among all processors.
+   */
+  bool useDistributedSampleData() { return _use_distributed_sample_data; };
+
+  /**
+   * Return the total number of samples.
+   */
+  virtual dof_id_type getTotalNumberOfSamples() = 0;
+
+  /**
+   * Return the total number of matrices.
+   */
+  virtual dof_id_type getTotalNumberOfMatrices() = 0;
+
 protected:
   /**
    * Get the next random number from the generator.
@@ -164,9 +179,8 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
 
   /**
    * Reinitialize the offsets and row counts.
-   * @param data Samples, as returned from getSamples(), to use for re-computing size information
    */
-  void reinit(const std::vector<DenseMatrix<Real>> & data);
+  void reinit();
 
   /// Map used to store the perturbed parameters and their corresponding distributions
   std::vector<Distribution const *> _distributions;
@@ -177,6 +191,9 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
   /// Sample names
   std::vector<std::string> _sample_names;
 
+  /// Sampler data is distributed among all processors
+  bool _use_distributed_sample_data;
+
 private:
   /// Random number generator, don't give users access we want to control it via the interface
   /// from this class.

framework/src/samplers/Sampler.C

@@ -30,6 +30,8 @@ validParams<Sampler>()
   params.addRequiredParam<std::vector<DistributionName>>(
       "distributions", "The names of distributions that you want to sample.");
   params.addParam<unsigned int>("seed", 0, "Random number generator initial seed");
+  params.addParam<bool>(
+      "use_distributed_sample_data", false, "Sample data is distributed across all processors.");
   params.registerBase("Sampler");
   return params;
 }
@@ -39,6 +41,7 @@ Sampler::Sampler(const InputParameters & parameters)
     SetupInterface(this),
     DistributionInterface(this),
     _distribution_names(getParam<std::vector<DistributionName>>("distributions")),
+    _use_distributed_sample_data(getParam<bool>("use_distributed_sample_data")),
     _seed(getParam<unsigned int>("seed")),
     _total_rows(0)
 {
@@ -51,23 +54,23 @@ void
 Sampler::execute()
 {
   // Get the samples then save the state so that subsequent calls to getSamples returns the same
-  // random numbers until this execute command is called again.
-  std::vector<DenseMatrix<Real>> data = getSamples();
   _generator.saveState();
-  reinit(data);
+  reinit();
 }
 
 void
-Sampler::reinit(const std::vector<DenseMatrix<Real>> & data)
+Sampler::reinit()
 {
   // Update offsets and total number of rows
   _total_rows = 0;
   _offsets.clear();
-  _offsets.reserve(data.size() + 1);
+  unsigned int num_matrices = getTotalNumberOfMatrices();
+  unsigned int num_samples = getTotalNumberOfSamples();
+  _offsets.reserve(num_matrices + 1);
   _offsets.push_back(_total_rows);
-  for (const DenseMatrix<Real> & mat : data)
+  for (unsigned int i = 0; i < num_matrices; i++)
   {
-    _total_rows += mat.m();
+    _total_rows += num_samples;
     _offsets.push_back(_total_rows);
   }
 
@@ -137,7 +140,7 @@ Sampler::Location
 Sampler::getLocation(dof_id_type global_index)
 {
   if (_offsets.empty())
-    reinit(getSamples());
+    reinit();
 
   mooseAssert(_offsets.size() > 1,
               "The getSamples method returned an empty vector, if you are seeing this you have "
@@ -157,7 +160,7 @@ dof_id_type
 Sampler::getTotalNumberOfRows()
 {
   if (_total_rows == 0)
-    reinit(getSamples());
+    reinit();
   return _total_rows;
 }
 
@@ -168,7 +171,7 @@ dof_id_type
 Sampler::getLocalNumerOfRows()
 {
   if (_total_rows == 0)
-    reinit(getSamples());
+    reinit();
   return _local_rows;
 }
 
@@ -179,7 +182,7 @@ dof_id_type
 Sampler::getLocalRowBegin()
 {
   if (_total_rows == 0)
-    reinit(getSamples());
+    reinit();
   return _local_row_begin;
 }
 
@@ -190,6 +193,6 @@ dof_id_type
 Sampler::getLocalRowEnd()
 {
   if (_total_rows == 0)
-    reinit(getSamples());
+    reinit();
   return _local_row_end;
 }

modules/stochastic_tools/include/samplers/MonteCarloSampler.h

@@ -26,6 +26,10 @@ class MonteCarloSampler : public Sampler
 protected:
   virtual std::vector<DenseMatrix<Real>> sample() override;
 
+  virtual dof_id_type getTotalNumberOfSamples() override { return _num_samples; };
+
+  virtual dof_id_type getTotalNumberOfMatrices() override { return _num_matrices; };
+
   /// Number of matrices
   const dof_id_type _num_matrices;
 

modules/stochastic_tools/include/samplers/SobolSampler.h

@@ -27,6 +27,8 @@ class SobolSampler : public Sampler
   virtual std::vector<DenseMatrix<Real>> sample() override;
   virtual void sampleSetUp() override;
   virtual void sampleTearDown() override;
+  virtual dof_id_type getTotalNumberOfSamples() override { return _num_samples; };
+  virtual dof_id_type getTotalNumberOfMatrices() override { return _distributions.size() + 2; };
 
   /// Number of Monte Carlo samples to create for each Sobol matrix
   const std::size_t _num_samples;
@@ -37,4 +39,3 @@ class SobolSampler : public Sampler
   DenseMatrix<Real> _b_matrix;
   ///@}
 };
-

modules/stochastic_tools/include/vectorpostprocessors/SamplerData.h

@@ -27,6 +27,7 @@ class SamplerData : public GeneralVectorPostprocessor, SamplerInterface
   SamplerData(const InputParameters & parameters);
   void virtual initialize() override;
   void virtual execute() override;
+  void virtual finalize() override;
 
 protected:
   /// Storage for declared vectors
@@ -38,4 +39,3 @@ class SamplerData : public GeneralVectorPostprocessor, SamplerInterface
   /// Whether to output the number of rows and columns in the first two rows of output
   const bool & _output_col_row_sizes;
 };
-

modules/stochastic_tools/src/multiapps/SamplerFullSolveMultiApp.C

@@ -160,11 +160,30 @@ SamplerFullSolveMultiApp::getCommandLineArgsParamHelper(unsigned int local_app)
         oss << ";";
 
       if (_mode == StochasticTools::MultiAppMode::BATCH_RESET)
-        oss << cli_args_name[col] << "="
-            << Moose::stringify(matrix(_local_batch_app_index + _sampler.getLocalRowBegin(), col));
+      {
+        if (_sampler.useDistributedSampleData())
+        {
+          oss << cli_args_name[col] << "=" << Moose::stringify(matrix(_local_batch_app_index, col));
+          if (_local_batch_app_index >= matrix.m())
+            oss.str("");
+        }
+        else
+        {
+          oss << cli_args_name[col] << "="
+              << Moose::stringify(
+                     matrix(_local_batch_app_index + _sampler.getLocalRowBegin(), col));
+          if ((_local_batch_app_index + _sampler.getLocalRowBegin()) >= matrix.m())
+            oss.str("");
+        }
+      }
       else
+      {
         oss << cli_args_name[col] << "="
             << Moose::stringify(matrix(local_app + _first_local_app, col));
+
+        if ((local_app + _first_local_app) >= matrix.m())
+          oss.str("");
+      }
     }
   }
 

modules/stochastic_tools/src/samplers/MonteCarloSampler.C

@@ -30,18 +30,51 @@ MonteCarloSampler::MonteCarloSampler(const InputParameters & parameters)
     _num_matrices(getParam<dof_id_type>("n_matrices")),
     _num_samples(getParam<dof_id_type>("n_samples"))
 {
+  if (_use_distributed_sample_data && n_processors() > _num_samples)
+    mooseError("In MonteCarloSampler, the number of MPI processes cannot be larger than the number "
+               "of samples(",
+               _num_samples,
+               ") when sample data is distributed across all processor.");
 }
 
 std::vector<DenseMatrix<Real>>
 MonteCarloSampler::sample()
 {
   std::vector<DenseMatrix<Real>> output(_num_matrices);
-  for (MooseIndex(_num_matrices) m = 0; m < _num_matrices; ++m)
+
+  if (_use_distributed_sample_data)
   {
-    output[m].resize(_num_samples, _distributions.size());
-    for (MooseIndex(_num_samples) i = 0; i < _num_samples; ++i)
-      for (MooseIndex(_distributions) j = 0; j < _distributions.size(); ++j)
-        output[m](i, j) = _distributions[j]->quantile(rand());
+    dof_id_type local_rows, local_row_begin, local_row_end;
+    MooseUtils::linearPartitionItems(
+        _num_samples, n_processors(), processor_id(), local_rows, local_row_begin, local_row_end);
+
+    for (MooseIndex(_num_matrices) m = 0; m < _num_matrices; ++m)
+    {
+      output[m].resize(local_rows, _distributions.size());
+      dof_id_type local_i = 0;
+      for (MooseIndex(_num_samples) i = 0; i < _num_samples; ++i)
+      {
+        for (MooseIndex(_distributions) j = 0; j < _distributions.size(); ++j)
+        {
+          Real quant = _distributions[j]->quantile(rand());
+          if (i >= local_row_begin && i < local_row_end)
+            output[m](local_i, j) = quant;
+        }
+        if (i >= local_row_begin && i < local_row_end)
+          local_i++;
+      }
+    }
   }
+  else
+  {
+    for (MooseIndex(_num_matrices) m = 0; m < _num_matrices; ++m)
+    {
+      output[m].resize(_num_samples, _distributions.size());
+      for (MooseIndex(_num_samples) i = 0; i < _num_samples; ++i)
+        for (MooseIndex(_distributions) j = 0; j < _distributions.size(); ++j)
+          output[m](i, j) = _distributions[j]->quantile(rand());
+    }
+  }
+
   return output;
 }

modules/stochastic_tools/src/transfers/SamplerPostprocessorTransfer.C

@@ -97,7 +97,8 @@ SamplerPostprocessorTransfer::finalizeFromMultiapp()
   if (processor_id() == 0)
   {
     _results->initialize();
-    const dof_id_type n = _sampler->getTotalNumberOfRows();
+    const dof_id_type n =
+        (_sampler->getTotalNumberOfMatrices()) * (_sampler->getTotalNumberOfSamples());
     for (MooseIndex(n) i = 0; i < n; i++)
     {
       Sampler::Location loc = _sampler->getLocation(i);

modules/stochastic_tools/src/transfers/SamplerTransfer.C

@@ -101,7 +101,6 @@ SamplerTransfer::initializeToMultiapp()
 void
 SamplerTransfer::executeToMultiapp()
 {
-
   SamplerReceiver * ptr = getReceiver(processor_id());
   std::vector<Real> row = getRow(_global_index);
   ptr->transfer(_parameter_names, row);

modules/stochastic_tools/src/vectorpostprocessors/SamplerData.C

@@ -78,3 +78,16 @@ SamplerData::execute()
               _sample_vectors[i]->begin() + offset);
   }
 }
+
+void
+SamplerData::finalize()
+{
+  if (_sampler.useDistributedSampleData())
+  {
+    std::vector<DenseMatrix<Real>> data = _sampler.getSamples();
+
+    auto n = data.size();
+    for (MooseIndex(data) i = 0; i < n; ++i)
+      _communicator.gather(0, *_sample_vectors[i]);
+  }
+}

modules/stochastic_tools/src/vectorpostprocessors/StochasticResults.C

@@ -38,9 +38,11 @@ StochasticResults::initialize()
 
   // Resize and zero vectors to the correct size, this allows the SamplerPostprocessorTransfer
   // to set values in the vector directly.
-  std::vector<DenseMatrix<Real>> data = _sampler->getSamples();
-  for (MooseIndex(data) i = 0; i < data.size(); ++i)
-    _sample_vectors[i]->resize(data[i].m(), 0);
+  unsigned int num_matrices = _sampler->getTotalNumberOfMatrices();
+  unsigned int num_samples = _sampler->getTotalNumberOfSamples();
+
+  for (unsigned i = 0; i < num_matrices; ++i)
+    _sample_vectors[i]->resize(num_samples, 0);
 }
 
 VectorPostprocessorValue &

modules/stochastic_tools/test/tests/multiapps/batch_commandline_control/gold/master_multiple_out_sampler_data_0002.csv

@@ -0,0 +1,7 @@
+mat_0
+8.970799818498
+5.7991231115664
+6.9539363296392
+9.6824257791673
+8.2825361623213
+7.7957496369701

modules/stochastic_tools/test/tests/multiapps/batch_commandline_control/gold/master_multiple_out_storage_0002.csv

@@ -1,4 +1,4 @@
 sample_0
-8.47877
-6.82795
-5.96624
+8.9708
+6.95394
+8.28254

modules/stochastic_tools/test/tests/multiapps/batch_commandline_control/gold/master_single_out_sampler_data_0002.csv

@@ -0,0 +1,4 @@
+mat_0
+8.970799818498
+5.7991231115664
+6.9539363296392

modules/stochastic_tools/test/tests/multiapps/batch_commandline_control/gold/master_single_out_storage_0002.csv

@@ -1,4 +1,4 @@
 sample_0
-9.68243
-8.28254
-7.79575
+8.9708
+5.79912
+6.95394

modules/stochastic_tools/test/tests/multiapps/batch_commandline_control/master_multiple.i

@@ -50,6 +50,11 @@
   [storage]
     type = StochasticResults
   []
+  [./sampler_data]
+    type = SamplerData
+    sampler = sample
+    execute_on = 'INITIAL'
+  [../]
 []
 
 [Controls]

modules/stochastic_tools/test/tests/multiapps/batch_commandline_control/master_single.i

@@ -51,6 +51,11 @@
   [storage]
     type = StochasticResults
   []
+  [./sampler_data]
+    type = SamplerData
+    sampler = sample
+    execute_on = 'INITIAL'
+  [../]
 []
 
 [Controls]