February 28, 2024
While many physics-based closure model forms have been posited for the SFS in LES, vast amounts of data available from DNS create opportunities to leverage data-driven modeling techniques. Albeit flexible, data-driven models still depend on the dataset and the functional form of the model chosen. Increased adoption of such models requires reliable uncertainty estimates both in the data-informed and out-of-distribution regimes. In this work, we employ BNNs to capture both epistemic and aleatoric uncertainties in a reacting flow model. In particular, we model the filtered progress variable scalar dissipation rate which plays a key role in the dynamics of turbulent premixed flames. We demonstrate that BNN models can provide unique insights about the structure of uncertainty of the data-driven closure models. We also propose a method for the incorporation of out-of-distribution information in a BNN, which can be used for out-of-distribution query detection. The efficacy of the model is demonstrated by a priori evaluation on a dataset consisting of a variety of flame conditions and fuels.
Large Eddy Simulation,Bayesian Neural Networks,Uncertainty Quantification,Progress variable dissipation rate
Numerical simulations of turbulent reacting flows are usually associated with high computational cost due to the large range of spatiotemporal scales that need to be resolved [1]. Concurrent with significant computational advances [2], several modeling approaches have since been devised to avoid resolving the smallest spatio-temporal scales. The two most prominent are RANS [3] and LES [4]. We will focus on LES in this work, but the approximations made in both approaches introduce the need for closure modeling [5], [6]. In LES, one only resolves the largest fluid flow scales. This is achieved by applying a low-pass filter to the original governing equations, which in turn requires to model the sub-filter scales. The modeling of this closure term is the primary focus of this paper. While several physics-based modeling strategies for closure terms have emerged, they tend to make strict assumptions about the unresolved scales [7], [8]. While physics models increasingly avoid these assumptions [9], data-driven strategies are becoming ever more popular, given their flexibility [10], data availability, and the development of SciML frameworks [11]–[16]. Recent SciML developments have focused on enforcing known physics constraints, such as with PINNs [17], [18]. This approach can be useful if a deconvolution procedure is used to compute closure terms [19], [20]. However, when directly approximating the closure term, no physics-law are typically available, and a purely data-driven approach is more appropriate.
There exists a considerable body of literature applying data-driven techniques for closure modeling broadly [21]. For turbulent combustion applications, multiple data-driven strategies have been shown to be accurate approximators of closure terms, such as artificial neural networks [10], [22]–[26], convolutional neural networks [27], [28], and neural ordinary differential equations [29]. While flexible, data-driven closure modeling strategies still depend on the model form and the training procedure, both of which are typically optimized via a hyper-parameter search. In addition, the choice of the dataset itself can induce significant variability on the closure model [30]. These choices can result in a significant uncertainty thereby motivating the need to equip these ML methods with reliable uncertainty estimates that can be propagated through the governing equations. From a practical perspective, providing objective uncertainty estimates of data-driven models is necessary before they can be confidently adopted or not by the engineering community. Uncertainty estimates are also needed for decision-making tasks, such as design optimization, where one relies on numerical simulations to optimize specific quantities of interest, e.g. ignition time [31], [32], maximum temperature [33], or mechanical structure [34], [35].
A variety of machine learning methods have been developed to model uncertainty [36]–[38]. Popular methods to quantify uncertainty in neural networks include the dropout method [39], [40] and its Bayesian interpretation [41]. While inexpensive and easy to implement, dropout generates probabilistic predictions at inference time but does not quantify the model parameter uncertainty [42]. In computational science and engineering, Gaussian processes have long been applied to uncertainty quantification tasks [43] owing to their flexibility and interpretability. The primary drawback of Gaussian processes is the \(\mathcal{O}(n^3)\) computation expense to train and \(\mathcal{O}(n^2)\) cost to evaluate due to the required matrix inversion and multiplication, where \(n\) is the data size. Despite efforts, scalability for large datasets remains a challenge [44]–[48]. New approaches are needed to handle massive datasets generated by DNS, especially in the context of rapid online evaluation within LES codes.
BNNs are an attractive method to estimate and predict modeling uncertainty due to their ability to ingest large amounts of data, relatively fast inference cost (compared to Gaussian processes), rigorous characterization of uncertainty, and expressivity. Recent advancements employing variational inference [40] have made the training of BNNs tractable for large models and amenable to large datasets [49], [50]. BNNs reformulate deterministic deep learning models as point-estimators and emulate the construction of an ensemble of neural nets by assigning a probability distribution to each network parameter [51]. Thus, they generate a predictive distribution by sampling the parameter distributions and collecting the resulting distribution of point estimates.
The central contribution of this paper is the development and application of a Bayesian neural network model to modeling the sub-filter progress variable dissipation rate of premixed turbulent flames [52]. In turbulent combustion modeling, uncertainty quantification tasks primarily target kinetic rates [53]–[56], operating and boundary conditions [33] and model coefficient in closure models [57], but do not address uncertainties of data-driven closure models. Recent applications of BNNs for combustion modeling have mostly focused on the prediction of macroscopic quantities such as ignition delay [58], fuel properties [59], or to accelerate data assimilation [60]. To our knowledge, this work is the first to explore the use of BNNs for quantifying both epistemic and aleatoric uncertainties in data-driven closure models [61]. While the framework developed and demonstrated in the context of progress variable dissipation rate, the techniques may be readily applied to other data-driven closure models.
Beyond the development of the model itself, we focus on two important applications. First, we utilize the learned epistemic uncertainty to derive physical insights about the potential weaknesses of the model. This knowledge can be used, for example, to inform what future data collection should be performed to improve the model’s predictive quality. Second, we provide a method for endowing a BNN model with meaningful predictive power in the out-of-distribution setting, i.e. when tasked with extrapolation in a regression setting. This allows one to strongly enforce a prior or detect scenarios where the data-driven model may fail and resort to another method or heuristic.
Like other data-driven closure models, BNNs can be introduced into existing reacting LES code in place of physics-based closure models, allowing for closure modeling prediction with uncertainty quantification. The integration of the BNN within a reacting flow solver will be the object of future work. The remainder of this paper is organized as follows. In Section 2, we formulate the target problem, including sources of uncertainty, dataset generation, and Bayesian neural network modeling. Numerical experiments and their results assessing the model performance are presented in Section 3. We discuss the utility of the model for downstream applications and its limitations in Section 4. Concluding remarks are given in Section 5.
In this section we first review the nature of aleatoric and epistemic uncertainty in Section 2.1. Next, we discuss generation of the training data in Section 2.2. We then provide an overview of Bayesian neural networks in Section 2.3 as well as practical considerations (Section 2.3.1 and Section 2.3.2). We conclude by providing the architecture for our experiments in Section 2.4.
While there are multiple sources of uncertainty arising from model inputs, numerical errors, or experiments, in this work we primarily consider the uncertainties driven by the data that the model trains on. In the machine learning context, the use of deep neural networks is typically justified by a universal approximation theorem [62] valid in the asymptotic limit of data availability. Thus, we concern ourselves with uncertainty in the learned parameters and resultant functional representation due to the empirical data distribution used for training.
Broadly speaking, uncertainties can be categorized as either aleatoric or epistemic [61], [63], [64]. As an illustrative example, consider the example shown in Fig. 1. The underlying generating function for the data is \[y = x ^3 + 0.1 (1.5+x)\varepsilon,\] where \(\varepsilon\sim\mathcal{N}(0,\sigma^2)\). Here, the noise in the data increases as the variable \(x\) increases, however only 50 data points are retained in the region on the left, compared to 500 data points in the region on the right. The left region has low aleatoric uncertainty due to the lack of noise but has high epistemic uncertainty due to the lack of data. On the other hand, any model on the right would be well informed by the quantity of the data and so would have low epistemic uncertainty, but would have a relatively high level of aleatoric uncertainty due to high noise level.
In the context of data-driven closure modeling, epistemic uncertainties arise from a lack of knowledge when too few training data are available. An ideal LES model in the sense of Langford et al. [65] obtained via a data-driven method, could be subject to epistemic uncertainty when there does not exist sufficient data to confidently determine the optimal estimator of the closure term. Epistemic uncertainty also presents itself in the extrapolation setting, which may be viewed as an extreme case of lacking data availability.
Aleatoric uncertainty is also known as statistical or stochastic uncertainty and is inherent to a problem or experiment. It is irreducible in the sense that it cannot be reduced by additional data collection. A clear example is the measurement tolerance of a sensor: additional data collection would not inform one to a greater precision than the instrument rating. In the data-driven closure modeling context, aleatoric uncertainty primarily presents itself as coarse–graining or filtering uncertainty. The act of filtering destroys information, as there are potentially many DNS realizations that, when filtered, would present identical filtered measurements. Inferring a mapping from filtered to unfiltered data is therefore subject to an irreducible uncertainty [20], [65]–[67]. Likewise, the selection of input features can also be seen as a coarse-graining step.
From a practical perspective in the context of closure modeling, a quantified epistemic uncertainty allows one to decide where to collect additional data samples before retraining a model. The aleatoric uncertainty may help reassess the choice of input features. Quantifying epistemic and aleatoric uncertainties is the first step to propagate them through a fluid simulation, and incorporate their effect into decision-making tasks.
The dataset used in this work is obtained from several DNS of turbulent premixed flames [68], [69]. The reader is referred to Yellapantula et al. [52] for a more detailed description of the data generation and preparation processes. Hereafter, only key details of the dataset are repeated. The DNS data was Favre filtered using Gaussian filter kernels with different filter sizes to prepare input features for a priori model evaluation that would match the LES features available in the a posteriori setting. The objective is to learn an optimal mapping between the filtered quantities and the quantity of interest, here the unresolved contribution to the filtered progress variable dissipation rate, \(\chi_{\rm SFS}\) defined as \[\chi_{\rm SFS} = \widetilde{\chi} - \chi_{\widetilde{C}},\] where \(\chi_{\widetilde{C}}\) is the resolved progress variable dissipation rate defined as \[\chi_{\widetilde{C}} = 2\widetilde{D}_C |\nabla \widetilde{C} |^2,\] and \(\widetilde{\chi}\) is the total progress variable dissipation rate defined as \[\widetilde{\chi} = \widetilde{2 D_C |\nabla C |^2},\] where \(D_C\) is the progress variable diffusivity, \(C\) is the progress variable, and \(\widetilde{.}\) denotes the Favre filtering operation.
The relevance of approximating closure models for the subfilter progress variable dissipation rate has been extensively documented elsewhere and is not repeated here [52], [70]. The input features are the same as those considered in [52] and are summarized in Table 1. The input features were selected based on a feature importance metric [71] and only contain filtered quantities as the model is envisioned to be deployed online within a larger LES model. Unlike [52], the output layer only predicts the unresolved contribution of the progress variable dissipation rate \(\chi_{\rm SFS}\). Other derived quantities include the principle rates of filtered strain rate. The principal rates of strain, \(\alpha > \beta > \gamma\) are computed using the strain rate tensor constructed from the filtered velocity field. Additionally, the alignments between the eigenvectors and local gradient of the filtered progress variable were added to the list of input features [52]. The input space \(\mathbf{x}\in\mathbb{R}^{d_i}\) is mapped to an output \(\mathbf{y}\in\mathbb{R}^{d_o}\), where \(d_i=10\) is the dimensionality of the input space and \(d_o\) is the dimensionality of the output space, in this case \(d_o=1\). The effect of the number of input parameters on the estimated uncertainties is further discussed in 8.
| Input Feature | Description |
|---|---|
| \(\widetilde{C}\) | Filtered progress variable |
| \(\widetilde{C^{\prime\prime 2}}\) | Filtered progress variable variance |
| \(2\widetilde{D}_C |\nabla \widetilde{C} |^2\) | Resolved progress variable dissipation rate |
| \(\widetilde{D}_{C}\) | Filtered progress variable diffusivity |
| \(\alpha\), \(\beta\), \(\gamma\) | First, second, and third principal rate of strain |
| \(e_{\alpha,\beta,\gamma} \cdot \nabla \widetilde{C} / |\nabla \widetilde{C} |\) | Alignment of local progress variable gradient with principal eigenvectors |
After the selection of the input features, a near-uniform in phase-space sampling [30] was performed by clustering the data with 40 clusters in the input feature space and uniformly selecting data from each cluster. This process allows for appropriately capturing the reaction zone of all the DNS cases in the dataset with a reduced number of data points. The final dataset contains \(7.88\times 10^6\) data points for training and \(2.63\times 10^6\) data points for testing, which is about \(10^3\) times less than the total number of DNS data points initially available.
Classical neural networks are only capable of deterministically generating point estimates and thus are not well suited for the task of assessing uncertainty. On the other hand, Bayesian neural networks model epistemic uncertainty estimation by placing a parametric distribution over the neural network parameters [50], [72]–[74]. By sampling the parameter distributions, one selects an ensemble of DNNs. Note that sampling multiple candidate neural networks to estimate uncertainty is consistent with peer-modeling approaches for uncertainty propagation [75]. In particular, we consider the parametric distribution for the weights to be a Gaussian [50]. The mean and variance of the weights become the parameters to be learned in the BNN representation. To delineate between the sampled parameters defining a DNN and the distribution placed over the weights in the BNN representation, we will set the notation that “weights” \(\mathbf{w}\) are sampled from a (parametric) distribution placed upon the weights \(q(\mathbf{w}|\theta)\), where \(\theta\) are the parameters of the BNN. Under this notation, \(\mathbf{w}\sim q(\mathbf{w}|\theta)\).
In the variational inference formulation of BNN training [76], one seeks to minimize the KL divergence between the weight distribution and the true Bayesian posterior conditioned on the dataset \(\boldsymbol{D}\), \[\theta^* = \operatorname*{arg\,min}_\theta \operatorname{KL}\left[ q(\mathbf{w}|\theta) || p(\mathbf{w}|\boldsymbol{D}) \right].\] Expansion and rearrangement of this formulation yields the ELBO objective function [50], \[\label{eq:elbo} \operatorname*{arg\,min}_\theta \underbrace{\operatorname{KL}\left[ q(\mathbf{w}|\theta) || p(\mathbf{w}) \right]}_{\text{prior-informed}} - \underbrace{\mathbb{E}_{q(\mathbf{w}|\theta)}\left[ \log p(\boldsymbol{D}|\mathbf{w}) \right]}_{\text{data-informed}}.\tag{1}\] The first term in the ELBO is the KL divergence between the learned distribution on the BNN parameters and the prior. The second term in the ELBO is a data misfit term given by the expected negative log-likelihood of the data over the distribution of plausible models encoded by the distribution of the weights. The consequences of this optimization formulation are explored in the following sections. Hereafter, this probabilistic representation of the model predictions is discussed for regression tasks.
Although the primary interest of modelers that deploy closure models in a reacting flow solver is the epistemic uncertainty, one may wish to evaluate the aleatoric uncertainty as well (see Sec. 2). One method for doing so is to further parameterize each output dimension so as to capture heteroskedastic uncertainty in the data [77]–[79]. We employ the modeling choice from Kendall and Gal [77] of a Gaussian form for the output random variable.
In this work, a dense feed-forward neural network parameterization of the closure term is adopted. That is, layers comprise a weight matrix \(W\) acting upon the previous layers’ outputs, a bias term \(b\), and a non-linear activation \(\sigma\), \[x_{\ell+1} = \sigma_\ell \left(W_\ell x_\ell + b_\ell\right).\] Here \(x_\ell\) denotes the output from the last layer, and \(\mathbf{x}\in\mathbb{R}^{d_i}\) again represents an input data point. The final layer is mapped to \((\mathbf{y}_{\boldsymbol{\mu}}, \mathbf{y}_{\boldsymbol{\sigma}})^T\in\mathbb{R}^{2d_o}\) in the case of a Kendall and Gal [77] style architecture. Here \(\mathbf{y}_{\boldsymbol{\mu}}, \mathbf{y}_{\boldsymbol{\sigma}}\) parameterize the distribution of the output variable as \(\mathbf{y}\in\mathbb{R}^{d_o}\sim\mathcal{N}\left(\mathbf{y}_{\boldsymbol{\mu}}, \operatorname{diag}(\mathbf{y}_{\boldsymbol{\sigma}})\right)\). This parameterization allow to model epistemic uncertainty (through the variance of \(\mathbf{y}_{\boldsymbol{\mu}}\)) the aleatoric uncertainty (through the variance of \(\mathbf{y}_{\boldsymbol{\sigma}}\)). The difference in the epistemic only architecture and the epistemic and aleatoric architectures is shown graphically for \(d_o=1\) in Fig. 2 (a) and Fig. 2 (b), respectively. In this work, we adopt the notation that the output \(\mathbf{y}\) is the result of sampling the output of a BNN parameterized by \(\theta\) at a datum \(\mathbf{x}\) after appropriate sampling of the weights \(\mathbf{w}\) and the resultant random variable.
In this framework, one can use the law of total variance [80] to decompose the variance (uncertainty) in the predictions into its epistemic and aleatoric components [79], [81], \[\label{eqn:uncertainty-decomposition} \operatorname{Var}(\mathbf{y}) = \underbrace{\mathbb{E}_{q(\mathbf{w}|\theta)} \left[\operatorname{Var}(\mathbf{y} \vert \mathbf{x})\right]}_{\text{aleatoric}} + \underbrace{\operatorname{Var}_{q(\mathbf{w}|\theta)}\left(\mathbb{E}\left[\mathbf{y}\vert \mathbf{x}\right]\right)}_{\text{epistemic}}.\tag{2}\] The predictive variance \(\operatorname{Var}(\mathbf{y})\) is decomposed into an aleatoric component \(\mathbb{E}_{q(\mathbf{w}|\theta)} \left[\operatorname{Var}(\mathbf{y} \vert \mathbf{x})\right]\), the mean variability in the estimates, and an epistemic component \(\operatorname{Var}_{q(\mathbf{w}|\theta)}(\mathbb{E}[\mathbf{y}\vert \mathbf{x}])\), the variability in the model’s mean predictions. For a given \(\mathbf{w}\sim q(\mathbf{w}|\theta)\), the predictive mean is given by the first network output, \(\mathbf{y}_{\boldsymbol{\mu}}\in\mathbb{R}^{d_o}\). Detailed algorithms for computation of the uncertainties are provided in 9. For a concrete example, the magnitudes of the uncertainties for the one-dimensional example (Fig. 1) are shown in Fig. 3. It can be observed that the regions with high and low epistemic uncertainty are appropriately characterized in Fig. 3 (a) and that the predictive envelope captures the entire data distribution in Fig. 3 (b).


