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