A Weno-Tvd Implementation for Solving Some Problems of Hyperbolic Conservation Laws

This work deals with a numerical implementation of a fifth order CENTRAL WENO-TVD (Weighted Essentially Non-Oscillatory-Total Variation Dimimishing) scheme [50] applied to the convective terms of some hyperbolic conservation laws problems, in a volume finite framework. The WENO-TVD scheme is used to solve the 1D advection and Burgers equations. For this case is implemented two different numerical fluxes: The Lax-Friedrichs and TVD fluxes. In the TVD fluxes the schemes applied are in flux-limiter form. The schemes implemented for this flux are: Van Albada-1 [44], van Albada-2 [22], van Leer [18] and MINMOD [19]. The WENO type schemes are characterized for their high order approximation, and do not produce spurious oscilations near discontinuities, shocks and higher gradients. A third order Runge-Kutta TVD [13] for the temporal variable is used. Qualitative and quantitative comparison are presented. The numerical solutions are computed with an in-house computer code developed in MATLAB software. In future works, it will develope a paralelization of computer code for solving systems of conservation laws, e.g. Euler equations of gas dynamics.


Introduction
The conservation laws are mathematical expressions of a basic principle that allows to describe the temporal evolution of an amount of interest (temperature, pressure of a fluid, among others).They refer to the physical laws that postulate that during the temporal evolution of an isolated system, where certain magnitudes have a constant value, allowing to preserve the properties of a certain quantity, which can be mass, energy or moment.. Particularly, the hyperbolic laws of conservation describe several physical problems in diverse dynamic areas such as of fluids, astrophysics, shallow water equations, prediction of the climate, compressible gas dymamics, among others.In the face of the need to solve problems in these areas and in turn the lack and/or difficulty of analytical processes for its solution, the numerical methods are used in order to obtain a true approximation to the physical phenomenon.
In order to approximate the conservation laws the finite volume method is used for fluid dynamics problems, with the purpose of studying the conservation of mass u, of some chemical present in a fluid transported by a tube or pipette, with x distance to the Length of the pipette, for a time t.The MVF subdivide the domain into very small finite parts, called cells or control volumes.This method is based on the fact that many physical laws are conservation laws, suggesting that the input flow of an interest amount, in a cell, is identical to the flow that leaves the adjacent cell.Following this idea, we proceed to formulate the ruling equations in flow conservation equations, defined in an integral way using the averages of cells [9].
To find the solution, we implement the high resolution method WENO (Weighted Essentially Non-Oscillatory) in the convective part, which is based on an interpolation of the average points of cells, standing out for its ability to achieve a high order of precision in smooth regions, it does not produce oscillatory solutions, maintains the shape in the transition due to contact discontinuities and guarantees convergence.
The first WENO scheme was developed by Liu, Chan and Osher in 1994 [29], which was a finite volume version of third order in a single dimension.In 1996, finite difference schemes of third and fifth order in multispace dimensions were constructed by Jiang and Shu in [21], introducing smoothness indicators and non-linear weights.Then, Shu in 1997, developed [36], ENO schemes (Essentially Non-Oscillatory) and WENO for hyperbolic conservation laws, in 1997.A year later, together with Hu, they developed a scheme WENO for triangular meshes [38].Harten in [16], introduced the notion of TVD schemes (Total Variation Disminishing) of second order.Later, Titarev and Toro, in 2005, designed [42], an ENO and WENO scheme of second order based on TVD upwind and central fluxs, which is then improved in [49], when using a third-order TVD flow, constructing a fifth-order WENO scheme, developed by Yousef Hashem Zahran in 2006.This work is framed in the area of numerical analysis and fluid mechanics, as it is intended to undertake a study of the WENO/WENO-TVD schemes of Fifth Order, by implementing a computational code in MATLAB, for the convective part of Hyperbolic conservation laws for the one-dimensional case.For this purpose, the analytical solutions of the working equations, such as the advection and Burgers equations (non-viscous), are compared qualitatively and quantitatively with the numerical solution obtained with WENO/WENO-TVD using the Lax-Friedrichs numerical flux and TVD of third order presented in [49], as well as flow limiters like van Leer [45], van Albada 1 [44], Minmod [32] and van Albada 2 [22].In order to achieve high-order precision in the temporal discretization, the third-order Runge-Kutta TVD method is used.This document is organized as follows.In section 2, conservation laws and work equations are studied.In section 3, the numerical methodology of finite volumes, and the types of numerical schemes are discussed in detail.In section 4, the WENO method is explained, based on the ENO reconstruction.Additionally, in section 5 presents the numerical flux schemes, the TVD property and the flux limiter schemes.Time discretization is treated with the Runge-Kutta method in section 6.Finally, in section 7 presents the results in tables and grahps that verify the efficiency of the method.

