{'fraction_non_alphanumeric': 0.03366297727421573, 'fraction_numerical': 0.013786102560743991, 'mean_word_length': 4.363364945556726, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhanced fifth order WENO Shock-Capturing Schemes with Deep Learning

In this paper we enhance the well-known fifth order WENO shock-capturing scheme by using deep learning techniques. This fine-tuning of an existing algorithm is implemented by training a rather small neural network to modify the smoothness indicators of the WENO scheme in order to improve the numerical results especially at discontinuities. In our approach no further post-processing is needed to ensure the consistency of the method, which simplifies the method and increases the effect of the neural network. Moreover, the convergence of the resulting scheme can be theoretically proven. We demonstrate our findings with the inviscid Burgers' equation, the Buckley-Leverett equation and the 1-D Euler equations of gas dynamics. Hereby we investigate the classical Sod problem and the Lax problem and show that our novel method outperforms the classical fifth order WENO schemes in simulations where the numerical solution is too diffusive or tends to overshoot at shocks.


Introduction
Typically, numerical fluid mechanics deals with nonlinear hyperbolic partial differential equations (PDEs). In its simplest one-dimensional form, these equations can be represented as where x represents space, t denotes time, u(x, t) is conserved quantity and f (u(x, t)) is its flux. Discontinuous initial data develops after finite time a discontinuity ('shock wave') or a rarefaction wave (regardless of the smoothness of the initial or boundary data). Hence, suitable numerical methods for these problems must be designed to properly handle especially the accurate solution of the discontinuities. First in 1980, Crandall and Majda [1] proposed the class of monotone schemes that are nonlinearly stable in the L 1 norm and satisfy certain entropy conditions. It can be proven that the corresponding solutions converge to bounded variation entropy solutions including error estimates. However, these schemes are only first order accurate and by the fundamental Godunov theorem [2] it is known that one has to consider nonlinear non-oscillatory schemes to overcome this accuracy barrier. In this direction so-called shock-capturing schemes were developed that were able to resolve sharply a shock or a sharp gradient region without introducing too much diffusion or overshoot behaviour [3]. Additionally, to remedy the above mentioned drawback, at regions with smooth flow these schemes exhibit a rather high order of convergence. The well-known representative of this class of methods are the essentially non-oscillatory (ENO) schemes [4] with high order accuracy in smooth regions and sharply resolving shocks in an essentially non-oscillatory way using a smoothness indicator function, see e.g. [5]. Later on Jiang and Shu [6] further improved these schemes and proposed a weighted ENO (in the sequel abbreviated with WENO-JS) scheme, that is still regarded as a state-of-the-art solution approach.
Recently, machine learning was widely used to compute the solution of PDEs. We refer to [17,18,19], where the neural network algorithm is used to approximate a solution of a particular PDE problem. Following that approach, the solution of a particular PDE is a result of a neural network training procedure. Another idea is to improve a specific numerical scheme using neural networks. The training of a neural network is made offline and results in a new numerical scheme applicable to a wider class of problems. This idea was recently used by Beck et al. [20] for discontinuous Galerkin methods or in [21] for learning iterative PDE solvers and we also follow this approach.
The recent work of Stevens and Colonius [22] introduces new WENO-NN scheme based on neural network algorithm. In their work, the finite-volume coefficients of the WENO-JS scheme are perturbed, while maintaining the original smoothness indicators and nonlinear weights. However, the resulting scheme presented in their paper has only first order of convergence. Another neural network based WENO scheme was developed by Liu and Wen [23], where the new smoothness indicators are an output of the neural network algorithm.
We implement in our work another WENO extension based on deep learning. This approach is used to improve the classic WENO-JS and WENO-Z [7] scheme in this paper, but could be efficiently applied also to other WENO modifications. For this purpose we will train a rather small neural network to perturb the smoothness indicator functions of the WENO-JS scheme. As we do not develop any new smoothness indicators as in [23], but only their multiplicative perturbations, we are able to prove the convergence of the resulting scheme. We call this new scheme WENO-DS (Deep Smoothness), as we modify the smoothness indicators using deep neural networks. This scheme has less diffusion and less overshoot in shocks than the WENO-JS and the WENO-Z scheme, while maintaining high order accuracy in smooth regions. Unlike the recent work of Stevens and Colonius [22], we do not need any post-processing steps, which simplifies the procedure and also increases the effect of the deep learning part.
Finally, let us note that the use of WENO methods is not limited to hyperbolic PDEs, see e.g. [24] for an application in finance.
The paper is structured as follows. In Section 2 we introduce the WENO-JS and WENO-Z schemes under consideration in detail. Next, in Section 3 we introduce our deep learning approach, where neural networks are used to further improve WENO methods without any post-processing. The corresponding convergence proofs for two WENO schemes are given in Section 4. Then we present our numerical results in Section 5, which illustrate the improvements of our proposed method. Finally, we conclude our work in Section 6.

The WENO Scheme
Let {I i } be the partition of a spatial domain with the i-th cell ]. We consider a uniform grid defined by the points x i = x 0 + i∆x, i = 0, . . . , N , which are the centers of the cells with cell boundaries defined by The spatial discretization of one-dimensional conservation laws (1) yields a system of ordinary differential equations ('method of lines') and the resulting semi-discrete scheme is wheref is a numerical approximation of the flux function f . Following [6], if we define a function h implicitly by then (2) is approximated by where h i± 1 2 approximates the numerical fluxf ± 1 2 with fifth order of accuracŷ This results in a conservative numerical scheme.
To guarantee the stability of the method and avoid entropy violating solutions, the flux splitting method is applied, thus The numerical fluxf i± 1 2 is then represented byf i± 1 and the final scheme is formed as Next, only the construction off + i± 1 2 is considered. The negative part can be then obtained using symmetry (see e.g. [25]).