Figure 2: (a) A BNN in the style of [50] that captures epistemic uncertainty only and (b) a BNN in the style of [77] that captures both epistemic and aleatoric uncertainty..


Figure 3: Comparison of two Bayesian neural networks trained on the dataset presented in Fig. 1, (a) models only the epistemic (model-form) uncertainty, while (b) models both epistemic and aleatoric uncertainty..
Bayesian inference provides a way to integrate empirical data with prior information about phenomenon being modeled. In practice, the prior drives the first term in the ELBO objective function (Eq. 1 ), and thus the selection of the prior has a direct impact on the quality of the resultant model [82]. In extreme cases, prior misspecification can only be overcome in the asymptotic limit of data, with the resultant model strictly adhering to the prior belief.
There are two ways of viewing BNN priors: the weight-space view and the functional-space view. Since BNNs operate on a parameterized distribution of the weights, \(q(\mathbf{w}\vert\theta)\), a prior is inherently defined on the BNN parameters. However, the effect of a single weight on the output is difficult to gauge due to interactions with other weights and nonlinear layers. In other words, the lack of interpretability of the neural network parameter prevents one from meaningfully assigning a prior to the weights. In this case, using a non-informative prior such as an isotropic Gaussian could be envisioned, however, this has been shown to be sub-optimal [83]. Often one does have intuition about the outputs and many methods have been derived to “tune” priors empirically. These include warm-start methods [84], initial fitting of the model to Gaussian process realizations [85], and initialization from a trained deterministic model [86]. One may be tempted to instead ignore the contribution of the prior altogether, however, it has been shown that direct minimization of the negative log-likelihood is inappropriate for BNNs [87]. In this work, we employ a “trainable” prior that updates the prior mean every epoch while while keeping the variance fixed. This choice changes the interpretation of the KL term to that of a data-driven regularization and has connections to the KL annealing method in [88].
The aforementioned strategies provide a way to quantify epistemic and aleatoric uncertainty where data is available. Deploying data-driven models in place of physics-based models requires that the phase space spanned by the training data appropriately encompasses the phase space over which the model is queried. In general, this requirement is difficult to guarantee, which has traditionally led to adaptive methods [89] or data-driven methods complemented with synthetic data [90]. In the absence of training data, one would like for extrapolatory predictions to closely match those coming from physics models. Solely relying on the extrapolatory nature of neural network models is unrealistic, especially since it is still not fully understood how neural networks extrapolate [91], [92]. In the present case, if training data is not available, one may prefer that the sub-filter progress variable dissipation rate should fall back to the predictions of a physics-based model. One way to ensure this behavior would be to detect that a model is called outside of the range of the training data, in which case a separate closure model may be used. However, detecting whether or not an individual input is OOD is a challenging problem. An alternative approach is to include data obtained from the physics-based model in the training dataset. To ensure that the synthetic data does not “pollute” the original dataset, the data is generated only in the OOD region feature space that is separate from the true data distribution, \(\boldsymbol{D}\). The synthetic data can be labeled with the prediction of the physics model and an additive noise in case the physics model is subject to uncertainty.
A demonstration of these principles is presented in Fig. 4. Here, the underlying true data is the same as for Fig. 1. On the other hand, the functional form of the synthetic data is chosen to be \[\hat{y} = x^2 + \eta,\] where \(\eta\sim\mathcal{N}(0,0.25^2)\). In Fig. 4 (a), it is shown that training on the true data pairs \((\mathbf{x}_i,y_i)_i\) after an initial training phase to the “prior” data \((\mathbf{x}_i,\hat{y}_i)_i\) demonstrates an inability to reproduce the prior. However, when the two datasets are concatenated as in Fig. 4 (b), the model is able to more faithfully reproduce the prior. Clearly, there will be a trade-off between the amount of synthetic, low-fidelity, or “prior” data included as well as its distance from the true data pairs necessary to reasonably enforce the prior without spoiling the representation learned from the actual data. These trade-offs are explored in Section 3.4.