Results
The finite volume method is implemented with a fifth order WENO reconstruction, for the spatial variable, while we use the third order Runge-Kutta TVD method for the time discretization of the problem.
The problems to solve are defined in a spatial domain Ω ⊂ R, considering periodic boundary conditions (PBC) in the given interval, for a determined final time T. It denotes N as the number of partitions in the size space ∆x, M is the number of time partitions with step ∆t, and c = a ∆x ∆t is the Courant-Friedrichs-Lewy number (CFL), with a wave speed.
In addition, the flux limiter is denoted by TVD-3.
In this problem, the rate of convergence is checked for long times.In table 1 the results obtained from problem 2.1 are presented.In it, it can be seen that the WENO-5-LF and WENO-5-TVD schemes reach an order of precision of the fourth order, in the standards L 1 and L ∞ , even after a long integration time.In addition, it is noted that a higher order is achieved with the first-order flux of Lax-Friedrichs than with the third-order TVD flux, and that better results are obtained with the norm L ∞ .However, both schemes maintain a minimum of third order.
In the figure 2 it is shown how the numerical solution adjusts to the exact solution, to the point of not distinguishing itself with the naked eye, due to the small errors obtained in the numerical results.Precisely, an extension of this graph is presented in the figure 2, in which the solutions are seen at the point x ≈ 0.5, allowing to observe that the WENO-5-LF is closer to the solution Exactly the WENO-5-TVD scheme, corroborating that the Lax-Friedrichs flux works better than the TVD flux.
It is important to note that in order to obtain a "visible difference " in the graphic comparisons of the schemes it is necessary to make enlargements in the graphs considering a scale of up to 10 −5 , ratifying the closeness of the solutions, as indeed is shown in the figure 2.  On the other hand, the graph 3 compares the flux limiters of van Leer, van Albada 1, van Albada 2, and Minmod together with the third-order TVD flux used in the WENO-5-TVD scheme.Regardless of whether the limiter satisfies the TVD properties or symmetry, good results are obtained when solving hyperbolic problems with smooth initial condition.An extension of the graph 3 is presented in the figure 4 at the point x ≈ 0.5, in which it can be seen that the flux limiter that works best is van Albada 2, and that in addition the second order limiters present a very similar precision, but still better than the third order TVD limiter.
Precisely, a remarkable "difference " between the flux limiters considered can be seen in the figure 5, where the solution at the point x ≈ 0 is taken as reference, noting that The best resolution flux limiter is the van Albada 2, followed by van Leer, van Albada 1, Minmod and finally the TVD-3.
In this problem the initial condition is modified by a function with greater physical oscillations, in such a way that possible deteriorations of precision can be detected due to the strong oscillations in the parameters that determine the grid.According to the results presented in table ??, it is seen that more irregular orders are obtained, since orders of precision are obtained higher than those of problem 1, but at the same time it is also possible to find second order for some values of N, obtaining errors up to E-7 when in problem 1 the errors were up to E-8.With respect to the numerical flows, it can be said that there is a subtle improvement in favor of the Lax-Friedrichs flux, even though the TVD flux achieves fifth order of precision for N = 160 with the standard L ∞ .As can be seen in the figure 6, there is no deterioration in the accuracy of the solution despite strong oscillations due to the change in the initial condition.It is clear to see, that the numerical solution fits very well to the exact solution both in the valleys and in the crests of the wave.

