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             : #pragma once
      10             : 
      11             : #include "Calculators.h"
      12             : 
      13             : #include "libmesh/dense_matrix.h"
      14             : 
      15             : namespace StochasticTools
      16             : {
      17             : /**
      18             :  * Calculator for computing Sobol sensitivity indices according to the paper by Saltelli (2002)
      19             :  * https://doi.org/10.1016/S0010-4655(02)00280-1
      20             :  *
      21             :  * The data provided is stacked vectors provided by the SobolSampler. Example use of this object
      22             :  * is also available in the stochastic_tools unit testing.
      23             :  */
      24             : template <typename InType, typename OutType>
      25             : class SobolCalculator : public Calculator<std::vector<InType>, std::vector<OutType>>
      26             : {
      27             : public:
      28             :   SobolCalculator(const libMesh::ParallelObject & other, const std::string & name, bool resample);
      29             : 
      30             : protected:
      31             :   virtual void initialize() override;
      32             :   virtual void update(const InType & data) override;
      33             :   virtual void finalize(bool is_distributed) override;
      34      192146 :   virtual std::vector<OutType> get() const override { return _sobol; }
      35             : 
      36             : private:
      37             :   /// Set to true if the resampling matrix exists for computing second-order indices
      38             :   const bool _resample;
      39             :   /// Matrix containing dot products of data
      40             :   DenseMatrix<OutType> _amat;
      41             :   /// The returned sobol indices
      42             :   std::vector<OutType> _sobol;
      43             : };
      44             : 
      45             : template <typename InType, typename OutType>
      46             : class SobolCalculator<std::vector<InType>, std::vector<OutType>>
      47             :   : public Calculator<std::vector<std::vector<InType>>, std::vector<std::vector<OutType>>>
      48             : {
      49             : public:
      50             :   SobolCalculator(const libMesh::ParallelObject & other, const std::string & name, bool resample);
      51             : 
      52             : protected:
      53             :   virtual void initialize() override;
      54             :   virtual void update(const std::vector<InType> & data) override;
      55             :   virtual void finalize(bool is_distributed) override;
      56        7007 :   virtual std::vector<std::vector<OutType>> get() const override { return _values; }
      57             : 
      58             : private:
      59             :   /// Set to true if the resampling matrix exists for computing second-order indices
      60             :   const bool _resample;
      61             :   /// Sobol calculator for each entry in vector data
      62             :   std::vector<SobolCalculator<InType, OutType>> _sobol_calcs;
      63             :   /// Sobol indices
      64             :   std::vector<std::vector<OutType>> _values;
      65             : };
      66             : 
      67             : template <typename InType, typename OutType>
      68           7 : SobolCalculator<std::vector<InType>, std::vector<OutType>>::SobolCalculator(
      69             :     const libMesh::ParallelObject & other, const std::string & name, bool resample)
      70             :   : Calculator<std::vector<std::vector<InType>>, std::vector<std::vector<OutType>>>(other, name),
      71           7 :     _resample(resample)
      72             : {
      73           7 : }
      74             : 
      75             : template <typename InType, typename OutType>
      76             : void
      77        7007 : SobolCalculator<std::vector<InType>, std::vector<OutType>>::initialize()
      78             : {
      79        7007 :   _sobol_calcs.clear();
      80        7007 :   _values.clear();
      81        7007 : }
      82             : 
      83             : template <typename InType, typename OutType>
      84             : void
      85       50050 : SobolCalculator<std::vector<InType>, std::vector<OutType>>::update(const std::vector<InType> & data)
      86             : {
      87             :   // Need to transpose data
      88             :   std::vector<InType> data_update;
      89      750750 :   for (const auto & d : index_range(data))
      90     2102100 :     for (const auto & i : index_range(data[d]))
      91             :     {
      92     1401400 :       if (data_update.size() <= i)
      93      100100 :         data_update.emplace_back(data.size());
      94     1401400 :       data_update[i][d] = data[d][i];
      95             :     }
      96             : 
      97             :   // Build calculators and update
      98      150150 :   for (const auto & i : index_range(data_update))
      99             :   {
     100      100100 :     if (_sobol_calcs.size() <= i)
     101             :     {
     102       28028 :       _sobol_calcs.emplace_back(*this, "SOBOL_" + std::to_string(i), _resample);
     103       14014 :       _sobol_calcs[i].initializeCalculator();
     104             :     }
     105      100100 :     _sobol_calcs[i].updateCalculator(data_update[i]);
     106             :   }
     107       50050 : }
     108             : 
     109             : template <typename InType, typename OutType>
     110             : void
     111        7007 : SobolCalculator<std::vector<InType>, std::vector<OutType>>::finalize(bool is_distributed)
     112             : {
     113             :   // Need to create calculators here no data was added
     114        7007 :   if (is_distributed)
     115             :   {
     116        7007 :     auto ncalc = _sobol_calcs.size();
     117        7007 :     this->_communicator.max(ncalc);
     118        7007 :     for (const auto & i : make_range(_sobol_calcs.size(), ncalc))
     119             :     {
     120           0 :       _sobol_calcs.emplace_back(*this, "SOBOL_" + std::to_string(i), _resample);
     121           0 :       _sobol_calcs[i].initializeCalculator();
     122             :     }
     123             :   }
     124             : 
     125       35035 :   for (const auto & i : index_range(_sobol_calcs))
     126             :   {
     127       14014 :     _sobol_calcs[i].finalizeCalculator(is_distributed);
     128       14014 :     const auto val = _sobol_calcs[i].getValue();
     129             : 
     130      392392 :     for (const auto & ind : index_range(val))
     131             :     {
     132      378378 :       if (_values.size() <= ind)
     133      189189 :         _values.emplace_back(_sobol_calcs.size());
     134      378378 :       _values[ind][i] = val[ind];
     135             :     }
     136             :   }
     137        7007 : }
     138             : 
     139             : } // namespace