LCOV - code coverage report
Current view: top level - src/reporters - DirectPerturbationReporter.C (source / functions) Hit Total Coverage
Test: idaholab/moose stochastic_tools: f45d79 Lines: 77 86 89.5 %
Date: 2025-07-25 05:00:46 Functions: 11 11 100.0 %
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 "DirectPerturbationReporter.h"
      11             : #include "DirectPerturbationSampler.h"
      12             : 
      13             : registerMooseObject("StochasticToolsApp", DirectPerturbationReporter);
      14             : 
      15             : InputParameters
      16         160 : DirectPerturbationReporter::validParams()
      17             : {
      18         160 :   InputParameters params = GeneralReporter::validParams();
      19         160 :   params.addClassDescription("Compute local sensitivities using the direct perturbation method.");
      20             : 
      21         320 :   params.addRequiredParam<SamplerName>("sampler",
      22             :                                        "Direct Perturbation sampler used to generate samples.");
      23         320 :   params.addParam<std::vector<VectorPostprocessorName>>(
      24             :       "vectorpostprocessors",
      25             :       {},
      26             :       "List of VectorPostprocessor(s) to utilize for sensitivity computations.");
      27         320 :   params.addParam<std::vector<ReporterName>>(
      28             :       "reporters", {}, "List of Reporter values to utilize for sensitivity computations.");
      29             : 
      30         320 :   params.addParam<bool>("relative_sensitivity",
      31         320 :                         false,
      32             :                         "If the reporter should return relative or absolute sensitivities.");
      33             : 
      34         160 :   return params;
      35           0 : }
      36             : 
      37          80 : DirectPerturbationReporter::DirectPerturbationReporter(const InputParameters & parameters)
      38             :   : GeneralReporter(parameters),
      39          80 :     _sampler(getSampler<DirectPerturbationSampler>("sampler")),
      40         160 :     _relative_sensitivity(getParam<bool>("relative_sensitivity")),
      41          80 :     _initialized(false)
      42             : {
      43         240 :   if (getParam<std::vector<ReporterName>>("reporters").empty() &&
      44          80 :       getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors").empty())
      45           0 :     mooseError(
      46             :         "The 'vectorpostprocessors' and/or 'reporters' parameters must be defined and non-empty.");
      47          80 : }
      48             : 
      49             : void
      50          80 : DirectPerturbationReporter::initialize()
      51             : {
      52             :   // It is enough to do this once
      53          80 :   if (_initialized)
      54           0 :     return;
      55             : 
      56             :   // Sensitivities from Reporters
      57             :   std::vector<std::string> unsupported_types;
      58         480 :   for (const auto & r_name : getParam<std::vector<ReporterName>>("reporters"))
      59             :   {
      60         320 :     if (hasReporterValueByName<std::vector<Real>>(r_name))
      61         240 :       declareValueHelper<Real>(r_name);
      62          80 :     else if (hasReporterValueByName<std::vector<std::vector<Real>>>(r_name))
      63          80 :       declareValueHelper<std::vector<Real>>(r_name);
      64             :     else
      65           0 :       unsupported_types.push_back(r_name.getCombinedName());
      66             :   }
      67             : 
      68          80 :   if (!unsupported_types.empty())
      69           0 :     paramError("reporters",
      70             :                "The following reporter value(s) do not have a type supported by the "
      71             :                "DirectPerturbationReporter:\n",
      72           0 :                MooseUtils::join(unsupported_types, ", "));
      73             : 
      74             :   // Sensitivities from VPPs
      75          80 :   for (const auto & vpp_name :
      76         160 :        getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors"))
      77           0 :     for (const auto & vec_name :
      78           0 :          _fe_problem.getVectorPostprocessorObjectByName(vpp_name).getVectorNames())
      79           0 :       declareValueHelper<Real>(ReporterName(vpp_name, vec_name));
      80          80 :   _initialized = true;
      81          80 : }
      82             : 
      83             : template <typename DataType>
      84             : void
      85         320 : DirectPerturbationReporter::declareValueHelper(const ReporterName & r_name)
      86             : {
      87             :   // We create a reporter value for the sensitivities based on the input
      88         320 :   const auto & data = getReporterValueByName<std::vector<DataType>>(r_name);
      89         320 :   const std::string s_name = r_name.getObjectName() + "_" + r_name.getValueName();
      90         960 :   declareValueByName<std::vector<DataType>, DirectPerturbationReporterContext<DataType>>(
      91         320 :       s_name, REPORTER_MODE_ROOT, _sampler, data, _relative_sensitivity);
      92         320 : }
      93             : 
      94             : template <typename DataType>
      95         320 : DirectPerturbationReporterContext<DataType>::DirectPerturbationReporterContext(
      96             :     const libMesh::ParallelObject & other,
      97             :     const MooseObject & producer,
      98             :     ReporterState<std::vector<DataType>> & state,
      99             :     DirectPerturbationSampler & sampler,
     100             :     const std::vector<DataType> & data,
     101             :     const bool relative_sensitivity)
     102             :   : ReporterGeneralContext<std::vector<DataType>>(other, producer, state),
     103         320 :     _sampler(sampler),
     104         320 :     _data(data),
     105         320 :     _relative_sensitivity(relative_sensitivity)
     106             : {
     107         320 :   this->_state.value().resize(_sampler.getNumberOfCols());
     108         320 : }
     109             : 
     110             : template <typename DataType>
     111             : void
     112         320 : DirectPerturbationReporterContext<DataType>::finalize()
     113             : {
     114         320 :   const dof_id_type offset = _sampler.getLocalRowBegin();
     115         320 :   const dof_id_type num_columns = _sampler.getNumberOfCols();
     116             : 
     117             :   // If we are computing the relative sensitivity, we will need
     118             :   // the reference value. So the process that has that will communicate it
     119             :   // to everybody. We reuse the initialization function for the reference
     120             :   // value as well
     121         320 :   auto reference_value = initializeDataType(_data[0]);
     122         320 :   if (_relative_sensitivity)
     123             :   {
     124          64 :     if (_sampler.getLocalRowBegin() == 0)
     125          40 :       reference_value = _data[0];
     126             : 
     127          64 :     this->comm().sum(reference_value);
     128             :   }
     129             : 
     130        1280 :   for (const auto param_i : make_range(num_columns))
     131             :   {
     132             :     // If we need relative coefficients we need to normalize the difference
     133         960 :     const auto interval = _relative_sensitivity ? _sampler.getRelativeInterval(param_i)
     134         768 :                                                 : _sampler.getAbsoluteInterval(param_i);
     135             : 
     136             :     dof_id_type left_i;
     137             :     dof_id_type right_i;
     138             :     // Depending on which method we use, the indices will change
     139         960 :     if (_sampler.perturbationMethod() == "central_difference")
     140             :     {
     141         576 :       left_i = 2 * param_i + 1;
     142         576 :       right_i = 2 * param_i + 2;
     143             :     }
     144             :     else
     145             :     {
     146         384 :       left_i = param_i + 1;
     147             :       right_i = 0;
     148             :     }
     149             : 
     150             :     const bool left_i_in_owned_range =
     151         960 :         _sampler.getLocalRowBegin() <= left_i && left_i < _sampler.getLocalRowEnd();
     152             :     const bool right_i_in_owned_range =
     153         960 :         _sampler.getLocalRowBegin() <= right_i && right_i < _sampler.getLocalRowEnd();
     154             : 
     155         960 :     if (left_i_in_owned_range || right_i_in_owned_range)
     156             :     {
     157         192 :       const dof_id_type copy_i = left_i_in_owned_range ? left_i - offset : right_i - offset;
     158         768 :       this->_state.value()[param_i] = initializeDataType(_data[copy_i]);
     159             :     }
     160             : 
     161             :     // We add the contribution from one side
     162         768 :     if (left_i_in_owned_range)
     163         150 :       addSensitivityContribution(
     164         600 :           this->_state.value()[param_i], _data[left_i - offset], reference_value, interval);
     165             : 
     166             :     // We add the contribution from the other side
     167         960 :     if (right_i_in_owned_range)
     168         600 :       addSensitivityContribution(
     169         600 :           this->_state.value()[param_i], _data[right_i - offset], reference_value, -interval);
     170             :   }
     171             : 
     172             :   // We gather the contributions across all processors
     173         320 :   this->vectorSum();
     174         320 : }
     175             : 
     176             : template <typename DataType>
     177             : void
     178         300 : DirectPerturbationReporterContext<DataType>::addSensitivityContribution(
     179             :     DataType & add_to,
     180             :     const DataType & to_add,
     181             :     const DataType & reference_value,
     182             :     const Real interval) const
     183             : {
     184             :   // DataType is a numeric type that we can sum (excluding bool)
     185             :   if constexpr (std::is_arithmetic<DataType>::value && !std::is_same<DataType, bool>::value)
     186             :   {
     187         900 :     add_to += to_add / (_relative_sensitivity ? interval * reference_value : interval);
     188         900 :     return;
     189             :   }
     190             :   // DataType is a vector type
     191             :   else if constexpr (is_std_vector<DataType>::value)
     192             :   {
     193             :     // It can still be anything in the vector elements, so we will check this later
     194             :     using VectorValueType = typename DataType::value_type;
     195             : 
     196             :     mooseAssert(add_to.size() == to_add.size(), "The vectors for summation have different sizes!");
     197             : 
     198             :     // Check if the vector elements are of a numeric type
     199             :     if constexpr (std::is_arithmetic<VectorValueType>::value &&
     200             :                   !std::is_same<VectorValueType, bool>::value)
     201             :     {
     202        1200 :       for (const auto index : index_range(add_to))
     203         900 :         add_to[index] +=
     204         900 :             to_add[index] / (_relative_sensitivity ? interval * reference_value[index] : interval);
     205         300 :       return;
     206             :     }
     207             :     // Check if the vector elements are also vectors
     208             :     else if constexpr (is_std_vector<VectorValueType>::value)
     209             :     {
     210             :       // This is as deep as we will go for now
     211             :       using InnerValueType = typename VectorValueType::value_type;
     212             : 
     213             :       // Check if the inner vector elements are of a numeric type
     214             :       if constexpr (std::is_arithmetic<InnerValueType>::value &&
     215             :                     !std::is_same<InnerValueType, bool>::value)
     216             :       {
     217             :         // Iterate over each inner vector in the outer vector
     218             :         for (auto & inner_index : index_range(add_to))
     219             :         {
     220             :           mooseAssert(add_to[inner_index].size() == to_add[inner_index].size(),
     221             :                       "The vectors for summation have different sizes!");
     222             :           for (const auto index : index_range(add_to[inner_index]))
     223             :             add_to[inner_index][index] +=
     224             :                 to_add[inner_index][index] /
     225             :                 (_relative_sensitivity ? interval * reference_value[inner_index][index] : interval);
     226             :         }
     227             :         return;
     228             :       }
     229             :       else
     230             :         static_assert(Moose::always_false<DataType>,
     231             :                       "Sensitivity coefficient computation is not implemented for the given type!");
     232             :     }
     233             :   }
     234             :   else
     235             :     static_assert(Moose::always_false<DataType>,
     236             :                   "Sensitivity coefficient computation is not implemented for the given type!");
     237             : }
     238             : 
     239             : template <typename DataType>
     240             : DataType
     241         272 : DirectPerturbationReporterContext<DataType>::initializeDataType(
     242             :     const DataType & example_output) const
     243             : {
     244             :   // DataType is a numeric type so we just return 0
     245             :   if constexpr (std::is_arithmetic<DataType>::value && !std::is_same<DataType, bool>::value)
     246             :     return 0.0;
     247             :   // DataType is a vector type
     248             :   else if constexpr (is_std_vector<DataType>::value)
     249             :   {
     250             :     // It can still be anything in the vector elements, so we will check this later
     251             :     using VectorValueType = typename DataType::value_type;
     252             : 
     253             :     // Check if the vector elements are of a numeric type
     254             :     if constexpr (std::is_arithmetic<VectorValueType>::value &&
     255             :                   !std::is_same<VectorValueType, bool>::value)
     256         272 :       return std::vector<VectorValueType>(example_output.size(), 0);
     257             :     // Check if the vector elements are also vectors
     258             :     else if constexpr (is_std_vector<VectorValueType>::value)
     259             :     {
     260             :       // This is as deep as we will go for now
     261             :       using InnerValueType = typename VectorValueType::value_type;
     262             : 
     263             :       // Check if the inner vector elements are of a numeric type
     264             :       if constexpr (std::is_arithmetic<InnerValueType>::value &&
     265             :                     !std::is_same<InnerValueType, bool>::value)
     266             :       {
     267             :         auto return_vector = std::vector<VectorValueType>(example_output.size());
     268             :         // Iterate over each inner vector in the outer vector
     269             :         for (auto & inner_index : index_range(example_output))
     270             :           return_vector[inner_index].resize(example_output.size(), 0.0);
     271             :         return return_vector;
     272             :       }
     273             :     }
     274             :     else
     275             :       static_assert(Moose::always_false<DataType>,
     276             :                     "Sensitivity coefficient computation is not implemented for the given type!");
     277             :   }
     278             :   else
     279             :     static_assert(Moose::always_false<DataType>,
     280             :                   "Sensitivity coefficient initialization is not implemented for the given type!");
     281             : }

Generated by: LCOV version 1.14