Preprints
For instance, at the point x ≈ −0.5,The third-order TVD scheme works better, because the approximate value is closer to the exact value, as shown in the figure 7.
This problem, better known as the Zalesak problem [51], is used to observe the resolution properties of the scheme, such as the control of oscillations near discontinuities [43].In the figures 8 and 9 above it can be seen that the solution has no spurious oscillations, since the scheme controls oscillations near discontinuities.Even when the initial soft condition has been changed, by a function in sections, the resolution of the WENO and WENO-TVD scheme remains faithful to the exact solution of the problem, even when it is an initial condition with great variety of shape in its sections.In addition, you can see how the solution improves by taking N = 400 instead of N = 200, especially in the non-smooth parts such as the "peaks "and ""bars.In the first "peaks"of the solution, reached in x ≈ −0.7, a better approximation with the TVD flux is achieved than with the Lax-Friedrichs flux, as can be seen in the figure 10.
An analogous result is obtained at the point x ≈ −0.3, where a high order of precision is achieved using the TVD flux, in such a way that the numerical solution remains faithful to the exact solution considering a very small scale in the expansion of the graph in the figure 9 (see figure 11  In the figure 12, we can see once again how the WENO schema works well close to discontinuities and "peaks ".In this case, we observe the solutions at the point x ≈ 0.1, being able to conclude that the flow of Lax-Friedrichs continues to be less precise than the TVD flux, when it comes to this type of graphs.Finally, when observing the final section of the graph 9 presented in the figure 13, it can be said that it is very similar to the wave crests of the figure 1 of problem 1, and as in that case, the monotonous flux of Lax-Friedrichs presents better results in this region of the graph, as occurs for example at the point x ≈ 0.5.
T = 0.33 and T = 1.5, where f (u) = 1 2 u 2 .This problem is used to show the resolution capability of the schemes in the presence of a shock.Two integration times are considered, T = 0.33 to show the solution before the shock and T = 1.5 to illustrate the solution after the shock which occurs at T = 0.5.Clearly, the WENO/WENO-TVD schemes perform efficiently for nonlinear equations, as can be seen in the figure 14, in which it is observed that the schemes continue maintaining high precision (before the shock) in regions smooth as in the linear problems previously considered.In the figures 15 and 16 there are enlargements of the graph (14), which correspond to the solutions (exact and numerical) at the points x ≈ −0.5 and x ≈ 0, respectively.At both points, it is noted that the Lax-Friedrichs flux works better than the TVD flux, because its graph is closer to the exact solution.

Preprints
On the other hand, the figure 17 shows the resolution capacity of the WENO schemes after the shock, allowing observing how the numerical solutions are adjusted to the exact solution in regions with high precision in soft regions (see figure 18) and handling "small" precision errors in the "peaks".This is the case, of the intervals [0.1, 0.5] and [0.5, 0.8] in which we can see how the WENO scheme yields better results than the WENO-TVD scheme in the smooth regions of the solution.Otherwise, in the "peaks" and non-smooth regions of the solution, where the WENO-TVD scheme shows better resolution for the nonlinear problem after the shock, as shown in the figures 20.If now, the flux limiters considered above are compared, according to the figure 21 it can be said that the TVD-3 limiter is the best resolution, because it fits better to the exact solution of the problem, unlike the other limiters that tend to produce oscillations and move away from the exact solution in the "peaks" of it.Finally, of the second order flux limiters, it can be said that the Minmod works better than the others, except at the point x ≈ 0.45 in which the solution presents imprecision with respect to the exact solution.In addition, you can see that the lowest performance is obtained with the van Albada 2 limiter, followed by van Leer and van Albada 1.

Formatting of Mathematical Components
A conservation law is the mathematical expression of a basic principle that allows describing the time evolution of a certain amount of interest u (temperature, pressure of a fluid or the concentration of a chemical, among others) [47].
According [27], a conservation law can be expressed in abbreviated differential form as: where f is a function of u and its derivatives, and it's called the conserved flow of the variables.This project focuses on hyperbolic PDEs and therefore on problems modeled by this type of equations, called hyperbolic problems, particularly working with Hyperbolic Conservation Laws (HCL).

Work equations 1. Advection linear equation
The advection or transport equation describes the advection of the scalar field u(x, t) transported by a constant velocity flow a. Examples of this type of equation are the transport of an atmospheric property, such as humidity, by the effect of wind or the variation of the concentration of a chemical species within a fluid-dynamic current [7].The equation is given by where u is the advection rate and a is the wave propagation speed [23].

Non-viscous Burgers equation
This equation represents the propagation of weakly non-linear waves in which we can consider one-dimensionality in the spatial variable [7].The best simplification of the Navier-Stokes equations that preserves the characteristic of having shock waves is given by the hyperbolic EDP of non-viscous Burgers, which is given by: In this equation the speed of propagation depends on the solution itself u [23].

Finite Volume Method (FVM)
The finite volume method is a numerical method developed by Patankar and Spalding (1972), which allows the resolution of conservation laws.The FVM subdivide the domain into very small finite parts, called cells or control volumes.This method is based on the fact that many physical laws are conservation laws, suggesting that the input flow of an interest amount, in a cell, is identical to the flow that leaves the adjacent cell.Following this idea, we proceed to formulate the governing equations in flow conservation equations, defined using the integral form of a conservation law using the cell averages [9].
where f : R → R is a C 2 function called the physical flow of the problem.Then, each numeral is expanded.

Discretize of problem domain
Let a spatial discretization step ∆x > 0 and a time discretization ∆t > 0, the sequences are defined Let Ω ⊂ R the continuous domain of the semi-discrete problem (4).The process of discretization of the domain, consists of decomposing Ω in a finite number N (discrete) of subdomains I i (i = 1, dots, N), called control volumes (CVs).This is, For the one-dimensional case as shown in the figure 23, the CVs ] are subintervals of the problem interval, where the edge points of the subintervals xI ± 1 2 are called nodos.In this work, uniform grid are considered, that is, grid in which the size of the cell ∆x i = (x i+ 1 2 ) is equal to the spatial discretization step ∆x = max 1≤i≤N ∆x i , this is, ∆x i = ∆x.