Figure 4: Demonstration of “catastrophic forgetting” when a BNN model is (a) first trained to a prior-enforcing dataset and then the true dataset (the warm-start method) and (b) trained on combined dataset..
To generate OOD synthetic data, two methods are explored: a SBO approach [93] and a NF approach [94]. In the SBO approach, the OOD data is generated based on a distance metric computed with respect to all the available data. The method is particularly useful for small datasets and allows controlling the distance between the OOD region spanned and the training data. However, the distance computation can become computationally prohibitive when used with large datasets. In the NF approach, synthetic data is uniformly generated over a larger domain than the original training data. Here the data spans a hypercube whose edge size is 40% larger than the spanned amplitude of each feature. The synthetic data that overlaps with training data is subsequently discarded. To evaluate whether a synthetic data point overlaps with the original training data, the likelihood of the training dataset is estimated via a NF. Any synthetic data whose likelihood is higher than a given threshold is discarded. The method is summarized in Algorithm 5.
Unlike SBO, the NF approach does not allow to control the distance between the synthetic OOD dataset and the original dataset, however, it allows fast processing of very large datasets [30]. Here, the likelihood threshold \(p_0\) is the minimal likelihood predicted over the true training set, which ensures that the OOD region spanned by the synthetic data is strictly disjoint from the original dataset \(\boldsymbol{D}\).
To account for extrapolation uncertainty, the OOD dataset generated can be labeled with a label distribution that describes how much uncertainty can be expected if the model is queried outside the original dataset \(\boldsymbol{D}\). To simplify the discussion on extrapolation uncertainty, the OOD data is labeled as
\[\label{eq:oodlabel} \chi_{\rm SFS,OOD} \sim \mathcal{N}(\mu_{\rm OOD}, \sigma_{\rm OOD}),\tag{3}\]
where \(\mathcal{N}\) denotes a normal distribution with mean \(\mu_{\rm OOD}\) being arbitrarily set to 0. The uncertainty in the OOD is set to \(\sigma_{\rm OOD} = \operatorname{Var}_{\boldsymbol{D}}(\chi_{\rm SFS})^{1/2}\) which simplifies the implementation and the evaluation (Sec. 3.4). More complex strategies for generating OOD labels could be formulated, such as labeling the synthetic OOD data using a LRM. Here, an arbitrary marker is used which can primarily serve to detect OOD queries (see Sec. 4.2). In Sec. 3.4, the performance of both synthetic data generation strategies are evaluated with respect to the number of data generated and distance to the training set.
In this study, a feed-forward fully connected architecture is considered, similar to that in [52]. To account for variability in
performance due to model architecture and form, a hyper-parameter tuning study was performed. Specifically, a grid search over all 336 combinations of the following parameters: hidden dimensions \(n_h\in\{5,10,15,20\}\),
number of hidden layers \(n_l\in\{\)2, 3, 4\(\}\), batch size \(M\in\{\)256,512,1024,2048,4096,8192,16384\(\}\), and
learning rate \(\eta\in\{\) \(1\mathrm{e}{-03}\), \(1\mathrm{e}{-04}\), \(1\mathrm{e}{-05}\), \(1\mathrm{e}{-06}\) \(\}\) with each model trained for 1500 epochs. Hyperparameter realizations were generated with scikit-learn [95] and were distributed using the Texas Advanced Computing Center’s Launcher utility [96]. The most successful model architecture is summarized in Table 2. Training took approximately two hours per model on the Eagle computing
system at the National Renewable Energy Laboratory. Inference over the roughly two million data points in the validation dataset took on the order of one second. The model was implemented with TensorFlow Probability [97], [98] and the code is available in a companion repository
https://github.com/NREL/mluq-prop.
| Symbol | Name | Value |
|---|---|---|
| \(n_h\) | Number of units in hidden layer | 20 |
| \(n_l\) | Number of hidden layers | 4 |
| \(M\) | Batch size | 2048 |
| \(\eta\) | Learning rate | \(1\mathrm{e}{-04}\) |
| \(\sigma\) | Nonlinear activation function | Sigmoid |
In this section, the performance of the BNN model is demonstrated through a priori testing. As noted in Sec. 2.2, we use the dataset from Ref. [52] where input parameters are generated from the filtered DNS data. Model predictions using these input features are compared against filtered DNS quantities. We assess the BNN model quality on the validation dataset in Section 3.1. The predicted aleatoric and epistemic uncertainties and their utility for dataset augmentation are discussed in Section 3.2. Flame contours are modeled and studied in Section 3.3. Finally, the extrapolatory behavior of the model is evaluated in Section 3.4.
The predictive distribution of a BNN model trained only on the dataset presented in Sec. 2.2, with no additional synthetic data is shown in Fig. 6 and Fig. 7. The mean of the model predictions on the unseen validation data are compared to the true progress variable progress variable dissipation rate as computed by the DNS in Fig. 6. An excellent agreement can be observed between the model and data, especially in regions with copious amounts of data. The vast majority of the data clustered along the one-to-one line between the mean prediction values and the true values. At high dissipation rate, the model performance is degraded primarily because the choice of input features does not uniquely map to \(\chi_{\rm SFS}\). This is confirmed by Fig. 7 (b) which shows elevated aleatoric uncertainty especially when \(\chi_{\rm SFS}/\chi_{\rm lam}\) (where \(\chi_{\rm lam}\) is the maximum of the progress variable dissipation rate in the freely propagating laminar flame) exceeds 2. The trained BNN model is compared with a DNN model trained with the same architecture as the BNN, as well as the physics-based LRM [7]. The mean-squared error on the test dataset for these models is reported in Table 3. In the case of the BNN, the error is computed using the average predictions of the model. We observe that both machine learning approaches outperform the LRM (98.5% error reduction for the BNN and 99.4% error reduction for the DNN). The DNN slightly outperforms the BNN which is unsurprising given that it is trained to only minimize the mean squared error, while the BNN also accounts for the KL divergence term (Eq. 1 ). Overall, the BNN allows for obtaining uncertainty estimates while minimally compromising on accuracy.


