A machine learning approach to enhance the SUPG stabilization method for advection-dominated differential problems

1.   Introduction

The Galerkin-finite element (FE) method applied to partial differential equations (PDEs) with advection terms dominating over diffusion ones may suffer of numerical instability ^[22,34]. These numerical instabilities cause the numerical solution to exhibit oscillations that increase in amplitude with a local increment of the transport dominance over diffusion, i.e., as soon as the local Péclet number becomes larger than one. In order to eliminate (or at least mitigate) numerical instabilities, a general employed strategy is to consider the generalized Galerkin method, i.e., the Galerkin method with additional stabilization terms ^[34]. Examples of stabilization methods to reduce numerical oscillations in the advection dominated regimes are for instance the upwind and streamline-diffusion method, both methods that are not strongly consistent. Instead, examples of strongly consistent methods are the Galerkin-least-squares, the streamline upwind Petrov-Galerkin (SUPG), and the Douglas-Wang methods ^[6,34,41]. In particular, for these strongly consistent stabilization methods, the formulation depends on the residual of the PDE (in strong formulation), other than on a parameter – called stabilization parameter – whose definition is crucial for the success of the stabilization strategies. Determining the values of such stabilization parameter is not straightforward, especially for 2D and 3D problems. In addition, a universal formulation of such parameter is lacking, especially as it is strongly dependent on data of the model and the numerical setting used for the FE approximation of the PDEs, e.g., low and high order FE methods, spectral element methods, isogeometric methods, etc. ^[7,12,42]. Some of the formulations for the stabilization parameter have been derived from analytical considerations on the differential problem in 1D, others are suitable only for a specific numerical approximation method, while the most comprehensive ones incorporate some (empirical) dependence on numerical setting, as e.g., the order of the FE ^{[5,6,8,9,10,16,18,38,43,44,45]}.

In this work, we propose using machine learning (ML) ^[19,30,50] to make computers learn autonomously the values of the stabilization parameter for the FE approximation of advection dominated PDEs. Artificial neural networks (ANNs) ^[27,46] are widely popular in ML and Deep Learning for a wide array of applications: they are indeed very versatile tools that are increasingly finding their way in scientific computing ^[29,31], especially in the context of numerical approximation of PDEs. For instance, as substitute to standard numerical methods, ANNs can be employed as meshless methods in physics informed neural networks (PINNs) to directly approximate the solution of the PDE as it is trained by minimizing the (strong) residual of the PDE ^[35,36,37]. ANNs are also largely employed in a data-driven fashion in the context of model order reduction for parametric PDEs ^{[17,20,21,40,51]} and to enhance the stability properties of numerical methods for PDEs ^[13]. In fluid dynamics modelling, ANNs are massively adopted for flow features extractions, modelling, optimization and control. For flow features extractions, ANNs are used through clustering and classification to classify wake topologies ^[11]; for modeling fluid dynamics by reconstructing specific flows such as the near wall field in a turbulent flow ^[28] and for flow optimization ^[32], and control for aerodynamics applications ^[4]. ANNs are also used in a data-driven framework as a manner for providing alternative closure models for stress tensor ^[14] in Reynolds-Average Navier-Stokes (RANS) equations, for sub-grid scale models in large eddy simulation (LES) turbulence models ^{[23,24,48,49,52]}, or for model learning input-output relationships in complex physical processes ^[39].

In this work, we use ANNs to learn the optimal stabilization parameter in advection dominated PDEs that are discretized by means of the FE method: the goal is to enhance the accuracy of the SUPG FE method by optimally selecting the stabilization parameter under different data of the PDE and numerical settings of the FE method. We found that the proposed ANN-enhanced stabilization method allows to improve accuracy and stabilization properties of the numerical solution compared to those results obtained by analytical expressions of the stabilization parameter. The numerical results obtained shed light also on the possibility to apply the presented strategy to learn closure laws for stabilization and turbulence models of fluid dynamics, as for instance to learn the stabilization parameters in the Variational Multiscale–LES model to model transitional and turbulent flows ^[3,15,54].

