WENO3-NN: A maximum-order three-point data-driven weighted essentially non-oscillatory scheme

https://doi.org/10.1016/j.jcp.2021.110920Get rights and content

Highlights

  • •
    A low-dissipation data-driven third-order WENO scheme (WENO3-NN) is proposed.
  • •
    Additional loss on the reconstruction weights yields maximum-order convergence.
  • •
    The network smoothness measure facilitates interpretation and generalization.
  • •
    Very good performance on canonical test cases, including strong shock interactions.
  • •
    Error magnitudes and wavenumber resolution comparable to or better than WENO5-JS.

Abstract

Neural networks have become more and more relevant for computational fluid dynamics. In recent works, neural network based weighted essentially non-oscillatory schemes have been developed. Challenges faced with such schemes are to ensure maximum-order convergence on narrow stencils and the ENO property. In this work, we use a neural network as a weighting function in the WENO scheme and address these shortcomings. Based on the input stencil, the neural network calculates a convex combination of local interpolation polynomials. We use a Galilean invariant embedding in the input layer and introduce an additional loss on the reconstruction weights, such that the WENO scheme inherently recognizes a smooth input function and achieves maximum-order convergence. The performance of the WENO3-NN scheme is demonstrated for one- and two-dimensional test cases, including strong shocks and shock-density wave interactions. The WENO3-NN scheme shows very good generalizability across all benchmark cases and different resolutions, and exhibits a performance similar to or better than the classical WENO5-JS scheme. By analyzing the approximate dispersion relation of the WENO3-NN scheme, we find that the neural network scheme learns a highly non-trivial dispersion-dissipation relation. Especially, data-driven schemes may introduce vanishing dissipation near the cutoff wavenumber which is counterintuitive to classical discretization-design principles.

Keywords

WENO
WENO-JS
WENO-Z
Neural network
Machine learning
Euler equations

1. Introduction

The solution of hyperbolic partial differential equations, e.g. the time-dependent Euler equations for compressible flows, requires numerical schemes capable of shock capturing and, at the same time, capable of resolving small-scale flow features. These seemingly contradictory demands, i.e. sufficient numerical dissipation around shock discontinuities and low-dissipation in smooth regions of the flow field, render the development of numerical schemes a long-standing challenge in computational fluid mechanics.
Weighted essentially non-oscillatory (WENO) [1] schemes are among the most popular choices to address aforementioned issues. The full stencil of a -order WENO scheme can be subdivided into r stencils, each of which forms a local interpolation polynomial. WENO schemes use a solution-adaptive nonlinear convex combination of these low-order polynomials. The contribution of the low-order polynomials is determined by nonlinear coefficients in the convex combination, so called weights. These weights are calculated from local smoothness measures of the solution such that in smooth parts of the flow field a high-order approximation is achieved while interpolation across discontinuities is avoided and sufficient numerical dissipation is introduced to allow for sharp shock capturing. The low-order polynomial of a stencil which contains a discontinuity receives essentially zero weight (ENO property) so that oscillations are suppressed. In smooth parts of the solution, the low-order polynomials are combined such that the background central upwinding scheme of order is recovered.
As shown by Henrick et al. [2], the classical WENO-JS scheme looses its full-order of accuracy at critical points. The proposed WENO-M scheme maps the classical WENO-JS weights such that the overall scheme recovers optimal-order of convergence. Borges et al. [3] proposed an improved WENO-Z scheme which assigns higher weights on less smooth stencils, for the first time, taking into account a global smoothness measure on the full stencil. In order to further decrease numerical dissipation, Hu et al. [4] propose an adaptive central-upwind 6-th-order WENO-CU6 scheme. WENO-CU6 has been further developed into an implicit LES model [5]. The family of targeted ENO (TENO) schemes by Fu et al. [6] is based on low-order polynomials with incrementally increasing stencil sizes. This allows for better treatment of multiple discontinuities while further minimizing the numerical dissipation.
In recent years, machine learning has become increasingly present in the field of fluid mechanics. Hybrid numerical schemes try to enhance classical numerical schemes by data-driven components. In the finite-volume method, increased interest has been placed on enhancing the cell face reconstruction process with data-driven algorithms, e.g. [7]. Stevens and Colonius [8] have proposed a WENO5-NN scheme in which a neural network maps the classical WENO5-JS weights in order to improve the solution. However, subsequent postprocessing of the new weights is necessary and the overall accuracy of the method is decreased to first-order.
In this work, we follow a different approach which is based on learning a complete smoothness measure that weights the standard ENO polynomials. This approach is similar to our previous work in [7], where we have used Harten polynomials as the basis of a data-driven finite-volume scheme for nonclassical shocks. By using ENO polynomials and a Galilean invariant input embedding for the neural network we incorporate prior knowledge and facilitate the machine learning task. This allows for very good generalization to unseen configurations and yields a highly interpretable neural network scheme. The resulting data-driven WENO scheme is trained offline on a set of canonical functions and later is applied to hyperbolic conservation laws. We show that by adaptive penalization of the deviation of the neural network output from the ideal weights (i.e. the weights that recover the background linear scheme), we can train a WENO-NN scheme which inherently achieves maximum-order convergence in smooth regions of the solution while maintaining the essentially non-oscillatory property at discontinuities. We focus on a third-order scheme due to its narrow stencil which is advantageous for practical applications. Comparing the weights of the proposed WENO3-NN scheme with the classical WENO3-JS scheme, the WENO3-NN maintains near-ideal weights over a large part of the solution field. One- and two-dimensional test cases show that WENO3-NN has very low dissipation and achieves a performance similar to the WENO5-JS scheme. An a posteriori analysis of the dispersion and dissipation characteristics of WENO3-NN schemes reveals that data-driven reconstruction schemes can generate a complex dispersion-dissipation relation which is counterintuitive to classical discretization-design principles. E.g. one variant of the WENO3-NN scheme introduces vanishing dissipation near the cutoff wavenumber.
The remainder of the paper is structured as follows. In Sec. 2, we review the classical WENO-JS scheme. In Sec. 3, we introduce the WENO3-NN method and describes the training procedure. Convergence behavior and characteristic behavior of the WENO3-NN scheme are discussed in detail. Section 4 presents applications to several 1D and 2D model problems, including the linear advection and Euler equations. Finally, we draw a conclusion and present an outlook in Sec. 5.

2. Review of WENO schemes

