# Exploring the interactions between experiments and simulations using Bayesian polynomials

The LS89 turbine geometry is a well-known standard test case for turbomachinery flow analysis. Fluid dynamicists from the von Karman Institute have published reports on experimental results, including one from Arts et al. in the 1990’s[1].

Fast forward to the 2010’s, the developers of the open source computational fluid dynamics (CFD) code SU2 have picked this test geometry up as a calibration test case [2], and various other studies have performed further numerical studies on the same geometry (such as shape optimisation in [3]).

The availability of both simulation and experimental studies makes this the ideal example to explore the trade-off between the two paradigms prevalent in engineering. Given a certain supply of experimental and simulation data, how can one construct the best model? While exprimental data measures the quantities of interest in a real-world setting, the setup of experiments can be costly. Meanwhile, setting up simulations is comparatively quick and straightforward. However, due to the complexity of flow phenomena, there is inevitably some amount of approximations and assumptions one needs to make in a CFD simulation (one exmaple being the modelling of turbulence). In this scenario, an approach that would fuse together these two sources of data, while noting their respective limitations, is valuable.

So, where does Bayesian polynomial come into this picture? In recent work by our group (preprint available at [2203.03508] Prior-informed Uncertainty Modelling with Bayesian Polynomial Approximations), we build upon classical ideas in polynomial approximations and extend these formalisms in a Bayesian perspective, and find that the Bayesian approach offers a quantitative way of measuring the uncertainty of the data we are given in a polynomial approximation. As a consequence, natural methods to fuse together different sources of knowledge can be derived. For example, one can provide information to a polynomial approximation based on experimental data using simulation data, and this improves the accuracy of the fit. In this post, we explore an example of this approach by analysing the experimental and simulation data of the LS89 geometry, using code from equadratures.

## Bayesian polynomial approximations

(Follow along this section with the notebook here: )

Let’s start by constructing a Bayesian polynomial with the experimental data for the LS89. In Arts et al., 21 data points are provided for a test case measuring the wall heat transfer (H_{int}), as a function of three variable boundary conditions: the turbulence intensity (Ti), exit Mach number (Ma) and exit Reynolds number (Re). Reserving 15 data points for training and 6 data points to test the accuracy of the model, we fit a cubic Bayesian polynomial to the data. The following code loads in the data…

import equadratures as eq
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

no_of_training_samples = 15
my_rng = np.random.default_rng(0)

random_order = my_rng.permutation(len(vki_data['experimental_outputs']))
train_inds = random_order[:no_of_training_samples]
test_inds = random_order[no_of_training_samples:]
input_training_data = vki_data['experimental_inputs'][train_inds]
input_testing_data = vki_data['experimental_inputs'][test_inds]
output_training_data = vki_data['experimental_outputs'][train_inds]
output_testing_data = vki_data['experimental_outputs'][test_inds]


and the following defines the polynomial using Parameters and Basis.

Ti = eq.Parameter(distribution='uniform', lower=1, upper=6, order=3,endpoints='both')
Ma = eq.Parameter(distribution='uniform', lower=0.7, upper=1.1, order=3,endpoints='both')
Re = eq.Parameter(distribution='uniform', lower=5e5, upper=2e6, order=3,endpoints='both')

my_params = [Ti, Ma, Re]
my_basis = eq.Basis('total-order', orders=[3,3,3])
print('Number of polynomial coefficients: %d' % my_basis.get_cardinality())


It is worth noting that we have more polynomial coefficients than number of training data points here. In a Bayesian polynomial approximation, one first specifies a prior probability distribution on the coefficients. Then, with the training data, one can obtain the likelihood function, which can be used to deduce the posterior distribution on the coefficients via Bayes’ rule:

\text{Posterior} \propto \text{Likelihood} \cdot \text{Prior}.

Using the posterior distribution on coefficients, predictions on test data can be obtained. In this example, we assume that the prior, likelihood and posterior distributions can be specified as Gaussian distributions, so they can be fully specified by providing the mean and covariance parameters (Note that this is not a necessary assumption, as is explored in our paper at [2203.03508] Prior-informed Uncertainty Modelling with Bayesian Polynomial Approximations). In EQ, the priors can be specified during the construction of a Polybayes instance:

basis_cardinality = my_basis.get_cardinality()
prior_coefficients_mean = np.zeros(basis_cardinality)
prior_coefficients_cov = np.eye(basis_cardinality)
data_std = 3.0

