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             : #include "SobolCalculators.h"
      10             : 
      11             : namespace StochasticTools
      12             : {
      13             : template <typename InType, typename OutType>
      14       32166 : SobolCalculator<InType, OutType>::SobolCalculator(const libMesh::ParallelObject & other,
      15             :                                                   const std::string & name,
      16             :                                                   bool resample)
      17       32166 :   : Calculator<std::vector<InType>, std::vector<OutType>>(other, name), _resample(resample)
      18             : {
      19       32166 : }
      20             : 
      21             : template <typename InType, typename OutType>
      22             : void
      23      426326 : SobolCalculator<InType, OutType>::initialize()
      24             : {
      25      426326 :   _amat.resize(0, 0);
      26      426326 : }
      27             : 
      28             : template <typename InType, typename OutType>
      29             : void
      30    37253444 : SobolCalculator<InType, OutType>::update(const InType & data)
      31             : {
      32    37253444 :   if (_amat.n() == 0)
      33             :   {
      34      426182 :     _amat.resize(data.size(), data.size());
      35      426182 :     if (_amat.n() < 2 || (_resample && _amat.n() % 2 != 0))
      36           0 :       mooseError("Size of input data is inconsistent with Sobol resampling scheme.");
      37             :   }
      38             :   mooseAssert(_amat.n() == data.size(), "Size of data updating SobolCalculator has changed.");
      39             : 
      40   546794260 :   for (const auto & c : index_range(data))
      41  4319066536 :     for (const auto & r : make_range(c + 1))
      42  3809525720 :       _amat(r, c) += data[r] * data[c];
      43    37253444 : }
      44             : 
      45             : template <typename InType, typename OutType>
      46             : void
      47      426326 : SobolCalculator<InType, OutType>::finalize(bool is_distributed)
      48             : {
      49             :   // The size of the A matrix
      50      426326 :   std::size_t s = _amat.n();
      51             : 
      52      426326 :   if (is_distributed)
      53             :   {
      54      416320 :     this->_communicator.max(s);
      55      416320 :     if (_amat.n() == 0)
      56         144 :       _amat.resize(s, s);
      57             :     mooseAssert(_amat.n() == s, "Size of data updating SobolCalculator has changed.");
      58      416320 :     this->_communicator.sum(_amat.get_values());
      59             :   }
      60             : 
      61             :   // Number of independent parameters
      62      426326 :   const std::size_t n = (s - 2) / (_resample ? 2 : 1);
      63             : 
      64             :   // The data to output
      65      426326 :   DenseMatrix<OutType> S(n, n);
      66      426326 :   std::vector<OutType> ST(n);
      67             : 
      68             :   // First order
      69             :   {
      70      426326 :     auto E = _amat(0, s - 1);
      71      426326 :     auto V = _amat(s - 1, s - 1) - E;
      72     2344218 :     for (std::size_t i = 0; i < n; ++i)
      73             :     {
      74     1917892 :       auto U0 = _amat(i + 1, s - 1);
      75     1917892 :       if (_resample)
      76             :       {
      77     1916832 :         auto U1 = _amat(0, i + n + 1);
      78     1916832 :         S(i, i) = (((U0 + U1) / 2) - E) / V;
      79             :       }
      80             :       else
      81        1060 :         S(i, i) = (U0 - E) / V;
      82             :     }
      83             :   }
      84             : 
      85             :   // Total-effect
      86             :   {
      87      426326 :     auto E = _amat(0, s - 1);
      88      426326 :     auto V = _amat(0, 0) - E;
      89     2344218 :     for (std::size_t i = 0; i < n; ++i)
      90             :     {
      91     1917892 :       auto U0 = _amat(0, i + 1);
      92     1917892 :       if (_resample)
      93             :       {
      94     1916832 :         auto U1 = _amat(i + n + 1, s - 1);
      95     1916832 :         ST[i] = 1 - (((U0 + U1) / 2) - E) / V;
      96             :       }
      97             :       else
      98        1060 :         ST[i] = 1 - (U0 - E) / V;
      99             :     }
     100             :   }
     101             : 
     102             :   // Second-order
     103      426326 :   if (_resample)
     104             :   {
     105     2342981 :     for (std::size_t i = 0; i < n; ++i)
     106             :     {
     107     1916832 :       auto E = _amat(i + 1, i + n + 1);
     108     5428796 :       for (std::size_t j = 0; j < i; ++j)
     109             :       {
     110     3511964 :         auto V = _amat(j + n + 1, j + n + 1) - E;
     111     3511964 :         auto U0 = _amat(i + 1, j + n + 1);
     112     3511964 :         auto U1 = _amat(j + 1, i + n + 1);
     113     3511964 :         auto Sc = (((U0 + U1) / 2) - E) / V;
     114     3511964 :         auto S2 = Sc - S(i, i) - S(j, j);
     115     3511964 :         S(j, i) = S2;
     116             :       }
     117             :     }
     118             :   }
     119             : 
     120             :   // Output the data
     121      426326 :   _sobol.clear();
     122      426326 :   if (_resample)
     123      426149 :     _sobol.reserve(n * (n + 3) / 2);
     124             :   else
     125         177 :     _sobol.reserve(2 * n);
     126             : 
     127             :   // First-order
     128     2344218 :   for (std::size_t i = 0; i < n; ++i)
     129     1917892 :     _sobol.push_back(S(i, i));
     130             : 
     131             :   // Total-effect
     132     2344218 :   for (std::size_t i = 0; i < n; ++i)
     133     1917892 :     _sobol.push_back(ST[i]);
     134             : 
     135             :   // Second-order
     136      426326 :   if (_resample)
     137             :   {
     138     2342981 :     for (std::size_t i = 0; i < n; ++i)
     139     5428796 :       for (std::size_t j = i + 1; j < n; ++j)
     140     3511964 :         _sobol.push_back(S(i, j));
     141             :   }
     142      426326 : }
     143             : 
     144             : template class SobolCalculator<std::vector<Real>, Real>;
     145             : template class SobolCalculator<std::vector<std::vector<Real>>, std::vector<Real>>;
     146             : } // namespace