Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme

For a model of nonlinear elastodynamics, we construct a finite volume scheme which is able to capture nonclassical shocks (also called undercompressive shocks). Those shocks verify an entropy inequality but are not admissible in the sense of Liu. They verify a kinetic relation which describes the jump, and keeps an information on the equilibrium between a vanishing dispersion and a vanishing diffusion. The scheme pre-sented here is by construction exact when the initial data is an isolated nonclassical shock. In general, it does not introduce any diffusion near shocks, and hence nonclas-sical solutions are correctly approximated. The method is fully conservative and does not use any shock-tracking mesh. This approach is tested and validated on several test cases. In particular, as the nonclassical shocks are not diffused at all, it is possible to obtain large time asymptotics.


Introduction
Hyperbolic systems of conservation laws often arise as limits of systems of partial differential equations including small scales effects, like diffusion and dispersion, when the parameters driving the small scales effects vanish.For example, the compressible Euler equations are the (formal) limit of the compressible Navier-Stokes equations, when the diffusion parameter tends to zero.However, the solutions of the hyperbolic system do not inherit all the properties of the augmented system, and it particular they are not always smooth.More importantly, in the class of weak solutions, the Cauchy problem for the hyperbolic system may admit infinitely many solutions, plenty of them being non relevant toward the augmented system.Thus it is necessary to add some criterions that select the solutions that are indeed limits, as the parameters vanish, of the solutions of the augmented diffusivedispersive system.One of this criterion yields to the selection of classical solutions, which roughly speaking corresponds to the case where diffusion is dominant.But when two small scale effects are of comparable strength (say diffusion and dispersion) nonclassical solutions appear.They contain shocks that do not verify the classical criterion but do arise as limits of the diffusive-dispersive system.
Let us sketch a general framework for nonclassical solutions of hyperbolic systems.Consider the Cauchy problem for a system of conservation laws where t is positive, x belongs to R, U is the vector regrouping the N ≥ 1 unknowns belonging to a open convex set Ω of R N , and f = R N → R N is the flux function, which we assume to be smooth.We suppose that System (1) is strictly hyperbolic, i.e. that the Jacobian matrix of f is diagonalizable with simple eigenvalues λ 1 (U ) < • • • < λ n (U ).We denote by r i (U ) a right eigenvector associated to the eigenvalue λ i (U ).
Solutions of hyperbolic systems are in general not smooth, even when the initial data is smooth.Discontinuities appear in finite time and it is necessary to consider weak solutions.Doing so, uniqueness is lost, and an additional criterion has to be imposed to select a unique solution.For example, the vector of unknown U is asked to verify, in addition to the initial System (1), a so called entropy inequality (2) The entropy U and the entropy flux W usually come from considerations on an augmented version of the system under study, where small physical terms like diffusion and dispersion are not neglected.Such an augmented system has additional properties, for example the total energy is conserved.Inequality (2) states that at the limit, this energy is dissipated.
In general, from the augmented system only one entropy inequality is deduced, thus the solutions of (1) are not expected to verify other entropy inequalities than (2).When the system is genuinely nonlinear, i.e. when and when the entropy U is stricly convex, the entropy inequality (2) contains enough information to select only one weak solution of (1).In the sequel we suppose that the entropy U is strictly convex.In the general case where (3) does not hold, the addition of a single entropy inequality is not sufficient to select a unique weak solution.
The question of uniqueness of the weak solution is strongly related to the choice of a criterion to decide wether a shock is admissible or not.Two states U L and U R of Ω are linked by a shock if there exists a real s, called the speed of the shock, so that the Rankine-Hugoniot relations hold.The shock is entropy satisfying if it verifies (2) in a weak sense, i.e. if (5) When the system is genuinely nonlinear (i.e. when (3) holds), ( 5) is equivalent to the Lax inequalities [Lax57] ∃i which themselves are an extension to systems (N > 1) of Oleinik's criterion [Ole57].
In [Liu74], Liu generalized this criterion to the non genuinely nonlinear case where ∇λ i • r i may vanish.Liu proved that his criterion selects a unique selfsimilar solution to the Riemann problem (1-2).We recall that a Riemann problem is the case of a piecewise constant initial data For systems that are not genuinely nonlinear, the Liu criterion is stronger than (5) (whereas it is equivalent for genuinely nonlinear systems).An admissible shock in the sense Liu always verifies (5), but some shocks may verify (5) without satisfying the Liu criterion.These shocks are called nonclassical.Moreover, unlike in the genuinely nonlinear case, the addition of the single entropy inequality (2) is not enough to regain uniqueness of the weak solution, because too many shocks are allowed.In other words, the inequality entropy (2) does not contain enough information on the small scale effects in the augmented system.Uniqueness can be regain by strongly constraining the nonclassical shocks to verify an algebraic relation of the form The function Φ ♭ is called a kinetic relation.
It has been proved that the addition of such a kinetic relation selects a unique weak solution to the Riemann problem for several models, including models of phase transitions in solids in [AK91], [LT01] and [LT02] or in liquids [LT00], [HL00], [SY95] and a model of magnetohydrodynamics [HL00].The link between kinetic relations and augmented diffusive dispersive equations is explored in [JMS95], [LeF02] and [HL97] in the scalar case and in [Sle89] for gas dynamics with a Van der Waals pressure law.The literature on the subject is so wide that we only cited a few key references.The interested reader will find a more comprehensive bibliography on this topic in [LeF02].
The numerical approximation of nonclassical solutions is challenging, because they are by nature very sensitive to the equilibrium between the small scale effects.The stake is to preserve this equilibrium at the numerical level.Roughly speaking, usual finite volume schemes introduce a numerical viscosity which destroys the equilibrium between the small scales effects.The diffusion becomes dominant and the schemes converge toward the classical solution, even though they are based on nonclassical Riemann solvers.
To overcome this difficulty, two types of schemes have been proposed.On the one hand, it is possible to use high order schemes consistent with the augmented system (see for example [HL98], [LMR02], [LR00] and [KR10]).As outlined in [HL98], it is necessary to discretize the flux of the hyperbolic part with a high enough order, once again to keep the exact balance between the small scales.On the other hand, some schemes relies on the tracking of nonclassical shocks (see for example [CL03], [CG08], [BCLL08], [Per11], [CCER12] and [CDMG14] ).In that case the kinetic relation is taken into account the scheme, typically through the use of an exact nonclassical Riemann solver.
The conservative finite volume scheme presented here belongs to this last category.It is built to be exact when the initial data is a nonclassical shock.In general, it does not introduce any numerical diffusion near nonclassical shocks (a particular class of nonclassical solutions) and turns out to capture correctly nonclassical solutions.This scheme extends to hyperbolic systems the discontinuous reconstruction scheme introduced in the scalar case in [BCLL08] and recently used in [CDMG14].
In the first section, we present the hyperbolic system admitting nonclassical solution for which we construct the scheme, namely the model of nonlinear elasticity studied in [LT01] and [HL00].The second section is devoted to the construction of the scheme itself.We present a way to select cells in which a nonclassical shock is reconstructed, explain how this shock is reconstructed and give the numerical flux.Eventually, we prove that the scheme is exact when the initial data is a nonclassical shock.In the third and last section, we propose several test cases involving nonclassical shocks.It appears that the proposed scheme capture nonclassical shocks very sharply.Let us outline that the scheme is fully conservative (unlike Glimm type schemes [CL03]) and uses a fixed grid, which allows the computation of solutions containing several interacting nonclassical shocks.
Acknowledgment The author warmly thanks Frédéric Lagoutière for his support and advice.

