Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "SobolStatistics.h"
      11             : #include "SobolSampler.h"
      12             : #include "SobolCalculators.h"
      13             : #include "BootstrapCalculators.h"
      14             : 
      15             : // MOOSE includes
      16             : #include "MooseVariable.h"
      17             : #include "ThreadedElementLoopBase.h"
      18             : #include "ThreadedNodeLoop.h"
      19             : 
      20             : #include "libmesh/quadrature.h"
      21             : 
      22             : #include <numeric>
      23             : 
      24             : registerMooseObjectDeprecated("StochasticToolsApp", SobolStatistics, "11/03/2021 12:00");
      25             : 
      26             : InputParameters
      27          64 : SobolStatistics::validParams()
      28             : {
      29          64 :   InputParameters params = GeneralVectorPostprocessor::validParams();
      30          64 :   params.addClassDescription(
      31             :       "Compute SOBOL statistics values of a given VectorPostprocessor objects and vectors.");
      32         128 :   params.addParam<SamplerName>("sampler", "SobolSampler object.");
      33         128 :   params.addParam<VectorPostprocessorName>(
      34             :       "results", "StochasticResults object containing data to use for calculation.");
      35             : 
      36         128 :   params.addParam<std::vector<Real>>(
      37             :       "ci_levels",
      38          64 :       std::vector<Real>(),
      39             :       "A vector of confidence levels to consider, values must be in (0, 1).");
      40         128 :   params.addParam<unsigned int>(
      41             :       "ci_replicates",
      42         128 :       10000,
      43             :       "The number of replicates to use when computing confidence level intervals.");
      44         128 :   params.addParam<unsigned int>("ci_seed",
      45         128 :                                 1,
      46             :                                 "The random number generator seed used for creating replicates "
      47             :                                 "while computing confidence level intervals.");
      48          64 :   return params;
      49           0 : }
      50             : 
      51          32 : SobolStatistics::SobolStatistics(const InputParameters & parameters)
      52             :   : GeneralVectorPostprocessor(parameters),
      53          32 :     _sobol_sampler(getSampler<SobolSampler>("sampler")),
      54          64 :     _ci_levels(getParam<std::vector<Real>>("ci_levels")),
      55          64 :     _ci_replicates(getParam<unsigned int>("ci_replicates")),
      56          96 :     _ci_seed(getParam<unsigned int>("ci_seed"))
      57             : {
      58          32 : }
      59             : 
      60             : void
      61          32 : SobolStatistics::initialSetup()
      62             : {
      63          32 :   const VectorPostprocessorName & vpp_name = getParam<VectorPostprocessorName>("results");
      64          32 :   const VectorPostprocessor & vpp_object = _fe_problem.getVectorPostprocessorObjectByName(vpp_name);
      65          32 :   const std::set<std::string> & vpp_vectors = vpp_object.getVectorNames();
      66          32 :   _sobol_ci_vectors.resize(_ci_levels.size());
      67          64 :   for (const auto & vec_name : vpp_vectors)
      68             :   {
      69             :     ReporterName rname(vpp_name, vec_name);
      70          32 :     _result_vectors.push_back(std::make_pair(
      71          32 :         &getReporterValueByName<VectorPostprocessorValue>(rname), vpp_object.isDistributed()));
      72          32 :     _sobol_stat_vectors.push_back(&declareVector(vpp_name + "_" + vec_name));
      73             : 
      74         176 :     for (const auto & l : index_range(_ci_levels))
      75             :     {
      76         144 :       std::stringstream vname; /// Vectors computed by this object
      77         288 :       vname << vpp_name << "_" << vec_name << "_" << _ci_levels[l] * 100 << "CI";
      78         144 :       _sobol_ci_vectors[l].push_back(&declareVector(vname.str()));
      79         144 :     }
      80             :   }
      81          32 : }
      82             : 
      83             : void
      84          32 : SobolStatistics::execute()
      85             : {
      86          64 :   TIME_SECTION("execute", 3, "Executing Sobol Statistics");
      87             : 
      88             :   StochasticTools::SobolCalculator<std::vector<Real>, Real> calc(
      89          32 :       *this, "SOBOL", _sobol_sampler.resample());
      90          32 :   auto boot_calc = _ci_levels.empty()
      91          32 :                        ? nullptr
      92          64 :                        : makeBootstrapCalculator(MooseEnum("percentile", "percentile"),
      93             :                                                  *this,
      94             :                                                  _ci_levels,
      95          16 :                                                  _ci_replicates,
      96          16 :                                                  _ci_seed,
      97          16 :                                                  calc);
      98          64 :   for (std::size_t i = 0; i < _result_vectors.size(); ++i)
      99             :   {
     100          32 :     const std::size_t ncol = _sobol_sampler.resample() ? 2 * _sobol_sampler.getNumberOfCols() + 2
     101           0 :                                                        : _sobol_sampler.getNumberOfCols() + 2;
     102             :     const std::vector<std::vector<Real>> data =
     103          32 :         StochasticTools::reshapeVector(*(_result_vectors[i].first), ncol, /*row_major =*/true);
     104          64 :     (*_sobol_stat_vectors[i]) = calc.compute(data, _result_vectors[i].second);
     105             : 
     106          32 :     if (boot_calc)
     107             :     {
     108          16 :       std::vector<std::vector<Real>> sobol_ci = boot_calc->compute(data, _result_vectors[i].second);
     109         106 :       for (const auto & l : index_range(sobol_ci))
     110          90 :         (*_sobol_ci_vectors[l][i]) = std::move(sobol_ci[l]);
     111          16 :     }
     112          32 :   }
     113          64 : }