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 578 : PolynomialChaosReporter::validParams()
19 : {
20 578 : InputParameters params = GeneralReporter::validParams();
21 578 : params.addClassDescription(
22 : "Tool for extracting data from PolynomialChaos surrogates and computing statistics.");
23 1156 : params.addRequiredParam<std::vector<UserObjectName>>(
24 : "pc_name", "Name(s) of PolynomialChaos surrogate object(s).");
25 1156 : params.addParam<bool>("include_data",
26 1156 : false,
27 : "True to output information on the polynomial chaos model, including "
28 : "polynomial types, orders, and coefficients.");
29 578 : MultiMooseEnum stats("mean=1 stddev=2 skewness=3 kurtosis=4", "");
30 1156 : params.addParam<MultiMooseEnum>("statistics", stats, "Statistics to compute.");
31 1156 : params.addParam<bool>("include_sobol", false, "True to compute Sobol indices.");
32 1156 : params.addParam<std::vector<std::vector<Real>>>(
33 : "local_sensitivity_points",
34 1156 : std::vector<std::vector<Real>>(0),
35 : "Points for each polynomial chaos surrogate specifying desired location of sensitivity "
36 : "measurement.");
37 1156 : params.addParam<std::vector<SamplerName>>(
38 : "local_sensitivity_sampler",
39 1156 : std::vector<SamplerName>(0),
40 : "Sampler for each polynomial chaos surrogate specifying desired location of sensitivity "
41 : "measurement.");
42 578 : return params;
43 578 : }
44 :
45 289 : PolynomialChaosReporter::PolynomialChaosReporter(const InputParameters & parameters)
46 : : GeneralReporter(parameters),
47 : SurrogateModelInterface(this),
48 578 : _loc_point(getParam<std::vector<std::vector<Real>>>("local_sensitivity_points"))
49 : {
50 612 : for (const auto & sn : getParam<std::vector<SamplerName>>("local_sensitivity_sampler"))
51 34 : _loc_sampler.push_back(&getSamplerByName(sn));
52 :
53 935 : for (const auto & nm : getParam<std::vector<UserObjectName>>("pc_name"))
54 : {
55 : // Gather polynomial chaos surrogates
56 357 : _pc.push_back(&getSurrogateModelByName<PolynomialChaos>(nm));
57 : // PolynomialChaos data reporter
58 714 : if (getParam<bool>("include_data"))
59 : {
60 408 : auto & pc_ptr = declareValueByName<const PolynomialChaos *>(nm, REPORTER_MODE_ROOT);
61 136 : pc_ptr = _pc.back();
62 : }
63 : // Statistics reporter
64 1054 : for (const auto & stat : getParam<MultiMooseEnum>("statistics"))
65 1020 : declareValueByName<std::pair<Real, std::vector<Real>>, PCStatisticsContext<Real>>(
66 680 : nm + "_" + stat.name(), REPORTER_MODE_ROOT, *_pc.back(), stat);
67 : // Sobol reporter
68 714 : if (getParam<bool>("include_sobol"))
69 : declareValueByName<std::pair<std::vector<Real>, std::vector<std::vector<Real>>>,
70 476 : PCSobolContext<Real>>(nm + "_SOBOL", REPORTER_MODE_ROOT, *_pc.back());
71 : // Local sensitivity points reporter
72 357 : if (_loc_point.size() > 0)
73 : {
74 102 : 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 306 : _loc_point_sense.push_back(&declareValueByName<std::vector<std::vector<Real>>>(
78 204 : nm + "_POINT_SENSITIVITY", REPORTER_MODE_ROOT));
79 : }
80 : // Local sensitivity sampler reporter
81 357 : if (_loc_sampler.size() > 0)
82 : {
83 34 : 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 34 : _loc_sampler_sense.push_back(
87 34 : &declareValueByName<std::vector<std::vector<Real>>,
88 102 : StochasticReporterContext<std::vector<Real>>>(
89 68 : nm + "_SAMPLE_SENSITIVITY", REPORTER_MODE_ROOT, *_loc_sampler[_pc.size() - 1]));
90 : }
91 : }
92 289 : }
93 :
94 : void
95 289 : PolynomialChaosReporter::execute()
96 : {
97 289 : if ((_loc_point.size() + _loc_sampler.size()) == 0)
98 : return;
99 :
100 170 : for (const auto & i : index_range(_pc))
101 : {
102 102 : const auto & pc = *_pc[i];
103 : const auto nparam = pc.getNumberOfParameters();
104 102 : if (_loc_point.size() > 0)
105 : {
106 102 : 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 102 : const auto npoint = _loc_point[i].size() / nparam;
112 102 : _loc_point_sense[i]->resize(npoint);
113 306 : for (const auto & p : make_range(npoint))
114 : {
115 204 : const std::vector<Real> data(_loc_point[i].begin() + p * nparam,
116 204 : _loc_point[i].begin() + (p + 1) * nparam);
117 204 : (*_loc_point_sense[i])[p] = computeLocalSensitivity(pc, data);
118 204 : }
119 :
120 102 : if (_loc_sampler.size() > 0)
121 : {
122 34 : 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 : pc.name(),
129 : ".");
130 2234 : for (const auto & p : make_range(_loc_sampler[i]->getNumberOfLocalRows()))
131 2200 : (*_loc_sampler_sense[i])[p] =
132 4400 : computeLocalSensitivity(pc, _loc_sampler[i]->getNextLocalRow());
133 : }
134 : }
135 : }
136 : }
137 :
138 : std::vector<Real>
139 2404 : PolynomialChaosReporter::computeLocalSensitivity(const PolynomialChaos & pc,
140 : const std::vector<Real> & data)
141 : {
142 2404 : std::vector<Real> sense(data.size());
143 2404 : const auto val = pc.evaluate(data);
144 7348 : for (const auto & d : index_range(data))
145 4944 : sense[d] = data[d] / val * pc.computeDerivative(d, data);
146 2404 : return sense;
147 0 : }
148 :
149 : void
150 88 : to_json(nlohmann::json & json, const PolynomialChaos * const & pc)
151 : {
152 88 : pc->store(json);
153 88 : }
154 :
155 : template <typename OutType>
156 340 : 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 340 : _pc(pc),
164 340 : _stat(stat)
165 : {
166 340 : }
167 :
168 : template <typename OutType>
169 : void
170 340 : PCStatisticsContext<OutType>::finalize()
171 : {
172 : // Compute standard deviation
173 : OutType sig = OutType();
174 340 : if (_stat == "stddev" || _stat == "skewness" || _stat == "kurtosis")
175 204 : sig = _pc.computeStandardDeviation();
176 :
177 340 : OutType & val = this->_state.value().first;
178 340 : val = OutType();
179 : // Mean
180 340 : if (_stat == "mean" && this->processor_id() == 0)
181 88 : val = _pc.computeMean();
182 : // Standard Deviation
183 252 : else if (_stat == "stddev" && this->processor_id() == 0)
184 88 : val = sig;
185 : // Skewness
186 164 : else if (_stat == "skewness")
187 34 : val = _pc.powerExpectation(3) / (sig * sig * sig);
188 : // Kurtosis
189 130 : else if (_stat == "kurtosis")
190 34 : val = _pc.powerExpectation(4) / (sig * sig * sig * sig);
191 340 : this->_communicator.sum(val);
192 :
193 340 : ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>::finalize();
194 340 : }
195 :
196 : template <typename OutType>
197 : void
198 220 : PCStatisticsContext<OutType>::storeInfo(nlohmann::json & json) const
199 : {
200 220 : ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>::storeInfo(json);
201 220 : json["stat"] = _stat.name();
202 220 : }
203 :
204 : template <typename OutType>
205 119 : 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 119 : _pc(pc)
213 : {
214 0 : }
215 :
216 : template <typename OutType>
217 : void
218 119 : PCSobolContext<OutType>::finalize()
219 : {
220 119 : const unsigned int nparam = _pc.getNumberOfParameters();
221 119 : std::vector<OutType> & val = this->_state.value().first;
222 119 : val.clear();
223 :
224 : // Compute variance
225 119 : auto var = _pc.computeStandardDeviation();
226 119 : var *= var;
227 :
228 : // First order
229 697 : for (const auto & i : make_range(nparam))
230 1156 : val.push_back(_pc.computeSobolIndex({i}) / var);
231 : // Total
232 697 : for (const auto & i : make_range(nparam))
233 578 : val.push_back(_pc.computeSobolTotal(i) / var);
234 : // Second order
235 697 : for (const auto & i : make_range(nparam))
236 1751 : for (const auto & j : make_range(i + 1, nparam))
237 2346 : val.push_back(_pc.computeSobolIndex({i, j}) / var);
238 119 : }
239 :
240 : template <typename OutType>
241 : void
242 77 : PCSobolContext<OutType>::store(nlohmann::json & json) const
243 : {
244 77 : SobolReporterContext<std::vector<OutType>, OutType>::storeSobol(
245 77 : json, this->_state.value(), _pc.getNumberOfParameters());
246 77 : }
247 :
248 : template class PCStatisticsContext<Real>;
249 : template class PCSobolContext<Real>;
|