Figure 7: (a) Scatter plot showing mean predictions for individual data points shaded by the BNN model prediction of epistemic uncertainty of \(\chi_{\rm SFS}/\chi_{\rm lam}\) and (b) shaded by the BNN model prediction of aleatoric uncertainty of \(\chi_{\rm SFS}/\chi_{\rm lam}\)..
| Model | Mean Squared Error |
|---|---|
| Linear Relaxation Model [7] | \(3.98 \times 10^{-1}\) |
| Deterministic Neural Network [52] | \(2.52 \times 10^{-3}\) |
| Bayesian Neural Network | \(5.69 \times 10^{-3}\) |
For further evaluation of the model performance, the conditional structure of the progress variable dissipation rate with respect to the input features is inspected. In particular, the performance with respect to three important physical parameters is investigated: the filtered progress variable \(\widetilde{C}\), the sub-filter variance of the progress variable \(\widetilde{C^{\prime\prime 2}}\), and the resolved progress variable diffusivity \(\widetilde{D}_C\). Conditional means were computed with \(250\) bins discretizing the input feature domain. Credible intervals are constructed by computing the BNN model output across the dataset for 250 realizations of \(\mathbf{w}\sim q(\mathbf{w}|\theta)\) and retaining the \((1-\alpha) \times 100\%\) interquantile range of predicted trajectories. In this case, \(\alpha=0.1\) corresponding to the \(90^{\text{th}}\) percentile. For further discussion of the nuances of predictive, credible, and confidence intervals the reader is referred to [63].