my_polybayes = eq.Polybayes(my_params, my_basis, prior_mean=prior_coefficients_mean, prior_cov=prior_coefficients_cov, sigma_data=data_std)


The parameter data_std refers to the uncertainty in the training data we observe expressed as a standard deviation. The posterior coefficients can then be computed by supplying our Polybayes object with the training data.

my_polybayes.compute_posterior_coefficients(input_training_data, output_training_data)


The posterior coefficient distribution parameters are now stored inside this object. We can then evaluate posterior predictions on some input data.

predicted_mean_train, predicted_cov_train = my_polybayes.get_posterior_fit(input_training_data)
predicted_std_train = np.sqrt(np.diag(predicted_cov_train))
predicted_mean_test, predicted_cov_test = my_polybayes.get_posterior_fit(input_testing_data)
predicted_std_test = np.sqrt(np.diag(predicted_cov_test))


The posterior distribution of the predicted outputs are Gaussian, with a certain mean vector and covariance matrix. Pointwise standard deviations can be evaluated by taking the square root of the diagonal of the covariance matrix, and can be used to construct credible intervals (uncertainty estimates of the predictions). Let’s now plot the data and evaluate the root mean squared prediction error per point.

plt.figure(figsize=(6, 5))
plt.errorbar(output_training_data, predicted_mean_train, yerr=predicted_std_train, color='limegreen',ecolor='limegreen', fmt='o', \
solid_capstyle='projecting', capsize=5, markersize=10, alpha=0.5, label='Training')

plt.errorbar(output_testing_data, predicted_mean_test, yerr=predicted_std_test, color='orangered', ecolor='orangered', fmt='o', \
solid_capstyle='projecting', capsize=5, markersize=10, alpha=0.5, label='Testing')
plt.plot(np.linspace(0, 130, 10), np.linspace(0, 130, 10), 'k--')
plt.legend()
plt.xlabel('Experimental data, $H_{int}$', fontsize=14)
plt.ylabel('Posterior polynomial fit, $H_{int}$', fontsize=14)
plt.tight_layout()
plt.show()

rms_err_train = np.linalg.norm(predicted_mean_train - output_training_data) / len(output_training_data)
rms_err_test = np.linalg.norm(predicted_mean_test - output_testing_data) / len(output_testing_data)
print('Training RMS error: %f' % rms_err_train)
print('Test RMS error: %f' % rms_err_test)


As we would expect, the error on unseen test data is higher than the training error. Can we improve this fit using extra knowledge of the task?

## Conditioning on output mean

One source of knowledge we can incorporate into the model is an estimate of the output mean. This need not come from data, but can be informed from expert knowledge about the physics of the problem. Provided with this estimate, the space of posterior distributions on the coefficients can be constrained, which improves the prediction accuracy. Let’s demonstrate this by incorporating the output mean estimated from the experimental data of the LS89. In code, this amounts to an extra argument supplied to the prediction method.

predicted_mean_train, predicted_cov_train = my_polybayes.get_posterior_fit(input_training_data, estimated_mean=np.mean(vki_data['experimental_outputs']))
predicted_std_train = np.sqrt(np.diag(predicted_cov_train))
predicted_mean_test, predicted_cov_test = my_polybayes.get_posterior_fit(input_testing_data, estimated_mean=np.mean(vki_data['experimental_outputs']))
predicted_std_test = np.sqrt(np.diag(predicted_cov_test))


Now, we evaluate the fit in a similar manner as above.

plt.figure(figsize=(6, 5))
plt.errorbar(output_training_data, predicted_mean_train, yerr=predicted_std_train, color='limegreen',ecolor='limegreen', fmt='o', \
solid_capstyle='projecting', capsize=5, markersize=10, alpha=0.5, label='Training')

plt.errorbar(output_testing_data, predicted_mean_test, yerr=predicted_std_test, color='orangered', ecolor='orangered', fmt='o', \
solid_capstyle='projecting', capsize=5, markersize=10, alpha=0.5, label='Testing')
plt.plot(np.linspace(0, 130, 10), np.linspace(0, 130, 10), 'k--')
plt.legend()
plt.xlabel('Experimental data, $H_{int}$', fontsize=14)
plt.ylabel('Posterior polynomial fit, $H_{int}$', fontsize=14)
plt.tight_layout()
plt.show()

rms_err_train = np.linalg.norm(predicted_mean_train - output_training_data) / len(output_training_data)
rms_err_test = np.linalg.norm(predicted_mean_test - output_testing_data) / len(output_testing_data)
print('Training RMS error: %f' % rms_err_train)
print('Test RMS error: %f' % rms_err_test)