In this section we briefly review the classical third-order weighted essentially non-oscillatory scheme by Jiang and Shu (WENO3-JS) [1], [9] applied to hyperbolic conservation laws in a finite-volume framework. Without loss of generality, we consider the one dimensional scalar hyperbolic conservation law(1) Application to multiple dimensions or systems of equations is handled via the so called dimension-by-dimension technique. We discretize the spatial domain into finite-volumes of uniform size Δx and denote by the cell center and by the cell faces of the i-th cell. The semi-discrete finite-volume formulation of Eq. (1) reads(2) where(3) is the cell average of in cell i. According to the method of lines, Eq. (2) yields a system of ordinary differential equations (ODEs) which can be integrated in time by any ODE-solver, e.g. TVD Runge-Kutta schemes. The physical flux in Eq. (2) is approximated by the numerical flux function , so that in practice we solve(4) where the numerical flux is a two-argument function of left and right cell face values(5) The unknown cell face values have to be reconstructed from the known cell average values . In the following, we will only present the calculation of . The computation of is analogous. We will drop the superscript + for simplicity.
The classical third-order WENO3-JS scheme builds as a convex combination of the two-point substencils and . Each substencil, and , implies a second-order approximation of the cell face value , so that the convex combination reads(6) For third order, the interpolants are given as(7) The classical WENO3-JS scheme finds the normalized nonlinear weights in Eq. (6) as(8) are unnormalized weights and are the smoothness measures. are the ideal weights which generate the third-order upwind-biased scheme for the values . The parameter p enhances a scale separation of the smoothness parameters. Here we set , and avoids division by zero. The smoothness indicators are calculated as(9) The evaluation of the smoothness indicators in Eq. (9) gives(10) We summarize two central ideas behind the WENO methodology: On the one hand, in smooth flow regions maximal accuracy (here, third order) is desirable. Therefore, in smooth parts of the solution the smoothness indicators are similar to each other, thus the weights approach the ideal weights . On the other hand, the influence of a stencil containing a discontinuity should be diminished in the interpolation process Eq. (6). The smoothness measure of a discontinuous stencil is , so that the corresponding weight is relatively small (nearly zero) compared to smoother stencils. This guarantees that the interpolation process is essentially non-oscillatory (ENO property).

3. WENO3-NN

3.1. Neural network basics

Neural networks are parameterizable nonlinear compound functions that map any input x to an output , where θ are free and learnable parameters. Deep neural networks (DNNs) consist of multiple hidden layers of units (so called neurons) between in- and output layer. They perform successive elementary nonlinear transformations to map x to y. The numerical values in each layer are called hidden-state activations.
In multilayer perceptrons (MLPs), neurons in adjacent layers are densely connected. The vector of activations in layer l is computed from the activations of the previous layer by first applying an affine linear transformation, followed by an element-wise nonlinearity . The activation of the i-th neuron in layer l denoted as is calculated as(11) where indicates the weight matrix linking layers (l−1) and l, bil−1 is the bias vector, and Nl−1 indicates the number of neurons in layer (l−1). The error between the network prediction y and the true output yˆ is calculated by a suitable loss function L(y,yˆ). Training a neural network means finding a set of parameters θ that approximately minimizes the selected loss function. Typically, the loss function is minimized via mini-batch gradient descent or the popular Adam optimizer [10].

3.2. WENO3-NN Architecture