Fifth order WENO scheme
For a construction off i+ 1 2 , the fifth order WENO method uses a 5-point stencil divided into three candidate substencils, which are given by To form the numerical flux over the entire 5-point stencil, the numerical flux for each of these substencilsf m is calculated. These fluxes are then averaged in such a way, that fifth order convergence is ensured in the smooth regions. In regions with discontinuities, the weights should partly remove the contribution of these stencils so that the solution near the shock can be approximated in more stable manner.
Letf m (x) be the polynomial approximation of h(x) on each of the substencils (9). Then, evaluated at i + 1 2 we obtain where c m,j are the Lagrangian interpolation coefficients, dependent on m (see [6]). They take an explicit form f 0 and the numerical fluxesf m i− 1 2 can be obtained by shifting each index by −1. Using the Taylor series expansion it can be shown that: Then, the convex combination of the interpolated valuesf m (x i± 1 2 ) given byf yields the WENO approximation of the value h i± 1 2 , where ω m are the nonlinear weights defined as, cf. [6] The scheme using these nonlinear weights is denoted as the WENO-JS scheme. The parameter guarantees that the denominator does not become zero, and the coefficients d 0 , d 1 and d 2 are called ideal weights, which would form the upstream fifth order central scheme for the 5-point stencil and satisfy (5). Their values are: The parameter β m is called the smoothness indicator and is analyzed in the next section.

Smoothness Indicators
The role of smoothness indicators is to measure the regularity of the polynomial approximationf m (x) in each of three substencils. As developed in [6], they are defined as: Corresponding to the flux approximationf i+ 1 2 they take an explicit form The convergence analysis was performed in [6] and it was shown that if with D being a non-zero constant independent of m, the condition (23) is satisfied and the scheme has the expected fifth order accuracy. However, it was shown in [13] that at the critical points where the first derivative of f vanishes, the convergence order of the scheme from [6] decreases to the third order. Moreover, if the second derivative also vanishes, the order is further reduced to the second order. For a further explanation of this problem we refer to [13].

The WENO-Z scheme
In this paper we consider the modified WENO scheme of Borges et al. [7] with a new global smoothness indicator, which is characterized by It is easy to see from the equations (18) that The new WENO-Z weights are then defined by Borges et al. [7] have shown that when using these nonlinear weights, fifth order convergence is preserved, even at the critical points where f (u) = 0.