It can be observed that the predictions are indeed improved. One caveat with the results here is that we are using information from the test data in our predictions (because we estimated the mean using the entire dataset). In practice, the information would come from other sources.

## Using a physics-informed prior

Another source of information can come from the simulation data. Let’s load in another dataset obtained from CFD simulations of the LS89 by solving the RANS equations (i.e. a relatively low fidelity CFD, which can be done with a low computational budget). We form a least squares polynomial model using this data over the same polynomial basis as before.

cfd_inputs = vki_data['CFD_inputs']

cfd_outputs = vki_data['CFD_outputs'].reshape(64, 1)

cfd_poly = eq.Poly(parameters=[s1,s2,s3], basis=my_basis, method='least-squares', \
sampling_args={'sample-points': cfd_inputs,
'sample-outputs': cfd_outputs})
cfd_poly.set_model()


Now, set the prior coefficients of the experimental polynomial model to the least squares coefficients of the CFD polynomial model, and solve for the posterior coefficients and predictions.

basis_cardinality = my_basis.get_cardinality()
prior_coefficients_mean = cfd_poly.get_coefficients().reshape(-1)
prior_coefficients_cov = np.eye(basis_cardinality)
data_std = 3.0

my_polybayes = eq.Polybayes(my_params, my_basis, prior_mean=prior_coefficients_mean, prior_cov=prior_coefficients_cov, sigma_data=data_std)
my_polybayes.compute_posterior_coefficients(input_training_data, output_training_data)

predicted_mean_train, predicted_cov_train = my_polybayes.get_posterior_fit(input_training_data)
predicted_std_train = np.sqrt(np.diag(predicted_cov_train))
predicted_mean_test, predicted_cov_test = my_polybayes.get_posterior_fit(input_testing_data)
predicted_std_test = np.sqrt(np.diag(predicted_cov_test))


Evaluating this fit:

plt.figure(figsize=(6, 5))
plt.errorbar(output_training_data, predicted_mean_train, yerr=predicted_std_train, color='limegreen',ecolor='limegreen', fmt='o', \
solid_capstyle='projecting', capsize=5, markersize=10, alpha=0.5, label='Training')

plt.errorbar(output_testing_data, predicted_mean_test, yerr=predicted_std_test, color='orangered', ecolor='orangered', fmt='o', \
solid_capstyle='projecting', capsize=5, markersize=10, alpha=0.5, label='Testing')
plt.plot(np.linspace(0, 130, 10), np.linspace(0, 130, 10), 'k--')
plt.legend()
plt.xlabel('Experimental data, $H_{int}$', fontsize=14)
plt.ylabel('Posterior polynomial fit, $H_{int}$', fontsize=14)
plt.tight_layout()
plt.show()

rms_err_train = np.linalg.norm(predicted_mean_train - output_training_data) / len(output_training_data)
rms_err_test = np.linalg.norm(predicted_mean_test - output_testing_data) / len(output_testing_data)
print('Training RMS error: %f' % rms_err_train)
print('Test RMS error: %f' % rms_err_test)


Again, an improvement in predictive accuracy can be observed.

## Conclusions

As demonstrated by the test case here, adopting a Bayesian approach to polynomial approximations does not only allow us to quantify additional sources of uncertainty. The Bayesian approach leads to methods that can incorporate multiple sources of information in a quantitative manner, which can improve predictive accuracy. In our paper, we show more ways of using Bayesian polynomials, including structural priors and coregional models that leverage knowledge about the sparsity pattern of the quantity of interest and data from similar correlated models respectively. Experiment with the new Polybayes module in this commit:

[1]: T. Arts, M. L. de Rouvroit, Aero-thermal performance of a two-dimensional highly loaded transonic turbine nozzle guide vane: A test case for inviscid and viscous flow computations, Journal of Turbomachinery 114 (1) (1992) 147-154.

[2]: S. Vitale, G. Gori, M. Pini, A. Guardone, T. D. Economon, F. Palacios, J. J. Alonso and P. Colonna. “Extension of the SU2 open source CFD code to the simulation of turbulent flows of fuids modelled with complex thermophysical laws,” AIAA 2015-2760. 22nd AIAA Computational Fluid Dynamics Conference. June 2015.

[3]: I. S. Torreguitart, T. Verstraete., L. Mueller. Optimization of the LS89 Axial Turbine Profile Using a CAD and Adjoint Based Approach †. Int. J. Turbomach. Propuls. Power 2018, 3, 20.

2 Likes