The Riemann problem for a nonlinear elasticity model
This paper is devoted to the numerical approximation of the solutions of the system of conservation laws where the stress σ is twice differentiable, and verifies We are interested in weak solutions of (7) that are morally limits when ǫ tends to 0 + of the augmented system Figure 1: Notations used for solving the Riemann problem.Φ −♮ (w) is the first Liu shock reachable from w.The entropy dissipation between w and Φ ♭ ∞ is zero.The states in red can be linked to w by a Liu shock, the blue ones by a nonclassical shock.
The parameter α is positive.From (9) we deduce the following entropy inequality for System (7) (see [LeF02]): In this section we briefly recall the results of Thanh and LeFloch [LT01] on the Riemann problem for System (7), i.e. the case where We adopt the notation of this paper.It is recalled on Figure 1.If the left state (v L , w L ) and the right state (v R , w R ) are linked by a 1-shock, its speed is equal to −s(w L , w R ), where and v R is given by If they are linked by a 2-shocks, its speed is +s(w L , w R ), and Graphically, the speed of a shock s(w L , w R ) corresponds to the slope of the segment joining (w L , σ(w L )) to (w R , σ(w R )).
Let us now focus on when a 1-shock is admissible or not, either in the sense of the entropy dissipation (10) or in the sense of Liu.The results for 2-shocks are the same, with the roles of w L and w R reversed in what follows.A careful study of Equation (10) applied to the shock yields that it verifies the entropy inequality if and only if where Φ ♭ ∞ (w L ) is the real having the opposite sign than w L for which the entropy dissipation (5) vanishes.The function Φ ♭ ∞ is its own inverse.For System (7), a shock verifies the Liu criterion if and only if for all w between w L and w R , Geometrically, it means that the segment joining The fact that System (7) is not genuinely nonlinear yields Thus, if w R lies in between Φ ♭ ∞ (w L ) and Φ −♮ (w L ), the shock between (v L , w L ) and (v R , w R ) dissipates the entropy (10) but is not admissible in the sense of Liu.Such shocks are called nonclassical shocks.
Many Riemann problems for System (7-10) admit an infinity of solutions (see [LT01], Section 3).Uniqueness can be obtained by imposing a kinetic relation that strongly constrain the nonclassical shocks by imposing where but are not necessarily equal.We now state the main result of [LT01].
The exact Riemann solver associated to Theorem 1.1 is the foundation of the scheme presented in the next section.Moreover, we will use the inner structure of the Riemann solution.Thus for the sake of completeness, we describe below the forward 1-wave and the backward 2-wave.
Proposition 1.2.[Thanh, LeFloch] Let (v L , w L ) be a fixed left state, and (v M , w M ) be a state linked to (v L , w L ) by a 1-wave.Then the nature of this 1-wave is the following: otherwise, it is just a 1-classical shock (in which w changes sign).
Proposition 1.3.[Thanh, LeFloch] Let (v R , w R ) be a fixed right state, and (v M , w M ) be a state such that (v M , w M ) and (v R , w R ) are linked by a 2-wave.Then the nature of this 2-wave is the following: otherwise, it is just a 2-classical shock (in which w changes sign).
The distinction between the last two cases in the last point of the enumerations of Proposition 1.2 and 1.3 can be reexpressed geometrically, as illustrated on Figure 2.
Proposition 1.4.The exists a function ϕ ♯ such that for i ∈ {1, 2}, for all real w and for all real w M such that ww M < wϕ −♭,i (w), the quantity |s(w, ϕ Figure 2: The solution contains a 2-nonclassical shock (in blue) followed by 2-classical shock (in red) if their speeds are correctly ordered, i.e if the slope of the blue segment is smaller than the red one (right).This depends on the relative position of Φ ♭,2 (w M ) and Φ ♯ (w M , w L ).

The scheme
We now describe the scheme.The time discretization is denoted by and at each time t n , the space is discretized with cells having all the same size ∆x.We denote by x n j−1/2 their extremities and by x n j there centers.The centers of the cells will move at each time step but every cell will always be of size ∆x.The speed of the mesh is defined by Integrating the first line of (7) over the set we obtain the integral formula The formula is identical for the variable w with the flux The principle of finite volume methods is based on this integral formula.A finite volume scheme writes where v n j plays the role of the mean value and f n,v j+1/2 plays the role of the flux on the edge (and similarly for w n j and f n,w j+1/2 ).The scheme is initialized with the exact average of the initial data A particular finite volume scheme is characterized by a particular choice of a formula expressing the numerical fluxes f n,v j+1/2 and f n,w j+1/2 in terms of the (v n k ) k∈Z and (w n k ) k∈Z .Our choice is described in the next sections.
Recall that our aim is to derive a scheme which is exact when the initial data is a nonclassical shock.In that case, the initial sampling (16) introduces a small amount of dissipation (unless if by miracle the shock initially falls on an interface).The numerical initial data contains an intermediate value that does not correspond to any pointwise value of (v 0 , w 0 ).More generally, at the end of every time step, finite volume schemes contains a L 1projection on the mesh.Details of the solution are lost in that step, and intermediate values are created in the shocks profile.The main idea of the discontinuous reconstruction scheme is to rebuilt, from the mean values (v n j , w n j ) j∈Z , an initial data that contains more details.This idea is also the foundation of well-known other schemes, like the MUSCL scheme [vL97] and its central version [NT90].Our reconstruction consists in adding nonclassical shocks in some cells in which we are able to detect that a nonclassical shock where lying inside the cell before the L 1 -projection.Thus the reconstruction will be very precise near nonclassical shocks, and not on the smooth parts of the solution, taking the opposite strategy of [vL97] and [NT90] Following this idea, our scheme can be decomposed in three elementary steps.
• First, we detect the special cells which are associated with nonclassical shocks.This detection step is described in Section 2.1.
• Then, we reconstruct nonclassical shocks inside those particular cells, in the sense that we replace the mean values (v n j , w n j ) by piecewise constant functions having the form of a nonclassical shock.This procedure is explained in Section 2.2.
• Eventually, the numerical fluxes are computed by letting the reconstructed nonclassical shocks evolve during the time step t n+1 − t n .The use of a moving mesh makes this computation easy, see Section 2.3.
In the last section, we extend the scheme to reconstruct classical shocks as well.

Detection of nonclassical shocks
The key idea of the scheme is to see, whenever it is possible, each mean value (v n j , w n j ) produced by the finite volumes scheme as the average of some nonclassical shock located somewhere inside the cell.Of course, where (v n j−1 , w n j−1 ) and (v n j+1 , w n j+1 ) are linked by a nonclassical shock, it is the one that should be reconstructed.In general, (v n j−1 , w n j−1 ) and (v n j+1 , w n j+1 ) are not linked by a single nonclassical shock, but by the full pattern of two waves, each of them likely to contain a succession of a classical shock or a rarefaction followed by a nonclassical shock (see Propositions 1.2 and 1.3).In that case, we reconstruct one of the nonclassical shock appearing in the solution of the Riemann problem.It is chosen thanks to the following Lemma.
If they are linked by a 2-nonclassical shock, then Proof.The condition w L w R < 0 always holds through a nonclassical shock.The second one is a straightforward consequence of the Rankine-Hugoniot relations (4).In the case of nonlinear elastodynamics (7), they write with i = 1 for 1-shocks and i = 2 for 2-shocks.The positive real s(w L , w R ) is defined in Formula (12).Eventually, the last condition is a reminder of Proposition 1.4.
Definition 2.2.For all integer j and all positive integer n, we denote by (v n j,⋆ , w n j,⋆ ) the intermediate state appearing in the Riemann problem with left state (v n j−1 , w n j−1 ) and right state (v n j+1 , w n j+1 ).The left and right desired reconstructed states, denoted respectively by (v n j,L , wn j,L ) and (v n j,R , wn j,R ), are defined as follows.
• If ) and w n j−1 w n j,⋆ < 0, the solution contains a 1-shock in which w changes sign.
-If this shock is nonclassical, it is preceded by a classical wave and we set and no mass is loss when replacing the mean value by v n rec inside the j-th cell.Reasoning similarly for the w variable, we replace w n j by We now state one trivial but crucial property of this reconstruction procedure.
Proposition 2.4.If (v n j−1 , w n j−1 ) and (v n j+1 , w n j+1 ) are linked by a single nonclassical shock, and if it exists α in (0, 1) such that In particular if the initial data is a nonclassical shock, we have v 0 rec = v 0 and w 0 rec = w 0 : the numerical diffusion introduced by the initial sampling ( 16) is cancelled by the reconstruction.
In general the two distances d n,v j and d n,w j are different, and it is possible that at least one of these distances does not belong to (0, ∆x).In that case we cancel the reconstruction, considering that seeing (v n j , w n j ) as the mean value of the detected nonclassical shock is not relevant.
Definition 2.5.The left and right reconstructed states are defined by: According to the Rankine-Hugoniot relations (4), the nonclassical shock reconstructed in cell j has velocity and we set arbitrary s n j to σ ′ (w j ) when no reconstruction is performed (which reads w n j,L = w n j,R = w n j ).