In this section, we introduce the WENO3-NN architecture. The polynomial weights ωk from above are functions of the cell values in the full input stencil S3=∪k=0r−1Sk with r=2, i.e.(12)ωk=ωk(u¯i−1,u¯i,u¯i+1), subject to ∑kωk=1 and 0≤ωk≤1∀k. Therefore, any function f:R3→R2 subject to aforementioned restrictions can be regarded as a third-order WENO-type weighting function. Naturally, any artificial neural network with suitable in- and output space represents a parameterizable WENO weighting function(13)ωkNN=NN(u¯i−1,u¯i,u¯i+1;θ).The proposed WENO3-NN architecture is displayed in Fig. 1. First, the input values from the 3-point stencil S3={xi−1,xi,xi+1} are passed to the so called Delta layer which calculates input features for the neural network. In this work, we use the following features(14)Δ˜1=|u¯i−u¯i−1|,Δ˜2=|u¯i+1−u¯i|,Δ˜3=|u¯i+1−u¯i−1|,Δ˜4=|u¯i+1−2u¯i+u¯i−1|,(15)Δj=Δ˜j/max⁡(Δ˜1,Δ˜2,ϵ)j∈{1,2,3,4}, where ϵ=10−15. The four features {Δ1,Δ2,Δ3,Δ4} are essentially the normalized amplitudes of first and second-order derivatives of the local input field, and thus are expected to give a good measure of regularity of the underlying function. In line with classical WENO weighting functions, the features are designed to be Galilean invariant such that the overall weighting function is Galilean invariant as well. We note that the proper design of the input features injects appropriate prior knowledge into the neural network and eases generalization. {Δ1,Δ2,Δ3,Δ4} are then passed through a multilayer perceptron. Here, the network consists of an input layer, 3 hidden layers with 16 nodes each, and a softmax output layer. The softmax function is a natural choice for the output activation as it implies ∑kωkNN=1 and 0≤ωkNN≤1∀k. We use the swish activation function [11] in the hidden layers, swish(x)=x⋅σ(x), where σ(x)=(1+exp⁡(−x))−1 is the sigmoid function. Upon proper training, we expect the softmax activation function to assign low weights to low-order polynomials that contain discontinuities. During numerical experimentation, we found that it is difficult for the network to output essentially zero weights due to the saturation of the softmax activation function. This is mainly a problem for test cases including very strong shocks such that spurious oscillations are no longer suppressed and negative densities or pressures might be reconstructed, e.g. the interacting blast wave test case in Sec. 4. Therefore, at test time, we pass the ωkNN weights through a so called ENO layer which restores the ENO property. This is essentially a sharp cutoff function followed by a renormalization (similarly to [6]), i.e.(16)ω˜kNN=ωkNNψk∑jωjNNψjwithψj={0,ifωjNN<ceno,1,otherwise, where ceno is the cutoff threshold. We use ceno=2⋅10−4, see A.4 for more details.
Fig. 1
  1. Download: Download high-res image (168KB)
  2. Download: Download full-size image

Fig. 1. Schematic of the WENO3-NN architecture. The cell-averaged input stencil S3={xi−1,xi,xi+1} is passed to the Delta layer which extracts the Galilean invariant features Δj. The Δj are the input of the neural network which calculates ωkNN. During training ωkNN are used for reconstructing the cell-face values. At test time, the final WENO weights ω˜kNN are obtained after application of the ENO layer. For legibility, the neural network is depicted with two hidden layers.

Upon proper training, the neural network inherently learns a smoothness measure and a weight function. Restricting the neural network as a weight function has several advantages: Firstly, we make use of the well defined ENO interpolants. Secondly, the network outputs are interpretable. Finally, we can identify a posteriori the analytical relation that the network has learned. Thus machine learning enables a new avenue of finding improved smoothness measures.
As mentioned above, suitable WENO weights should be close to the ideal weights in smooth flow regions while discontinuous stencils should be assigned effectively zero weights. In previous works, WENO neural networks schemes have been only trained on a reconstruction error, i.e. finding a suitable mapping such that the cell face values are reconstructed in a best possible way. Unsurprisingly, such schemes were not able to achieve full-order convergence upon mesh refinement. Ideally, we want the network to fall back to the ideal weights when the input stencil is smooth enough.
Therefore, we introduce the following loss function(17)L=Lr+βdLd+βW‖W‖22,(18)Lr=1Nb∑s=1Nb(γs)α(ui+1/2NN,s−ui+1/2s)2,(19)Ld=1Nb∑s=1Nb(1−(γs)α)∑k=0r−1(ωkNN,s−dk)2. The summation of the loss functions is over a mini-batch consisting of Nb samples. The total loss L is composed of three individual parts: the cell face reconstruction loss Lr, the deviation from the ideal weights Ld, and an ℓ2-regularization loss on the neural network weights ‖W‖22 to prevent overfitting. γ is a data dependent parameter that adaptively weights the two loss components Lr and Ld for each sample. γ∈[0,1] measures the well-resolvedness of the function in the given stencil. For a completely resolved function (e.g. a linear function) γ=0, and we only require that the network tries to approximate the ideal reconstruction weights dk. For a discontinuous function γ=1, and the network has no a priori bias w.r.t. the ideal weights, only the reconstruction loss is penalized. For a function that is not well resolved on the given stencil, 0<γ<1, so that the network has to strike a balance between reconstruction loss and deviation from the ideal reconstruction weights. This a priori bias incentivizes the scheme to output the ideal weights dk while providing the neural network with the liberty to deviate from dk if the reconstruction loss is significantly reduced. The α exponent creates a scale separation mechanism. Increasing α or βd pushes the neural network towards the ideal weights dk. We highlight the influence of α in A.3. Note that α=0 penalizes only the reconstruction, α→∞ would recover the background linear scheme.
In order to find a suitable γ function for the WENO3-NN scheme, we expand u¯N in terms of a Taylor series around the cell center xi(20)u¯N(x)=u¯N+∂u¯N∂xΔx+12∂2u¯N∂x2Δx2+O(Δx3), where the right hand side is to be evaluated at the cell center xi. We consider a function to be smooth in the neighborhood of xi when the contribution of high-order derivatives is small compared to the contribution of their low-order counterparts. E.g., when only information on first and second-order derivatives is available, we consider a function as smooth if(21)|12∂2u¯N∂x2Δx2||∂u¯N∂xΔx|<<1. Therefore, for the third-order WENO3-NN scheme, we use(22)γ=|u¯i−1−2u¯i+u¯i+1||u¯i−u¯i−1|+|u¯i+1−u¯i|+ϵγ, where ϵγ=10−15 is a small number to avoid division by zero.

3.3. Training

The training dataset is composed of canonical functions that mimic local solution features of hyperbolic conservation laws. Table 1 summarizes the training dataset. We use polynomials up to degree 3, jump discontinuities, sawtooth functions, and trigonometric functions. Polynomials and the tanh functions are evaluated on the domain [−1,1], while all other functions are evaluated on [0,1]. We use a discretization of Δx=0.01 for the dataset. For jumps and sawtooth functions, only stencils that include a discontinuity are included in the dataset. During network training, we split the data between training and validation set with a validation split of 0.1. We train the network with the Adam optimizer [10] for 100 epochs. The learning rate is fixed at η=10−3, and we use a mini-batch size of Nb=100. During our hyperparameter study, we find two optional WENO3-NN schemes, one with α=0.1,βd=0.1 (in the following denoted as WENO3-NN1) and one with α=0.03,βd=0.1 (in the following denoted as WENO3-NN2), that perform very well over a number of different test cases. In the following, we will focus on these two WENO3-NN variants. Further details on the datasets and the training process are provided in A.1 and A.2, respectively.

Table 1. Training dataset. U represents a uniform distribution, B represents the Bernoulli distribution.

Function f(x)ParametersNumber of samples
∑k=0nakxkak∈U(−1,1)∀k4000 for each k
ul(x < 0.5)+ur(x > 0.5)ul,ur∈U(−10,10)8000
(−1)ax + δ(x > 0.5)a∈B(0.5),δ∈U(0.5,2)4000
sin⁡(kπx)k∈U(2,20)4000
tanh⁡(kx)k∈U(5,30)4000

3.4. Convergence behavior

In the following, we show that including an additional loss term Ld on the deviation from the ideal polynomial weights dk restores full-order accuracy. To this end, we compare the convergence behavior of the WENO3-NN with (α=0.1) and without (α=0) the ideal weights penalty in Fig. 2. As test functions, we use tanh⁡(2x) and sin3⁡(πx) on the domain [−1,1]. When no penalty on the deviation from the ideal weights is applied (α=0), the WENO3-NN scheme yields second-order convergence. The neural network has not learned to adapt the reconstruction in the asymptotic limit to recover the third-order central-upstream scheme. The overall accuracy degenerates to the accuracy of the low-order polynomials (in this case second-order polynomials). However, when a deviation from the ideal weights is penalized (α=0.1), the neural network learns to assemble the upstream linear background scheme for smooth inputs, and the WENO3-NN scheme achieves third-order convergence for the tanh test case. For the more complex sin3⁡(πx) that involves smooth extrema and critical points, the WENO3-NN scheme achieves a convergence order of O(2.5) which is in accordance with the WENO3-JS and WENO3-Z scheme. Note, that both WENO3-NN schemes show a drop in the L2 error in the coarse regime for sin3⁡(πx). This corroborates that the WENO-NN schemes presented here comply with our intention: they have the potential to improve on classical WENO methods for coarsely resolved flow features while maintaining full-order convergence upon mesh refinement.
Fig. 2
  1. Download: Download high-res image (203KB)
  2. Download: Download full-size image

Fig. 2. Convergence behavior of the WENO3-NN compared with other third-order WENO schemes.

3.5. Behavior at smooth extreme points and discontinuities

We analyze the behavior of WENO3-NN1 (α=0.1,βd=0.1) and WENO3-NN2 (α=0.03,βd=0.1) in smooth regions and near discontinuities by computing the weights ω0 and ω1 for the function(23)u(x)={sin⁡(2πx)if0≤x≤0.5,1−sin⁡(2πx)if0.5<x≤1.0. Fig. 3 shows the weight ω0 at all cell faces xi+1/2 (ω1=1−ω0 is omitted for legibility). The black line represents the ideal weight d0=1/3 for smooth regions. Both WENO3-NN schemes are able to follow the ideal weight for most parts of the domain, and show smaller deviations than the WENO3-N scheme [12]. At both extreme points as well as at the jump discontinuity the WENO3-NN scheme degenerates and gives ω0=0 or ω0=1, respectively. The analysis of the WENO3-NN weighting function underscores the high interpretability of the proposed framework. Compared to WENO3-NN1, the WENO3-NN2 scheme with α=0.03 deviates slightly stronger from the ideal weights as a smaller α prioritizes the reconstruction error over the ideal weights error. We note that the WENO3-NN scheme has learned a weighting function which shows some qualitative similarities to the WENO-F3+ scheme [13], although the deviation from the ideal weight is much smaller for WENO3-NN.
Fig. 3
  1. Download: Download high-res image (213KB)
  2. Download: Download full-size image

Fig. 3. Comparison of the smoothness measures of selected WENO schemes for Eq. (23).

Additionally, we are interested in the reconstruction behavior of the WENO3-NN scheme around saddle points. Fig. 4 shows the weight distribution ω0 for the function(24)u(x)=sin3⁡(πx),−1≤x≤1. While WENO3-JS and WENO3-Z degenerate at the saddle point (x=0), the WENO3-NN schemes are able to distinguish the saddle point, similar to the WENO-F3+ scheme, and adapt the reconstruction accordingly.
Fig. 4
  1. Download: Download high-res image (267KB)
  2. Download: Download full-size image

Fig. 4. Comparison of the smoothness measures of selected WENO schemes for Eq. (24).

4. Results

In the following, we illustrate the WENO3-NN scheme by solving the linear advection equation and the Euler equations. The one-dimensional linear advection equation is given by(25)∂u∂t+∂u∂x=0. The three-dimensional, compressible Euler equations are(26)(ρρuρvρwE)t+(ρuρu2+pρuvρuw(E+p)u)x+(ρvρvuρv2+pρvw(E+p)v)y+(ρwρwuρwvρw2+p(E+p)w)z=0. E=ρe+12ρ(u2+v2+w2) is the total energy per unit volume, where e is the internal energy per unit mass. Here, the internal energy is closed by the equation of state for an ideal gas, e=p/((γ−1)ρ) with γ=1.4 if not stated otherwise.
We use Local Lax-Friedrichs (LLF) flux splitting (also called Rusanov method). The 3rd-order TVD Runge Kutta scheme is used to integrate the equations in time. Unless stated otherwise, all computations are carried out with CFL=0.6 and performed on a uniform mesh. For the Euler equations, we apply the WENO-reconstruction on the characteristic variables [14]. Throughout this section, we compare the WENO3-NN scheme α=0.1,βd=0.1 (denoted as WENO3-NN1) and the WENO3-NN scheme with α=0.03,βd=0.1 (denoted as WENO3-NN2) with the third and fifth-order classical WENO-JS schemes (WENO3-JS and WENO5-JS, respectively), the third-order (improved) WENO3-Z scheme [15], the fifth-order WENO5-Z scheme [3], and the WENO3-N scheme [12]. Although it is well known that the classical WENO5-JS scheme is strongly dissipative and that many improved WENO5 schemes have been developed [2], [3], [16], the approximation quality of the fifth-order WENO5-JS scheme is still a good reference for any improved third-order WENO scheme due to its wider five-point stencil and the therein contained information on higher order derivatives. For all classical WENO schemes, we use the parameter ϵ=10−15 to avoid division by zero. Reference solutions labeled as “Exact” are obtained by the WENO5-JS scheme on a uniform mesh with N=4000 points.

4.1. Linear advection

We test the WENO3-NN scheme on the linear advection equation (25). We consider a discontinuous initial condition consisting of a Gaussian, a square, a triangle, and a semi-ellipse (GSTE), given by(27)u0(x)={16(G(x,β,z−δ)+G(x,β,z+δ)+4G(x,β,z))if−0.8<x<−0.6,1if−0.4<x<−0.2,1−|10(x−0.1)|if0.0<x<0.2,16(F(x,β,a−δ)+F(x,β,a+δ)+4F(x,β,a))if0.4<x<0.6,0otherwise, where G(x,β,z)=eβ(x−z)2, F(x,α,a)=max⁡(1−α2(x−a2),0). The constants are z=−0.7,δ=0.005,β=log⁡236δ2,a=0.5, and α=10. The domain is [−1,1] and we apply periodic boundary conditions. The resolution is N=200. The initial condition is transported up to t=10.0.
Fig. 5 compares the results of the proposed WENO3-NN scheme with WENO3-JS/Z/N and WENO5-JS/Z. The results for WENO3-NN outperform the WENO3-JS/Z/N schemes, and are much closer to WENO5-JS. The comparatively low dissipation of the WENO3-NN scheme yields better results especially near spikes and discontinuities compared to WENO3-N. While WENO3-NN is not able to achieve the same performance as WENO5-JS/Z for functions with a local maximum (i.e. the Gaussian, the triangle, and the semi-ellipse), it is interesting to note that the diffusion of the square discontinuity is similar for WENO3-NN and WENO5-JS/Z. WENO3-NN2 achieves slightly better results than the WENO3-NN1 scheme.
Fig. 5
  1. Download: Download high-res image (302KB)
  2. Download: Download full-size image

Fig. 5. Linear advection of GSTE at t = 10.0. Resolution N = 200.

Another standard test for WENO reconstruction schemes is the advection of a nonlinear discontinuity defined by the initial condition(28)u0(x)={−sin⁡(πx)−12x3if−1≤x<0,−sin⁡(πx)−12x3+1if0≤x≤1. We advance the solution up to t=10.0. Fig. 6 shows the numerical solutions of the different WENO schemes. Compared with WENO3-N, the WENO3-NN schemes yield substantially better results and are able to capture the discontinuity better. Note that the WENO3-NN schemes also give better results in the smooth part of the flow field, especially near the smooth extreme points where they are nearly identical with WENO5-JS.
Fig. 6
  1. Download: Download high-res image (186KB)
  2. Download: Download full-size image

Fig. 6. Linear advection of the discontinuous initial condition in Eq. (28) at t = 10.0. Resolution N = 200.

Finally, we consider the advection of the initial distribution u0(x)=sin3⁡(πx) up to t=10.0. Fig. 7 shows the results of the WENO3-NN scheme at a resolution of N=100 points. Note that this is a coarser resolution than the training dataset which is generated at Δx=0.01. In this test case, the behavior at the extreme points and the saddle point are of interest. As discussed earlier, the proposed WENO3-NN schemes are able to detect the saddle point and adapt the reconstruction accordingly. The results are remarkably better than other third-order WENO schemes. The low dissipation allows a good reconstruction at the extreme points.
Fig. 7
  1. Download: Download high-res image (195KB)
  2. Download: Download full-size image

Fig. 7. Linear advection of sin3⁡(πx) at t = 10.0. Resolution N = 100.

Additionally, we provide the pointwise error distributions and integral L1 error norms for the linear advection tests in Appendix A.6.

4.2. Shock-tube problems

We test the WENO3-NN schemes on the shock-tube test problems: namely, the Sod problem [17], the Lax problem [18], and the 123 problem [19]. We use a resolution of N=200 grid points for all three test problems. The initial condition of the Sod problem is(29)(ρ,u,p)={(1,0,1)if0≤x<0.5,(0.125,0,0.1)if0.5≤x<1.0, and the final simulation time is t=0.2. Fig. 8 shows the solution at the final simulation time. The WENO3-NN schemes capture the contact discontinuity and the right-moving shock much sharper than WENO3-JS/Z/N, and give a nearly identical performance to WENO5-JS.
Fig. 8
  1. Download: Download high-res image (162KB)
  2. Download: Download full-size image

Fig. 8. Sod shock tube at t = 0.2. Resolution N = 200.

The initial condition of the Lax problem is(30)(ρ,u,p)={(0.445,0.698,3.528)if0≤x<0.5,(0.5,0,0.5710)if0.5≤x<1.0, and the final simulation time is t=0.14. Fig. 9 shows density and velocity distributions at the final time. The WENO3-NN schemes do not introduce any oscillations at the shock wave. WENO3-NN captures the contact discontinuity and shock very sharply. The results of the WENO3-NN schemes are better than the WENO5-JS scheme and are closer to the performance of WENO5-Z. The plateau value is exceptionally well captured by both WENO3-NN schemes. The WENO3-NN2 slightly outperforms WENO3-NN1. Due to reduced dissipation, a slight overshoot, similar in strength to WENO5-JS, in the velocity profile is visible for WENO3-NN2. Compared to WENO3-JS/Z/N, WENO3-NN needs considerably less points to resolve the discontinuities.
Fig. 9
  1. Download: Download high-res image (164KB)
  2. Download: Download full-size image

Fig. 9. Lax shock tube at t = 0.14. Resolution N = 200.

The initial condition of the 123 problem is(31)(ρ,u,p)={(1,−2,0.4)if0≤x<0.5,(1,2,0.4)if0.5≤x<1.0, and the final simulation is t=0.15. The 123 problem consists of two strong rarefaction waves and poses a demanding test for numerical schemes as the pressure of the intermediate state is close to zero. Fig. 10 provides a comparison between the established WENO schemes and WENO3-NN. WENO3-NN is able to handle the strong rarefaction-waves very well without introducing invalid pressure values. A detailed view of the head section of the rarefaction wave indicates that WENO3-NN clearly outperforms WENO3-N and WENO5-JS.
Fig. 10
  1. Download: Download high-res image (163KB)
  2. Download: Download full-size image

Fig. 10. 123 problem at t = 0.15. Resolution N = 200.

We provide the pointwise error distributions and integral L1 error norms for the shocktube tests in A.6.

4.3. Interacting blast waves

We consider the interacting blast wave test case from Woodward and Colella [20]. The initial condition is(32)(ρ,u,p)={(1,0,1000)if0≤x<0.1,(1,0,0.01)if0.1≤x<0.9,(1,0,100)if0.9≤x<1. The computational domain is [0,1] and reflective boundary conditions are applied at x=0 and x=1. The final simulation time is t=0.038. Fig. 11 shows the simulation results on a uniform mesh with N=400 (top) and N=800 (bottom) cells. From all three-point stencils, WENO3-NN provides the best approximation of the density profile. The difference is especially pronounced near the valley and at the right density peak. The quality of the WENO3-NN approximation at 400 points is similar to the WENO3-JS scheme at twice the resolution. Both WENO3-NN give very similar results over most of the domain, however WENO3-NN2 approximates the right density peak slightly better.
Fig. 11
  1. Download: Download high-res image (318KB)
  2. Download: Download full-size image

Fig. 11. Interacting blast waves at t = 0.038. Top: N = 400, bottom: N = 800.

4.4. Shock-density wave interaction

We consider the shock-density interaction test case by Shu and Osher [21]. The initial condition is a Mach 3 shock running into a perturbed density field,(33)(ρ,u,p)={(3.857,2.629,10.333)if0≤x<1,(1+0.2sin⁡(5x),0,1)if1≤x<10. The computational domain is [0,10]. The final simulation time is t=1.8. Fig. 12 shows the simulation results on a uniform mesh with N=400 (top) and N=800 (bottom) cells. WENO3-NN manages to resolve the density waves much better compared to the smeared solution of WENO3-JS/Z/N, indicating that smooth flow features are very well captured by WENO3-NN. The shock wave is well captured by all schemes. WENO5-Z provides the best performance among all investigated schemes. However, the WENO3-NN results are sharper compared to WENO3-JS and less points are needed to resolve the discontinuity. The solution of WENO3-NN at N=400 is comparable to the WENO3-JS solution at twice the resolution. For N=800, we observe the WENO3-NN result approaching WENO5-JS.
Fig. 12
  1. Download: Download high-res image (373KB)
  2. Download: Download full-size image

Fig. 12. Shock-density wave interaction at t = 1.8. Top: N = 400, bottom: N = 800.

4.5. Gresho vortex advection

We consider the unsteady Gresho vortex. A uniform flow with(34)(ρ0,u0,v0,p0)=(1,M0c0,0,1) is superposed with a rotating vortex of radius R placed at the center (0.5,0.5) of the computational domain [0,1]×[0,1]. We use a uniform mesh with a resolution of 96×96 cells. c0=γp0/ρ0 is the speed of sound and M0 is the Mach number. Periodic boundary conditions are applied. The numerical solution is calculated on a 96×96 grid and the final simulation time is t=1/u0=84.515. The initial conditions of the vortex are given in terms of the radial distance from the vortex core r=(x−0.5)2+(y−0.5)2. The tangential velocity uϕ and the pressure are given by(35)uϕ=u0{2r/Rif0≤r<R/2,2(1−r/R)ifR/2≤r<R,0ifR≤r, and(36)p=p0+u02{2r2/R2+2−ln⁡(16)if0≤r<R/2,2r2/R2−8r/R+4ln⁡(r/R)+6ifR/2≤r<R,0ifR≤r. We choose R=0.4 and M0=0.01. The Gresho vortex is used to test the low Mach number behavior of numerical schemes. Here, we are mainly interested in the dissipative properties of the reconstruction schemes. Fig. 13 shows the pressure field at the final simulation time. The relatively strong dissipation of WENO3-JS and WENO3-Z have led to considerable deformations of the vortex. WENO3-N retains the overall vortex shape and provides an acceptable approximation, while the WENO5-JS solution approximates the vortex very well apart from small spurious disturbances. For the WENO3-NN schemes, the vortex structure looks reasonable and is only slightly smeared out when compared to the WENO5-JS solution. Compared to WENO3-NN2, the vortex in the WENO3-NN1 solution is approximated sharper. At the same time, the solution features more spurious noise due to low dissipation.
Fig. 13
  1. Download: Download high-res image (441KB)
  2. Download: Download full-size image

Fig. 13. Gresho vortex at t = 84.515. Pressure contours are shown from lowest (blue) to highest (red). Resolution 96 × 96. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

4.6. Double Mach reflection of a strong shock

We consider the double Mach reflection test of a strong shock from Woodward and Colella [20]. A right-moving Mach 10 shock is reflected from a wall. The initial condition is(37)(ρ,u,v,p)={(1.4,0,0,1)ify≤1.732(x−0.1667),(8,7.145,−4.125,116.5)otherwise. Reflective boundary conditions are applied at the lower wall, at the top boundary the flow variables are prescribed as to follow the exact evolution of the moving shock. We use a mesh resolution of 1024×256 cells. The final simulation time is t=0.2. Fig. 14 shows density contours in a detailed view of the region [2,3]×[0,1].
Fig. 14
  1. Download: Download high-res image (363KB)
  2. Download: Download full-size image

Fig. 14. Double Mach reflection at t = 0.2. Density contours. Resolution 1024 × 256. The figure is drawn with 40 equally spaced density contours between 1.85 and 21.5.

Both WENO3-NN schemes are able to resolve considerably finer structures than the WENO3-JS/Z/N schemes and show an overall sharper density field. This becomes especially visible in the jet region and at the primary slip line which shows the development of wave structures. Due to excess dissipation, the slip line appears rather smooth for the classical WENO3 schemes. Here, WENO3-NN1 resolves ever so slightly finer structures than WENO3-NN2. The WENO3-NN schemes show a slightly more pronounced post shock instability behind the incident shock wave when compared to WENO3-N or WENO5-JS. The overall solution structure looks very much similar to the WENO5-JS result. However, the WENO3-NN schemes manage to suppress some of the spurious oscillations present in the WENO5-JS solution.

4.7. Rayleigh-Taylor instability

The initial condition of the Rayleigh-Taylor instability is(38)(ρ,u,v,p)={(2,0,−0.025acos⁡(8πx),1+2y)if0<y<0.5,(1,0,−0.025acos⁡(8πx),1.5+y)if0.5<y<1, with the speed of sound a=γp/ρ and the ratio of specific heats γ=5/3. The computational domain is [0,0.25]×[0,1]. Reflective boundary conditions are applied for left and right boundary. Dirichlet boundary conditions are used for top and bottom boundary, i.e. (ρ,u,v,p)=(1,0,0,2.5) and (ρ,u,v,p)=(2,0,0,1), respectively. The final simulation time is t=1.95.
Fig. 15 compares the solutions of all WENO schemes at a resolution of 128×512 corresponding to a grid spacing Δx=Δy=1/512. The WENO3-NN result has much finer structures than WENO3-JS and WENO3-Z. The slightly more spherical shapes in the WENO3-NN results indicate low dissipation. Between the two WENO3-NN schemes, WENO3-NN1 resolves very sharp and fine flow structures.
Fig. 15
  1. Download: Download high-res image (484KB)
  2. Download: Download full-size image

Fig. 15. Rayleigh-Taylor instability at t = 1.95. 20 equally spaced density contours between 0.85 and 2.25 are drawn. Resolution 128 × 512.

Fig. 16 compares the numerical solutions at a higher resolution of 256×1024 corresponding to a grid spacing Δx=Δy=1/1024. The vortical structures starting to form in the WENO3-NN solution are considerably smaller than the WENO3-JS/Z/N ones. While the WENO3-NN2 solution has qualitative similarity with the WENO3-Z solution, the density contours of WENO3-NN1 display very fine and detailed flow features. Interestingly, WENO3-NN1 has less prominent vortical structures compared to WENO3-NN2, and WENO5-JS does not show vortical structures at all at the current resolution level.
Fig. 16
  1. Download: Download high-res image (455KB)
  2. Download: Download full-size image

Fig. 16. Rayleigh-Taylor instability at t = 1.95. 20 equally spaced density contours between 0.85 and 2.25 are drawn. Resolution 256 × 1024.

4.8. Approximate dispersion relation

The dissipation-dispersion relation provides further insight into the properties of nonlinear shock-capturing schemes. Following Pirozzoli's approximate dispersion relation [22], we analyze the properties of the WENO3-NN schemes and compare them with the background linear central upwind-schemes of corresponding order (CU1, CU3, CU5) as well as with third and fifth-order WENO-JS methods. Fig. 17 shows the ADR for WENO3-NN1 and -NN2. ξ denotes the wavenumber. ξ˜I and ξ˜R are the imaginary and real part of the modified wavenumber, respectively. Although the performance of both WENO3-NN schemes in above tests was very similar, the dispersion relations of both schemes differ quite drastically. The WENO3-NN2 scheme shows a more classical dissipation-dispersion behavior. Over the whole wave number range WENO3-NN2 is considerably less dissipative compared to WENO3-JS. Interestingly, WENO3-NN2 is even slightly less dissipative than the WENO5-JS scheme for higher wave numbers. The corresponding dissipation curve has a monotonic behavior and does not show the distinct valleys of WENO3-JS or WENO5-JS. The dispersive behavior of WENO3-NN2 shows qualitative similarities to the linear central-upwind schemes. The ADR of WENO3-NN1 shows a very intersting behavior. The dissipation is kept to a minimum over the whole wavenumber range. Especially, at the cutoff wavenumber zero dissipation is introduced. Here, WENO3-NN1 strikes a complex balance between dissipation and dispersion. We interpret these results such that underresolved wave packets near the cutoff wavenumber are maintained by the WENO3-NN1 scheme, but do not interact with the resolved scales, akin to soliton solutions. These findings could explain the good numerical results of WENO3-NN1 for the Gresho vortex test case and the much sharper flow structures in the Double Mach reflection test. While WENO3-NN1 introduces near zero dissipation at the cutoff wavenumber, this distinct behavior might at the same time impede the development of flow structures near the cutoff wavenumber, e.g. the vortical structures in the Rayleigh-Taylor test case.
Fig. 17
  1. Download: Download high-res image (229KB)
  2. Download: Download full-size image

Fig. 17. Approximate dispersion relations for selected schemes.

4.9. Online computational cost

Despite the fact that WENO3-NN predominantly is intended to explore machine learning supported design of nonlinear discretization schemes we assess in this section its computational performance in comparison with established WENO schemes. For a realistic estimation of the online computational cost, we evaluate the computational efficiency of the proposed WENO3-NN schemes for the double Mach reflection test case from Sec. 4.6. We measure the average computational time required for advancing the state of one cell for one time step. Since we use the 3rd-order TVD Runge Kutta scheme for time integration, one time step corresponds to three flux evaluations and therefore to three cell face reconstructions. The simulations are performed at a resolution of 1024×256 points and are integrated over 3700 time steps. The algorithm is implemented in Jax [23]. We use a Nvidia GTX2080 GPU for all simulations. Table 2 shows the absolute and normalized average wall clock times per cell for a computational step.

Table 2. Average computational performance of selected WENO schemes.

SchemeGPU time per cell and time step (in ns)Normalized cost
WENO3-JS3891.00
WENO3-Z3961.02
WENO3-N4111.06
WENO3-NN116764.31
WENO3-NN217434.48
WENO5-JS5611.44
Neural network evaluations at each stencil render WENO3-NN as clearly more expensive than classical schemes with analytically derived smoothness measures. However, we emphasize that neural network evaluations have not been optimized for computational efficiency nor have we optimized the size of the network. Keeping in mind the strongly improved results of the WENO3-NN schemes, we believe that with the rapid development of dedicated computational hardware for fast neural network evaluations (e.g. tensor processing units), data-driven reconstruction schemes will become a viable alternative to classical approaches.

5. Conclusion

Inspired by the recent success of machine learning enhanced numerical methods for computational fluid dynamics, we have proposed a neural-network-weighted essentially-non-oscillatory scheme. The characteristics of the WENO3-NN scheme have been extensively analyzed and a plethora of one- and two-dimensional benchmark test cases, involving strong shocks and shock-density interactions, has been performed. We demonstrate that data-driven training can lead us to new schemes with unexpected properties and significantly improved performance.
The WENO3-NN scheme utilizes a neural network as a parameterizable weighting function of low-order local interpolation polynomials. We have shown that embedding a priori knowledge, such as a Galilean invariant input layer and the reconstruction via Harten polynomials, into the design of the WENO3-NN scheme yields an interpretable and generalizable data-driven scheme. A newly introduced loss term on the network weights allows the WENO3-NN scheme to reach full-order convergence whereas earlier data-driven WENO schemes required a rescaling of the network output in order to guarantee convergence. The WENO3-NN scheme, although trained on a handful of one-dimensional elementary functions at a single resolution level, has shown very good performance at different resolution levels and in multidimensional problems. In fact, the WENO3-NN scheme outperforms the classical WENO3-JS/Z/N schemes across all benchmark problems and, in some test cases, has shown similarities to the performance of WENO5-JS. The approximate dispersion relation of the WENO3-NN schemes has revealed a very interesting behavior. While the WENO3-NN2 scheme features a rather classical, low-dissipation behavior over the whole wavenumber range, the ADR of the WENO3-NN1 scheme reveals a rather complex and unexpected dispersion and dissipation behavior. Especially, the WENO3-NN1 scheme introduces zero dissipation at the wavenumber cutoff. Machine learning might provide the necessary tools to further investigate such non-trivial dispersion-dissipation relations of numerical schemes. More generally, data-driven ansatze might have the potential to push the boundaries of conventional and long-standing concepts in the design of numerical schemes for fluid mechanics. The combination of such an analysis with an extension of the proposed WENO-NN methodology to higher order schemes seems like a promising avenue for the development of data-driven numerical schemes and is the subject of ongoing research. Finally, we want to point out that the high interpretability of the WENO3-NN weighting function may help to find analytical smoothness measures/weighting functions for higher order WENO schemes. The proposed framework could enable a new line of development of WENO schemes in which firstly a deep WENO-NN scheme is trained on a rich dataset and, then, a functional relation of the learned weighting scheme is identified (e.g. using sparse regression). Aforementioned issues motivate future work.

CRediT authorship contribution statement

Deniz A. Bezgin: Conceptualization, Formal analysis, Investigation, Methodology, Software, Writing – original draft. Steffen J. Schmidt: Conceptualization, Supervision, Writing – review & editing. Nikolaus A. Adams: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 667483).

Appendix A.

A.1. Dataset

One has to take special care when generating finite-volume data (i.e. cell-averaged data) on a three-point stencil in order to avoid contradicting training samples. Especially, jump discontinuities and extreme points can introduce inconsistencies into the dataset. For example, consider the finite-volume representation of a jump discontinuity in the middle of cell i, see Fig. 18 (left). In the cell-averaged representation the discontinuous example and the smooth linear example look identical. However, the corresponding cell face values are different. The network, therefore, would get two different labels for the same input. We avoid this predicament in training samples by placing jumps on cell faces only.
Fig. 18
  1. Download: Download high-res image (97KB)
  2. Download: Download full-size image

Fig. 18. Exemplary three-point stencil in which the finite-volume representation (green points) of a discontinuous function (orange line) and a smooth function (blue line) becomes identical.

The same holds true for extreme points (e.g. of a sine function) and jump discontinuities, see Fig. 18 (right). The finite-volume representation of the sine function is identical with a jump placed at xi−1/2. To this end, for trigonometric functions, we exclude all stencils from the training dataset for which cell face values lie outside the cell averages in the corresponding stencil.

A.2. Model training

Table 3 shows the value ranges of the model hyperparameters. Each model is initialized with 5 sets of different random weights. For each set of model hyperparameters, we choose the model with lowest validation error that also passes all 1D benchmark problems. Fig. 19 shows the training and validation loss during the optimization process. The losses are fully converged at the end of the training process.

Table 3. Model hyperparameters.

ParameterValue range
α[0, 1E-2, 3E-2, 1E-1, 3E-1, 1]
βd[1E-2, 1E-1, 1]
βW1E-9
ceno2E-4
Fig. 19
  1. Download: Download high-res image (140KB)
  2. Download: Download full-size image

Fig. 19. Loss history for WENO3-NN1 (left) and WENO3-NN2 (right).

A.3. Influence of hyperparameters

The analysis of the hyperparameters on the model performance is very instructive. We evaluate the weight ω0 for sin3⁡(πx) while keeping either α or βd fixed and varying the other parameter in the loss function Eq. (17). The left part of Fig. 20 shows that for fixed βd=0.1 decreasing α provides the WENO3-NN model with more liberty to deviate further from the ideal weights dk. As mentioned earlier, α introduces a scale separation mechanism. α values closer to one smoothen out the input field and push the WENO3-NN outputs closer to the ideal weights over the whole domain. For α values closer to zero, network outputs deviate stronger from the ideal dk. The right part of Fig. 20 shows that increasing βd for fixed α pushes the network output more towards the ideal weights. For βd=1, ω0 deviates from the ideal weight only at the smooth extreme points. Both behaviors are in line with our understanding of the loss function Eq. (17).
Fig. 20
  1. Download: Download high-res image (141KB)
  2. Download: Download full-size image

Fig. 20. Influence of the hyperparameters on the NN weight function for sin3⁡(πx). Left: βd = 0.1, right: α = 0.1.

A.4. Cutoff threshold of the ENO layer

WENO3-NN schemes detect discontinuities in the flow field very well and adjust the reconstruction accordingly. However, during numerical experimentation we found that WENO3-NN schemes without an ENO layer assign small non-zero weights to discontinuous stencils. For example, the left of Fig. 21 shows the output weights of WENO3-NN1 without ENO layer at the cell face xi+1/2 evaluated for a single jump discontinuity placed at x=1.0, i.e.(A.1)u(x)={0if0≤x≤1,1if1<x≤2. We observe that WENO3-NN assigns a small, but non-zero weight of around 10−4 to the discontinuous stencils. During training the neural network learns to detect discontinuities and adapts the reconstruction accordingly. However, the neural network does not decrease the output weights below a certain threshold (here 10−4) since the loss at such discontinuous stencils is already small compared to other training samples. Additionally, the saturation of the softmax output activation makes it difficult to output exact zeros and ones. In practice, this does not pose a problem for many flow applications. However in the vicinity of very strong shocks, e.g. the interacting blast wave test by Woodward and Colella in Sec. 4.3, ωk≈10−4 might already lead to the reconstruction of negative pressures or densities. To prevent the reconstruction of such inadmissible states, we postprocess the network output at test time. We pass the output weights ωkNN through a simple cutoff function, the so called ENO layer, which restores the ENO property.(A.2)ω˜kNN=ωkNNψk∑jωjNNψjwithψj={0,ifωjNN<ceno,1,otherwise. As ωkNN≈10−4 for discontinuities, we choose ceno=2⋅10−4 as the threshold value for the ENO layer. Note, that the cutoff function can be easily implemented via a ReLU activation function with corresponding threshold, i.e.(A.3)ω˜kNN=ReLU(ωkNN−ceno)∑jReLU(ωjNN−ceno). The weights ωkNN with active ENO layer are visualized in the right part of Fig. 21. The ENO layer only affects the numerical solution near discontinuities.
Fig. 21
  1. Download: Download high-res image (143KB)
  2. Download: Download full-size image

Fig. 21. Evaluation of the output weights of WENO3-NN for Eq. (A.1) with and without ENO layer. Left: pure network output (without ENO layer), right: post-processed network output (with ENO layer).

A.5. Convergence for the linear advection equation

We provide further details on the convergence behavior of the WENO3-NN scheme for the linear advection equation. We apply the WENO3-NN scheme to the linear advection of u0(x)=sin3⁡(πx) on the domain [−1,1]. The initial condition is integrated up to t=2.0. We choose a small enough time step (CFL=0.05) to exclude any errors by the time stepping scheme. Fig. 22 shows the grid convergence in L1, L2, and L∞ errors for WENO3-JS/Z/N/NN1/NN2 and WENO5-JS/Z, respectively. WENO3-NN1 and -NN2 achieve lower absolute errors than WENO3-JS/Z/N and also show better convergence behavior.
Fig. 22
  1. Download: Download high-res image (325KB)
  2. Download: Download full-size image

Fig. 22. L1,L2, and L∞ errors of WENO3-JS/Z/N/NN1/NN2 and WENO5-JS/Z for the linear advection equation with u0(x)=sin3⁡(πx).

A.6. Pointwise error for linear advection and shocktube tests

In the following, we provide the pointwise errors for the linear advection test cases from Sec. 4.1 and the shocktube tests from Sec. 4.2. Fig. 23 shows the error distributions for linear advection of the multiwave Eq. (27), the discontinuity Eq. (28), and sin3⁡(πx). The error distributions show that WENO3-NN schemes consistently give a lower error level near discontinuities and in smooth regions when compared to WENO3-JS/Z/N. The multiwave test case indicates that initial discontinuities are smeared out similarly for WENO3-NN and WENO5-JS. The error level of WENO3-NN2 is slightly lower than WENO3-NN1 for all linear advection test cases.
Fig. 23
  1. Download: Download high-res image (519KB)
  2. Download: Download full-size image

Fig. 23. Pointwise error distributions for the linear advection of the multiwave Eq. (27) (top left), the discontinuity Eq. (28) (top right), and of sin3⁡(πx) (bottom center) at t = 10.0.

We also provide the pointwise density and velocity error plots for the Sod, Lax, and 123 shocktube problems, see Eqs. (29), (30), and (31). We compare the WENO approximations to the cell-averaged exact solution. Fig. 24 shows the error distributions for the Sod problem. For the density, WENO3-NN1 and -NN2 show considerably lower error levels than other three-point schemes. Around the head and tail of the rarefaction wave and also around the discontinuities, WENO3-NN schemes even outperform WENO5-JS. The error distributions for the Lax problem in Fig. 25 underline the improved shock-capturing capabilities of the WENO3-NN schemes. Near the contact discontinuity and the shock, WENO3-NN schemes outperform WENO3-JS/Z/N and WENO5-JS. Here, the error levels of WENO3-NN schemes are much closer to WENO5-Z. The 123 problem indicates that WENO3-NN schemes are able to resolve strong rarefactions very well, see Fig. 26. WENO3-NN schemes have the lowest overall error from all three-point stencil schemes. However, WENO3-NN schemes show a slightly larger error in the region between the two rarefaction waves. We note that at times, the pointwise errors plots are difficult to interpret. For example in the Sod or 123 tests, WENO3-JS shows lower error magnitudes than WENO5-JS in some regions but delivers much worse approximations in other regions. Therfore, we summarize the integral L1 errors in Table 4. The L1 error for the Euler equations is calculated as the sum of the L1 errors of density, velocity, and pressure, i.e. L1=L1ρ+L1u+L1p. In the linear advection tests, WENO3-NN schemes clearly outperform WENO3-JS/Z/N and are much closer to the WENO5-JS results. WENO5-Z yields the best performance across all linear advection tests. For the nonlinear shocktube problems, WENO3-NN schemes consistently deliver lower error magnitudes than WENO3-JS/Z/N and even outperform WENO5-JS. For the Sod and Lax problems, WENO5-Z achieves the lowest error magnitudes among all the WENO schemes considered in this work. However in the 123 problem, WENO3-NN schemes deliver better approximations than WENO5-Z.
Fig. 24
  1. Download: Download high-res image (254KB)
  2. Download: Download full-size image

Fig. 24. Absolute pointwise error distribution of density and velocity for the Sod test case defined by Eq. (29).

Fig. 25
  1. Download: Download high-res image (285KB)
  2. Download: Download full-size image

Fig. 25. Absolute pointwise error distribution of density and velocity for the Lax test case defined by Eq. (30).

Fig. 26
  1. Download: Download high-res image (297KB)
  2. Download: Download full-size image

Fig. 26. Absolute pointwise error distribution of density and velocity for the 123 test case defined by Eq. (31).

Table 4. L1 errors of various WENO schemes for the linear advection equation and shocktube problems of the Euler equation. Linear advection tests include sin3⁡(πx), the multiwave in Eq. (27), and the discontinuity in Eq. (28). The Sod, Lax, and 123 shocktube problems are defined by Eqs. (29), (30), and (31).

Empty Cellsin3⁡(πx)GSTEDisc.SodLax123
WENO3-JS1.180E-13.702E-18.260E-21.636E-23.949E-22.540E-2
WENO3-Z5.552E-22.343E-15.902E-21.272E-22.931E-22.317E-2
WENO3-N4.197E-21.969E-15.225E-21.188E-23.153E-22.255E-2
WENO3-NN12.661E-21.651E-14.355E-21.078E-22.527E-21.973E-2
WENO3-NN22.445E-21.609E-14.235E-21.035E-22.416E-21.967E-2
WENO5-JS2.572E-39.964E-22.634E-21.106E-22.590E-22.393E-2
WENO5-Z1.930E-37.731E-22.335E-28.113E-31.900E-22.282E- 2

References

Cited by (12)

  • JAX-Fluids 2.0: Towards HPC for differentiable CFD of compressible two-phase flows

    2025, Computer Physics Communications

  • A deep reinforcement learning framework for dynamic optimization of numerical schemes for compressible flow simulations

    2023, Journal of Computational Physics
    Citation Excerpt :

    Data-driven methods based on artificial neural networks (ANNs) or Gaussian processes (GPs) gain significance in CFD due to their ability to learn complex, nonlinear relation from data. They have been integrated with the numerical solvers to produce physically consistent solutions [20–26]. Due to different characteristics of subgrid scales, traditional optimization strategies, such as those based on spectral properties [7] or based on evolutionary [15] as well as Bayesian optimization [24,25], may not be effective in finding the optimal solutions.

  • JAX-Fluids: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows

    2023, Computer Physics Communications
    Citation Excerpt :

    Upon proper training, they are then plugged into an existing CFD solver for evaluation of down-stream tasks. Examples include training of explicit subgrid scale models in large eddy simulations [33], interface reconstruction in multiphase flows [34,35], and cell face reconstruction in shock-capturing schemes [36,37]. Although the offline training of ML models is relatively easy, there are several drawbacks to this approach.

  • Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations

    2024, Nature Machine Intelligence

  • Review of the High-Order TENO Schemes for Compressible Gas Dynamics and Turbulence

    2023, Archives of Computational Methods in Engineering

  • BAYESIAN OPTIMIZATION ON FIFTH-ORDER TARGETED ENO SCHEME FOR COMPRESSIBLE FLOWS

    2022, WCCM-APCOM 2022 - 15th World Congress on Computational Mechanics and 8th Asian Pacific Congress on Computational Mechanics: Pursuing the Infinite Potential of Computational Mechanics
View all citing articles on Scopus
View Abstract
© 2021 Elsevier Inc. All rights reserved.