Formulation of integral balance equations for each control volume.
For the mean value for integrals, the cell averages are approximate as follows, where ∆x = x i+1/2 − x i−1/2 is the length of the cell, considering a uniform grid for simplicity [27].3. Approximation of integrals by numerical integration.
It must be ensured that the numerical method is conservative, thus tracing the natural phenomenon [27].Then ∑ n i=1 U n i ∆x approximates the integral (5) over the entire interval.If the conservative form is used, the flows vary at the extremes of the interval (see figure 24), and so by the integral form of the conservation law [27], one has to

Approximation of function values and derivatives by interpolation with nodal values.
Let the cell averages U n i , in a time t n , we want to approximate U n+1 i in the next time t n+1 , for a time step ∆t = t n+1 − t n .Then, integrating (6) from t n to t n+1 have Dividing between ∆x and clearing we have 5. Assembling and solution of discrete algebraic system.
Since (5), the discrete scheme conservatively remains where F n i±1/2 is the approximation to the average flow in the points x = x i 1/2 , respectively, with: where F is a numerical flux function (see [27]) that satisfies the following conditions: • F is a continuous Lipschitz function in all its variables; • F is consistent with the physical flow f , this is, F (u, . . ., u) = f (u).
For example, the Lax-Friedrichs monotonous flux in the interval [a, b] given by where σ = max u | f (u)| is a constant [36].

WENO schemes
WENO schemes are high order schemes to approximate the convective part of hyperbolic conservation laws or other partial equations.They are based on high-order reconstruction in order to achieve high precision, not produce oscillatory solutions, preserve the transition form due to contact discontinuities and guarantee convergence [39].
Next, the WENO method is described for the one-dimensional case.