The Deep Learning Approach for WENO Schemes
To better capture discontinuities and avoid oscillations, we propose to apply deep learning to develop new smoothness indicators. We construct them as products of the original smoothness indicators β m and multipliers δ m which are outputs of a neural network algorithm. We refer to these new smoothness indicators as β DS m , where index DS corresponds to the new WENO-DS scheme: where C is a constant, which role is crucial for the proof of consistency and convergence and we will explain how to choose it in the Section 4. We emphasize that this formulation is very advantageous in a sense that the consistency and convergence properties of the original WENO method are preserved. In the case that the solution is smooth and the original smoothness indicator β m converges to zero, the smoothness indicator β DS m behaves in the same way. If the smoothness indicator β m is O(1), the multiplier δ m can change it so that the final scheme performs better.
In the original WENO method, the stencil (8)  we use (10) and shift each index by −1. In our approach we proceed as in the classical WENO method [5] and compute the smoothness indicators as described in (16) in the Section 2.2. We use them for a flux reconstructionf i+ and the values δ 0 , δ 1 , δ 2 are obtained by simple index shift so that it holds Finally, we obtain the flux approximations in the same way as in classical WENO schemes, but using the new smoothness indicators (29): which are used for approximating of a solution in a point x i . Superscripts p and n indicate the difference between the values at the same location smoothness indicators computation nonlinear weights computation nonlinear weights computation substencil flux computation Figure 1: The structure of the WENO method combined with the neural network algorithm. The white parts of the graph correspond to the original WENO method. The grey parts are added to this method so that the whole graph corresponds to the new method WENO-DS. 2k + 1 is the size of the receptive field of the whole CNN, x denotes the element-wise multiplication.
when we approximate the solution in points x i−1 and x i ). We present the whole algorithm of the method in the Figure 1.
As we mentioned before, the flux splitting technique (6) is used. Each part of a flux, f + and f − represents different type of input data to the neural network. Therefore we use two neural networks, for the positive and negative part of a flux with the input values f + (x i ) to the first neural network and f − (x i ) to the second neural network, i = 0, . . . , N . Each of the neural networks produces different outputs, multipliers corresponding to the positive and negative part of a flux. For simplicity we further drop the superscripts ± and when we talk about the inputs to the neural network we always mean both f + (x i ) and f − (x i ). We denote byf (x) the vector (f (x 0 ), f (x 1 ), . . . f (x N )) and formulate the neural network as a function F f (x) .
To ensure consistency we propose the use of a convolutional neural network (CNN). Firstly, this type of neural network is computationally efficient and secondly, it makes the multipliers independent of their position in the spatial grid so that the final numerical scheme is spatially invariant. To ensure the convergence of the method, we use the differentiable activation functions like the exponential linear unit (ELU) and the sigmoid function. If the layers of the neural network are differentiable functions, then their composition, the neural network function F (·), is also a differentiable function. The ELU activation function has the advantage that it does not cause the dying gradient problem, the sigmoid activation function ensures that the output of the neural network is between 0 and 1.
Let us note that we use for the implementation Python with the deep learning library PyTorch [26] (https://pytorch.org/), which is capable of GPU acceleration. Let us express the multipliers δ m,i for the smoothness indicators β m,i± 1 2 , m = 0, 1, 2 used in the node x i as the outputs of a neural network function. Following (30) and using the fact, that the layers of the CNN are spatially invariant differentiable functions, we can write where 2k + 1 is the size of the receptive field of the whole CNN and Φ is the function composition F •f . Then using (24) it holds Here we can proceed as in [13], but for the reader's convenience we repeat the steps of the proof: insert (35) into (14) and take = 0 This implies that where we used the fact that 2 m=0 d m = 1. Finally, substituting into (14) we where the superscript DS denotes the enhancement of the nonlinear weights (14) using our novel method. We see, that neither the condition (20), nor (23) is satisfied. However, as (38) holds, we can still guarantee convergence for the WENO-JS scheme with the smoothness indicators (28)  Let us now analyse the convergence of the scheme (27) with the new smoothness indicators (28). From (18) we see that the smoothness indicators β m,i± 1 2 are of the form and the global smoothness indicator (25) Then it holds We take = 0 and choose C such that Φ(x i ) + C > κ > 0, with κ fixed.
Then we see that in the non-critical points where f and the condition (20) is satisfied. Here we use the superscript DS denoting the enhancement of the nonlinear weights (27) using our novel method. Since we ensure P (x i ) > C > κ > 0, the multipliers P (x i ) do not introduce any further critical points. Therefore the analysis of the critical points with f x = 0 remains the same as in [7]. Thus we can guarantee the fifth order convergence of the scheme (27) with the smoothness indicators (28) also in the critical points.

Numerical Results
For the system of ordinary differential equations resulting from (2) we use a third-order total variation diminishing (TVD) Runge-Kutta method [27] given by ) and u n is the solution at the time step n. For (6) we consider in our examples the Lax-Friedrichs flux splitting where α = max u |f (u)|.

