PolynomialChaos

Computes and evaluates polynomial chaos surrogate model.

Overview

Polynomial chaos is a surrogate modeling technique where a quantity of interest (QoI) that is dependent on input parameters is expanded as a sum of orthogonal polynomials. Given a QoI dependent on a set of parameters , the polynomial chaos expansion (PCE) is:

(1)

where is the multidimensional polynomial order and are coefficients that are to be computed. These coefficients can be found using intrusive and non intrusive techniques. The intrusive technique is quite difficult to generalize and very computationally demanding. Since the polynomial basis is orthogonal, a non intrusive technique is developed where the coefficients are found by performing a Galerkin projection and integrating:

(2)

where,

The weight function () and bases () are typically products of one-dimensional functions:

The weighting functions are defined by the probability density function of the parameter and the polynomials are based on these distributions, Table 1 is a list of commonly used distributions and their corresponding orthogonal polynomials.

Table 1: Common probability density functions and their corresponding orthogonal polynomials

DistributionDensity Function ()Polynomial ()Support
NormalHermite
UniformLegendre
BetaJacobi
ExponentialLaguerre
GammaGeneralized Laguerre

Computing Coefficients

There are currently two techniques to compute the coefficients in Eq. (1), the first performs the integration described by Eq. (2). For Monte Carlo sampling, the integration is in the form

And for using a QuadratureSampler

The numerical quadrature method is typically much more efficient that than the Monte Carlo method and has the added benefit of exactly integrating the polynomial basis. However, the quadrature suffers from the curse of dimensionality. The naive approach uses a Cartesian product of one-dimensional quadratures, which results in quadrature points to be sampled. Sparse grids can help mitigate the curse of dimensionality significantly.

The other technique is using ordinary least-squares (OLS) regression, like in PolynomialRegressionTrainer. In the majority of cases, OLS is more accurate than integration for Monte Carlo sampling, as shown in the figure below.

Comparison between OLS regression and integration methods for computing polynomial chaos coefficients

Generating a Tuple

In polynomial chaos, a tuple describes the combination of polynomial orders representing the expansion basis (). Again, the naive approach would be to do a tensor product of highest polynomial order, but this is often wasteful since generating a complete monomial basis is usually optimal. Below demonstrates the difference between a tensor basis and a complete monomial basis:

The tuple is generated and stored as matrix in the userobject, below is an example of this matrix with and :

Computing Statistical Moments

Statistical moments are based on the expectation of a function of the quantity of interest:

The first four statistical moments, and the most common ones, are defined as:

Because of the orthogonality of the polynomials, mean and variance are trivial to compute:

where is known analytically. The higher order moments are significantly more taxing to compute since it does not take advantage of orthogonality:

where the polynomial norms are computed using one-dimensional quadrature. We see here the number of operations required to compute Kurtosis is approximately . If the number of coefficients is sufficiently high, these moments would probably be best computed inexactly by sampling the PCE surrogate.

Implementation

The PolynomialChaos user object takes in a list of distributions and constructs a polynomial class based on their type. Given a sampler and a vectorpostprocessor of results from sampling, it then loops through the MC or quadrature points to compute the coefficients. The statistical moments are then computed based on user preferences and the model can be evaluated using the evaluate function from any other moose object that has the reference. The algorithm uses the parallelization of the sampler to compute the coefficients, no other part of the algorithm is parallelized.

Example Input File Syntax

The example involves a homogeneous, one-dimensional diffusion-reaction problem, where the diffusion coefficient () and reaction coefficient () are uncertain with a uniform probability:

where the QoI is the integrated average of .

After the sub-app is set with the diffusion-reaction problem, distributions are set for each uncertain parameter:

[Distributions]
  [D_dist]
    type = Uniform
    lower_bound = 2.5
    upper_bound = 7.5
  []
  [S_dist]
    type = Uniform
    lower_bound = 2.5
    upper_bound = 7.5
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

A sampler is then defined, either using Monte Carlo,

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 100
    distributions = 'D_dist S_dist'
    execute_on = initial
  []
[]

[Trainers]
  [poly_chaos]
    type = PolynomialChaosTrainer
    execute_on = timestep_end
    order = 5
    distributions = 'D_dist S_dist'
    sampler = sample
    response = storage/data:avg:value
    regression_type = integration
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_mc.i)

or quadrature,

[Samplers]
  [quadrature]
    type = Quadrature
    distributions = 'D_dist S_dist'
    execute_on = INITIAL
    order = 5
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

It is important that the order in the quadrature sampler input matches the order in the PolynomialChaos input. The sampler is then used by the MultiApp and Transfers to sample the sub-app, the QoI from the app is then put in a reporter:

[Reporters]
  [storage]
    type = StochasticReporter
  []
  [pc_samp]
    type = EvaluateSurrogate
    model = poly_chaos
    sampler = sample
    parallel_type = ROOT
    execute_on = final
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

All this information is ready to be sent to the PolynomialChaos trainer:

[Surrogates]
  [poly_chaos]
    type = PolynomialChaos
    trainer = poly_chaos
  []
[]

[Trainers]
  [poly_chaos]
    type = PolynomialChaosTrainer
    execute_on = timestep_end
    order = 5
    distributions = 'D_dist S_dist'
    sampler = quadrature
    response = storage/data:avg:value
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

Input Parameters

  • filenameThe name of the file which will be associated with the saved/loaded data.

    C++ Type:FileName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the file which will be associated with the saved/loaded data.

  • trainerThe SurrogateTrainer object. If this is specified the trainer data is automatically gathered and available in this SurrogateModel object.

    C++ Type:UserObjectName

    Unit:(no unit assumed)

    Controllable:No

    Description:The SurrogateTrainer object. If this is specified the trainer data is automatically gathered and available in this SurrogateModel object.

Optional Parameters

  • control_tagsAdds user-defined labels for accessing object parameters via control logic.

    C++ Type:std::vector<std::string>

    Unit:(no unit assumed)

    Controllable:No

    Description:Adds user-defined labels for accessing object parameters via control logic.

  • enableTrueSet the enabled status of the MooseObject.

    Default:True

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Set the enabled status of the MooseObject.

Advanced Parameters