This work is organized as follows: in Section 2, we recall the SUPG stabilization method for advection-diffusion equations; in Section 3, we present our numerical strategy to compute an optimal SUPG stabilization parameter through a feed-forward fully connected ANN. In Section 4 we validate our method by comparing the ANN results with those obtained with the 1D advection-diffusion problem from which the expression of the theoretical stabilization parameter has been derived. In Section 5 we first show the ANN's training and we present our numerical results on the 2D advection-diffusion problem used for training and we finally generalize our findings using the ANN's prediction on a different advection-diffusion problem. Finally, in Section 7 we draw our conclusions highlighting possible future developments.

2.   The SUPG method for advection-diffusion problems

We briefly recall the advection-diffusion differential problem and the SUPG stabilization method for the advection dominated regime.

Let $\Omega \in \mathbb{R}^d, \; d = 1, 2, 3$ be the physical domain with $\partial \Omega$ being its boundary. We consider the following problem in the unknown function $u$:

where $\mu$, $\boldsymbol{\beta}$, and $f$ are assigned functions or constants, with $\mu \in L^\infty(\Omega)$, $\boldsymbol{ \beta} \in [L^\infty(\Omega)]^d$, with $\nabla \cdot \boldsymbol{ \beta} \in L^2(\Omega)$, and $f \in L^2(\Omega)$. The Dirichlet datum on the boundary is $g\in H^{1/2}(\partial \Omega)$. Let $V_g = \{ v \in H^1(\Omega): \, v|_{\partial \Omega} = g \}$ and $V_0 = H_0^1(\Omega)$; the weak formulation of Eq (2.1) reads

with the bilinear form $a(u, v)$ and the linear functional $F(v)$ respectively defined as:

We consider a family of function spaces $V_h \subset V$ (either for $V_g$ and $V_0$) dependent on a parameter $h$ such that $\dim(V_h) = N_h < \infty$. Let $X_r^h = \{ v^h \in C^0(\overline{ \Omega}):\, v^h|_K \in \mathbb{P}_r, \, \, \text{for all } \, K \in \mathcal T_h\}$ be the function space of the FE discretization with piecewise Lagrange polynomials of degree $r \geq 1$, $\mathcal T_h$ a triangulation of $\Omega$ and $h$ the characteristic size of the mesh, comprised of elements $K \in \mathcal T_h$. By setting $V_h = {X}_h^r \bigcap V_0$, the Galerkin FE method applied to Eq (2.2) reads

where $V_{g, h} = {X}_h^r \bigcap V_g$. The standard Galerkin-FE method, in Eq (2.3), can generate numerical oscillations on $u_h$ if the problem is dominated by the advection term. In particular, these numerical instabilities can arise if the local Péclet number is $\mathbb{P}\mathrm{e}_h > 1$, where

The generalized Galerkin method (^[33]) allows for eliminating or mitigating numerical oscillations by adding stabilization terms to the standard Galerkin formulation as

In the SUPG method, the additional stabilization term reads:

where $R(u_h)$ is the residual in strong formulation of Eq (2.1), which is defined as:

The term $\tau_K$, appearing in Eq (2.6), is the stabilization parameter, which is the focus of this work. The stabilization parameter ${\tau}_K$ is generally defined locally, i.e., mesh element by element. In this paper, we consider a uniform stabilization parameter, thus $\tau_K = \tau$ for all $K \in \mathcal T_h$.

A universal and optimal definition of $\tau$ in terms of the problem data and numerical settings like the mesh size and FE degree is lacking. An extensive review of stabilization parameters for the SUPG method is reported in ^[25]. The most common choices for $\tau$ come from an analytic derivation made for the advection-diffusion problem in 1D with $f = 0$ and approximated by means of linear finite elements, which reads:

where $\xi(\theta)$ is the upwind function:

