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 "PolynomialChaos.h"
11 : #include "Sampler.h"
12 : #include "CartesianProduct.h"
13 :
14 : registerMooseObject("StochasticToolsApp", PolynomialChaos);
15 :
16 : InputParameters
17 394 : PolynomialChaos::validParams()
18 : {
19 394 : InputParameters params = SurrogateModel::validParams();
20 394 : params.addClassDescription("Computes and evaluates polynomial chaos surrogate model.");
21 394 : return params;
22 0 : }
23 :
24 197 : PolynomialChaos::PolynomialChaos(const InputParameters & parameters)
25 : : SurrogateModel(parameters),
26 197 : _order(getModelData<unsigned int>("_order")),
27 394 : _ndim(getModelData<unsigned int>("_ndim")),
28 394 : _ncoeff(getModelData<std::size_t>("_ncoeff")),
29 394 : _tuple(getModelData<std::vector<std::vector<unsigned int>>>("_tuple")),
30 394 : _coeff(getModelData<std::vector<Real>>("_coeff")),
31 197 : _poly(
32 394 : getModelData<std::vector<std::unique_ptr<const PolynomialQuadrature::Polynomial>>>("_poly"))
33 : {
34 197 : }
35 :
36 : Real
37 3034 : PolynomialChaos::evaluate(const std::vector<Real> & x) const
38 : {
39 : mooseAssert(x.size() == _ndim, "Number of inputted parameters does not match PC model.");
40 :
41 3034 : DenseMatrix<Real> poly_val(_ndim, _order);
42 :
43 : // Evaluate polynomials to avoid duplication
44 9958 : for (unsigned int d = 0; d < _ndim; ++d)
45 35304 : for (unsigned int i = 0; i < _order; ++i)
46 28380 : poly_val(d, i) = _poly[d]->compute(i, x[d], /*normalize =*/false);
47 :
48 : Real val = 0;
49 61844 : for (std::size_t i = 0; i < _ncoeff; ++i)
50 : {
51 58810 : Real tmp = _coeff[i];
52 217270 : for (unsigned int d = 0; d < _ndim; ++d)
53 158460 : tmp *= poly_val(d, _tuple[i][d]);
54 58810 : val += tmp;
55 : }
56 :
57 3034 : return val;
58 3034 : }
59 :
60 : const std::vector<std::vector<unsigned int>> &
61 40 : PolynomialChaos::getPolynomialOrders() const
62 : {
63 40 : return _tuple;
64 : }
65 :
66 : unsigned
67 0 : PolynomialChaos::getPolynomialOrder(const unsigned int dim, const unsigned int i) const
68 : {
69 0 : return _tuple[i][dim];
70 : }
71 :
72 : const std::vector<Real> &
73 40 : PolynomialChaos::getCoefficients() const
74 : {
75 40 : return _coeff;
76 : }
77 :
78 : Real
79 40 : PolynomialChaos::computeMean() const
80 : {
81 : mooseAssert(_coeff.size() > 0, "The coefficient matrix is empty.");
82 40 : return _coeff[0];
83 : }
84 :
85 : Real
86 133 : PolynomialChaos::computeStandardDeviation() const
87 : {
88 : Real var = 0;
89 42448 : for (std::size_t i = 1; i < _ncoeff; ++i)
90 : {
91 : Real norm = 1.0;
92 213885 : for (std::size_t d = 0; d < _ndim; ++d)
93 171570 : norm *= _poly[d]->innerProduct(_tuple[i][d]);
94 42315 : var += _coeff[i] * _coeff[i] * norm;
95 : }
96 :
97 133 : return std::sqrt(var);
98 : }
99 :
100 : Real
101 28 : PolynomialChaos::powerExpectation(const unsigned int n) const
102 : {
103 : std::vector<StochasticTools::WeightedCartesianProduct<unsigned int, Real>> order;
104 28 : order.reserve(_ndim);
105 28 : std::vector<std::vector<Real>> c_1d(n, std::vector<Real>(_coeff.begin() + 1, _coeff.end()));
106 84 : for (unsigned int d = 0; d < _ndim; ++d)
107 : {
108 56 : std::vector<std::vector<unsigned int>> order_1d(n, std::vector<unsigned int>(_ncoeff - 1));
109 840 : for (std::size_t i = 1; i < _ncoeff; ++i)
110 3528 : for (unsigned int m = 0; m < n; ++m)
111 2744 : order_1d[m][i - 1] = _tuple[i][d];
112 56 : order.push_back(StochasticTools::WeightedCartesianProduct<unsigned int, Real>(order_1d, c_1d));
113 56 : }
114 :
115 : dof_id_type n_local, st_local, end_local;
116 28 : MooseUtils::linearPartitionItems(
117 : order[0].numRows(), n_processors(), processor_id(), n_local, st_local, end_local);
118 :
119 : Real val = 0;
120 411628 : for (dof_id_type i = st_local; i < end_local; ++i)
121 : {
122 411600 : Real tmp = order[0].computeWeight(i);
123 607725 : for (unsigned int d = 0; d < _ndim; ++d)
124 : {
125 566805 : if (MooseUtils::absoluteFuzzyEqual(tmp, 0.0))
126 : break;
127 196125 : std::vector<unsigned int> comb(n);
128 960305 : for (unsigned int m = 0; m < n; ++m)
129 764180 : comb[m] = order[d].computeValue(i, m);
130 196125 : tmp *= _poly[d]->productIntegral(comb);
131 196125 : }
132 411600 : val += tmp;
133 : }
134 :
135 28 : return val;
136 28 : }
137 :
138 : Real
139 2224 : PolynomialChaos::computeDerivative(const unsigned int dim, const std::vector<Real> & x) const
140 : {
141 2224 : return computePartialDerivative({dim}, x);
142 : }
143 :
144 : Real
145 2224 : PolynomialChaos::computePartialDerivative(const std::vector<unsigned int> & dim,
146 : const std::vector<Real> & x) const
147 : {
148 : mooseAssert(x.size() == _ndim, "Number of inputted parameters does not match PC model.");
149 :
150 2224 : std::vector<unsigned int> grad(_ndim);
151 4448 : for (const auto & d : dim)
152 : {
153 : mooseAssert(d < _ndim, "Specified dimension is greater than total number of parameters.");
154 2224 : grad[d]++;
155 : }
156 :
157 2224 : DenseMatrix<Real> poly_val(_ndim, _order);
158 :
159 : // Evaluate polynomials to avoid duplication
160 6896 : for (unsigned int d = 0; d < _ndim; ++d)
161 30272 : for (unsigned int i = 0; i < _order; ++i)
162 25600 : poly_val(d, i) = _poly[d]->computeDerivative(i, x[d], grad[d]);
163 :
164 : Real val = 0;
165 113984 : for (std::size_t i = 0; i < _ncoeff; ++i)
166 : {
167 111760 : Real tmp = _coeff[i];
168 495440 : for (unsigned int d = 0; d < _ndim; ++d)
169 383680 : tmp *= poly_val(d, _tuple[i][d]);
170 111760 : val += tmp;
171 : }
172 :
173 2224 : return val;
174 2224 : }
175 :
176 : Real
177 721 : PolynomialChaos::computeSobolIndex(const std::set<unsigned int> & ind) const
178 : {
179 721 : linearPartitionCoefficients();
180 :
181 : // If set is empty, compute mean
182 721 : if (ind.empty())
183 0 : return computeMean();
184 :
185 : // Do some sanity checks in debug
186 : mooseAssert(ind.size() <= _ndim, "Number of indices is greater than number of parameters.");
187 : mooseAssert(*ind.rbegin() < _ndim, "Maximum index provided exceeds number of parameters.");
188 :
189 : Real val = 0.0;
190 170181 : for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
191 : {
192 169460 : Real tmp = _coeff[i] * _coeff[i];
193 316600 : for (unsigned int d = 0; d < _ndim; ++d)
194 : {
195 310615 : if ((ind.find(d) != ind.end() && _tuple[i][d] > 0) ||
196 200720 : (ind.find(d) == ind.end() && _tuple[i][d] == 0))
197 : {
198 147140 : tmp *= _poly[d]->innerProduct(_tuple[i][d]);
199 : }
200 : else
201 : {
202 : tmp = 0.0;
203 : break;
204 : }
205 : }
206 169460 : val += tmp;
207 : }
208 :
209 : return val;
210 : }
211 :
212 : Real
213 238 : PolynomialChaos::computeSobolTotal(const unsigned int dim) const
214 : {
215 238 : linearPartitionCoefficients();
216 :
217 : // Do some sanity checks in debug
218 : mooseAssert(dim < _ndim, "Requested dimension is greater than number of parameters.");
219 :
220 : Real val = 0.0;
221 64998 : for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
222 64760 : if (_tuple[i][dim] > 0)
223 42120 : val += _coeff[i] * _coeff[i] * _poly[dim]->innerProduct(_tuple[i][dim]);
224 :
225 238 : return val;
226 : }
227 :
228 : void
229 959 : PolynomialChaos::linearPartitionCoefficients() const
230 : {
231 959 : if (_n_local_coeff == std::numeric_limits<dof_id_type>::max())
232 49 : MooseUtils::linearPartitionItems(_ncoeff,
233 : n_processors(),
234 : processor_id(),
235 49 : _n_local_coeff,
236 49 : _local_coeff_begin,
237 49 : _local_coeff_end);
238 959 : }
239 :
240 : void
241 40 : PolynomialChaos::store(nlohmann::json & json) const
242 : {
243 80 : json["order"] = _order;
244 80 : json["ndim"] = getNumberOfParameters();
245 40 : json["ncoeff"] = getNumberofCoefficients();
246 40 : json["tuple"] = getPolynomialOrders();
247 40 : json["coeff"] = getCoefficients();
248 160 : for (const auto & p : _poly)
249 : {
250 : nlohmann::json jsonp;
251 120 : p->store(jsonp);
252 120 : json["poly"].push_back(jsonp);
253 : }
254 40 : }
|