Figure 8: Comparison of the model predictive mean and predictive envelope to the validation data mean ((a), (c), (e)) and the distribution of the data conditioned on one dimension of phase space ((b), (d), (f)) for: the filtered progress variable ((a), (b)), the subfilter progress variable variance ((c), (d)), and the filtered progress variable diffusivity ((e), (f))..
The conditional mean profiles displayed in Fig. 8 demonstrate excellent agreement between the data and the model predictive mean. Moreover, any discrepancies, such as those in Fig. 8 (c) are clearly contained by the predictive envelope. The total predictive uncertainty is dominated by the aleatoric uncertainty as suggested by Fig. 7.
To estimate and visualize the distribution of epistemic and aleatoric uncertainties, 250 realizations of the model \(\mathbf{w}\sim q(\mathbf{w}|\theta)\) are generated and the sample variance estimate of the epistemic uncertainty \(\operatorname{Var}_{q(\mathbf{w}|\theta)}(\mathbb{E}[\mathbf{y}\vert \mathbf{x}])\), is computed. The conditional mean of this quantity is computed with respect to the input features and is presented in Fig. 9. The predicted aleatoric uncertainty is computed similarly as \(\mathbb{E}_{q_{c}(\mathbf{w}|\theta)}\left(\operatorname{Var}\left[\mathbf{y}\vert \mathbf{x}\right]\right)\) and is also shown in Fig. 9.
First, the epistemic uncertainty is about two orders of magnitude lower than the aleatoric uncertainty, which echos the finding that aleatoric uncertainty dominates the predictive uncertainty (Fig. 7). The epistemic uncertainty follows the same trends as the ones observed for the aleatoric uncertainty for the conditioning done with respect to the progress variable (Fig. 9 (a)) and the filtered diffusivity (Fig. 9 (c)). However, as the progress variable variance increases (Fig. 9 (b)), the distribution of epistemic and aleatoric uncertainty deviate, which is further discussed in Sec. 4.



Figure 9: Conditional average of the aleatoric and epistemic uncertainty (\(\times\) 100) with respect to (a) filtered progress variable, (b) subfilter progress variable variance, and (c) filtered diffusivity..
The predicted total filtered progress variable dissipation rates, \(\widetilde{\chi} = \chi_{\widetilde{C}} + \chi_{\rm SFS}\) and associated uncertainties (epistemic and aleatoric) are visualized for two different flames: an n-heptane flame with Karlovitz number \(7\) [69] which was not included in the dataset (referred to as \(B_{\rm Le}\)), and a n-heptane flame with Karlovitz number \(256\) which was included in the dataset (referred to as \(D_{\rm Le}\)). Since the training dataset can only be constructed with a limited set of operating conditions, testing the model outside of those operating conditions evaluates whether the model is only applicable over the conditions included in the training dataset. The prediction results and the ground truth data are shown for filter width 4 (Fig. 10) and 16 (Fig. 11).
Figure 10 shows that the BNN is in agreement with the filtered DNS data, including for the \(B_{Le}\) flame which was not included in the dataset. The BNN also appropriately captures the effect of the Karlovitz number in that the total filtered progress variable dissipation rate is about one order of magnitude higher for the \(D_{\rm Le}\) than for the \(B_{\rm Le}\) case. The effect of the filter width is also appropriately captured as shown in Fig. 11. At higher filter width, the total filtered progress variable dissipation rate peaks at lower values for both flames.


Figure 10: Top: \(\widetilde{\chi} / \chi_{\rm lam}\) at a midplane for the ground truth (left) and mean BNN prediction (right). Bottom: Epistemic (left) and aleatoric (right) uncertainties associated with the prediction of \(\chi_{\rm SFS} / \chi_{\rm lam}\). The bottom of the contours denotes burnt gas and the top of the contours denotes unburnt gas. Results are shown for a filterwidth of 4.. b — \(B_{\rm Le}\) case., d — \(D_{\rm Le}\) case.
In line with the results found in Sec. 3.2, the epistemic uncertainty is consistently about 2 orders of magnitude lower than the aleatoric uncertainty. The spatial distribution of epistemic and aleatoric uncertainties are similar as also noted in Sec. 3.2. However, the epistemic uncertainty is more localized, which was also noted earlier (Fig. 9 (b)). At larger filter widths, both the epistemic and the aleatoric uncertainty increase, which is the result of a higher amount of information being destroyed by the filtering operation. In the case of the \(D_{\rm Le}\) flame with a filter width of 16, \(\widetilde{\chi}\) prediction is not as smooth as the ground truth data, especially where the epistemic uncertainty is high. We attribute this to the dataset being insufficiently rich.