Advection of the reconstructed discontinuities
The fluxes are computed by letting the reconstructed nonclassical shocks evolve during the time step and by computing exactly what goes through the interfaces.However, two discontinuities reconstructed in adjacent cells can interact and the waves resulting from the interaction can meet the line x = x n j+1/2 within the time step.It follows that if we want to use a fixed grid (x n+1 j = x n j , or equivalently V n mesh = 0), the flux along the interface x = x n j+1/2 cannot be computed without resolving the wave interaction, which is obviously extremely costly.
This can been avoiding by using a moving mesh.Let us recall that by moving mesh, we mean that the centers of the cells are moving from time to time, but that their size remains constant equals to ∆x.Thus the numerical difficulties of handling cells with different and varying widths is avoided.The mesh speed is chosen such that it is larger than the maximum of the waves speed: and the time step such that a wave cannot cross more than an entire cell during the time step: These two hypothesis insure that any wave created at time t n inside the j-th cell can only cross the right interface x = x n j+1/2 + V n mesh t if V n mesh < 0 and the left interface It also imply that if two waves interact inside the jcell during the time step, the waves resulting from their interaction will not have time to catch up the left or the right interfaces.This is depicted on Figure 3. Thus, computing the flux along the j + 1/2-th interface x = x n j+1/2 + V n mesh t boils down to compute the time at which the nonclassical shock reconstructed in cell j crosses the interface if V n mesh < 0 (in cell j + 1 otherwise).On Figure 3, we also see that if no reconstruction is performed in cells j and j + 1, the flux is trivial to compute and is given by a simple left or right decentering (depending on the sign of V n mesh ), which exactly corresponds to the staggered Lax-Friedrichs flux.
We now compute the flux in details, when V n mesh is negative and verifies (23), like in Figure 3. Then the nonclassical shock reconstructed in cell j may only cross the right interface x = x n j+1/2 + V n mesh (t − t n ).The discontinuity is not located at the same place for the variables v and w, therefore we compute two different crossing times Once again, Condition (23) insures that the flux passing through the j + 1/2 interface only comes from waves arriving from the cell j.It can be clearly seen on Figure 3 that it is Remark 2.7.The idea of using a moving mesh to simplify the computation of the fluxes goes back to the Lax-Friedrichs scheme.The Rusanov scheme follows the same idea by using a local mesh speed V j,n mesh to make the scheme less diffusive.Higher order extensions based on piecewise polynomial reconstructions are possible, see for example [NT90].

