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