Figure 11: Same as Fig. 10 but with a filterwidth of 16.. b — \(B_{\rm Le}\) case., d — \(D_{\rm Le}\) case.
In this section, the performance of the model outside of the original dataset distribution \(\boldsymbol{D}\) is evaluated. The model’s learned representation in the OOD regime is mainly driven by the choice of synthetic dataset (Sec. 2.3.2). The role of the synthetic dataset is to enforce an OOD behavior without spoiling the performances in distribution. The objective of this section is to understand how should the synthetic dataset be generated. In particular, the choice of SBO or NF is evaluated, the balance between the size of the synthetic data and the original data is discussed, and the effect of the distance between the synthetic data and the original dataset is described where applicable. A baseline model is first trained with the architecture prescribed in Tab. 2 only using the filtered DNS dataset described in Sec. 2.2. The in-distribution performance of the baseline model was previously presented in Sec. 3.1 and Sec. 3.2. The variational posterior of this reference model is denoted by \(q_r(\mathbf{w}\vert\theta)\). Comparison models with identical architectures to the reference model are then trained with datasets augmented by various amounts of synthetic data generated as described in Sec. 2.3.2. The variational posterior of these candidate models is denoted \(q_c(\mathbf{w}\vert\theta)\).
The normalized \(L_2\) error for the first moment outside the data distribution is computed as
\[\label{eq:normErrmom1} \| \mu_{\rm OOD} - \mathbb{E}_{q_{c}(\mathbf{w}|\theta)}\left(\mathbb{E}\left[\mathbf{y}\vert \mathbf{x}\right]\right)\|_2/\sqrt{N},\tag{4}\]
where the \(L_2\) norm is computed over the set of synthetic data \(\boldsymbol{D_{\rm OOD}}\) and \(\mathbf{x}\) \(N\) is the number of synthetic data points in \(\boldsymbol{D_{\rm OOD}}\), and the normalization \(1/\sqrt{N}\) accounts for the different synthetic dataset sizes. Likewise, the normalized \(L_2\) error for the second moment is computed as
\[\label{eq:normErrmom2} \|\sigma_{\rm OOD}^2 - \mathbb{E}_{q_{c}(\mathbf{w}|\theta)}\left(\operatorname{Var}\left[\mathbf{y}\vert \mathbf{x}\right]\right)\|_2/\sqrt{N}.\tag{5}\]
Convergence plots for these metrics with respect to the amount of synthetic data are presented in Fig. 12. Unsurprisingly, the higher the amount of synthetic data, the lower the error in the space spanned by the synthetic dataset for both the first and the second moment. It can also be observed that for the same amount of data, generating the OOD data with the NF method consistently outperforms all the SBO methods. Furthermore, the distance between the synthetic and the in-distribution dataset that can be controlled with the \(d\) parameter does not bridge the gap with the NF method. This observation suggests that the difference in performance between NF and SBO may be instead due to synthetic dataset distribution. Unlike NF, the SBO method does not generate uniformly distributed data in the OOD region. In addition, the bounds of the OOD region are explicitly specified with the NF method but not in the SBO. The gap in performance between the NF and the SBO methods could be due to undersampled parts of the OOD region. Finally, although the error metric converges exponentially fast with respect to the amount of synthetic data, the convergence rates appear to decrease when the amount of synthetic data approaches the amount of in-distribution data.


Figure 12: Convergence of the normalized \(L_2\) error between (a) the model predictive means and the synthetic data in the extrapolation regime and (b) the aleatoric uncertainty estimate and the true aleatoric uncertainty associated with the synthetic data generation..
With the addition of increasing amounts of synthetic data, the model is tasked with representing a larger space and it is natural to expect that the performance should degrade on the original dataset. Figure 13 shows the effect of the amount of synthetic data on the relative \(L_2\) error on the original dataset. For a given model \(M\), the predictive mean is noted as \(P_M = \mathbb{E}_{q_{M, c}(\mathbf{w}|\theta)}\left[\mathbf{y}\vert \mathbf{x}\right]\). The error metric considered is a relative error between the mean predictions of a model trained without synthetic data \(P_{\rm ref}\) and a model trained with different amounts of synthetic data \(P_{\rm OOD}\), as
\[\label{eq:normErrmom1-inD} \| P_{\rm OOD} - P_{\rm ref} \|_2 / \| P_{\rm ref} \|_2.\tag{6}\]
The relative errors in the predictive mean are shown in Fig. 13. It is observed that the quality of the predictive mean is affected by the presence of synthetic data. However, the relative error increases at a slower rate than the error reduction in the OOD regime. The relative error increases at a faster rate once the number of synthetic data points approaches the size of the in-distribution dataset, especially for the SBO data generation. Once more, it appears that the NF tends to generate higher quality data.
In this section, the uncertainties estimated with the BNN are discussed to describe how they can be used to improve data-driven modeling of closure terms, either to improve the definition of input features, the data collection effort or detect out-of-distribution queries. We discuss the role of the uncertainty decomposition in Section 4.1. We discuss how the extrapolation behavior of the BNN can be used for out-of-distribution detection in Section 4.2. Finally, an outlook for a posteriori uncertainty propagation through high-fidelity forward simulations is discussed in Section 4.3.
For both the BNN (Fig. 7) and the DNN [52], higher errors were primarily observed for large values of \(\chi_{\rm SFS}/\chi_{\rm lam}\). There, large aleatoric uncertainties are consistently observed (Fig. 7 (b)) which suggests that the input features are insufficiently fine-grained for large filterwidth ratio, or for intense turbulence. In the future, it could be useful to design input features that specifically minimize the aleatoric uncertainty in this region.
The distribution of aleatoric uncertainty with respect to the progress variable (Fig. 8 (b)) suggests that most of the aleatoric uncertainty is located in the burning region of the flame, while reduced uncertainty can be observed either in the fully burnt or unburnt regions. The distribution of aleatoric uncertainty is also skewed towards the high end of the progress variable distribution which may suggest that small-scale variability is mostly controlled by the end of the burning process, as also noted in Ref. [99].
The subfilter variance of the progress variable \(\widetilde{C^{\prime\prime 2}}\) may be understood as a marker for the filter width of the LES. As the filter width increases, the aleatoric uncertainty rapidly increases which suggests the presence of small-scale variability that cannot be captured with the filtered features (Fig. 8 (d)). The aleatoric uncertainty also appears to plateau at high progress variable variance \(\widetilde{C^{\prime\prime 2}}\). Since \(\widetilde{C^{\prime\prime 2}}\) characterizes small scale variability, it is increasingly difficult for the filtered input features to capture \(\chi_{\rm SFS}\). This would result in monotonically increasing aleatoric uncertainty with respect to \(\widetilde{C^{\prime\prime 2}}\). However, the plateauing behavior instead suggests that the subgrid scale variability of \(\chi_{\rm SFS}\) that can be captured by the filtered input features saturates at medium values of \(\widetilde{C^{\prime\prime 2}}\).
Overall, the consistent distribution of epistemic and the aleatoric uncertainties showns in Fig. 9 could be a consequence of the near-uniform phase-space sampling done as a pre-processing step (Sec. 2.2). Aleatoric uncertainty is not uniformly distributed over features space, meaning that not all feature regions “need" the same amount of data. This suggests that the downsampling procedure should attempt to over-represent high-aleatoric uncertainty regimes, instead of uniformly distributing the dataset.
When conditioned on the progress variable sub-filter variance (Fig. 9 (b)), the epistemic and aleatoric uncertainties are initially similarly distributed. However, as the progress variable variance increases, aleatoric uncertainty plateaus, while epistemic uncertainty increases. This observation suggests that highest filter width variances were undersampled in the final dataset. This may be a consequence of the inadequacy of performing uniform-in-phase-space sampling with a clustering approach as discussed in Ref. [30].
Finally, it was consistently observed that the aleatoric uncertainty exceeded the epistemic uncertainty (Fig. 9, Fig. 10 and Fig. 11). Therefore, most of the errors may be attributed to the coarse-graining uncertainty stemming from the filtering operation and the feature selection rather than the lack of data. Albeit small, the epistemic uncertainty could also rapidly increase for a model deployed a posteriori in a chaotic system such as a turbulent combustion simulation [99], [100].
The extrapolative capabilities enforced with the synthetic dataset can serve two purposes: first, they can ensure that extrapolation is at least as accurate as some low-fidelity model; second, they can serve as a marker for OOD queries. In this work, the synthetic data labeling was set arbitrarily. Therefore the trained BNNs described in Sec. 3.4 best serve the second purpose: OOD query detection. The task of ensuring accuracy OOD is left for future work and would require labeling the synthetic dataset with the low-fidelity model prediction and is left for future work.
Using a BNN trained with synthetic OOD data, the OOD query detection method can be done as described in Algo. 14.
In Algo. 14, the distance \(d_{\rm OOD}\) to the arbitrary distribution of labels in the OOD region is used as a marker for whether or not an input \(\mathbf{x}\) is OOD. Clearly the choice of the threshold \(T\) impacts efficacy of the method. A systematic method to choose \(T\) would depend on the labels adopted for the synthetic OOD data and is out-of-scope of this work. In the following, \(T\) set to \(0.6 \times \sqrt{\mu_{\rm OOD}^2 + \sigma_{\rm OOD}^2}\) which was chosen through a sensitivity analysis (See 11).
The OOD query detection is applied to the \(D_{Le}\) and \(B_{Le}\) flame contours already used in Sec. 3.3 for validation of the BNN predictions. The BNN model used is the one trained with \(10^7\) synthetic OOD datapoints using the normalizing flow method (Algo. 5). Figure 15 (top) shows the in-distribution index that characterizes whether a query point (each point where the BNN is queried) is in or out-of-distribution. An in-distribution index is constructed with a value \(1\) if \(d_{\rm OOD} > T\) and \(0\) if \(d_{\rm OOD} < T\). Figure 15 (bottom) shows the distance between the query point and its nearest neighbor in the training and testing datasets. If this distance (hereafter referred to as the the nearest neighbor metric) is large, then the query point is far from any data point seen by the BNN during training, and the query point should be considered to be OOD. The metric is used to assess whether the BNN accurately labels OOD query points.