Detection of classical shocks
It is also possible to detect and reconstruct shocks in which w does not change sign.Those shocks are always classical.Their detection is based on the following proposition.
Proposition 2.8.If the states (v L , w L ) and (v R , w R ) are linked by a 1-shock in which w does not change sign, then either If they are linked by a 2-shock in which w does not change sign, then either There is no conflict with the detection of nonclassical shocks of Proposition 2.1, and we can straightforwardly extend Definitions 2.2 and 2.5 to take into account those shocks.

Exact approximation of isolated nonclassical shocks
The aim of this section is to prove that the scheme described above is exact when the initial data is an isolated nonclassical shock, i.e. (11) with left state (v L , w L ) and right state (v R , w R ) verifying the Rankine-Hugoniot relation (18) and constrained by the kinetic relation ( 14).We recall that in that case the exact solution is with i = 1 for a 1-shock and i = 2 for a 2-shock.
Theorem 3.1.The scheme described in the previous section is exact when the initial data is a single nonclassical shock.In other words, the numerical solution is the L 1 -projection on the mesh of the exact solution at time t n : for all n ≥ 0 and for all j ∈ Z, v exa (t n , x) dx and w n j = 1 ∆x Figure 4: Structure of w when a 1-phase transition is detected in cell −1.In that case, the Riemann solution (on the right) contains three shocks.
Proof.We prove the result for a 1-nonclassical shock such that w L < 0 < w R , the other cases being exactly similar.With the initial sampling ( 16), the property holds true at time t 0 .Suppose that is holds true at some time t n .At this time, we renumbered the cells in a way that the discontinuity lies inside the cell numbered 0, and we denote by δ its distance to x n −1/2 .We have and Note that as The detection step of Proposition 2.1 detects a 1-nonclassical shock in cell 0. By Proposition 2.4, the nonclassical shock is reconstructed in that cell, and it is be placed exactly at the right position in the reconstruction step: If no reconstruction is performed in the other cells, our scheme gives the correct solution at time t n+1 by construction.This is easy to check, but a little tedious, and we refer the reader to [Agu14] for a detailed proof.Suppose that w n 0 > 0. Then by Proposition 2.8, no classical shock is detected in cell 1, and by Proposition 2.1, a 1-nonclassical shock is detected in cell −1.Let us focus on the solution Riemann problem between the state (v L , w L ) in cell −2 and the state (v n 0 , w n 0 ) in cell 0. The notation are recalled on Figure 4.If wn −1,L ≥ w L , both wn −1,L and wn are larger than w L .Thus it is be impossible to have d n,w −1 ∈ (0, ∆x) and no reconstruction occurs in cell −1.Suppose now that w L is larger than wn −1,L .Then w L and wn −1,L are linked by a 1-classical shock, and wn −1 in (0, 1).Indeed, the kinetic function φ ♭,1 decreases, thus states that along this curve, v is an increasing function of w.Therefore vn −1,R is larger than v R .On the other hand, the Rankine-Hugoniot relations applied to the 2-shock yields that vn −1,R is smaller than v n −1 , which contradicts the fact that v n −1 < v R .The case where a 1-nonclassical shock is detected in cell 1 is simpler.We have w n 0 < 0 and nothing is is detected on cell −1.Moreover, wn 1,L ≤ wn 1,R , thus it is not possible to reconstruct w in a conservative manner if wn 1,R ≤ w R .On the other hand, if wn 1,R > w R , the solution contains a 2-classical shock and (4) yields that vn 1,R ≤ v R .The Rankine-Hugoniot applied to the 1-nonclassical shock implies that vn 1,L ≤ vn 1,R .Thus it is impossible that d n,v 1 belongs to (0, ∆x).
Remark 3.2.It is crucial to require that both v and w are reconstructed in a conservative manner.In the case where the nonclassical shock is detected on cell 1, it is clear that the proof collapses if we authorize d n,v 1 outside (0, ∆x).Moreover if we authorized d n,v 1 outside (0, ∆x), the solution of the Riemann problem between (v n 1 , w n 1 ) and (v R , w R ) might contain a 1-classical shock, a 1-nonclassical shock and 2-rarefaction wave, in which case and a reconstruction is performed in cell 1.
The case β = 1/2 corresponds to the classical solutions ; the choice β = 1 corresponds to the case where the entropy dissipation (10) is zero across nonclassical choice.It does not fall in the theory of [LeF02], but the Riemann problem can be solved (see [LT01]) and it is possible to explore that case numerically.
Test 1: Isolated nonclassical shock This test case illustrates Theorem 3.1.The initial data is the Riemann problem: With m = 1 and β = 2/3, the solution is an isolated nonclassical 1-shock.On the top of Figure 5, we plot the exact nonclassical solution and the solution given by the reconstruction scheme.As expected they are exactly the same.On the bottom of Figure 5, we plot the solution given by the Godunov scheme based on an exact nonclassical Riemann solver.It does not capture the nonclassical solution but the classical one, which in that case in a rarefaction followed by a shock.This is a general phenomenon: usual finite volume schemes are not able to capture nonclassical solutions.Thus in the sequel we used the Glimm scheme to compare our scheme with.The CFL number is set to 0.45, the final time is T = 0.038 and the space interval [−0.5, 0.5] is discretized with 200 cells.

