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