Lack of convergence in Sobol indices and standard deviation using Gauss quadrature in chaospy

[Masilokwa-bantu]

Hello Jonathan,

Firstly, thank you for developing chaospy and its great documentation.I’m using it for a sensitivity analysis of a CFD (finite volume) model coupled with an ODE. There are 3 uniformly distributed uncertain parameters shown below:

R1 = [47975880.0, 79959800.0] R2 = [300844875.0, 501408125.0] C = [3.07206e-09, 5.120099999999999e-09]

I used gauss quadrature to generate gauss nodes to evaluate the CFD for polynomial orders 1 to 3. The mean, standard deviation, and Sobol indices were computed using chaospy.E(), chaospy.Std(), and chaospy.Sens_m() respectively.

However, the Sobol indices and standard deviation do not show convergence or a consistent trend with increasing polynomial order. The figure attached shows that the indices vary between the polynomial order without a similar trend.

For reference the convergence metrics are shown below (the metric is the mean relative error between successive polynomial orders):

pressure mean convergence: [12.1479473   0.77818001] %
pressure std convergence: [83.37969258 54.11365791] %
pressure SI convergence: [4.79972721 3.39656272] %

A 3rd-order expansion was expected to be sufficient, as observed in similar studies (albeit with two uncertain parameters and finite element models). The same methodology for sensitivity analysis was used.

Below is the function used to conduct the analysis.

def gauss_sobol_analysis_data(distribution, samples_evals, coeff_order, output=1, rule="gaussian"):
    gauss_quads = chaospy.generate_quadrature(coeff_order, distribution, rule=rule)
    nodes, weights = gauss_quads
    gauss_evals = samples_evals
    gauss_expansion = chaospy.generate_expansion(coeff_order, distribution) 
    gauss_model_approx = chaospy.fit_quadrature(gauss_expansion, nodes, weights, gauss_evals[:, output, -51:-1])

    mean = chaospy.E(gauss_model_approx, distribution)
    std = chaospy.Std(gauss_model_approx, distribution)
    S1 = chaospy.Sens_m(gauss_model_approx, distribution)

    return mean, std, S1, gauss_model_approx

Any insights into possible causes for the lack of convergence or suggestions for improving stability would be greatly appreciated.

Thank you in advance!

[jonathf]

The scale of C is many magnitudes smaller than R1 and R2, making for bad condition number for the underlying methods.

Try switching to generalized polynomial chaos and see if that solves the issue.

[jonathf]

Use J(Uniform(-1, 1), Uniform(-1, 1), Uniform(-1, 1)) as your base.

[Masilokwa-bantu]

The scale of C is many magnitudes smaller than R1 and R2, making for bad condition number for the underlying methods.

Try switching to generalized polynomial chaos and see if that solves the issue.

I tried switching to the generalised polynomial chaos method. The function that I used to calculate the results shown in the figure (attached) and convergence metrics is shown below:

def general_PCE_sobol_analysis_data(distribution, gauss_evals, coeff_order, output=1, rule="gaussian"):
    distribution_q = distribution
    distribution_r = chaospy.J(chaospy.Uniform(-1, 1), chaospy.Uniform(-1, 1), chaospy.Uniform(-1, 1))

    nodes_r, weights_r = chaospy.generate_quadrature(coeff_order, distribution_r, rule=rule)
    gauss_expansion = chaospy.generate_expansion(coeff_order, distribution)
    gauss_model_approx = chaospy.fit_quadrature(gauss_expansion, nodes_r, weights_r, gauss_evals[:, output, -51:-1])

    mean = chaospy.E(gauss_model_approx, distribution_r)
    std = chaospy.Std(gauss_model_approx, distribution_r)
    S1 = chaospy.Sens_m(gauss_model_approx, distribution_r)
    return mean, std, S1, gauss_model_approx

Convergence for the mean and standard deviation increased significantly. While the sobol indices decreased. Though as seen in the figure (attached), the sobol indices are constant for the entire duration. With C fixed at 1 and R1 and R2 fixed at 0. Which ideally should not be the case.

Convergence metrics: pressure mean convergence: [59.79709013 43.36904433] % pressure std convergence: [87.75654617 66.52619814] % pressure SI convergence: [1.03260140e-15 2.99694806e-15] %

I left this out earlier but thought that it may help. When calculating the mean and std using chaospy.E() and chaospy.Std() the magnitude of the pressure changes. I have included plots of the mean and std calculated directly from the data. The overall shape of the profiles is maintained, though.

[Masilokwa-bantu]

Forgot to attach the image to my earlier comment. Above are the plots that I was referring to.

Thank you for your quick response earlier.

[jonathf]

My thoughts looking through what you have presented here:

• The gauss_evals seem detatched from gauss_quads. Can we assume that the former is generated using the latter? Is that also the case for the transformed case?
• Do you have the same issue with collocation method, or is it unique for quadrature?
• You have uniform distribution, does switching to Clenshaw-Curtis samples make a difference?

[Masilokwa-bantu]

• The gauss_evals are generated using the nodes from gauss_quads. The simulations are run separately using OpenFOAM by evaluating the CFD model on the set of gauss nodes for polynomial orders 1-3. For your suggestion to use J(Uniform(-1, 1), Uniform(-1, 1), Uniform(-1, 1)) (distribution_r) as my base, I checked that the transformed gauss nodes (distribution_q nodes_q) were the same as the ones I already had. So I presumed that I could use the gauss_evals I already had with distribution_r.

• The point collocation method from the attached figure seems to be worse than gauss quadrature, as the 2nd and 3rd order means and stds do not resemble the ones calculated from the data. The results from the figure are using the point collocation method without the generalised PCE. Using the generalised PCE for the point collocation method, I get the same figure as previously.

Added the regression parameter and if statement to switch between the point collocation and gauss quadrature.

def gauss_sobol_analysis_data(distribution, gauss_evals, coeff_order, output=1, rule="gaussian", regression=False):
    gauss_quads = chaospy.generate_quadrature(coeff_order, distribution, rule=rule)
    nodes, weights = gauss_quads
    gauss_expansion = chaospy.generate_expansion(coeff_order, distribution)
    if regression:
        gauss_model_approx = chaospy.fit_regression(gauss_expansion, nodes, gauss_evals[:, output, -51:-1])
    else:
        gauss_model_approx = chaospy.fit_quadrature(gauss_expansion, nodes, weights, gauss_evals[:, output, -51:-1])

    mean = chaospy.E(gauss_model_approx, distribution)
    std = chaospy.Std(gauss_model_approx, distribution)
    S1 = chaospy.Sens_m(gauss_model_approx, distribution)

    return mean, std, S1, gauss_model_approx

• I haven't tried using Clenshaw-Curtis samples. Would they be the same as the current gauss nodes? If they are different, it would mean I need to rerun the CFD model for the new sample set, which could take a week.

• Could it be the case that the PCE is not able to approximate the non-linearity of the results? Or is it too early to make that kind of judgement?

[jonathf]

Collocation method is generally more robust than quadrature, so I doubt Clenshaw-Curtis is going to help you now.

Your hypothesis about non-linearity is a much more likely candidate.

The best I can suggest is trying to increase the order to 4, but considering the cost of your evaluation, I understand if that is too high for to just test an hypothesis.

[Masilokwa-bantu]

Hello Jonathan,

Hope your week is off to a good start. Thank you for the suggestion. I tried running simulations for 3rd and 4th order for only 2 uncertain inputs using gauss quadrature. The results converged between the 3rd and 4th order. So I will give running a 4th order result for 3 uncertain inputs a try; hopefully the results will also converge. Otherwise, thank you for having made the time to answer my questions.