Riemann problems with non--local point constraints and capacity drop

In the present note we discuss in details the Riemann problem for a one--dimensional hyperbolic conservation law subject to a point constraint. We investigate how the regularity of the constraint operator impacts the well--posedness of the problem, namely in the case, relevant for numerical applications, of a discretized exit capacity. We devote particular attention to the case in which the constraint is given by a non--local operator depending on the solution itself. We provide several explicit examples. We also give the detailed proof of some results announced in the paper [Andreainov, Donadello, Rosini,"Crowd dynamics and conservation laws with non--local point constraints and capacity drop", which is devoted to existence and stability for a more general class of Cauchy problems subject to Lipschitz continuous non--local point constraints.


Point constraints in traffic modeling
Traffic modeling is an exciting and fast-developing field of research with plentiful applications to real life. While this subject was initially limited to the description and the management of vehicular traffic, we see a growing interest nowadays on different applications as crowd dynamics and bio-mathematics. This note is related to an extensive on-going research project concerning the theoretical and the numerical study of macroscopic models for which the definition of solution involves an artificial limitation of the flux in a finite number of points. From the modeling point of view, this may correspond to a narrow exit in crowd modeling, a toll gate in vehicular traffic, a cell membrane in bio-medical modeling.
In the pioneering paper [7], R. Colombo and P. Goatin introduced point constraints in the classical one-dimensional LWR road traffic model [17,18], with the goal to model the presence of obstacles on the road as toll gates and road lights. This model reads as f (ρ(t, 0±)) ≤ q(t) t ∈ R + (1b) where ρ = ρ(t, x) ∈ [0, R] is the unknown (mean) density at time t ∈ R + of vehicles moving along the road parameterized by the coordinate x ∈ R. Then, R ∈ R + is the maximal road density, f : [0, R] → R is the nonlinearity relating the flux in the direction of increasing x to the density, q : R + → R + is a given function prescribing the maximal flow allowed through the point x = 0, and ρ 0 : R → [0, R] is the initial (mean) density. Finally, ρ(t, 0−) denotes the left measure theoretic trace along the constraint implicitly defined by for all φ ∈ C ∞ c (R 2 ; R). The right measure theoretic trace, ρ(t, 0+), is defined analogously.
In the above setting, the authors of [7] were able to prove existence and wellposedness of solutions. Further theoretical considerations and numerics associated to this model have been developed in [4].
In view of the applications, however, we also need to consider the case in which the evolution of the constraint function q is not given beforehand, but depends on the solution ρ itself in a neighborhood of x = 0. In such situation we say that the point constraint is non-local. In this way we obtain crowd and cell membrane dynamics models described by coupled PDE-ODE systems for which the existence and well-posedness of solutions are not trivial matter. Nevertheless, this case has a practical interest. In crowd dynamics, as an example, the experimental observations by E. Cepolina in [6] prove that the irrational behavior of pedestrians at bottlenecks ends up by reducing the maximal possible outflow. This phenomenon, called capacity drop, is also related to other effects observed in crowd dynamics, such Faster Is Slower and the Braess' paradox.
In full generality, we may consider the constraint function q as follows where Q : C 0 [0, T ]; L 1 (R; [0, R]) → L 1 ([0, T ]; [0, R]). The minimal regularity properties to impose on Q in order to achieve well-posedness of solutions are not known at the moment, and they are the object of one of our current research projects.

An example of non-local point constraint
In the paper [3], B. Andreianov, C. Donadello and M. D. Rosini proposed a model which generalizes the one in [7] and consists of a Cauchy problem for a onedimensional hyperbolic conservation law as 1 subject to a non-local point constraint of the form Here p : R + → R + prescribes the maximal flow allowed through an exit placed in x = 0 as a function of the weighted average density of pedestrians, ρ, in a left neighborhood of the exit and w : R − → R + is the weight function used to average the density. The authors of [3] proved well-posedness of the Cauchy problem in L ∞ (R; [0, R]) for the model under the following assumptions on the regularity of f , w and p, see Figure 1, and it is a non-increasing map.
We recall that, in the previous literature, the only macroscopic model able to reproduce the capacity drop at bottlenecks is the CR model introduced by R. Colombo and M. D. Rosini in [9]. The Riemann solver for the model described in [9] is fairly intricate; in this note, our main goal is to describe in an exhaustive way the Riemann solver (or, rather, solvers) for the model of [3]. Notice that, in contrast to [9], a specific nonclassical Riemann solver has to be used only at the exit point x = 0 while the simple classical Riemann solver is used elsewhere.
The notion of solution we adopt is a natural extension of the one introduced in [7], for a constrained Cauchy problem of the form 1.
2. In addition q is linked to ρ by the relation 2.
If q is given a priori, then 3 is the definition of entropy weak solution to problem 1. We refer to Proposition 2.6 in [4] for a series of equivalent formulations of conditions 3.
The next proposition lists the basic properties of a entropy weak solution of 1, 2, for the case of a single valued p, the proof is given in [3].
be an entropy weak solution of 1, 2 in the sense of Definition 1.1. Then (1) It is also a weak solution of the Cauchy problem 1a, 1c.
(3) Any discontinuity away from the constraint is classical, i.e. satisfies the Lax entropy inequalities, see [12].
(4) Nonclassical discontinuities, see [16], may occur only at the constraint location x = 0, and in this case the flow at x = 0 is the maximal flow allowed by the constraint. Namely, if the solution contains a nonclassical discontinuity for all times t ∈ I, I open in R + , then for a.e. t in I If the constraint function p is multivalued the equalities in 2 and 4 should be interpreted as inclusions, and the result of the proposition remains true.
The existence result in [3] is achieved by a procedure which couples the operator splitting method [13], with the wave-front tracking algorithm, [11], see also [1] for a similar technique. This approach allows us to approximate our problem by a problem with "frozen" constraint, as 1, at each time step.
The regularity of p plays a central role in the well-posedness result. While existence still holds for the Cauchy problem when p is merely continuous, it is difficult to justify uniqueness in this case. Further, in [3], the authors give some basic examples illustrating that solutions of a Riemann problem for the case of a non-decreasing piecewise constant p, see (P2) below, may fail to be unique, L 1 loccontinuous and consistent.
The case in which p is piecewise constant is extremely important both for the theoretical study of the problem and its numerical applications. First, it is related to the construction of the Riemann solver, which is the basic building block for the wave-front tracking algorithm, a precious tool in the study of existence and stability of the solutions for the general Cauchy problem. Moreover, the piecewise constant case is essentially the only case in which solutions can be computed explicitly, which is an undeniable atout when looking for examples and applications. To this aim, it is relevant to provide a detailed study of the different pathological behaviors one may encounter. Remarkably, we show that these behaviors can be easily forecast and avoided when looking for explicit examples of solutions.
In the following section we develop a detailed proof of the fact that, when p is piecewise constant, the Riemann solver for 1 with a constraint of the form 2, is not unique and does not satisfy the minimal requirements needed to develop the classical wave-front tracking approach. Additionally, we compare the two extreme Riemann solvers: the one that minimizes the capacity drop, and the one that maximizes it. In particular, for any time T > 0 we estimate the distance between the profiles of the solutions produced by the two Riemann solvers starting from the same initial condition.

The constrained Riemann problem
In this section we study constrained Riemann problems of the form . The flux f and the weight function w satisfy (F) and (W), moreover, we adopt the following assumption on p (instead of (P1)) to allow an explicit construction of solutions to 5 is piecewise constant non-increasing map with a finite number of jumps. Unfortunately, the regularity of p required by (P2) is not enough to apply the well-posedness results in [3]. In particular, Example 2 in [3] illustrates the loss of uniqueness and stability of entropy weak solutions. In this section we present a systematic study of the possible pathological behaviors. We denote by R the classical Riemann solver. This means that the map is the unique entropy weak solution for the unconstrained problem 5a, 5c, see for example [5] for its construction. Whenever the classical weak solution given by R does not satisfy the constraint 5b, we replace it by a nonclassical weak solution, see [16] as a general reference, where the mapsρ,ρ : We stress that by Proposition 1 any nonclassical entropy weak solution is a classical entropy weak solution in the Kružkov sense, see [15], [12], in the half-planes R + ×R − and R + × R + . However, the jump at x = 0 is a nonclassical shock, in the sense that it does not satisfy the Lax entropy inequalities. First, we should notice that as soon as the constraint function t → q(t) is not constant, the solutions of the Riemann problem may not be self-similar.
When the traffic light is green, i.e. for t ∈ [2k, 2k + 1[, the flow at x = 0 is free from any constraint. Conversely, when the traffic light is red, i.e. for t ∈ [2k + 1, 2(k + 1)[, the admissible flow at x = 0 become zero. This means that the Riemann problem with initial condition at t = 0 given by ρ L = ρ R = 0 will not be self-similar, because the constant solution will not satisfy the constraint starting from t = 1.
The above example also shows that as soon as we use a nonclassical Riemann solver we lose the a priori BV bounds on the solution.
In the proof of Proposition 2 we show that any entropy weak solution of 5 is self-similar for sufficiently small times. Therefore, it makes sense to introduce nonclassical local Riemann solvers, see Definition 2.3. Then, the availability of a local Riemann solver allows us to construct a global solution to the Riemann problem 5 by a wave-front tracking algorithm in which the jumps in the map [t → q(t)] are interpreted as non-local interactions.
Aiming for a general construction of the solutions to 5, we allow p to be a multi-valued piecewise constant function, namely, see Fig. 1, right: . . , n, being ξ 0 = 0 and ξ n+1 = R.
In the following we will use the notationsρ i =ρ(p i ) andρ i =ρ(p i ).
As it will become clear in Proposition 2, the possible loss of uniqueness and stability can be easily forecast once the piecewise constant constraint p and the flux f are given. In particular, for some respective configurations of p and f the solution of the Riemann problem exists and is unique, locally in time, for any initial data where we say that (ρ L , ρ R ) satisfies condition (C) if it satisfies one of the following conditions: where we say that (ρ L , ρ R ) satisfies condition (N) if it satisfies one of the following conditions: Here, C stands for classical and N , for nonclassical, in relation with the nature of the shock appearing in the solution of 5 at x = 0. Observe that if the constraint function p is constant in a neighborhood of the state ρ L , then p(ρ L −) = p(ρ L +) and this simplifies the above conditions. Also a right or left continuity assumption on p would simplify the above definitions.
The next proposition says that uniqueness holds at least for small times if and only if the initial data are in C ∪ N . It is fundamental to remark that since nonuniqueness is possible only when p(ρ L −) = p(ρ L +) and p(ρ L −) ≥ f (ρ L ) ≥ p(ρ L +), non-uniqueness concerns at most a finite number of left states and therefore the region [0, R] 2 \ (C ∪ N ) is the union of a finite number of line segments.
is the unique entropy weak solution at least for t > 0 sufficiently small.
, then the corresponding constrained Riemann problem 5 admits more than one entropy weak solution.
The proof is deferred to Section 4.1 Remark 1. Once the function p is fixed, the time interval [0, τ ] on which the solution to the Riemann problem 5 is self-similar can be estimated from the initial datum (ρ L , ρ R ) provided it belongs to C ∪ N . On the contrary, when (ρ L , ρ R ) ∈ [0, R] 2 \ (C ∪ N ), we are not always able to forecast when the next "interaction with the constraint" will take place. In some situations, a whole one-parameter family of solutions exists, we refer to Example 2 in [3] for a detailed discussion of this case.

Riemann solvers
As the local solutions of the Riemann problem 5 are not unique in general, we are naturally led to question the existence of suitable selection criteria. All the solutions we introduce are solutions in the Kružkov sense in the open half-planes R + × R + and R + ×R − , so they satisfy the basic requirement of entropy dissipation. However, coming back to the real situations which our model aims to describe, we argue that the most interesting behaviors to track correspond to the extreme cases in which the flux at the exit is either the highest or the lowest possible from a given initial condition.
If, as an example, our model describes the evacuation of a narrow corridor, it is clear that the optimal solution corresponds to the highest admissible values of the flux at the exit. By opposition to the next case, we describe this situation as quiet behaviour. In analogy to the discussion in [14] we interpret all other possible solutions as consequences of an irrational behavior, which in literature is often described as panic. In particular, we can use the solution corresponding to the lowest admissible values of the flux at the exit to find an upper bound for the evacuation time.
From now on we restrict ourselves to the case in which p(ξ i ) can only take the values p i and p i+1 and not the intermediate values.
Definition 2.3. Two Riemann solvers R q and R p for 5 are defined as follows for t > 0 sufficiently small and x ∈ R: (NNN4), (NNN5) If (ρ L , ρ R ) satisfies one of these sets of conditions then R q [ρ L , ρ R ](t, x) takes the form 11 withp = p(ρ L −) and R p [ρ L , ρ R ](t, x) takes the form 11 with p = p(ρ L +).
In the next proposition we collect the main properties of the Riemann solvers R q and R p . In particular (R6) means that the Riemann solver R q is the one which allows for the fastest evacuation, while R p is associated to the slowest one.
Then, for = q, p: x)] is a weak solution to 5a, 5c.
(R6) R q [ρ L , ρ R ] maximizes the flux at the exit, in the sense that if E is the set of all entropy weak solutions of the Riemann problem 5, we have Analogously, R p [ρ L , ρ R ] minimizes the flux at the exit, in the sense that We recall that a Riemann solver is said to be consistent when the juxtaposition of the solutions of two Riemann problems with respective initial conditions (ρ L , ρ M ) and (ρ M , ρ R ) is the solution of the Riemann problem with datum (ρ L , ρ R ). Moreover the vice versa also holds true, in the sense that whenever the state ρ M is an intermediate state in the solution of the Riemann problem with initial condition (ρ L , ρ R ), then the solution consist of exactly the same states and waves which we would obtain by solving side by side the two Riemann problems with data (ρ L , ρ M ) and (ρ M , ρ R ).
The proof of Proposition 3 is deferred to Section 4.2.

Remark 2.
It is important to observe that even if p(ξ i ) can only take the two values p i and p i+1 , this is not enough to rule out the existence of infinitely many different solutions as the ones described in Example 2 of [3], in the case However, each of the extremes Riemann solvers R , = p, q, selects one of them because it sticks to the constant level of constraint prescribed by Definition 2.3, the level p i for R q and the level p i+1 for R p , until a non-local interaction takes place.

On the use of Riemann solvers R , = q, p
Although the Riemann solvers R are not L 1 loc -continuous, an existence result for the Cauchy problem 1 can be obtained from a wave-front tracking algorithm based on R , see for instance [10], [19]. Such approach using R does not require the operator splitting method. However, the non-local nature of the approximating problems prevents us from a direct application of the Riemann solvers R . In fact, even in a arbitrary small neighborhood of x = 0, to prolong the approximating solution ρ n beyond a time t =t > 0 it is not sufficient to know the traces ρ n (t, 0−), ρ n (t, 0+), but also the value 0 −iw w(x) ρ n (t, x) dx is needed. Roughly speaking, because of the non-local character of the constraint one cannot merely juxtapose the solution to the Riemann problem associated to the values of the traces at x = 0 with the solution to the Riemann problems away from the constraint. Finally, also jumps in [t → p (ξ(t))] have to be considered as (non-local) interactions. Therefore, the approach using R is considerably heavier and more technical than the one presented in [3].

On the comparison between the Riemann solvers R p and R q
In this section we aim to compare the solutions obtained by the two Riemann solvers introduced above, starting from the same initial condition (ρ L , ρ R ). It is clear from Definition 2.3 that as soon as (ρ L , ρ R ) belongs to C ∪ N the solutions obtained by the two Riemann solvers coincide. As a preliminary remark we stress that adapting the proof of Theorem 3.1 of [3] to the case in which p is discontinuous, one can get a rough upper bound, exponential in time, for the L 1 distance between two solutions ρ 1 , ρ 2 obtained from the same initial condition (not necessarily of Riemann type). Indeed, instead of the bound valid when p is Lipschitz continuous, in the discontinuous case one can use the bound where h is the maximal size of jump in p(·) and N is a constant (observe that if p(·) is a discretization of some Lipschitz function, see Section 4.1 of [3], then N can be taken independent of h). From the fundamental stability estimate of Proposition 2.10 of [4], using 8 and the Gronwall inequality one easily gets the bound whenever ρ 1 (0, ·) = ρ 2 (0, ·). This rough estimate is enough to show that as h goes to zero, the discrepancy between different solutions vanishes and this argument applies to any initial datum, not necessarily of the Riemann type. However, the exponential growth with respect to t of the upper bound 9 is clearly not optimal when we aim to compare the solutions of a Riemann problem, since it does not take into account the specific self-similar structure of solutions valid at least on a small interval of time [0, T ].
Let us demonstrate that in the case where different Riemann solvers co-exist, the L 1 distance of the associated solution grows at most linearly both in h and in t ∈ [0, T ] (see also the numerical experiment on Figure 6).
In order to keep our presentation as light as possible, we focus on only one of the possible cases in which R p and R q differ. All other cases can be handed in a similar way. Assume that p(ξ) = p 1 χ [0,ξ] (ξ) + p 2 χ ]ξ,R] (ξ), wheref > p 1 > f (ξ) > p 2 > 0, and ρ L =ξ, ρ R =ρ, see Figure 3. We only consider solutions in a small interval of time [0, T ] in which they are self-similar. Then we get where the valuesρ i andρ i , for i = 1, 2, are implicitly defined by the relations f (ρ i ) = f (ρ i ) = p i andρ i ≤ρ ≤ρ i . More explicitly, we can say that the solution corresponding to R p consists of a shock of negative speed λ p between ρ L andρ 2 , a stationary nonclassical shock betweenρ 2 andρ 2 and a shock of positive speed µ p betweenρ 2 andρ. The solution corresponding to R q consists of a rarefaction wave between ρ L andρ 1 , a stationary nonclassical shock betweenρ 1 andρ 1 and a shock of positive speed µ q betweenρ 1 andρ, see Figure 4.
As the characteristics of this problem propagate with finite speed, we expect the solutions associated to the two solvers coincide outside a bounded interval. The geometry of the problem, see Figure 3, implies that µ p > µ q and that λ p is smaller than all the propagation speeds in the rarefaction wave between ρ L andρ 1 . Therefore, at time t ∈ [0, T ] fixed, the two solutions coincide outside the interval [λ p t, µ p t].
The value of the distance R p [ρ L , ρ R ](t) − R q [ρ L , ρ R ](t) L 1 (R;R) corresponds, loosely speaking, to the value of the area between the profiles of solutions, see Figure 4.
Following the same technique as in [5], Chapter 7, we can estimate the distance between the solutions profiles at a fixed time t ∈ [0, T ]. For the reader convenience, we recall that the propagation speed of the shock discontinuity between states the ρ a and ρ b is given by the Rankine-Hugoniot condition and that the propagation speed of the characteristics in a rarefaction wave joining the states ρ a and ρ b varies between the values f (ρ a ) and f (ρ b ). Also, by definition Figure 4: The solutions corresponding to R p , left, and to R q , center, and the comparison between their profiles at fixed time, right, as described in Section 2.3.
This means that whenever the piecewise constant function p we consider is the discretization of a smooth function we can bound a priori the size of the error due to the lack of uniqueness and we can make it smaller and smaller as h tends to 0.

Numerical results
We present here some numerical experiments in order to illustrate the results of the above section. The scheme used for the simulations combines the ideas of [4] with an explicitly updated constraint computed from weighted space averages of the discrete solution at previous time step. We will justify in the future work [2] convergence of this scheme to an entropy solution of the nonlocally constrainted problem in the sense of Definition 1.1, where the constraint function p(·) must be taken multi-valued. While it is delicate or even impossible to identify a unique Riemann solver to which the scheme would converge, we can use the simulations on Figure 5 to illustrate the fact that non-uniqueness for the Riemann problem results as unstable behavior in a vicinity of some specific data.  Figure 5 is shown the computed solutions ρ 1 and ρ 2 corresponding respectively to the following initial states Note that we have considered ρ 1 L and ρ 2 L such that ρ 1 L − ρ 2 L ∆x/8, where the space step ∆x = 0.025. Finally we assume the time step ∆t = ∆x/10. Nonetheless, as shown in Figure 6, in practice the instability is limited to a behavior of kind where h is the maximal size of jump in p(·), α is close to 1 and C > 0 is a constant. Indeed we reported in Figure 6 the computation of the L 1 -discrete norms of the difference ρ 1 −ρ 2 when we take, p 2 = 0.05, 0.075, 0.1, 0.125 and 0.15 in the definition of the constraint function. Using logarithmic scales, we deduce that the distance between the two solutions is approximatively proportional to |p 1 − p 2 | 0.9 .

Proof of Proposition 2
First, we introduce the notation Therefore ξ(0) = ρ L and the map [t → ξ(t)] is continuous. We stress that any nonclassical entropy weak solution in the sense of Definition 1.1 is also a classical entropy weak solution in the Kružkov sense in the halfplanes R + × R − and R + × R + . Therefore, at least for t > 0 sufficiently small, by Proposition 1, assumption (P2) and the continuity of the map [t → ξ(t)], any nonclassical entropy weak solution of 5 must have the form, see Fig. 2,  Figure 6: The norm ρ 1 − ρ 2 L 1 with respect to |p 1 − p 2 | in log/log scale.
Observe that 11a is uniquely identified once we knowp which, by 4, satisfies We recall that 11b means in particular that the Rankine-Hugoniot jump condition is satisfied at x = 0 even when the solution to the Riemann problem is nonclassical. As a consequence of 11b, of assumption (P2) and of the continuity of [t → ξ(t)], we have thatp This implies that if p(ρ L +) = p(ρ L −), then p(ξ) is constant in a neighborhood of ρ L and, since the solution is in C 0 R + ; L 1 loc (R; [0, R]) , uniqueness is ensured by the results in [7]. However, the continuity of p at ρ L is not a necessary condition for uniqueness as we show in the following section.

Cases in which uniqueness holds
In this section we prove that: If (ρ L , ρ R ) ∈ C the corresponding classical solution satisfies 5 for all t > 0 sufficiently small and it is not possible to construct a different solution. If (ρ L , ρ R ) ∈ N the corresponding classical solution does not satisfy 5b, and there exists a unique nonclassical solution that satisfies 5. We list here two basic properties which will be of great help in the following case by case analysis.
Second, when the solution is nonclassical, due to the finite speed of propagation of the waves, the assumption (P2) and properties bp1 and bp2, we have Now we start the description of the possible cases and we proceed as follows. First, we show that for any initial datum satisfying (Ci), i = 1, . . . , 5, the problem actually has a unique solution and that the solution is classical. Second, we take into consideration the corresponding case (Ni), for which we prove that the classical solution is not suitable and that there exists a unique nonclassical solution.
In general the solutions to the constrained Riemann problem 5 are not selfsimilar. All the cases listed below describe self-similar solutions because we let the solutions evolve only on a small interval of time.
(N1) In this case [(t, x) → R[ρ L , ρ R ](x/t)] does not satisfy 5b because f (ρ R ) > p(ρ L +), see case (C1). Therefore, there does not exist any classical solution and we can consider only nonclassical solutions of the form 11. If p is continuous in ρ L , then by np4 we have thatp = p(ρ L ). If p experiences a jump at ρ L then, one may wonder which value in [p(ρ L +), p(ρ L −)] has to be chosen asp. As in the case (C1), the assumptions imply thatρ < ρ R and then thatρ(p) < ρ R andp ≤ f (ρ R ). Thenp is strictly smaller than f (ρ L ) andρ(p) > ρ L . As a consequence, property np2 forces us to choose the unique possible value ofp, which is p(ρ L +).
performs a shock with non negative speed σ(ρ L , ρ R ) and it satisfies 5b because f (ρ L ) ≤ p(ρ L +). Assume that there exists a nonclassical solution of the form 11. Observe that the assumptions ρ L < ρ R and f (ρ R ) ≥ f (ρ L ) together imply thatρ > ρ L . Then ρ(p) ≥ρ > ρ L and R[ρ L ,ρ(p)] is given by a shock with non positive speed if and only ifp ≤ f (ρ L ). Thusp ≤ f (ρ L ) ≤ p(ρ L +) and this implies by 11c thatp = f (ρ L ) = p(ρ L +) and that the nonclassical solution coincides with the classical one.
(C3) In this case [(t, x) → R[ρ L , ρ R ](x/t)] performs a possible null rarefaction on the right of the constraint and it satisfies 5b because f (ρ L ) ≤ p(ρ L +).
Assume that there exists a nonclassical solution of the form 11. Since ρ L ≤ ρ ≤ρ(p), R[ρ L ,ρ(p)] is given by a shock that has non positive speed if and only ifp ≤ f (ρ L ). Thereforep ≤ f (ρ L ) ≤ p(ρ L +) and this by 11c implies thatp = f (ρ L ) = p(ρ L +) and that the nonclassical solution coincides with the classical one.
does not satisfy 5b because f (ρ L ) > p(ρ L −), see case (C3). Therefore, there does not exist any classical solution and we can consider only nonclassical solutions of the form 11. By hypothesis and 11c we have f (ρ L ) > p(ρ L −) ≥p. Therefore ρ L ≤ρ <ρ(p) and by np2 we havep = p(ρ L +).
does not satisfy 5b because f (ρ) > p(ρ L −), see case (C4). Therefore, there does not exist any classical solution and we can consider only nonclassical solutions of the form 11.

Cases in which uniqueness is violated
Now we list the "pathological" cases, where we have more than one admissible solution. We stress once again that a necessary condition for non-uniqueness is . This is the situation considered in the Example 2 in [3]. Observe that such solutions are distinct as far as they correspond to distinct constraint levelsp, and that in any case there exist at least two distinct nonclassical solutions. (CNN5) Ifρ < ρ R ≤ ρ L , f (ρ R ) ≤ p(ρ L −), p(ρ L −) = p(ρ L +) and p(ρ L +) ≤ f (ρ L ), then the classical solution [(t, x) → R[ρ L , ρ R ](x/t)], which consists of a possible null rarefaction on the left of the constraint, as well as the nonclassical solutions of the form 11 corresponding top ∈ {p(ρ L +), f (ρ L )} satisfy 5. Observe that the two nonclassical solutions are distinct as far as they correspond to distinct constraint levelsp, and that in any case there exist at least two distinct solutions, one classical and one nonclassical. (NNN5) Ifρ < ρ R < ρ L , f (ρ R ) > p(ρ L −) ≥ f (ρ L ) ≥ p(ρ L +) and p(ρ L −) = p(ρ L +), then the nonclassical solutions of the form 11 corresponding top ∈ {p(ρ L +), f (ρ L ), p(ρ L −)} satisfy 5. Observe that such solutions are distinct as far as they correspond to distinct constraint levelsp, and that in any case there exist at least two distinct nonclassical solutions.

Proof of Proposition 3
(R1) Any solution given by R coincides on each side of the constraint with a solution given by the classical Riemann solver R. Therefore it satisfies the Rankine-Hugoniot jump condition along any of its discontinuities away from the constraint. Finally, by definition ofρ andρ, it satisfies the Rankine-Hugoniot jump condition also along the constraint.
(R2) It is clear by the proof of Proposition 2.
(R3) It is proved as in (R1) since any classical solution is in BV.
(R4) As was proved in [7], R is continuous on C ∪ N . If (ρ L , ρ R ) is not in C ∪ N then p experiences a jump at ξ = ρ L . Therefore, the local in time solutions of the Riemann problem for the initial conditions (ρ L + ε, ρ R ) and (ρ L − ε, ρ R ) are different and only one of the two converges to R [ρ L , ρ R ] as ε > 0 goes to zero.
(R5) We first stress once again that we can discuss the consistency property of our Riemann solvers only locally in time because, in general, the solutions may be not even self-similar globally in time. However, locally in time, the efficiency of the exit can be assumed to be constant and it is thus sufficient to proceed as in [7].
(R6) It is clear by the proof of Proposition 2.

Aknowledgements
All the authors are supported by French ANR JCJC grant CoToCoLa and Polonium 2014 (French-Polish cooperation program) No.331460NC. The second author is also supported by the Université de Frache-Comté, soutien aux EC 2014. The last author is also supported by ICM, University of Warsaw, and by Narodowe Centrum Nauki, grant 4140.