Figure 15: Top: in-distribution index predicted by Algo. 14 for filterwidth of 4 and 16. The isocontour denotes the flame contour based on the progress variable value. Bottom: distance between the query point and the nearest neighbor of the query point in the training and testing dataset.. b — \(B_{\rm Le}\) case., d — \(D_{\rm Le}\) case.
It is observed that the in-distribution index plots and the distance to nearest neighbor metric are in agreement. In particular, the BNN successfully labeled query points to be OOD when the distance to nearest neighbor is the largest. The OOD detection with the BNN can be done at the cost of BNN inference to compute the empirical averages of \(\mathbb{E} \left[\mathbf{y}\vert \mathbf{x}\right]\) (Algo. 18) and \(\mathbb{E} \left[\mathbf{y}\vert \mathbf{x}\right]\) (Algo. 19) and takes on the order of \(\mathcal{O}(10^{-3})\) seconds on a single CPU. The computation of nearest neighbor metric required 108 CPUh (using Intel Xeon Gold Skylake CPUs). Therefore, the OOD query detection using a BNN is computationally advantageous compared with more naive methods.
Figure 15 (b) (top) shows that the none of query points from the \(B_{Le}\) flame case could be considered out-of-distribution. This is a surprising result given that the \(B_{Le}\) flame was not included in the training dataset. Nevertheless, the phase-space spanned by the \(B_{Le}\) flame is similar to the training and testing dataset. This result highlights that good predictive performances of data-based closure models on cases not included in the training dataset does not test the extrapolative abilities of data-based closure models. Instead, it may test whether the dataset was sufficiently large so that the data-based model can be used on reacting flows not included in the training dataset.
In turn, Fig. 15 (d) (top) suggests that the case with filterwidth of 4 for the \(D_{Le}\) flame contains OOD query points. This is, again, a surprising result given that the \(D_{Le}\) flame was included in the training dataset. The nearest neighbor metric (Fig. 15 (d), bottom) confirms that this finding is reasonable since the OOD query points are indeed far from any training or testing data points. This effect could be explained by the fact that the training data was preprocessed with an approximate uniform-in-phase-space sampling to eliminate redundant datapoints. However the method used in Ref. [52] has been found to be inaccurate since it tends to remove too many rare datapoints [30]. Fortunately, the OOD query is more prominent for the filterwidth 4 which is also where the contribution of \(\chi_{\rm SFS}\) is negligible compared to \(\chi_{\widetilde{C}}\). This explains why the OOD queries did not affect the reconstruction of \(\widetilde{\chi}\) in this work (Sec. 3.3) or in Ref. [52].
For the purpose of decision-making, predictive simulations need to be paired with appropriate uncertainty estimates that reflect various uncertainty forms. Several efforts have shown that it was possible to estimate uncertainty due to boundary conditions [33], [57] and the closure model adopted for LES [57], [75] and RANS [101]. Model form uncertainty is typically represented by the variation a few (at most 3 in [57]) model parameters, the ones that appear in physics-based closure models [102]–[104] or that are used to superimpose physics-based models [75]. Here, a key complication is that the number of uncertain parameters is the number of weights (typically \(\mathcal{O}(10^3)\)) in the neural network. To enable uncertainty propagation, a key first step would be to reduce the number of uncertain weights. This task is left for future work.
This work presented the first demonstration of a Bayesian neural network approach for data-driven closure modeling equipped with aleatoric and epistemic uncertainty estimates. The main focus was on modeling subfilter progress variable dissipation rate, but the methods could be applied in general to closure modeling. a priori tests showed a good mean prediction of the subfilter progress variable dissipation rate, which suggests that including uncertainty estimates during training does not adversely affect the model accuracy.
Overall, the aleatoric uncertainty was found to outweigh the epistemic uncertainty. Epistemic and aleatoric uncertainty were found to similarly vary over phase-space, with the exception of large progress variable regions. This is likely a consequence of the approximate uniform-in-phase-space data sampling procedure. These findings motivate the use of non-uniform data selection which could compensate for local pockets of high aleatoric uncertainty in phase space.
A strategy for enforcing a specific OOD behavior was proposed using synthetic data generation. It was found that using a uniform in-phase space data generation led to the best performances and that the performance in-distribution did not degrade until the synthetic dataset size approached that of the original dataset. The OOD behavior can be used as a marker that efficiently tests whether the inputs passed to the BNN are far or not from the data-distribution. It was found that the strategy proposed is computationally efficient for identifying OOD queries.
Future work will include propagation of epistemic uncertainty through LES codes and analysis of the a posteriori uncertainty as well as dimension reduction strategies to efficiently propagate the uncertainty without sampling the high dimensional weight-space. Additionally, future work will explore the application of active learning or OED to inform future data collection and both a priori and a posteriori analysis of the model before and after the informed data collection.
GP: Methodology, Investigation, Software, Writing - original draft. MH: Conceptualization, Methodology, Investigation, Software, Supervision, Visualization, Writing - review & editing. SY: Conceptualization, Data curation, Funding acquisition, Supervision, Writing - review & editing.
The authors would like to thank Bruce Perry and Michael Mueller for fruitful discussions.
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0021110. This work was authored in part by the National Renewable Energy Laboratory (NREL), operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. This work was supported as part of DEGREES, an Energy Earthshot Research Center (EERC) funded by the U.S. Department of Energy, and by DOE’s Advanced Scientific Computing Research (ASCR) program. SY was supported by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Vehicle Technologies Office (VTO). The research was performed using computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.
To illustrate the effect of choice of input features on the uncertainties estimated, a BNN that uses 12 dimensions (referred to as 12D BNN) was trained with the following two additional features: 1) \(\widetilde{\dot{\omega_C}}\), the Favre-filtered reaction source term of the progress variable; and 2) \(\nabla \widetilde{T} \cdot \nabla \widetilde{C}\) where \(\widetilde{C}\) is the Favre-filtered temperature. Figure 16 shows the effect of expanding the feature space on the average aleatoric and epistemic uncertainties over the test dataset, between the 12D BNN and the original BNN trained with 10 input features (referred to as the 10D BNN). None of the datasets are augmented with synthetic OOD data.
When augmenting the input feature space dimension, the aleatoric uncertainty decreases. This can be explained by the fact that the additional input features reduce the loss of information in the filtered dataset. However, using a higher dimensional input space also increases the epistemic uncertainty. This phenomenon can be explained by the fact that data density is higher in the low-dimensional input space, thereby mitigating statistical errors that result from the lack of data.
In this appendix, algorithms for the computation of a predictive distribution (Alg. 17), the epistemic uncertainty (Alg. 18) and the aleatoric uncertainty (Alg. 19) are presented. Overall, the BNN allows for fast evaluation of the samples of weights, but appropriate averaging is needed to differentiate between epistemic and aleatoric uncertainties.
To compute the predictive distribution, we must both sample realizations of the weights (accounting for epistemic uncertainty) and sample from the resultant random variable (accounting for the aleatoric uncertainty). The can be achieved as a double for loop, with a third, outer loop describing the process of generating predictions across a set of input values. Naively, both inner sampling loops are achieved with Monte Carlo sampling and are subject to the standard \(\mathcal{O}(1/\sqrt{n})\) convergence rates, where \(n\) is the number of samples.
Computation of the epistemic uncertainty relies on the parameterization chosen for BNN output. In the case of a Kendall and Gal [77] style architecture, the model mean is directly parameterized as one of the model outputs. This fact may be exploited to avoid a sample average computation of the average. A sample variance calculation may then be performed on the collection of mean predictions. The Monte Carlo sampling and sample variance calculation is subject to the error prescribed by the Central Limit Theorem, again \(\mathcal{O}(1/\sqrt{N_e})\)
Computation of the aleatoric uncertainty similarly relies on the parameterization of the model outputs. With the variance, or standard deviation, as one of the model outputs, one may directly obtain this from a model evaluation after sampling the weights \(\mathbf{w}\sim q(\mathbf{w}\vert\theta)\). A sample average calculation of this collection then returns an estimate of the aleatoric uncertainty. Once more, the sample mean calculation is subject to the convergence properties of the Central Limit Theorem, \(\mathcal{O}(1/\sqrt{N_a})\).
To evaluate the impact of the prior selection as outlined in Sec. 2.3.1, a model was trained using an isotropic Gaussian prior for the weights, i.e. \(p(\mathbf{w})=\mathcal{N}(0,I)\). The conditional averages of the aleatoric and epistemic uncertainties for this model are presented in Figure 20. Compared to those presented in Fig. 9, we observe that the resultant epistemic uncertainty is larger, which we attribute to the prior term’s bias towards a model form that generates predictions with a mean of zero. There also appears to be a slight bias to the profile in Fig. 20 (a) as compared to Fig. 9 (a). Given the quantity of the data and the pre-processing steps to center the data, one might expect that an isotropic prior may be appropriate, and this indeed seems to be the case, however the results presented in the main text seem to have better agreement. As mentioned in Sec. 2.3.1, specification of the prior is an ongoing area of research.