Interlude: The Glimm scheme
In the sequel we use the random sampling method of Glimm [Gli65] to compute reference solution.The Glimm scheme is built as follow.Let (r n ) n∈N be a sequence of i.i.d.random variables uniformly distributed over [0, ∆x].Suppose that at time t n a piecewise constant approximation of the solution is given, and denote by (t, x) → U n exa (t, x) the exact solution at time t ≥ 0 for that initial data.The numerical solution of Glimm at time t n+1 = t n +∆t n is given by ∀j It ∆t n is small enough, i.e. if it is smaller than the maximum of the waves speed appearing in the Riemann problems multiplied by the cell size ∆x, U n exa is just the juxtaposition of Riemann solutions at each interface.
The main feature of the Glimm scheme is that is does not introduce any numerical diffusion, in the sense that the shock profiles are not smeared out at all.This is why this scheme has been used in [CL03] to approximate nonclassical solutions.The two drawbacks are that the scheme is not conservative and that the exact computation of the solution is costly.
Test 2: A Riemann problem with two nonclassical shocks The initial data is now The exact solution consists in a 1-classical shock, a 1-nonclassical shock, a 2-nonclassical shock and a 2-rarefaction wave.On Figure 6, we plot the exact solution and the numerical solution given by the reconstruction schemes.We compare the reconstruction where we detect nonclassical shocks only (cf Proposition 2.1), referred to in the legends as RecNC, and the reconstruction scheme were the classical shocks are also detected (cf Proposition 2.8), referred to in the legends as RecNC+C.Both of them capture very sharply the nonclassical shocks; the 1-classical shock is much more diffused when it is not reconstructed, and both scheme behaves in the same way in the rarefaction wave.The CFL number is 0.45, the space interval [−1, 1] contains 200 cells.
Test 3: Perturbation of a classical shock This test case is taken from [CL03].We still have β = 2/3 but now m = 2.The initial data is with ǫ ∈ {0, 0.05, 0.1} When ǫ = 0, the solution is a 1-classical shock.When ǫ > 0, the classical shock is split into a classical shock followed by a nonclassical shock.Their speeds are different but close to each other.We plot the solutions at time T = 0.4 on Figure 7, using 600 cells per unit interval and a CFL number of 0.45.The solutions are very well approached when ǫ > 0. When ǫ = 0, a spike appears.This spike corresponds to the state linked to (v R , w R ) by a nonclassical shock (on Figure 7 we see that it has the exact same height as the nonclassical shock appearing for ǫ > 0).Remark that ǫ = 0 is exactly the limit between the two cases in the last point of the enumeration of Proposition 1.2, for which the speed of the nonclassical and the classical shock coincide (see also Figure 2).Numerically the mechanism is the following.After one iteration in time, an intermediate value the cell just before, and this time the reconstruction succeeds.Thus the scheme is not exact in that case.Note that the proof of Theorem 3.1 uses the kinetic relation and thus collapses when the reconstructed shock does not verify the kinetic relation (and is hence classical).This phenomena does not prevent the scheme from converging.The numerical solution has the same shape than in the case ǫ > 0: a classical shock followed by a nonclassical shock, but this time they have the same speed, thus they remain at the same position.The spike is only two cells wide when the classical shocks are reconstructed, and hence not diffused.If they are not reconstructed, the spurious classical shock is diffused.
In conclusion, this test case shows a limitation of scheme, namely its incapability to approach exactly classical shock in which w changes sign, but also show its ability to capture finely nonclassical solutions, the gap between the shocks being very thin here.In particular, the apparition of the intermediate state is immediate, while it is linked to the ratio of the width of the gap at time T and of the cell size ∆x when using a random sampling based method like the Glimm scheme.
Test 4: Apparition of nonclassical waves from a smooth initial data We now focus on the case of a smooth initial data v 0 (x) = 3 sin(2πx), w 0 (x) = 1 + 3 cos(8πx), with β = 0.95 and m = 1 and periodic boundary condition.Nonclassical shocks appear around time t = 0.011 and then propagate in the solution.On Figures 8, 9 and 10 we compare the solution given by the reconstruction schemes (with or without reconstruction of classical shocks) at time t = 0.015, t = 0.06 and t = 1.The space interval contains 1 000 cells per unit interval and CFL number is set to 0.45.The reference solution is given by the Glimm scheme with 8 000 cells.We see that the reconstruction schemes capture accurately the nonclassical shocks.The result are poorer in smooth areas, because in those regions the reconstruction schemes behave as the Lax-Friedrichs scheme which is quite diffusive.This can be improved by using another scheme on interfaces where no reconstruction is performed.The use of a moving mesh encourages to chose a central scheme like the Nessyahu and Tadmor scheme [NT90].This scheme is second order accurate in smooth regions of the solutions.It is both easy to implement and fast, as it does not use any information on the Riemann problems.We reproduce here the test case of [CL03] for which β = 1 and the final time is large.The case β = 1 corresponds to the case where the entropy dissipation is zero through nonclassical shocks.It is a limit case in the theoretical framework of [LeF02].As the Riemann solver is perfectly defined for β = 1, it can be explored numerically.Its mean value is null.The parameter m is fixed to 0.05.On Figure 11, we plot the solution at time t = 40 with 2 000 and 8 000 points per unit interval.We can see that at that time w changes sign three times, instead of two in the initial solution, and does not converge to zero.The position of the nonclassical shocks are different with the two schemes.Let  us recall that the Glimm scheme is based on a random sampling of the solution, thus two realizations of the same test give different results.On Figure 12, we plot the histogram of the first nonclassical shock position at time t = 20 for 100 independent realizations of the Glimm scheme with 2 000 cells (bottom) and the comparison of the two schemes at the same time (with 8 000 points).Moreover, the structure with 2 nonclassical shocks very close to each other only appears in 31 out of those 100 realizations.Its indicates that the width of this structure is small (and probably smaller than the one appearing in the realization of [CL03]).

Figure 5 :
Figure5: The nonclassical Godunov scheme is unable to capture nonclassical shocks.On the contrary, the reconstruction scheme captures isolated nonclassical shocks exactly.

Figure 6 :
Figure 6: A Riemann problem with two nonclassical shocks at time 0.15.

Figure 11 :
Figure 11: Solution of test 5 at time t = 40 (the shape of v is similar)

Figure 12 :
Figure 12: Solution of test 5 at time t = 40 (the shape of v is similar)