LCOV - code coverage report
Current view: top level - src/reporters - PolynomialChaosReporter.C (source / functions) Hit Total Coverage
Test: idaholab/moose stochastic_tools: f45d79 Lines: 108 115 93.9 %
Date: 2025-07-25 05:00:46 Functions: 10 11 90.9 %
Legend: Lines: hit not hit

          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 "PolynomialChaosReporter.h"
      11             : 
      12             : #include "SobolReporter.h"
      13             : #include "StochasticReporter.h"
      14             : 
      15             : registerMooseObject("StochasticToolsApp", PolynomialChaosReporter);
      16             : 
      17             : InputParameters
      18         544 : PolynomialChaosReporter::validParams()
      19             : {
      20         544 :   InputParameters params = GeneralReporter::validParams();
      21         544 :   params.addClassDescription(
      22             :       "Tool for extracting data from PolynomialChaos surrogates and computing statistics.");
      23        1088 :   params.addRequiredParam<std::vector<UserObjectName>>(
      24             :       "pc_name", "Name(s) of PolynomialChaos surrogate object(s).");
      25        1088 :   params.addParam<bool>("include_data",
      26        1088 :                         false,
      27             :                         "True to output information on the polynomial chaos model, including "
      28             :                         "polynomial types, orders, and coefficients.");
      29         544 :   MultiMooseEnum stats("mean=1 stddev=2 skewness=3 kurtosis=4", "");
      30        1088 :   params.addParam<MultiMooseEnum>("statistics", stats, "Statistics to compute.");
      31        1088 :   params.addParam<bool>("include_sobol", false, "True to compute Sobol indices.");
      32        1088 :   params.addParam<std::vector<std::vector<Real>>>(
      33             :       "local_sensitivity_points",
      34        1088 :       std::vector<std::vector<Real>>(0),
      35             :       "Points for each polynomial chaos surrogate specifying desired location of sensitivity "
      36             :       "measurement.");
      37        1088 :   params.addParam<std::vector<SamplerName>>(
      38             :       "local_sensitivity_sampler",
      39        1088 :       std::vector<SamplerName>(0),
      40             :       "Sampler for each polynomial chaos surrogate specifying desired location of sensitivity "
      41             :       "measurement.");
      42         544 :   return params;
      43         544 : }
      44             : 
      45         272 : PolynomialChaosReporter::PolynomialChaosReporter(const InputParameters & parameters)
      46             :   : GeneralReporter(parameters),
      47             :     SurrogateModelInterface(this),
      48         544 :     _loc_point(getParam<std::vector<std::vector<Real>>>("local_sensitivity_points"))
      49             : {
      50         576 :   for (const auto & sn : getParam<std::vector<SamplerName>>("local_sensitivity_sampler"))
      51          32 :     _loc_sampler.push_back(&getSamplerByName(sn));
      52             : 
      53         880 :   for (const auto & nm : getParam<std::vector<UserObjectName>>("pc_name"))
      54             :   {
      55             :     // Gather polynomial chaos surrogates
      56         336 :     _pc.push_back(&getSurrogateModelByName<PolynomialChaos>(nm));
      57             :     // PolynomialChaos data reporter
      58         672 :     if (getParam<bool>("include_data"))
      59             :     {
      60         384 :       auto & pc_ptr = declareValueByName<const PolynomialChaos *>(nm, REPORTER_MODE_ROOT);
      61         128 :       pc_ptr = _pc.back();
      62             :     }
      63             :     // Statistics reporter
      64         992 :     for (const auto & stat : getParam<MultiMooseEnum>("statistics"))
      65         960 :       declareValueByName<std::pair<Real, std::vector<Real>>, PCStatisticsContext<Real>>(
      66         640 :           nm + "_" + stat.name(), REPORTER_MODE_ROOT, *_pc.back(), stat);
      67             :     // Sobol reporter
      68         672 :     if (getParam<bool>("include_sobol"))
      69             :       declareValueByName<std::pair<std::vector<Real>, std::vector<std::vector<Real>>>,
      70         448 :                          PCSobolContext<Real>>(nm + "_SOBOL", REPORTER_MODE_ROOT, *_pc.back());
      71             :     // Local sensitivity points reporter
      72         336 :     if (_loc_point.size() > 0)
      73             :     {
      74          96 :       if (_loc_point.size() < _pc.size())
      75           0 :         paramError("local_sensitivity_points",
      76             :                    "There must be a set of points for each inputted polynomial chaos model.");
      77         288 :       _loc_point_sense.push_back(&declareValueByName<std::vector<std::vector<Real>>>(
      78         192 :           nm + "_POINT_SENSITIVITY", REPORTER_MODE_ROOT));
      79             :     }
      80             :     // Local sensitivity sampler reporter
      81         336 :     if (_loc_sampler.size() > 0)
      82             :     {
      83          32 :       if (_loc_sampler.size() < _pc.size())
      84           0 :         paramError("local_sensitivity_sampler",
      85             :                    "There must be a sampler for each inputted polynomial chaos model.");
      86             :       _loc_sampler_sense.push_back(
      87          32 :           &declareValueByName<std::vector<std::vector<Real>>,
      88          96 :                               StochasticReporterContext<std::vector<Real>>>(
      89          64 :               nm + "_SAMPLE_SENSITIVITY", REPORTER_MODE_ROOT, *_loc_sampler[_pc.size() - 1]));
      90             :     }
      91             :   }
      92         272 : }
      93             : 
      94             : void
      95         272 : PolynomialChaosReporter::execute()
      96             : {
      97         272 :   if ((_loc_point.size() + _loc_sampler.size()) == 0)
      98             :     return;
      99             : 
     100         160 :   for (const auto & i : index_range(_pc))
     101             :   {
     102          96 :     const auto & pc = *_pc[i];
     103             :     const auto nparam = pc.getNumberOfParameters();
     104          96 :     if (_loc_point.size() > 0)
     105             :     {
     106          96 :       if (_loc_point[i].size() % pc.getNumberOfParameters() != 0)
     107           0 :         paramError("local_sensitivity_points",
     108             :                    "Number of values must be divisible by number of parameters in "
     109             :                    "Polynomial Chaos model.");
     110             : 
     111          96 :       const auto npoint = _loc_point[i].size() / nparam;
     112          96 :       _loc_point_sense[i]->resize(npoint);
     113         288 :       for (const auto & p : make_range(npoint))
     114             :       {
     115         192 :         const std::vector<Real> data(_loc_point[i].begin() + p * nparam,
     116         192 :                                      _loc_point[i].begin() + (p + 1) * nparam);
     117         384 :         (*_loc_point_sense[i])[p] = computeLocalSensitivity(pc, data);
     118             :       }
     119             : 
     120          96 :       if (_loc_sampler.size() > 0)
     121             :       {
     122          32 :         if (_loc_sampler[i]->getNumberOfCols() != nparam)
     123           0 :           paramError(
     124             :               "local_sensitivity_sampler",
     125             :               "Sampler ",
     126           0 :               _loc_sampler[i]->name(),
     127             :               " does not have the same number of columns as the number of dimensions in model ",
     128           0 :               pc.name(),
     129             :               ".");
     130        2032 :         for (const auto & p : make_range(_loc_sampler[i]->getNumberOfLocalRows()))
     131        2000 :           (*_loc_sampler_sense[i])[p] =
     132        4000 :               computeLocalSensitivity(pc, _loc_sampler[i]->getNextLocalRow());
     133             :       }
     134             :     }
     135             :   }
     136             : }
     137             : 
     138             : std::vector<Real>
     139        2192 : PolynomialChaosReporter::computeLocalSensitivity(const PolynomialChaos & pc,
     140             :                                                  const std::vector<Real> & data)
     141             : {
     142        2192 :   std::vector<Real> sense(data.size());
     143        2192 :   const auto val = pc.evaluate(data);
     144        6704 :   for (const auto & d : index_range(data))
     145        4512 :     sense[d] = data[d] / val * pc.computeDerivative(d, data);
     146        2192 :   return sense;
     147             : }
     148             : 
     149             : void
     150          80 : to_json(nlohmann::json & json, const PolynomialChaos * const & pc)
     151             : {
     152          80 :   pc->store(json);
     153          80 : }
     154             : 
     155             : template <typename OutType>
     156         320 : PCStatisticsContext<OutType>::PCStatisticsContext(
     157             :     const libMesh::ParallelObject & other,
     158             :     const MooseObject & producer,
     159             :     ReporterState<std::pair<OutType, std::vector<OutType>>> & state,
     160             :     const PolynomialChaos & pc,
     161             :     const MooseEnumItem & stat)
     162             :   : ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>(other, producer, state),
     163         320 :     _pc(pc),
     164         320 :     _stat(stat)
     165             : {
     166         320 : }
     167             : 
     168             : template <typename OutType>
     169             : void
     170         320 : PCStatisticsContext<OutType>::finalize()
     171             : {
     172             :   // Compute standard deviation
     173             :   OutType sig = OutType();
     174         320 :   if (_stat == "stddev" || _stat == "skewness" || _stat == "kurtosis")
     175         192 :     sig = _pc.computeStandardDeviation();
     176             : 
     177         320 :   OutType & val = this->_state.value().first;
     178         320 :   val = OutType();
     179             :   // Mean
     180         320 :   if (_stat == "mean" && this->processor_id() == 0)
     181          80 :     val = _pc.computeMean();
     182             :   // Standard Deviation
     183         240 :   else if (_stat == "stddev" && this->processor_id() == 0)
     184          80 :     val = sig;
     185             :   // Skewness
     186         160 :   else if (_stat == "skewness")
     187          32 :     val = _pc.powerExpectation(3) / (sig * sig * sig);
     188             :   // Kurtosis
     189         128 :   else if (_stat == "kurtosis")
     190          32 :     val = _pc.powerExpectation(4) / (sig * sig * sig * sig);
     191         320 :   this->_communicator.sum(val);
     192             : 
     193         320 :   ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>::finalize();
     194         320 : }
     195             : 
     196             : template <typename OutType>
     197             : void
     198         200 : PCStatisticsContext<OutType>::storeInfo(nlohmann::json & json) const
     199             : {
     200         200 :   ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>::storeInfo(json);
     201         200 :   json["stat"] = _stat.name();
     202         200 : }
     203             : 
     204             : template <typename OutType>
     205         112 : PCSobolContext<OutType>::PCSobolContext(
     206             :     const libMesh::ParallelObject & other,
     207             :     const MooseObject & producer,
     208             :     ReporterState<std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>>> & state,
     209             :     const PolynomialChaos & pc)
     210             :   : ReporterGeneralContext<std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>>>(
     211             :         other, producer, state),
     212         112 :     _pc(pc)
     213             : {
     214           0 : }
     215             : 
     216             : template <typename OutType>
     217             : void
     218         112 : PCSobolContext<OutType>::finalize()
     219             : {
     220         112 :   const unsigned int nparam = _pc.getNumberOfParameters();
     221         112 :   std::vector<OutType> & val = this->_state.value().first;
     222             :   val.clear();
     223             : 
     224             :   // Compute variance
     225         112 :   auto var = _pc.computeStandardDeviation();
     226         112 :   var *= var;
     227             : 
     228             :   // First order
     229         656 :   for (const auto & i : make_range(nparam))
     230        1088 :     val.push_back(_pc.computeSobolIndex({i}) / var);
     231             :   // Total
     232         656 :   for (const auto & i : make_range(nparam))
     233         544 :     val.push_back(_pc.computeSobolTotal(i) / var);
     234             :   // Second order
     235         656 :   for (const auto & i : make_range(nparam))
     236        1648 :     for (const auto & j : make_range(i + 1, nparam))
     237        2208 :       val.push_back(_pc.computeSobolIndex({i, j}) / var);
     238         112 : }
     239             : 
     240             : template <typename OutType>
     241             : void
     242          70 : PCSobolContext<OutType>::store(nlohmann::json & json) const
     243             : {
     244          70 :   SobolReporterContext<std::vector<OutType>, OutType>::storeSobol(
     245          70 :       json, this->_state.value(), _pc.getNumberOfParameters());
     246          70 : }
     247             : 
     248             : template class PCStatisticsContext<Real>;
     249             : template class PCSobolContext<Real>;

Generated by: LCOV version 1.14