The Neural Network Structure
The proposed neural network algorithm can be generally applied to any type of conservation laws. For the equations where discontinuities or shocks are present, we propose to train a neural network separately for each equation class. Then we can better adjust the size of a neural network and its structure as well as the loss function, which leads to better results.
As we mentioned earlier, the inputs to the CNN are the values f + (x i ) and f − (x i ), i = 0, . . . , N and we train two neural networks for a positive and negative part of a flux. (The superscripts ± will be further dropped.) The first layer of the neural network is not learned, but represents a preprocessing of the solution from the last time step into a set of features that we assume to be suitable inputs for the following learned layers. Since our goal is to improve the smoothness indicators, we first compute the first and second central finite differences of f (x i ), i = 0, . . . , N . These parameters give us information about the smoothness of the solution and can facilitate and speed up the training of the CNN, and we can use a rather small CNN that still remains powerful. So we have the following values as input for the first learned hidden layer: The values f diff1 , f diff2 computed from f + from (46) represent the input values for the first neural network and the values f diff1 , f diff2 computed from f − represent the input values for the second neural network. Next, we use a fixed number of hidden layers, each with a specific kernel size and number of channels. We set these CNN parameters separately for each of the equation classes and experimentally find the best setting for each equation, keeping the size of proposed CNN small. We move the kernel by one space step so the stride is set to 1, and we use an ELU activation function in all hidden layers except the last one where we use sigmoid. In all our experiments, we set C = 0.1 in (28), which we experimentally found to be efficient. Let us note, that due to subsequent normalization of β DS m during the computation of nonlinear weights, using large value of C would decrease the effect of the trained multipliers. On the other hand, for C close to zero the experimental order of convergence could be smaller on rough grids (but still achieved for ∆x → 0). We use the nonlinear weights as defined in (27), replacing β m with β DS m . The value of is set to 10 −13 . As the first choice of the loss function we use the mean square error where u i is a numerical approximation of u(x i ) and u ref i is the corresponding reference solution. An advantage of this L 2 -norm based loss function in contrast to the L 1 -norm based loss function are stronger gradients with respect to u i resulting in faster training.