If a uniform mesh is used, as we consider in this paper, then the value of $h$ is uniform over $\mathcal T_h$; if in addition $\boldsymbol{\beta}$ and $\mu$ are constant, then this implies that $\tau$ is uniform over $\mathcal T_h$. The choice of $\tau$ made in Eq (2.7) represents an optimal choice of the stabilization parameter as it yields a nodally exact numerical solution for the 1D advection-diffusion problem if $f = 0$ and the FE polynomial order is $r = 1$ ^[18,34]. Thus, the stabilization parameter of Eq (2.7) may not be fully effective to provide the optimal stabilization for a general advection-diffusion problem, e.g., to guarantee a nodally exact solution in 2D/3D or when using FE degrees larger than $r = 1$. A commonly used generalization of the formula of Eq (2.7) to higher FE degrees ($r > 1$) is presented in ^[18] and reads:

Differently from Eq (2.7), the stabilization parameter $\widetilde{\tau}_r$ in Eq (2.8) takes into account the contribution of higher polynomial degrees $r$. Still, this formula is not optimal as it does not guarantee a nodally exact numerical solution. Our goal is to find a general and optimal expression of the stabilization parameter holding for advection-diffusion problems and FE approximations of degree $r\geq 1$ to be dependent on: the dimension $d$, the FE degree $r$, the mesh size $h$, the forcing term $f$, the diffusion coefficient $\mu$ and the transport coefficient $\beta$.

3.   ANN-based approach for determining the optimal stabilization parameter

We present our approach to determine the optimal stabilization parameter $\tau$ for the SUPG stabilization method by using an ANN. We consider as features (inputs) of our ANN: the FE degree $r$, the mesh size $h$ and the global Péclet number $\mathbb Pe_g : = | \boldsymbol{\beta}| L / (2\mu)$ of the advection-diffusion problem, where $L$ is the characteristic length of the problem, $\boldsymbol \beta \in \mathbb R^d$ and $\mu \in \mathbb R$. As we consider $\boldsymbol{\beta}$ to be uniform and fixed, using $\mathbb{P}\mathrm{e}_g$ as input corresponds to varying the value of the diffusion coefficient $\mu$. The features (input) of the ANN read:

the target (output) of interest is the optimal SUPG stabilization parameter that we denote with $\tau^*$:

The first step consist in generating the dataset, i.e., pairs of inputs and outputs to be used for training the ANN (data generation step). This consists in choosing repeatedly and randomly values of the features $\mathbf x^{(i)}$ in given ranges. These features are used to feed an optimization problem that, through an optimizer and a suitable error measure, provide an optimal stabilization parameter $\mathbf y^{(i)}$, i.e., the target of the given feature. Such optimization problem considers as error measure $E(\tau)$: the mismatch between the numerical solution $u_h$ and the exact one $u_\mathrm{ex}$ over the nodes of the FE mesh. $E(\tau)$ reads as:

being ${\mathit{\boldsymbol{x}}}_k$ the $k$-th node of the FE mesh $\mathcal T_h$ and $K_h$ the number of nodes. $E(\tau)$ is an approximation of the $L_1(\Omega)$ norm. To find the minimum of $E(\tau)$ for different problem configurations we solve, for each instance of the input parameters $\mathbf x^{(i)}$, the following optimization problem:

Thus, the dataset generated consists of $m$ pairs of inputs-outputs $(\mathbf x^{(i)}, \mathbf y^{(i)})$, with $i = 1, \dots, m$. The latter is used for training the ANN. The latter will be then used to predict the optimal stabilization parameter $\mathbf y^{(j)} = \tau_\mathrm{ANN}$ to be used for the FE approximation of the advection-diffusion problem in the settings provided by $\mathbf x^{(j)}$.

We report in Figure 1 a sketch of the procedure used to build the ANN for the prediction of the optimal stabilization parameter for any new feature $\mathbf x^{(j)}$.

4.   Numerical validation

In this section, we validate the proposed numerical strategy by means of a 1D advection-diffusion problem from which the expression of $\widetilde{\tau}_1$ in Eq (2.7) has been derived. In particular, we consider the advection-diffusion problem in Eq (2.1) with $\Omega = (0, 1)$, $f = 0$ and with Dirichlet BCs $u = 0$ on $x = 0$ and $u = 1$ on $x = 1$. It admits the following exact solution:

