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 "PolynomialChaos.h"
      11             : #include "Sampler.h"
      12             : #include "CartesianProduct.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", PolynomialChaos);
      15             : 
      16             : InputParameters
      17         918 : PolynomialChaos::validParams()
      18             : {
      19         918 :   InputParameters params = SurrogateModel::validParams();
      20         918 :   params.addClassDescription("Computes and evaluates polynomial chaos surrogate model.");
      21         918 :   return params;
      22           0 : }
      23             : 
      24         459 : PolynomialChaos::PolynomialChaos(const InputParameters & parameters)
      25             :   : SurrogateModel(parameters),
      26         459 :     _order(getModelData<unsigned int>("_order")),
      27         918 :     _ndim(getModelData<unsigned int>("_ndim")),
      28         918 :     _ncoeff(getModelData<std::size_t>("_ncoeff")),
      29         918 :     _tuple(getModelData<std::vector<std::vector<unsigned int>>>("_tuple")),
      30         918 :     _coeff(getModelData<std::vector<Real>>("_coeff")),
      31         459 :     _poly(
      32         918 :         getModelData<std::vector<std::unique_ptr<const PolynomialQuadrature::Polynomial>>>("_poly"))
      33             : {
      34         459 : }
      35             : 
      36             : Real
      37        6364 : PolynomialChaos::evaluate(const std::vector<Real> & x) const
      38             : {
      39             :   mooseAssert(x.size() == _ndim, "Number of inputted parameters does not match PC model.");
      40             : 
      41        6364 :   DenseMatrix<Real> poly_val(_ndim, _order);
      42             : 
      43             :   // Evaluate polynomials to avoid duplication
      44       20328 :   for (unsigned int d = 0; d < _ndim; ++d)
      45       75464 :     for (unsigned int i = 0; i < _order; ++i)
      46       61500 :       poly_val(d, i) = _poly[d]->compute(i, x[d], /*normalize =*/false);
      47             : 
      48             :   Real val = 0;
      49      140184 :   for (std::size_t i = 0; i < _ncoeff; ++i)
      50             :   {
      51      133820 :     Real tmp = _coeff[i];
      52      499800 :     for (unsigned int d = 0; d < _ndim; ++d)
      53      365980 :       tmp *= poly_val(d, _tuple[i][d]);
      54      133820 :     val += tmp;
      55             :   }
      56             : 
      57        6364 :   return val;
      58        6364 : }
      59             : 
      60             : const std::vector<std::vector<unsigned int>> &
      61          88 : PolynomialChaos::getPolynomialOrders() const
      62             : {
      63          88 :   return _tuple;
      64             : }
      65             : 
      66             : unsigned
      67           0 : PolynomialChaos::getPolynomialOrder(const unsigned int dim, const unsigned int i) const
      68             : {
      69           0 :   return _tuple[i][dim];
      70             : }
      71             : 
      72             : const std::vector<Real> &
      73          88 : PolynomialChaos::getCoefficients() const
      74             : {
      75          88 :   return _coeff;
      76             : }
      77             : 
      78             : Real
      79          88 : PolynomialChaos::computeMean() const
      80             : {
      81             :   mooseAssert(_coeff.size() > 0, "The coefficient matrix is empty.");
      82          88 :   return _coeff[0];
      83             : }
      84             : 
      85             : Real
      86         323 : PolynomialChaos::computeStandardDeviation() const
      87             : {
      88             :   Real var = 0;
      89      103088 :   for (std::size_t i = 1; i < _ncoeff; ++i)
      90             :   {
      91             :     Real norm = 1.0;
      92      519435 :     for (std::size_t d = 0; d < _ndim; ++d)
      93      416670 :       norm *= _poly[d]->innerProduct(_tuple[i][d]);
      94      102765 :     var += _coeff[i] * _coeff[i] * norm;
      95             :   }
      96             : 
      97         323 :   return std::sqrt(var);
      98             : }
      99             : 
     100             : Real
     101          68 : PolynomialChaos::powerExpectation(const unsigned int n) const
     102             : {
     103             :   std::vector<StochasticTools::WeightedCartesianProduct<unsigned int, Real>> order;
     104          68 :   order.reserve(_ndim);
     105          68 :   std::vector<std::vector<Real>> c_1d(n, std::vector<Real>(_coeff.begin() + 1, _coeff.end()));
     106         204 :   for (unsigned int d = 0; d < _ndim; ++d)
     107             :   {
     108         136 :     std::vector<std::vector<unsigned int>> order_1d(n, std::vector<unsigned int>(_ncoeff - 1));
     109        2040 :     for (std::size_t i = 1; i < _ncoeff; ++i)
     110        8568 :       for (unsigned int m = 0; m < n; ++m)
     111        6664 :         order_1d[m][i - 1] = _tuple[i][d];
     112         136 :     order.push_back(StochasticTools::WeightedCartesianProduct<unsigned int, Real>(order_1d, c_1d));
     113         136 :   }
     114             : 
     115             :   dof_id_type n_local, st_local, end_local;
     116          68 :   MooseUtils::linearPartitionItems(
     117             :       order[0].numRows(), n_processors(), processor_id(), n_local, st_local, end_local);
     118             : 
     119             :   Real val = 0;
     120      905588 :   for (dof_id_type i = st_local; i < end_local; ++i)
     121             :   {
     122      905520 :     Real tmp = order[0].computeWeight(i);
     123     1336995 :     for (unsigned int d = 0; d < _ndim; ++d)
     124             :     {
     125     1246971 :       if (MooseUtils::absoluteFuzzyEqual(tmp, 0.0))
     126             :         break;
     127      431475 :       std::vector<unsigned int> comb(n);
     128     2112671 :       for (unsigned int m = 0; m < n; ++m)
     129     1681196 :         comb[m] = order[d].computeValue(i, m);
     130      431475 :       tmp *= _poly[d]->productIntegral(comb);
     131      431475 :     }
     132      905520 :     val += tmp;
     133             :   }
     134             : 
     135          68 :   return val;
     136          68 : }
     137             : 
     138             : Real
     139        4944 : PolynomialChaos::computeDerivative(const unsigned int dim, const std::vector<Real> & x) const
     140             : {
     141        4944 :   return computePartialDerivative({dim}, x);
     142             : }
     143             : 
     144             : Real
     145        4944 : PolynomialChaos::computePartialDerivative(const std::vector<unsigned int> & dim,
     146             :                                           const std::vector<Real> & x) const
     147             : {
     148             :   mooseAssert(x.size() == _ndim, "Number of inputted parameters does not match PC model.");
     149             : 
     150        4944 :   std::vector<unsigned int> grad(_ndim);
     151        9888 :   for (const auto & d : dim)
     152             :   {
     153             :     mooseAssert(d < _ndim, "Specified dimension is greater than total number of parameters.");
     154        4944 :     grad[d]++;
     155             :   }
     156             : 
     157        4944 :   DenseMatrix<Real> poly_val(_ndim, _order);
     158             : 
     159             :   // Evaluate polynomials to avoid duplication
     160       15376 :   for (unsigned int d = 0; d < _ndim; ++d)
     161       68032 :     for (unsigned int i = 0; i < _order; ++i)
     162       57600 :       poly_val(d, i) = _poly[d]->computeDerivative(i, x[d], grad[d]);
     163             : 
     164             :   Real val = 0;
     165      269504 :   for (std::size_t i = 0; i < _ncoeff; ++i)
     166             :   {
     167      264560 :     Real tmp = _coeff[i];
     168     1182640 :     for (unsigned int d = 0; d < _ndim; ++d)
     169      918080 :       tmp *= poly_val(d, _tuple[i][d]);
     170      264560 :     val += tmp;
     171             :   }
     172             : 
     173        4944 :   return val;
     174        4944 : }
     175             : 
     176             : Real
     177        1751 : PolynomialChaos::computeSobolIndex(const std::set<unsigned int> & ind) const
     178             : {
     179        1751 :   linearPartitionCoefficients();
     180             : 
     181             :   // If set is empty, compute mean
     182        1751 :   if (ind.empty())
     183           0 :     return computeMean();
     184             : 
     185             :   // Do some sanity checks in debug
     186             :   mooseAssert(ind.size() <= _ndim, "Number of indices is greater than number of parameters.");
     187             :   mooseAssert(*ind.rbegin() < _ndim, "Maximum index provided exceeds number of parameters.");
     188             : 
     189             :   Real val = 0.0;
     190      374563 :   for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
     191             :   {
     192      372812 :     Real tmp = _coeff[i] * _coeff[i];
     193      696520 :     for (unsigned int d = 0; d < _ndim; ++d)
     194             :     {
     195      683353 :       if ((ind.find(d) != ind.end() && _tuple[i][d] > 0) ||
     196      441584 :           (ind.find(d) == ind.end() && _tuple[i][d] == 0))
     197             :       {
     198      323708 :         tmp *= _poly[d]->innerProduct(_tuple[i][d]);
     199             :       }
     200             :       else
     201             :       {
     202             :         tmp = 0.0;
     203             :         break;
     204             :       }
     205             :     }
     206      372812 :     val += tmp;
     207             :   }
     208             : 
     209             :   return val;
     210             : }
     211             : 
     212             : Real
     213         578 : PolynomialChaos::computeSobolTotal(const unsigned int dim) const
     214             : {
     215         578 :   linearPartitionCoefficients();
     216             : 
     217             :   // Do some sanity checks in debug
     218             :   mooseAssert(dim < _ndim, "Requested dimension is greater than number of parameters.");
     219             : 
     220             :   Real val = 0.0;
     221      143050 :   for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
     222      142472 :     if (_tuple[i][dim] > 0)
     223       92664 :       val += _coeff[i] * _coeff[i] * _poly[dim]->innerProduct(_tuple[i][dim]);
     224             : 
     225         578 :   return val;
     226             : }
     227             : 
     228             : void
     229        2329 : PolynomialChaos::linearPartitionCoefficients() const
     230             : {
     231        2329 :   if (_n_local_coeff == std::numeric_limits<dof_id_type>::max())
     232         119 :     MooseUtils::linearPartitionItems(_ncoeff,
     233             :                                      n_processors(),
     234             :                                      processor_id(),
     235         119 :                                      _n_local_coeff,
     236         119 :                                      _local_coeff_begin,
     237         119 :                                      _local_coeff_end);
     238        2329 : }
     239             : 
     240             : void
     241          88 : PolynomialChaos::store(nlohmann::json & json) const
     242             : {
     243         176 :   json["order"] = _order;
     244         176 :   json["ndim"] = getNumberOfParameters();
     245          88 :   json["ncoeff"] = getNumberofCoefficients();
     246          88 :   json["tuple"] = getPolynomialOrders();
     247          88 :   json["coeff"] = getCoefficients();
     248         352 :   for (const auto & p : _poly)
     249             :   {
     250             :     nlohmann::json jsonp;
     251         264 :     p->store(jsonp);
     252         264 :     json["poly"].push_back(jsonp);
     253             :   }
     254          88 : }