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 864 : PolynomialChaos::validParams()
18 : {
19 864 : InputParameters params = SurrogateModel::validParams();
20 864 : params.addClassDescription("Computes and evaluates polynomial chaos surrogate model.");
21 864 : return params;
22 0 : }
23 :
24 432 : PolynomialChaos::PolynomialChaos(const InputParameters & parameters)
25 : : SurrogateModel(parameters),
26 432 : _order(getModelData<unsigned int>("_order")),
27 864 : _ndim(getModelData<unsigned int>("_ndim")),
28 864 : _ncoeff(getModelData<std::size_t>("_ncoeff")),
29 864 : _tuple(getModelData<std::vector<std::vector<unsigned int>>>("_tuple")),
30 864 : _coeff(getModelData<std::vector<Real>>("_coeff")),
31 432 : _poly(
32 864 : getModelData<std::vector<std::unique_ptr<const PolynomialQuadrature::Polynomial>>>("_poly"))
33 : {
34 432 : }
35 :
36 : Real
37 5792 : 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 5792 : DenseMatrix<Real> poly_val(_ndim, _order);
42 :
43 : // Evaluate polynomials to avoid duplication
44 18504 : for (unsigned int d = 0; d < _ndim; ++d)
45 68752 : for (unsigned int i = 0; i < _order; ++i)
46 56040 : poly_val(d, i) = _poly[d]->compute(i, x[d], /*normalize =*/false);
47 :
48 : Real val = 0;
49 129072 : for (std::size_t i = 0; i < _ncoeff; ++i)
50 : {
51 123280 : Real tmp = _coeff[i];
52 462360 : for (unsigned int d = 0; d < _ndim; ++d)
53 339080 : tmp *= poly_val(d, _tuple[i][d]);
54 123280 : val += tmp;
55 : }
56 :
57 5792 : return val;
58 5792 : }
59 :
60 : const std::vector<std::vector<unsigned int>> &
61 80 : PolynomialChaos::getPolynomialOrders() const
62 : {
63 80 : 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 80 : PolynomialChaos::getCoefficients() const
74 : {
75 80 : return _coeff;
76 : }
77 :
78 : Real
79 80 : PolynomialChaos::computeMean() const
80 : {
81 : mooseAssert(_coeff.size() > 0, "The coefficient matrix is empty.");
82 80 : return _coeff[0];
83 : }
84 :
85 : Real
86 304 : PolynomialChaos::computeStandardDeviation() const
87 : {
88 : Real var = 0;
89 97024 : for (std::size_t i = 1; i < _ncoeff; ++i)
90 : {
91 : Real norm = 1.0;
92 488880 : for (std::size_t d = 0; d < _ndim; ++d)
93 392160 : norm *= _poly[d]->innerProduct(_tuple[i][d]);
94 96720 : var += _coeff[i] * _coeff[i] * norm;
95 : }
96 :
97 304 : return std::sqrt(var);
98 : }
99 :
100 : Real
101 64 : PolynomialChaos::powerExpectation(const unsigned int n) const
102 : {
103 : std::vector<StochasticTools::WeightedCartesianProduct<unsigned int, Real>> order;
104 64 : order.reserve(_ndim);
105 64 : std::vector<std::vector<Real>> c_1d(n, std::vector<Real>(_coeff.begin() + 1, _coeff.end()));
106 192 : for (unsigned int d = 0; d < _ndim; ++d)
107 : {
108 128 : std::vector<std::vector<unsigned int>> order_1d(n, std::vector<unsigned int>(_ncoeff - 1));
109 1920 : for (std::size_t i = 1; i < _ncoeff; ++i)
110 8064 : for (unsigned int m = 0; m < n; ++m)
111 6272 : order_1d[m][i - 1] = _tuple[i][d];
112 128 : order.push_back(StochasticTools::WeightedCartesianProduct<unsigned int, Real>(order_1d, c_1d));
113 128 : }
114 :
115 : dof_id_type n_local, st_local, end_local;
116 64 : MooseUtils::linearPartitionItems(
117 : order[0].numRows(), n_processors(), processor_id(), n_local, st_local, end_local);
118 :
119 : Real val = 0;
120 823264 : for (dof_id_type i = st_local; i < end_local; ++i)
121 : {
122 823200 : Real tmp = order[0].computeWeight(i);
123 1215450 : for (unsigned int d = 0; d < _ndim; ++d)
124 : {
125 1133610 : if (MooseUtils::absoluteFuzzyEqual(tmp, 0.0))
126 : break;
127 392250 : std::vector<unsigned int> comb(n);
128 1920610 : for (unsigned int m = 0; m < n; ++m)
129 1528360 : comb[m] = order[d].computeValue(i, m);
130 392250 : tmp *= _poly[d]->productIntegral(comb);
131 : }
132 823200 : val += tmp;
133 : }
134 :
135 64 : return val;
136 64 : }
137 :
138 : Real
139 4512 : PolynomialChaos::computeDerivative(const unsigned int dim, const std::vector<Real> & x) const
140 : {
141 9024 : return computePartialDerivative({dim}, x);
142 : }
143 :
144 : Real
145 4512 : 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 4512 : std::vector<unsigned int> grad(_ndim);
151 9024 : for (const auto & d : dim)
152 : {
153 : mooseAssert(d < _ndim, "Specified dimension is greater than total number of parameters.");
154 4512 : grad[d]++;
155 : }
156 :
157 4512 : DenseMatrix<Real> poly_val(_ndim, _order);
158 :
159 : // Evaluate polynomials to avoid duplication
160 14048 : for (unsigned int d = 0; d < _ndim; ++d)
161 62336 : for (unsigned int i = 0; i < _order; ++i)
162 52800 : poly_val(d, i) = _poly[d]->computeDerivative(i, x[d], grad[d]);
163 :
164 : Real val = 0;
165 251392 : for (std::size_t i = 0; i < _ncoeff; ++i)
166 : {
167 246880 : Real tmp = _coeff[i];
168 1106720 : for (unsigned int d = 0; d < _ndim; ++d)
169 859840 : tmp *= poly_val(d, _tuple[i][d]);
170 246880 : val += tmp;
171 : }
172 :
173 4512 : return val;
174 4512 : }
175 :
176 : Real
177 1648 : PolynomialChaos::computeSobolIndex(const std::set<unsigned int> & ind) const
178 : {
179 1648 : linearPartitionCoefficients();
180 :
181 : // If set is empty, compute mean
182 1648 : 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 340568 : for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
191 : {
192 338920 : Real tmp = _coeff[i] * _coeff[i];
193 633200 : for (unsigned int d = 0; d < _ndim; ++d)
194 : {
195 621230 : if ((ind.find(d) != ind.end() && _tuple[i][d] > 0) ||
196 401440 : (ind.find(d) == ind.end() && _tuple[i][d] == 0))
197 : {
198 294280 : tmp *= _poly[d]->innerProduct(_tuple[i][d]);
199 : }
200 : else
201 : {
202 : tmp = 0.0;
203 : break;
204 : }
205 : }
206 338920 : val += tmp;
207 : }
208 :
209 : return val;
210 : }
211 :
212 : Real
213 544 : PolynomialChaos::computeSobolTotal(const unsigned int dim) const
214 : {
215 544 : 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 130064 : for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
222 129520 : if (_tuple[i][dim] > 0)
223 84240 : val += _coeff[i] * _coeff[i] * _poly[dim]->innerProduct(_tuple[i][dim]);
224 :
225 544 : return val;
226 : }
227 :
228 : void
229 2192 : PolynomialChaos::linearPartitionCoefficients() const
230 : {
231 2192 : if (_n_local_coeff == std::numeric_limits<dof_id_type>::max())
232 112 : MooseUtils::linearPartitionItems(_ncoeff,
233 : n_processors(),
234 : processor_id(),
235 112 : _n_local_coeff,
236 112 : _local_coeff_begin,
237 112 : _local_coeff_end);
238 2192 : }
239 :
240 : void
241 80 : PolynomialChaos::store(nlohmann::json & json) const
242 : {
243 160 : json["order"] = _order;
244 160 : json["ndim"] = getNumberOfParameters();
245 80 : json["ncoeff"] = getNumberofCoefficients();
246 80 : json["tuple"] = getPolynomialOrders();
247 80 : json["coeff"] = getCoefficients();
248 320 : for (const auto & p : _poly)
249 : {
250 : nlohmann::json jsonp;
251 240 : p->store(jsonp);
252 240 : json["poly"].push_back(jsonp);
253 : }
254 80 : }
|