WENO 1D reconstruction procedure
Consider the one-dimensional hyperbolic equation ( 1) with smooth initial condition.Assuming that the grid is uniform, then solving (1) using a conservative approach to the spatial derivative, we have to where is the numerical flux and U i (t) are the cell averages in the points (x i , t), which are given by A reconstruction procedure consists of obtaining an approximation, for example of a value u(x i+ 1 2 ), by means of an interpolation.Initially, consecutive cells should be selected near x i+ 1 2 , which should include I i o I i+1 .The collection of these i-cells is called a reconstruction stencil.Then, the task is summarized in finding a single polynomial p(x), at most k − 1 degree, such as the cell average over each cell I i in the stencil, coincides with the averages of given cells U i , in the points x i in the "small" stencil for r = 0, 1, . . ., k − 1 with r + s + 1 = k, where r, s ≥ 0 are the number of cells to the left and right, respectively, used in the construction of S r .Proposition 1 ([36]).Given the cell averages of a function u(x) : Then there is a unique plynomial p r (x), of degree at most k − 1, for each cell I i , such that it is a k order accurate approximation to the function u(x) inside I i : ), When considering the problem (4) with a smooth initial condition u 0 (x) = u(x 0 ), and U i the cell averages of the smooth function u(x) for each cell I i = [x i−1/2 , x i+1/2 ] for a given order of precision k, the WENO reconstruction is based on the ENO interpolation process, which consists of interpolating the function u, choosing the softer base of k bases considered as defined in (12).Then, a WENO approximation of order 2k − 1 to the function u(x) in the bounded cells (u(x i±1/2 )), is obtained according to [36], following the next steps: of k order of precision, based in the grids S r for r = 0, . . ., k − 1.
Clearly by the proposition (2.1), there exists a unique polynomial p r (x) of degree at the most k − 1, which interpolates the function u(x), this is, p r (x j ) = u j , over the points x j in the "small" stencil S r = {x i−r , x i−r+1 , . . ., x i+s } for r = 0, 1, . . .k − 1, with r + s + 1 = k.Soon, using u (r) i+1/2 = p r (x i+1/2 ) like a approximation to the value u(x i+1/2 ), y u and also because smoothness of u on S r , we have that: Now, since the U j are linear in the base S r , then they are also linear in the points u i+1/2 and u i−1/2 , therefore there are constants C rj and Crj , dependent on r on S r , such that Proposition 2 ([36]).The constants C rj in (14) are given by if the grid is uniform, for The table ?? presents values for the constants C rj to k = 6 [36].For instance, when k = 3, and r = 1 we have the "small" stencil

Preprints
as shown in the figure 25.Then, by the proposition (1) there is a unique polynomial p r (x) such that then by taking u These approximations are of third order of precision if the function u(x) is smooth in the stencil S r with r = 0, 1, 2.
The WENO reconstruction procedure takes a convex combination of the low-order approximations, to achieve a high order approximation at the boundary points, based on a "big" stencil S = {x i−k+1 , . . ., x i+k−1 }, which is the union of all the "small" stencils S r considered above.The objective is to find a unique polynomial of degree at the most 2k − 2, denoted by p(x), interpolating the function u(x) in the points of the "big" stencil.Suppose u i+1/2 = p(x i+1/2 ) as the approximation to the value u(x i+1/2 ), we have to by virtue of the smoothness of the function u(x) in the grid S. 2. Find the constants d r and dr , which are called the linear weights, such that are true, with d r , dr ≥ 0 y ∑ k−1 r=0 d r = 1 = ∑ k−1 r=0 dr , because of the consistency [36].For example, A definition for β r which meets these requirements is given in [36] by:

Preprints
For example, computing for k = 3, ( 22) get 5. Finally, we get the WENO reconstruction of 2k − 1 order, of u(x) in the boundary points given by: where the super-index (L) indicates that the base S is skewed to the left, meaning that more grid points are used to the left of x i+ 1 2 to the right, similarly the (R) indicates that the base S is skewed to the right, by using more mesh points to the right of x i− 1 2 to the left.

Numerical flux schemes
We intend to calculate the approximations of the numerical flow F i+1/2 and F i−1/2 , in the border points of the cell I i = [x i−1/2 , x i+1/2 ], for a certain numeric flux F .The schemes to work are:

• Lax-Friedrichs flux
It is a monotonous central flux of first order, probably the best known and used due to its simplicity.It is given by [37,49]: where σ = max u | f (u)|.

• Third order TVD flux
It is a TVD flux of third order presented in [49], which is given for the linear scalar case by: where with where donde θ i is called the local flow parameter and is called the upwind-downward flow parameter and η is defined by Observation 3 ([49]).For the non-linear case, in the scheme (26), the velocity is a dependent function of u, that is a = a(u), and therefore Roe rates are used [36], defined by: and we change a by a i+1/2 .In this way, the local flux parameter θ i is redefined by problem, sometimes it is not enough, and it is necessary to use flux limiters, which together with an appropriate high resolution scheme, causes the solutions to decrease the total variation (see [18,19,27]).

