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>;
|