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 "SobolReporter.h"
11 : #include "SobolSampler.h"
12 : #include "SobolCalculators.h"
13 : #include "BootstrapCalculators.h"
14 :
15 : // MOOSE includes
16 : #include "MooseVariable.h"
17 : #include "ThreadedElementLoopBase.h"
18 : #include "ThreadedNodeLoop.h"
19 :
20 : #include "libmesh/quadrature.h"
21 :
22 : #include <numeric>
23 :
24 : registerMooseObject("StochasticToolsApp", SobolReporter);
25 :
26 : InputParameters
27 56 : SobolReporter::validParams()
28 : {
29 56 : InputParameters params = GeneralReporter::validParams();
30 56 : params.addClassDescription("Compute SOBOL statistics values of a given VectorPostprocessor or "
31 : "Reporter objects and vectors.");
32 112 : params.addParam<SamplerName>("sampler", "SobolSampler object.");
33 :
34 112 : params.addParam<std::vector<VectorPostprocessorName>>(
35 : "vectorpostprocessors",
36 : "List of VectorPostprocessor(s) to utilized for statistic computations.");
37 :
38 112 : params.addParam<std::vector<ReporterName>>(
39 : "reporters", {}, "List of Reporter values to utilized for statistic computations.");
40 :
41 56 : params.addParam<std::vector<Real>>(
42 : "ci_levels",
43 56 : std::vector<Real>(),
44 : "A vector of confidence levels to consider, values must be in (0, 1).");
45 112 : params.addParam<unsigned int>(
46 : "ci_replicates",
47 112 : 10000,
48 : "The number of replicates to use when computing confidence level intervals.");
49 112 : params.addParam<unsigned int>("ci_seed",
50 112 : 1,
51 : "The random number generator seed used for creating replicates "
52 : "while computing confidence level intervals.");
53 56 : return params;
54 0 : }
55 :
56 28 : SobolReporter::SobolReporter(const InputParameters & parameters)
57 : : GeneralReporter(parameters),
58 28 : _sobol_sampler(getSampler<SobolSampler>("sampler")),
59 56 : _ci_levels(getParam<std::vector<Real>>("ci_levels")),
60 56 : _ci_replicates(getParam<unsigned int>("ci_replicates")),
61 56 : _ci_seed(getParam<unsigned int>("ci_seed")),
62 28 : _initialized(false)
63 : {
64 : // CI levels error checking
65 28 : if (!_ci_levels.empty())
66 : {
67 28 : if (*std::min_element(_ci_levels.begin(), _ci_levels.end()) <= 0)
68 0 : paramError("ci_levels", "The supplied levels must be greater than zero.");
69 28 : else if (*std::max_element(_ci_levels.begin(), _ci_levels.end()) >= 1)
70 0 : paramError("ci_levels", "The supplied levels must be less than 1.0");
71 : }
72 :
73 112 : if ((!isParamValid("reporters") && !isParamValid("vectorpostprocessors")) ||
74 70 : (getParam<std::vector<ReporterName>>("reporters").empty() &&
75 56 : getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors").empty()))
76 0 : mooseError(
77 : "The 'vectorpostprocessors' and/or 'reporters' parameters must be defined and non-empty.");
78 28 : }
79 :
80 : void
81 28 : SobolReporter::initialize()
82 : {
83 28 : if (_initialized)
84 : return;
85 :
86 : // Stats for Reporters
87 56 : if (isParamValid("reporters"))
88 : {
89 : std::vector<std::string> unsupported_types;
90 56 : const auto & reporter_names = getParam<std::vector<ReporterName>>("reporters");
91 63 : for (const auto & r_name : reporter_names)
92 : {
93 35 : if (hasReporterValueByName<std::vector<Real>>(r_name))
94 28 : declareValueHelper<Real>(r_name);
95 7 : else if (hasReporterValueByName<std::vector<std::vector<Real>>>(r_name))
96 7 : declareValueHelper<std::vector<Real>>(r_name);
97 : else
98 0 : unsupported_types.emplace_back(r_name.getCombinedName());
99 : }
100 :
101 28 : if (!unsupported_types.empty())
102 0 : paramError("reporters",
103 : "The following reporter value(s) do not have a type supported by the "
104 : "StatisticsReporter:\n",
105 0 : MooseUtils::join(unsupported_types, ", "));
106 28 : }
107 :
108 : // Stats for VPP
109 56 : if (isParamValid("vectorpostprocessors"))
110 : {
111 28 : const auto & vpp_names = getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors");
112 28 : for (const auto & vpp_name : vpp_names)
113 : {
114 : const VectorPostprocessor & vpp_object =
115 14 : _fe_problem.getVectorPostprocessorObjectByName(vpp_name);
116 14 : const std::set<std::string> & vpp_vectors = vpp_object.getVectorNames();
117 28 : for (const auto & vec_name : vpp_vectors)
118 : {
119 14 : ReporterName r_name(vpp_name, vec_name);
120 14 : declareValueHelper<Real>(r_name);
121 : }
122 : }
123 : }
124 :
125 28 : _initialized = true;
126 : }
127 :
128 : void
129 20 : SobolReporter::store(nlohmann::json & json) const
130 : {
131 20 : Reporter::store(json);
132 20 : if (!_ci_levels.empty())
133 20 : json["confidence_intervals"] = {{"method", "percentile"},
134 : {"levels", _ci_levels},
135 : {"replicates", _ci_replicates},
136 260 : {"seed", _ci_seed}};
137 :
138 40 : json["num_params"] = _sobol_sampler.getNumberOfCols();
139 20 : json["indices"].push_back("FIRST_ORDER");
140 20 : json["indices"].push_back("TOTAL");
141 20 : if (_sobol_sampler.resample())
142 30 : json["indices"].push_back("SECOND_ORDER");
143 280 : }
144 :
145 : template <typename T>
146 : void
147 49 : SobolReporter::declareValueHelper(const ReporterName & r_name)
148 : {
149 49 : const auto & mode = _fe_problem.getReporterData().getReporterMode(r_name);
150 49 : const auto & data = getReporterValueByName<std::vector<T>>(r_name);
151 49 : const std::string s_name = r_name.getObjectName() + "_" + r_name.getValueName();
152 49 : if (!_ci_levels.empty())
153 147 : declareValueByName<SobolState<T>, SobolReporterContext<std::vector<T>, T>>(
154 : s_name,
155 : REPORTER_MODE_ROOT,
156 : data,
157 : mode,
158 : _sobol_sampler,
159 147 : MooseEnum("percentile", "percentile"),
160 : _ci_levels,
161 : _ci_replicates,
162 : _ci_seed);
163 : else
164 0 : declareValueByName<SobolState<T>, SobolReporterContext<std::vector<T>, T>>(
165 : s_name, REPORTER_MODE_ROOT, data, mode, _sobol_sampler);
166 49 : }
167 :
168 : template <typename InType, typename OutType>
169 49 : SobolReporterContext<InType, OutType>::SobolReporterContext(
170 : const libMesh::ParallelObject & other,
171 : const MooseObject & producer,
172 : ReporterState<SobolState<OutType>> & state,
173 : const InType & data,
174 : const ReporterProducerEnum & mode,
175 : const SobolSampler & sampler)
176 : : ReporterGeneralContext<SobolState<OutType>>(other, producer, state),
177 49 : _data(data),
178 49 : _data_mode(mode),
179 49 : _sampler(sampler),
180 98 : _calc(other, "SOBOL", _sampler.resample())
181 : {
182 49 : }
183 :
184 : template <typename InType, typename OutType>
185 49 : SobolReporterContext<InType, OutType>::SobolReporterContext(
186 : const libMesh::ParallelObject & other,
187 : const MooseObject & producer,
188 : ReporterState<SobolState<OutType>> & state,
189 : const InType & data,
190 : const ReporterProducerEnum & mode,
191 : const SobolSampler & sampler,
192 : const MooseEnum & ci_method,
193 : const std::vector<Real> & ci_levels,
194 : unsigned int ci_replicates,
195 : unsigned int ci_seed)
196 49 : : SobolReporterContext<InType, OutType>(other, producer, state, data, mode, sampler)
197 : {
198 98 : _ci_calc_ptr =
199 : StochasticTools::makeBootstrapCalculator<std::vector<InType>, std::vector<OutType>>(
200 : ci_method, other, ci_levels, ci_replicates, ci_seed, _calc);
201 49 : }
202 :
203 : template <typename InType, typename OutType>
204 : void
205 49 : SobolReporterContext<InType, OutType>::finalize()
206 : {
207 49 : auto & val = this->_state.value(); // Convenience
208 49 : val.first.clear();
209 49 : val.second.clear();
210 :
211 49 : const bool is_dist = _data_mode == REPORTER_MODE_DISTRIBUTED;
212 49 : if (is_dist || this->processor_id() == 0)
213 : {
214 : const std::size_t ncol =
215 49 : _sampler.resample() ? 2 * _sampler.getNumberOfCols() + 2 : _sampler.getNumberOfCols() + 2;
216 49 : const std::vector<InType> data_reshape =
217 : StochasticTools::reshapeVector(_data, ncol, /*row_major =*/true);
218 :
219 98 : val.first = _calc.compute(data_reshape, is_dist);
220 49 : if (_ci_calc_ptr)
221 98 : val.second = _ci_calc_ptr->compute(data_reshape, is_dist);
222 49 : }
223 :
224 49 : ReporterGeneralContext<SobolState<OutType>>::finalize();
225 49 : }
226 :
227 : template <typename InType, typename OutType>
228 : void
229 35 : SobolReporterContext<InType, OutType>::store(nlohmann::json & json) const
230 : {
231 35 : storeSobol(json, this->_state.value(), _sampler.getNumberOfCols());
232 35 : }
233 :
234 : template <typename InType, typename OutType>
235 : void
236 70 : SobolReporterContext<InType, OutType>::storeSobol(nlohmann::json & json,
237 : const SobolState<OutType> & val,
238 : unsigned int nparam)
239 : {
240 70 : if (val.first.empty())
241 0 : return;
242 :
243 : // Convenience
244 70 : const unsigned int nlevels = val.second.size();
245 :
246 : std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>> first_order;
247 : std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>> total;
248 :
249 70 : first_order.first.assign(val.first.begin(), val.first.begin() + nparam);
250 70 : total.first.assign(val.first.begin() + nparam, val.first.begin() + nparam + nparam);
251 70 : if (nlevels > 0)
252 : {
253 35 : first_order.second.resize(nlevels);
254 35 : total.second.resize(nlevels);
255 140 : for (const auto & l : make_range(nlevels))
256 : {
257 105 : first_order.second[l].assign(val.second[l].begin(), val.second[l].begin() + nparam);
258 105 : total.second[l].assign(val.second[l].begin() + nparam,
259 : val.second[l].begin() + nparam + nparam);
260 : }
261 : }
262 70 : json["FIRST_ORDER"] = first_order;
263 140 : json["TOTAL"] = total;
264 :
265 70 : if (val.first.size() > 2 * nparam)
266 : {
267 : std::pair<std::vector<std::vector<OutType>>, std::vector<std::vector<std::vector<OutType>>>>
268 : second_order;
269 :
270 65 : second_order.first.resize(nparam, std::vector<OutType>(nparam));
271 395 : for (const auto & i : make_range(nparam))
272 330 : second_order.first[i][i] = val.first[i];
273 : unsigned int ind = 2 * nparam;
274 395 : for (const auto & i : make_range(nparam))
275 1035 : for (const auto & j : make_range(i + 1, nparam))
276 : {
277 705 : second_order.first[i][j] = val.first[ind++];
278 705 : second_order.first[j][i] = second_order.first[i][j];
279 : }
280 :
281 65 : if (nlevels > 0)
282 : {
283 30 : second_order.second.resize(nlevels, second_order.first);
284 :
285 125 : for (const auto & l : make_range(nlevels))
286 : {
287 625 : for (const auto & i : make_range(nparam))
288 530 : second_order.second[l][i][i] = val.second[l][i];
289 : ind = 2 * nparam;
290 625 : for (const auto & i : make_range(nparam))
291 1775 : for (const auto & j : make_range(i + 1, nparam))
292 : {
293 1245 : second_order.second[l][i][j] = val.second[l][ind++];
294 1245 : second_order.second[l][j][i] = second_order.second[l][i][j];
295 : }
296 : }
297 : }
298 130 : json["SECOND_ORDER"] = second_order;
299 : }
300 : }
301 :
302 : template void SobolReporter::declareValueHelper<Real>(const ReporterName & r_name);
303 : template class SobolReporterContext<std::vector<Real>, Real>;
304 : template void SobolReporter::declareValueHelper<std::vector<Real>>(const ReporterName & r_name);
305 : template class SobolReporterContext<std::vector<std::vector<Real>>, std::vector<Real>>;
|