Total Variation Disminishing property (TVD)
To measure the oscillations in a solution, we define the total variation TV(U n ) of a function U as where i = 1 and i = N are the boundary to the left and the right of the numerical grid, respectively.
[50].Therefore, a numerical scheme is TVD (Total Variation Disminishing) if: for all n and ∆t such that 0 ≤ n Deltat ≤ T. That is, the total variations do not grow over time, so that TV(U n ), for any n is bounded by TV(U 0 ) of the initial condition [13].

Flux limiter
The idea is to represent the numerical flux as an additive decomposition of a flux of under order F L (usually some monotone flux) and a flux of high order F H , such so that F is reduced to F H in smooth regions, or reduced to F L close to discontinuities [20].This is, where, φ is the limiter flux function, dependent on local flux parameter θ (or r in some cases), defined by: which represents the radius of the successive gradients in the given cell and measures the smoothness of the solution [15,49].
It is important to note that φ(θ i ) := φ i ≥ 0. When φ i = 0 (peak gradient, opposite slopes or zero gradient), the flux is represented by low order flux.Analogously, when φ i = 1 (smooth solution), the flux obtained is of high order (second order) [15].

Time discretization
In this section, time discretization is emphasized, since once the spatial part is worked for an initial time, the idea is to extend this procedure for a following time, which is implemented using a class of high Runge-Kutta methods.order [19], whose objective is to maintain the property TVD (32), when achieving a high order of precision in time [13].
This method is used to solve EDO systems in the same way where L(u) it is an approximation to the derivative −( f (u) x ) in the equation (1).
Analogous to the driscretization in space, we denote the time as t n = n • ∆t, so that we use U n i to denote the approximation to the exact solution u(x, t) in the mesh points (x i , t n ).
The Runge-Kutta methods that comply with the TVD property for some restriction in the CFL, are called Runge-Kutta SSP methods (Strong Stability Preserving) [14], in which the method is assumed from Euler forward is strongly stable, that is, ||U n i − ∆tL(U n i )|| ≤ ||U n i ||, and it is considered a restriction for the CFL, so that it is satisfied that ||U n+1 i || ≤ ||U n i || (see [13,50]).
According to [13], the third-order Runge-Kutta method TVD is given by It is important to remake that even when considering a good second-order spatial discretization TVD, if the temporal discretization is linearly stable but not TVD, then we can obtain oscillatory results (see [13]).Therefore it is suggested to use Runge-Kutta TVD/SSP methods for hyperbolic problems [14].

Discussion
In this work a WENO/WENO-TVD implementation was performed using the Lax-Friedrichs and TVD flux in the convective part for problems of hyperbolic conservation laws such as advection equations and Burgers (non-viscous).According to the results obtained, the following was concluded: First, that for initial conditions or smooth regions it is suggested to use the numerical flux of Lax-Friedrichs, while for non-linear problems with initial conditions in sections or in the presence of peaks and extremes, the TVD flux has better resolution.
According to the convergence tables, it was concluded that WENO schemes achieve a fourth order of precision, both when considering long times and when the initial condition has strong oscillations.
In addition, in the problems in which there are peaks like in the Burgers equation for T = 1.5 (after the shock) and in the Zalesak problem in which there are sharp, smooth and extreme peaks, it was observed that the best performance had was WENO-5-TVD.
Regarding the flux limiters, it was observed that in vanishing regions the Van Albada 2 limiter has better resolution, while the peaks and ends have a higher accuracy than the third order TVD limiter, guaranteeing greater resolution capacity than the other flux limiters.
Future work is intended to study the numerical solutions of the Euler equations of gas dynamics and shallow water, using variations of the WENO scheme.

Figure 23 .
Figure 23.Finite volume in one dimension.

Figure 24 .
Figure 24.Cell solution for the cell average U n i .

2 =
p r (x) as an approximation of u(x i+1/2 ), we can write such an approximation as a combination of the cell averages, then by virtue of the proposition (2) there are constants C rj such that Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 1 October 2018 doi:10.20944/preprints201810.0018.v1

Figure 25 .
Figure 25."Small" stencils for the reconstruction in x i+1/2 using three cells in each stencil.

Table 3 .
The constants C rj