We recall that the stabilization parameter $\widetilde{\tau}_1$ of Eq (2.7) yields a nodally exact numerical solution if Eq (2.1) is solved by linear FE ($r = 1$).

We plot in Figure 2 the error measure $E(\tau)$ against the value of the stabilization parameter $\tau$, with $h = 1/20$, $\mathbb{P}\mathrm{e}_h = 12.5$ and by using $r = 1$ (Figure 2a) and $r = 3$ (Figure 2b). We observe that $E(\tau)$ shows a minimum in $\tau^*$, thus suggesting the possibility to use an optimization algorithm to solve the problem. Specifically, we employ the L-BFGS-B optimization algorithm from $\tt{ SciPy}$, an open source library for $\tt{ Python}$ ^[47]. Instead, the advection-diffusion problem has been solved using the FE open source library $\tt{ FEniCS}$ ^[2], by applying the SUPG method on a uniform mesh. Particularly, the optimal value $\tau^*$ found for the case $r = 1$ corresponds to the one provided by the theory in Eq (2.7). For the case $r = 3$, we observe a optimal value of $\tau^*$ for which the error is minimized: in this case, it does not corresponds to the stabilization parameter $\widetilde{\tau}_r$ provided by the theory in Eq (2.8): as a matter of fact, the latter arises from empirical considerations to extend the parameter $\widetilde{\tau}_1$ in Eq (2.7) to polynomials degrees $r > 1$. However, this does not ensure a nodally exact numerical solution.

We solve the optimization problem for different values of the local Péclet number $1 \leq \mathbb{P}\mathrm{e}_h \leq 250$ and we plot in Figure 3 the optimal stabilization parameter $\tau^*$ against the local Péclet number for $r = 1, 2$, and $3$.

Figure 3a shows that, by employing a linear FE space $r = 1$, the optimal stabilization parameter $\tau^*$ obtained by minimizing $E(\tau)$ exactly matches the theoretical one $\widetilde{\tau}_1$ for every value of $\mathbb{P}\mathrm{e}_h$, thus confirming the findings of Figure 2a. This also serves as validation of the optimization procedure that we proposed. By recalling that $\widetilde{\tau}_1$ of Eq (2.7) guarantees the numerical solution to be exact at nodes for this particular advection-diffusion problem and $r = 1$, we can infer that the optimization procedure is meaningful and it can be therefore exploited further, for example for $r > 1$. On the other hand, Figure 3b and 3c show instead a mismatch between the theoretical $\widetilde{\tau}_r$ of Eq (2.8) and the optimal one $\tau^*$, thus confirming the findings reported in Figure 2b.

In order to better appreciate the differences in the numerical solutions due to the choice of the stabilization parameters, we report, in Figure 4, a comparison of the error $E(\tau)$ obtained by means of the optimal $\tau^*$ and theoretical $\widetilde{\tau}_r$ stabilization parameters for $r = 2$ and $3$. Moreover, we report in Figure 5 a comparison among the exact solution $u$, the SUPG-stabilized numerical solution $u^*_h$ obtained with the optimal stabilization parameter $\tau^*$, and the numerical solution $\widetilde{u}_h$ obtained with theoretical one $\widetilde{\tau}_r$. Specifically, we consider FE of degree $r = 3$ and two different values of $\mathbb{P}\mathrm{e}_h$. In both these cases, the optimal parameter $\tau^*$ leads to a more accurate solution with respect to using the theoretical parameter $\widetilde{\tau}_r$. In particular, when $\mathbb{P}\mathrm{e}_h$ is "small", $\widetilde{\tau}_r$ leads to overshooting in the numerical solution, while if $\mathbb{P}\mathrm{e}_h$ is "large" the theoretical stabilization parameter leads to undershooting. Conversely, the optimal parameter $\tau^*$ accurately intercepts the exact solution at the nodes.

5.   Numerical results

We first introduce the training set of our problem and we detail the setup of the ANN that we use in this work. Then, we show the prediction of the stabilization parameter by means of the ANN on different advection-diffusion problems.