Figure 20: Conditional average of the aleatoric and epistemic uncertainty (\(\times\) 100) with respect to (a) filtered progress variable, (b) filtered variance, and (c) filtered diffusion..
In Algo. 14, a key variable for the OOD of query points is the threshold \(T\) which discriminates between query points in the OOD regime and in-distribution regime. The value of \(T\) depends on the magnitude of the distance metric \(d_{\rm OOD}\), which itself depends on the magnitude of \(\mu_{\rm OOD}\) and \(\sigma_{\rm OOD}\). In Sec. 4.2, the value of \(T\) was therefore defined as \(\alpha \sqrt{\mu_{\rm OOD}^2 + \sigma_{\rm OOD}^2}\), where it was chosen \(\alpha = 0.6\). To define a reasonable value of \(\alpha\), a sensitivity analysis was performed: the value of \(\alpha\) was varied until the fields of in-distribution index became insensitive to \(\alpha\).




Figure 21: In-distribution index predicted by Algo. 14 for filterwidth of 4 and 16. Top: \(\alpha = 0.4\). Middle: \(\alpha=0.6\). Bottom: \(\alpha=0.8\). The isocontour denotes the flame contour based on the progress variable value.. c — \(B_{\rm Le}\) case., f — \(D_{\rm Le}\) case.
Figure 21 shows that the value of \(\alpha\) affects the in-distribution index fields. However, the in-distribution index fields is significantly less sensitive to \(\alpha\) if \(\alpha<0.6\) than if \(\alpha>0.6\), which explains the choice adopted for \(\alpha\) in Sec. 4.2. For all values of \(\alpha\), OOD queries detected are always clustered near the stretched portions of the flame front, which is indicative that the OOD queries coincide with inputs likely unseen in the training data. To systematically decide on the \(\alpha\) threshold, one could first define OOD data based on the nearest neighbor metric and choose the threshold such that the probability of a false positive is lower than a prescribed amount.