The Buckley-Leverett equation
In the first example, we apply our neural network algorithm to the Buckley-Leverett equation, which was also considered, for example, in [27,28,6]. It is a typical example with a non-convex flux function modeling a two-phase fluid flow in a porous medium [29]. The flux is given by where a < 1 is a constant indicating the ratio of the viscosities of the two fluids. The initial condition is set as and we use periodic boundary condition. In our implementation, we use the CNN with 3 hidden layers with the structure described in the Figure 2. The inputs to the learned hidden layers are the features (47). First, we create the dataset for which we compute the reference solution for the equation (1) with the flux (49) and the initial condition (50). We randomly generate the parameter a from a uniformly  For training, we proceed as follows. At the beginning, we randomly choose a problem and its reference solution from our dataset. The weights of the CNN are randomly initialized and we train our model on a solution where our computational domain is divided into 128×140 steps and successively compute the entire solution until the final time T . Using the solution at the time step n, we compute the solution at the time step n + 1 and during this computation the CNN is used to predict the multipliers of the smoothness indicators. After each of these time steps, we compute the loss and its gradient with respect to the weights of the CNN using backpropagation algorithm. Then we use this gradient to update the weights, using the well-known Adam optimizer [30] with learning rate 0.0001. After the last time step at time T , we test a model on a validation set and repeat the above steps. Then we select the model with the best performance on the validation set as our final model. For both training and comparing the performance of the models, we use the loss function defined as where LOSS MSE (u) is defined in (48) and represents the sum of the overflows of the solution above the maximum and below the minimum value of u, in our case u max = 1 and u min = 0. By adding this term to our loss function, we want to avoid the undesirable oscillations that occur especially in the first time steps of the solution. The Figure 3 shows how the value of the loss function for the problems from the validation set (which are not present in the training set) changes with increasing number of training cycles. As training cycle we denote a sequence of training steps performed on a solution for a single randomly chosen parameter a until the final time T . The loss is then evaluated at this final time T . We apparently see two optima for different values of a. If there are more than 50 training cycles, the loss begins to increase significantly for some problems, indicating that further training is not efficient. This may be caused by overfitting or divergence, so an early stopping algorithm could be implemented efficiently. Since we want only one final numerical scheme as output, we choose the model obtained after the 46th training cycle as the final model and present the result computed with it.
We compare the L ∞ -and L 2 -error in the Table 1 for the solution of the conservation law (1) with (49), a ∈ {0.25, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. Let us note that these parameters were neither in the training, nor in the validation set. We highlight the best performing WENO method using bold. In column denoted 'ratio' we divide the minimum error of WENO-JS and WENO-Z with an error of WENO-DS, showing how well our novel method works compared to the better one from the mentioned standard methods. WENO-DS outperforms the standard WENO methods in most cases. For a = 0.8, the error of WENO-DS is larger than the errors of the other two methods. How-   ever, this may be caused by the fact that the standard WENO methods perform disproportionately well for a = 0.8 compared to other values of the parameter a.
In the Figure 4, we present the solution of the Buckley-Leverett equation for the test problems with a = 0.25 and a = 0.5. We see that the WENO-DS gives a better solution quality than the WENO-JS or WENO-Z. Finally, we verify the theoretically proven fifth-order convergence of the WENO-DS scheme for a transport equation with a smooth solution given as with periodic boundary conditions. The results can be found in Table 2.
There is a great improvement when we compare our scheme with the WENO-NN scheme of Stevens and Colonius [22], where the resulting scheme has only first-order convergence.

The inviscid Burgers' equation
In the next example we consider the inviscid Burgers' equation, where the flux function in (1) is given by We consider following initial conditions where Using these initial conditions, we cover problems with both continuous and discontinuous initial conditions, and we simulate the shocks and discontinuities very well. We use periodic boundary condition. We first create the data set for training, in which we compute the reference solution of the Burgers' equation with the initial conditions (55)-(57). The computational domain is divided into 1024 space steps and 6400 time steps and the solution is computed up to time T = 0.3 using the WENO-Z scheme. We use the CNN with 3 hidden layers with the structure described in the Figure 5. For training, we proceed in the same way as in the previous example. The only differences are that the computational domain is divided into 128×100 steps and the learning rate used by the Adam optimizer is 0.001. In this example, we use the mean square error loss function (48) for training and validation. As the training on Burgers' equation exhibits much higher variance than in the Buckley-Leverett case, we performed 3 trainings each with 90 training cycles and finally selected the model showing the best performance on the validation set.    We compare the errors on the problems from the test set in Table 3 and 4. These were not in the training or validation set and the parameters were randomly generated. We observe rather small or no improvement for problems with the initial condition (56), but the improvement is significant for the solution with the discontinuous initial condition (55) as well as with    the initial condition (57). We conclude that the WENO-DS significantly outperforms the classical WENO methods. It should be noted that although our training set was created with the parameters sampled from uniform distribution as specified in (58), the method can also generalise for parameter values outside of these intervals, as can be seen in Table 4. Especially, we highlight the last two problems from Table 4, where we see a great improvement.
In the Figure 6 we show the solution of the Burgers' equation with the initial condition (55) for z 1 = 2.41, (56) for z 2 = 29.08, (57) for z 3 = 1.6 and (57) for z 3 = 2.12. We observe that WENO-DS captures shocks and discontinuities very well and gives us a better solution compared to WENO-JS and WENO-Z.

The one-dimensional Euler equations
We now investigate how WENO-DS behaves when applied to the onedimensional Euler system, which is considered a classical benchmark problem for methods for conservation laws. It has the form ∂ρ ∂t where ρ is the density, u is the velocity, p is the pressure and E is a total energy given by We take γ = 1.4, which is the ratio of the specific heats. To compute the fluxes, we use the characteristic decomposition of the system according to the steps in [5]. We use the Roe scheme to obtain the eigenvectors and eigenvalues [31] and the Lax-Friedrichs flux splitting to obtain the corresponding component of the flux. We take the solution based on [32] as the reference solution.
One of the most common benchmark problems is the Sod problem [33], where the initial condition is specified as , where u i is the local velocity and c i the local speed of sound. The solution consists of the left rarefaction wave, the right travelling contact wave and the right shock wave. We want to imitate this behavior of the solution, so we construct our data set as described in Appendix A.
We use the CNN with 3 hidden layers with the structure described in Figure 7. After projecting the flux and the solution on the characteristic fields using the left eigenvectors, we use Lax-Friedrichs flux splitting for each component of characteristic variables. From these values we compute the features (47), which are the inputs to the learned hidden layers. In this example we repeat the training procedure of the previous examples with some small modifications. To begin, we randomly generate the initial state from the dataset described earlier. We divide the spatial domain into 64 steps and compute the solution for the given initial state up to the time T = 0.1. After each time step we compute loss using the reference solution  from [32]. We use the gradient to update the weights, using Adam optimizer with learning rate 0.001. At the last time step we test the model and repeat the procedure with the new initial parameters (ρ, u, p). We use the loss function for training and validation. We show in the Figure 8 how the values of LOSS MSE (·) from (63) develop with the increasing number of training cycles. These values were obtained by testing the method on the validation problem, which was the Sod problem with the initial condition (61).
As we can see, loss decreases with the increasing number of training cycles. However, the final model of WENO-DS should be chosen carefully. In the Figure 9 and 10 we present the solution of the modified Sod problem with the initial condition (64) for ρ, u and p using 64 space points and using the models obtained after the 72th training cycle and after the 483th training cycle, where the minimal value of LOSS(ρ, u, p) is obtained. We compare the corresponding error values in Table 5. As we can see, the second model of WENO-DS leads to much smaller errors, compared to WENO-JS and WENO-Z. However, we observe an oscillation on the solution for p and u. On the other hand, using the first model obtained after the 72th training cycle, we see qualitatively very good solution, although the improvement on error is not so large. We also applied the trained models to the Lax problem [34] with an initial condition (ρ, u, p) = (0.445, 0.698, 3.528) 0 ≤ x ≤ 0.5, (0.5, 0, 0.571) 0.5 < x ≤ 1.
The solution also has a left rarefaction wave, a right traveling contact wave, and a right shock wave. In the Figure 11 and 12 we show the solution produced with both WENO-DS models mentioned before and see the similar behaviour of the solution. As indicated in Table 6, we observe even greater        improvement on errors than in the Sod problem using both models, especially using the second model. We can see, that the shape of the density profile is different, when compared to the Sod problem. Finally, let us note, that although we have not trained the presented models on the parameters that would lead to such a solution, the models are robust enough and can reliably detect the shocks.

Conclusion
In this work, we have improved the fifth-order WENO shock-capturing scheme by using deep learning techniques. To do this, we trained a relatively small neural network to obtain modified smoothness indicators of the WENO scheme. This was done in a way that avoided post-processing of the coefficients to ensure consistency. We applied our enhancement to the WENO-Z scheme, where the fifth-order convergence on the smooth solutions can be theoretically proven. Our new method, the WENO-DS scheme, is quite easy to use and significantly improves the numerical results, especially in the presence of discontinuities, even for cases that have not been trained before. We have demonstrated our results with the inviscid Burgers' equation, the Buckley-Leverett equation, and the 1-D Euler equations of gas dynamics. Finally, let us note, that this paper can be seen as a proof of concept, that neural networks can be efficiently combined with an existing numerical scheme, preserving its convergence order.