5.1.   Training the ANN for a 2D advection-diffusion problem

We apply the strategy presented so far to a 2D advection-diffusion problem to generate the dataset for the training of the ANN. Specifically, we consider in Eq (2.1): $\Omega = (0, 1)^2$, $f = 0$, and $\boldsymbol{\beta} = (1, 1)$. We prescribe the following exact solution on the whole boundary $\partial \Omega$:

We generate in $\Omega$ a structured mesh $\mathcal T_h$ of triangles with $\tt{ FEniCS}$ ^[2], as shown in Figure 6. We generate the dataset by repeatedly solving the optimization problem of Eq (3.4) for varying set of features as described in Section 3; specifically, we choose the features as reported in Eq (3.1). The complete dataset contains $m = \, 900$ examples with $r = \left \{1, 2, 3 \right \}$, $h = \left \{ \frac{\sqrt{2}}{10}, \frac{\sqrt{2}}{20}, \frac{\sqrt{2}}{40} \right \}$ and values of $\mathbb{P}\mathrm{e}_g$ randomly chosend in the range $[7, 70'710]$ (uniform distribution), which yields $\mu \in [10^{-5}, 10^{-1}]$. We summarize the details for generating the dataset in Table 1. Figure 7 provides an overview of the dataset used for the training of the ANN, specifically the features $\mathbf x^{(i)} = [r, \, h, \, \mathbb Pe_g]$ and the target $\mathbf y^{(i)} = [\tau^*]$. In this case, the characteristic length $L$ is fixed and set equal to $L = 1$.

We train a fully-connected feed-forward ANN on the generated dataset by using the open source library $\tt{ Keras}$ ^[26] built on top of $\tt{ TensorFlow}$ ^[1]. We divide the dataset into two parts: a training dataset that takes $80\%$ of the examples to be used for the ANN training and a validation dataset that takes the remaining $20\%$. We choose the loss function as the mean squared error that measures, for each training feature, the squared mismatch between the prediction of the ANN $\widehat{y}^{(j)}$ and the actual target $y^{(j)}$. Specifically, the loss function is defined as

We normalize the features by subtracting their sample mean and dividing by their sample standard deviation in order to help the weights to better adapt to the different scales of the features. Moreover, the targets and the feature $\mathbb{P}\mathrm{e}_g$ need special care as they are distributed over a wide range of values. Thus, we normalize them using by applying a base $10$ logarithm. The normalized features and targets read:

where $\overline{r}$, $\overline{h}$, $\overline{\log_{10}(\mathbb{P}\mathrm{e}_g)}$ are the sample mean of the training features $r$, $h$ and the logarithm of $\mathbb{P}\mathrm{e}_g$ respectively, while $\sigma_r$, $\sigma_h$ and $\sigma_{\log_{10}(\mathbb{P}\mathrm{e}_g)}$ are their sample standard deviations.

In order to find the best performing ANN architecture, we carried out a study by testing different numbers of hidden layers, nodes per layer, optimization algorithm, its learning rate, and batch size. A comparison of the loss function with different architectures is given in Figure 8. Using 3 hidden layers is beneficial in terms of validation error drops, while using more than 64 nodes per layer does not bring to considerable advantages. Finally, we choose an ANN with $3$ hidden layers, $64$ nodes per layer and an output layer with a single node, all using a rectified linear unit (ReLU) activation function. We display the ANN architecture in Figure 9. We trained the ANN with the SGD optimization algorithm, a constant learning rate of $0.01$, and mini-batch size of $32$ samples. Moreover, we employed a momentum of $0.9$ in the optimization algorithm to update the weights. The trained ANN is available in the $\tt{ GitLab}$ repository ^[53].

Regarding the computational efficiency of the proposed strategy, we stress the clear distinction between the offline phase (dataset generation and ANN's training) and the online phase, where we use the ANN to predict a new stabilization parameter – alongside the use of the FE solver – for unseen input parameter values. The most demanding part of the strategy is the dataset generation, requiring approximately 2 h on a standard laptop to repeatedly solve the optimization problem 900 times. Moreover, the ANN's training phase required approximately 10'. The testing (online) phase, which is the one that is performed for application purposes, is considerably inexpensive, requiring only the real-time evaluation of a composition of linear functions, comparable with the evaluation of the empirical relation that brings to the theoretical values (few milliseconds).

5.2.   Predictions of the stabilization parameter by ANN

Now, we compare the predictions of the ANN with the theoretical stabilization parameter $\widetilde{\tau}_r$ of Eq (2.8) ^[18,34]. We show in Figure 10 (left) the stabilization parameter $\tau_{\mathrm{ANN}}$ predicted by the ANN by varying mesh size $h$ and global Péclet number $\mathbb Pe_g$. For comparison, we report in Figure 10 (right) the corresponding theoretical stabilization parameter $\widetilde{\tau}_r$. We observe that the overall behavior of the trained ANN's predictions are qualitatively similar of the theoretical stabilization parameter $\widetilde{\tau}_r$. Nevertheless, it can be inferred that the values of the stabilization parameters are almost completely matched with linear FE ($r = 1$), while they quantitatively differ for $r = 2$ and $r = 3$. This was expected as $\widetilde{\tau}_r$ for $r > 1$ is an empirical extension of the formula for the case $r = 1$.

The ANN allows to make predictions with features outside the range of values for which it has been trained, that is for unseen values of such features. With this aim, we report in Figure 11 the comparison of the ANN's predictions $\tau_\mathrm{ANN}$ with $\widetilde{\tau}_r$ for the FE degree $r = 4$. This comparison shows a clear difference between the theoretical and ANN's $\tau$ for $r = 4$ even though the general trend is maintained similar.

5.2.1.   Test 1: predictions for the problem used in the ANN's training

We compare the numerical solutions $u^\mathrm{ANN}_h$ and $\widetilde u_h$ obtained by means of the SUPG stabilization method with the parameter $\tau_\mathrm{ANN}$ predicted by the ANN and the theoretical one $\widetilde{\tau}_r$, respectively; the comparison also involves the exact solution $u$ (5.1) of the 2D advection-diffusion problem used for the training of the ANN. In particular, we display the comparison of the former 2D solutions in Figure 12 along the line $(1-h, y)$ for any $y \in [0, 1]$ with: (a) $\mathbb{P}\mathrm{e}_h = 2$, $r = 1$, and $h = \sqrt{2}/10$ (Figure 12a); (b) $\mathbb{P}\mathrm{e}_h = 500$, $r = 3$, and and $h = \sqrt{2}/20$ (Figure 12b). We notice that in both the cases, the numerical solution $u_h^\mathrm{ANN}$ involving the stabilization parameter $\tau_\mathrm{ANN}$ provides more accurate results than with the theoretical stabilization paramter $\widetilde{\tau}_r$. In particular, in the case (a), $\widetilde{u}_h$ involves a much smoother boundary layer and overshoots the exact solutions $u$, conversely to the nearly nodally exact numerical solution $u_h^\mathrm{ANN}$. In the case (b), $u_h^\mathrm{ANN}$ provides a much better representation of the bundary layer, without the undershooting of the solution $u$ exhibited by $\widetilde u_h$. Furthermore, we report in Table 2 absolute errors between the numerical solutions ($u_h^\mathrm{ANN}$ and $\widetilde{u}_h$) and the exact one $u$. The errors, computed in $L^2(\Omega)$ and $H^1(\Omega)$ norms, show that the ANN-based SUPG stabilization method very often produces more accurate results than the ones obtained with theoretical stabilization parameter $\widetilde{\tau}_r$, especially for $r = 3$.

Moreover, to assess the ability of the ANN to predict the stabilization parameter out of the training range, we compare in Figure 13 the numerical solutions for high Péclet in the case $r = 4$. We observe that $u_h^\mathrm{ANN}$ is more accurate than $\widetilde{u}_h$ even with FE degree out of the training range.

5.2.2.   Test 2: predictions for an unseen problem with constant forcing term

We check the model generalization of the ANN trained for the 2D advection-diffusion problem with exact solution $u$ of Eq (5.1) by predicting $\tau_\mathrm{ANN}$ for an unseen 2D advection-diffusion problem. In particular, we consider the advection-diffusion problem of Eq (2.1) in $\Omega = (0, 1)^2$, with $f = 1$ and $\boldsymbol{\beta} = (1, 1)$. We prescribe the following exact solution on the whole boundary $\partial \Omega$:

We compare in Figure 14 the numerical solutions $u^\mathrm{ANN}_h$ and $\widetilde u_h$ with the novel exact solution $u$ for different Péclet and mesh sizes. As for the previous numerical tests, the stabilization parameter $\tau_\mathrm{ANN}$ predicted by the network provides more accurate numerical solutions $u_h^\mathrm{ANN}$ than with the theoretical stabilization parameter $\widetilde{\tau}_r$. In particular, the boundary layers and overall solution behaviours are better represented.

We better assess the former qualitative consideration, by computing the error $E(\tau)$ associated with the numerical solutions with SUPG stabilization for the parameters $\tau_\mathrm{ANN}$ and $\widetilde{\tau}_r$; different values of $\mathbb{P}\mathrm{e}_g$ and $r$ are considered. In particular, Figure 15 compares $E(\tau)$ against $\mathbb{P}\mathrm{e}_g$ (in logarithmic scale) for different values of $r$: the error obtained by using $\tau_\mathrm{ANN}$ is always lower than the one achieved with $\widetilde{\tau}_r$. These results indicate that $\tau_\mathrm{ANN}$, although trained on a specific advection-diffusion problem, can be used to make predictions of the optimal $\tau$ for an unseen advection-diffusion problem in place of the theoretical stabilization parameter.

5.2.3.   Test 3: predictions for and unseen problem with constant forcing term without exact solution

In this section, we consider an advection-diffusion problem with constant forcing term $f = 1$, advection coefficient $\boldsymbol{\beta} = (1, 1)$, and boundary conditions $u = 0$ on $\partial \Omega$. Although there is not an exact solution to the given problem, we consider as "ground truth" (reference) solution a numerical solution obtained on a much finer grid ($h = \sqrt 2/ 400$) without stabilization. In Figure 16, we compare the SUPG numerical solutions $\widetilde{u}_h$, $u_h^\mathrm{ANN}$ against our "ground truth" solution for two different Péclet numbers. The numerical solution with the stabilization parameter $\tau_\mathrm{ANN}$ provides more accurate results with respect to the the solution obtained with the theoretical parameter. In particular, we highlight the largest discrepancies observed in the corners of the extracted solution.

5.2.4.   Test 4: prediction for an unseen problem with a non-constant forcing term

In this section, we cope the case of a non costant forcing term by considering advection-diffusion problem of Eq (2.1) in $\Omega = (0, 1)^2$ with $\boldsymbol \beta = (1, 1)$, and we prescribe the following exact solution on the whole boundary $\partial \Omega$:

We display the exact solution of the considered problem in Figure 17.

In particular, Figures 18 and 19 show $E(\tau)$ against $\mathbb{P}\mathrm{e}_g$ for different values of $r$ with two different meshes with $h = \sqrt{2}/10$ and $h = \sqrt{2}/20$, respectively. The error obtained by using $\tau_\mathrm{ANN}$ is always comparable to the one achieved with $\widetilde{\tau}_r$, for both mesh levels and for all the FE degrees considered. These results suggest that $\tau_\mathrm{ANN}$, although trained on a specific advection-diffusion problem with $f = 0$, is robust with respect to unseen parameters and data in the advection-diffusion problem, including non constant $f$. Nevertheless, further studies can be conducted to better assess the role of a non-constant forcing term into the ANN's training phase, possibly encompassing the accuracy of the theoretical stabilization parameter in this scenario too.

6.   Direction of the advection field and stabilization parameters

We investigate the effect of the advection direction on the value of the stabilization parameter. We define $\theta$ as the angle between the advection velocity and the $x$-axis: $\boldsymbol{\beta} = (|\boldsymbol{\beta}| \cos \theta, \; |\boldsymbol{\beta}| \sin \theta)^T$. We carry out the optimization strategy introduced in Section 3 including also $\theta$ in the input parameter $\mathbf x^{(i)}$. We apply the optimization to the advection-diffusion problem with exact solution in Eq (5.3) varying $\theta$ in $[\pi / 12, \pi/2]$ and we report a comparison between the optimal, theoretical and ANN's stabilization parameter in Figure 20 (left). We set $|\boldsymbol \beta | = \sqrt 2$, $\mathbb{P}\mathrm{e}_g = 7'071$, $h = \sqrt{2}/20$, $r = 3$. We recall that, differently from the optimization strategy, in this plot we are still employing the ANN trained without accounting for the advection direction. We can observe that the optimal $\tau$ is affected by $\theta$, a feature that is not accounted by the theoretical stabilization parameter. Furthermore, in Figure 20 (right), we compare the mean errors obtained with the optimal, theoretical and ANN's stabilization parameter. It can be observed that the ANN still provides an advantage with respect to the theoretical stabilization. However, we believe that by including the advection direction as an additional feature in the training phase of the ANN, a higher accuracy can be obtained that the one achieved in this paper, where a single direction ($\theta = \frac{\pi}{4}$) has been considered.

7.   Conclusions

In this work, we presented an approach based on machine learning and ANN to compute the optimal stabilization parameter to be used in the SUPG FE approximation of advection-diffusion problems. Indeed, albeit the expression of the stabilization parameter is available for 1D problems and FE of degree $r = 1$, its extension to more general advection-diffusion problems (2D and 3D) and FE degrees $r > 1$ is still lacking.

We validated our approach against the 1D case, for which the ANN-stabilization parameter matches the already optimal, theoretical one for $r = 1$, leading to nodally exact numerical solutions Instead, for higher polynomials degrees $r > 1$, remarkable differences are observed among the theoretical and optimal stabilization parameter: we observed better accuracy and stabilization properties of the numerical solution with the optimal stabilization parameter with respect to the theoretical one.

Then, we generated the dataset on a 2D advection-diffusion problem, and we used it to train a fully-connected feed-forward ANN. We applied the predictions of the network on the original advection-diffusion problem for unseen input values and we also apply it to an unseen advection-diffusion problem to check for model generalization. Our numerical results showed that the proposed ANN-based approach provides more accurate numerical solutions than using the theoretical stabilization parameter for the SUPG method.

This work represents a step towards the enhancements of stabilization methods for the FE approximation of advection-dominated differential problems. In particular, this work is limited to few features for the ANN-based stabilization parameter and provides as output a single optimal stabilization parameter meant for the whole mesh FE $\mathcal T_h$. Natural extensions of this work will therefore involve local stabilization parameters, i.e., element by element over the mesh $\mathcal T_h$, to better account for non-uniform meshes, differential problems with varying coefficients, and capturing the local behaviour of the solution. In addition, to enhance the robustness of the ANN-based approach, possible additional inputs of the network are the angle of the advection velocity and the nodal value of the forcing term. A future development consists in the extension of the proposed approach to the 3D case for which, as for the 2D case, a universal definition of the stabilization parameter does not exist. The extension to the 3D case surely would present a larger cost in the training phase; however, once trained, the ANN could be used online real-time.

Finally, we envision extending the proposed machine learning approach to SUPG and variational multiscale stabilization methods for the FE approximation of the Navier-Stokes equations.

Acknowledgments

The authors acknowledge Prof. C. Canuto, Politecnico di Torino, for the fruitful discussions and useful suggestions. L. D. and A. Z. are members of the INdAM Research group GNCS. L. D. has been partially funded by the research project PRIN 2020, n.20204LN5N5, by MIUR. The authors ackwnowledge the anonymous reviewers for their insightful comments and suggestions.

Conflict of interest

The authors declare no conflict of interest.