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 : #pragma once
11 :
12 : #include "libmesh/utility.h"
13 :
14 : #include "DataIO.h"
15 : #include "Uniform.h"
16 : #include "Normal.h"
17 : #include "CartesianProduct.h"
18 : #include "nlohmann/json.h"
19 :
20 : /**
21 : * Polynomials and quadratures based on defined distributions for Polynomial Chaos
22 : */
23 : namespace PolynomialQuadrature
24 : {
25 : class Polynomial;
26 :
27 : std::unique_ptr<const Polynomial> makePolynomial(const Distribution * dist);
28 :
29 : /**
30 : * General polynomial class with function for evaluating a polynomial of a given
31 : * order at a given point
32 : */
33 : class Polynomial
34 : {
35 : public:
36 3598 : Polynomial() {}
37 : virtual ~Polynomial() = default;
38 : virtual void store(std::ostream & stream, void * context) const;
39 : virtual void store(nlohmann::json & json) const;
40 : virtual Real compute(const unsigned int order, const Real x, const bool normalize = true) const;
41 : /// Computes the mth derivative of polynomial: d^mP_n/dx^m
42 : virtual Real
43 : computeDerivative(const unsigned int order, const Real x, const unsigned int m = 1) const;
44 : virtual Real innerProduct(const unsigned int order) const = 0;
45 :
46 : virtual void gaussQuadrature(const unsigned int order,
47 : std::vector<Real> & points,
48 : std::vector<Real> & weights) const = 0;
49 : virtual void clenshawQuadrature(const unsigned int order,
50 : std::vector<Real> & points,
51 : std::vector<Real> & weights) const;
52 : Real productIntegral(const std::vector<unsigned int> order) const;
53 : };
54 :
55 : /**
56 : * Uniform distributions use Legendre polynomials
57 : */
58 : class Legendre : public Polynomial
59 : {
60 : public:
61 : Legendre(const Real lower_bound, const Real upper_bound);
62 : virtual void store(std::ostream & stream, void * context) const override;
63 : virtual void store(nlohmann::json & json) const override;
64 :
65 : /// Legendre polynomial using static function then scales by <P_n^2> = 1 / (2n+1)
66 : virtual Real
67 : compute(const unsigned int order, const Real x, const bool normalize = true) const override;
68 : /**
69 : * Compute derivative of Legendre polynomial
70 : * Recursive algorithm: d^m/dx^m P_n = (n + m - 1)d^(m-1)/dx^(m-1) P_{n-1} + xd^m/dx^m P_{n-1}
71 : */
72 : virtual Real computeDerivative(const unsigned int order,
73 : const Real x,
74 : const unsigned int m = 1) const override;
75 : Real computeDerivativeRef(const unsigned int order, const Real x, const unsigned int m = 1) const;
76 : virtual Real innerProduct(const unsigned int order) const override;
77 :
78 : /// Gauss-Legendre quadrature: sum(weights) = 2
79 : virtual void gaussQuadrature(const unsigned int order,
80 : std::vector<Real> & points,
81 : std::vector<Real> & weights) const override;
82 :
83 : /// Just normal Clenshaw Curtis with shifted points
84 : virtual void clenshawQuadrature(const unsigned int order,
85 : std::vector<Real> & points,
86 : std::vector<Real> & weights) const override;
87 :
88 : private:
89 : const Real _lower_bound;
90 : const Real _upper_bound;
91 : };
92 :
93 : /**
94 : * Normal distributions use Hermite polynomials
95 : */
96 : class Hermite : public Polynomial
97 : {
98 : public:
99 : Hermite(const Real mu, const Real sig);
100 : virtual void store(std::ostream & stream, void * context) const override;
101 : virtual void store(nlohmann::json & json) const override;
102 :
103 : /// Hermite polynomial using static function then scales by <P_n^2> = n!
104 : virtual Real
105 : compute(const unsigned int order, const Real x, const bool normalize = true) const override;
106 :
107 : /// Compute derivative of Hermite polynomial P'_n = nP_{n-1}
108 : virtual Real computeDerivative(const unsigned int order,
109 : const Real x,
110 : const unsigned int m = 1) const override;
111 : virtual Real innerProduct(const unsigned int order) const override;
112 :
113 : /// Gauss-Hermite quadrature: sum(weights) = sqrt(2\pi)
114 : virtual void gaussQuadrature(const unsigned int order,
115 : std::vector<Real> & points,
116 : std::vector<Real> & weights) const override;
117 :
118 : private:
119 : const Real _mu;
120 : const Real _sig;
121 : };
122 :
123 : /**
124 : * General multidimensional quadrature class
125 : */
126 : class Quadrature
127 : {
128 : public:
129 416 : Quadrature() {}
130 : virtual ~Quadrature() = default;
131 :
132 : /// Resulting number of quadrature points in grid
133 : virtual unsigned int nPoints() const = 0;
134 : /// Inputted number of dimensions
135 : virtual unsigned int nDim() const = 0;
136 :
137 : /// Quadrature point n
138 : virtual std::vector<Real> quadraturePoint(const unsigned int n) const = 0;
139 : /// Quadrature point n for dimension dim
140 : virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const = 0;
141 : /// Weight for quadrature point n
142 : virtual Real quadratureWeight(const unsigned int n) const = 0;
143 : };
144 :
145 : /**
146 : * Full tensor product of 1D quadratures
147 : */
148 : class TensorGrid : public Quadrature
149 : {
150 : public:
151 : TensorGrid(const std::vector<unsigned int> & npoints,
152 : std::vector<std::unique_ptr<const Polynomial>> & poly);
153 372 : TensorGrid(const unsigned int npoints, std::vector<std::unique_ptr<const Polynomial>> & poly)
154 372 : : TensorGrid(std::vector<unsigned int>(poly.size(), npoints), poly)
155 : {
156 372 : }
157 :
158 372 : virtual unsigned int nPoints() const override { return _quad->numRows(); };
159 372 : virtual unsigned int nDim() const override { return _quad->numCols(); };
160 :
161 0 : virtual std::vector<Real> quadraturePoint(const unsigned int n) const override
162 : {
163 0 : return _quad->computeRow(n);
164 : };
165 1913080 : virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const override
166 : {
167 1913080 : return _quad->computeValue(n, dim);
168 : };
169 163570 : virtual Real quadratureWeight(const unsigned int n) const override
170 : {
171 163570 : return _quad->computeWeight(n);
172 : };
173 :
174 : private:
175 : std::unique_ptr<const StochasticTools::WeightedCartesianProduct<Real, Real>> _quad = nullptr;
176 : };
177 :
178 : /**
179 : * Smolyak sparse grid
180 : */
181 : class SmolyakGrid : public Quadrature
182 : {
183 : public:
184 : SmolyakGrid(const unsigned int max_order, std::vector<std::unique_ptr<const Polynomial>> & poly);
185 :
186 22 : virtual unsigned int nPoints() const override { return _npoints.back(); }
187 22 : virtual unsigned int nDim() const override { return _ndim; }
188 :
189 : virtual std::vector<Real> quadraturePoint(const unsigned int n) const override;
190 : virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const override;
191 : virtual Real quadratureWeight(const unsigned int n) const override;
192 :
193 : private:
194 : /// Helper function to find which quadrature product to use
195 : unsigned int gridIndex(const unsigned int n) const;
196 : const unsigned int _ndim;
197 : /// Cumulative number of points for each quadrature product
198 : std::vector<std::size_t> _npoints;
199 : /// Modification of quadrature weight based on polynomial order
200 : std::vector<int> _coeff;
201 : /// Container for all the combinations of quadrature products
202 : std::vector<std::unique_ptr<const StochasticTools::WeightedCartesianProduct<Real, Real>>> _quad;
203 : };
204 :
205 : /**
206 : * Clenshaw-Curtis sparse grid
207 : */
208 : class ClenshawCurtisGrid : public Quadrature
209 : {
210 : public:
211 : ClenshawCurtisGrid(const unsigned int max_order,
212 : std::vector<std::unique_ptr<const Polynomial>> & poly);
213 :
214 22 : virtual unsigned int nPoints() const override { return _quadrature.size(); }
215 22 : virtual unsigned int nDim() const override { return _ndim; }
216 :
217 0 : virtual std::vector<Real> quadraturePoint(const unsigned int n) const override
218 : {
219 0 : return _quadrature[n].first;
220 : }
221 84720 : virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const override
222 : {
223 84720 : return _quadrature[n].first[dim];
224 : }
225 7060 : virtual Real quadratureWeight(const unsigned int n) const override
226 : {
227 7060 : return _quadrature[n].second;
228 : }
229 :
230 : private:
231 : const unsigned int _ndim;
232 : std::vector<std::pair<std::vector<Real>, Real>> _quadrature;
233 : };
234 :
235 : /**
236 : * Legendre polynomial of specified order. Uses boost if available
237 : */
238 : Real legendre(const unsigned int order,
239 : const Real x,
240 : const Real lower_bound = -1.0,
241 : const Real upper_bound = 1.0);
242 :
243 : /**
244 : * Hermite polynomial of specified order. Uses boost if available.
245 : */
246 : Real hermite(const unsigned int order, const Real x, const Real mu = 0.0, const Real sig = 1.0);
247 :
248 : /**
249 : * Generalized formula for any polynomial order. Resulting number of points is then:
250 : * N = order + 1
251 : * The sum of weights is 2
252 : * Uses the Golub-Welsch algorithm:
253 : * - Golub, Gene & Welsch, John ``Calculation of Gauss Quadrature Rules''.
254 : * https://web.stanford.edu/class/cme335/S0025-5718-69-99647-1.pdf
255 : */
256 : void gauss_legendre(const unsigned int order,
257 : std::vector<Real> & points,
258 : std::vector<Real> & weights,
259 : const Real lower_bound,
260 : const Real upper_bound);
261 :
262 : /**
263 : * Generalized formula for any polynomial order. Resulting number of points is then:
264 : * N = order + 1
265 : * The sum of weights is sqrt(2*\pi)
266 : * Uses the Golub-Welsch algorithm:
267 : * - Golub, Gene & Welsch, John ``Calculation of Gauss Quadrature Rules''.
268 : * https://web.stanford.edu/class/cme335/S0025-5718-69-99647-1.pdf
269 : */
270 : void gauss_hermite(const unsigned int order,
271 : std::vector<Real> & points,
272 : std::vector<Real> & weights,
273 : const Real mu,
274 : const Real sig);
275 :
276 : void
277 : clenshaw_curtis(const unsigned int order, std::vector<Real> & points, std::vector<Real> & weights);
278 : }
279 :
280 : template <>
281 : void dataStore(std::ostream & stream,
282 : std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
283 : void * context);
284 : template <>
285 : void dataLoad(std::istream & stream,
286 : std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
287 : void * context);
|
