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 : }
|