A Class of Well-Balanced Algorithms for Relativistic Fluids on a Schwarzschild Background

For the evolution of a compressible fluid in spherical symmetry on a Schwarzschild curved background, we design a class of well-balanced numerical algorithms up to third-order accuracy. We treat both the relativistic Burgers–Schwarzschild model and the relativistic Euler–Schwarzschild model and take advantage of the explicit or implicit forms available for the stationary solutions of these models. Our schemes follow the finite volume methodology and preserve the stationary solutions. Importantly, they allow us to investigate the global asymptotic behavior of such flows and determine the asymptotic behavior of the mass density and velocity field of the fluid.


Introduction
We are interested in the numerical approximation and the long time behavior of relativistic compressible fluid flows on a Schwarzschild black hole background. The flow is assumed to enjoy spherical symmetry and therefore we deal with nonlinear hyperbolic systems of partial differential equations (PDEs) in one space variable. This paper is twofold: on the one hand, designing and testing numerically finite volume algorithms that are well-balanced; on the other to perform a thorough investigation of the behavior of the solutions and numerically infer definite conclusions about the long-time behavior of such flows. Our study should provide first and useful insights for, on the one hand, further development concerning the mathematical analysis of the models and, on the other hand, further investigations to the same problem in higher dimensions without symmetry restriction.
Two models are of interest in the present paper. We treat first the relativistic Burgers-Schwarzschild model (as it is called in [15,16]): (1.1a) where v = v(t, r ) ∈ [−1, 1] is the unknown function and the flux and source terms read while the constant M > 0 represents the mass of the black hole. Obviously, the speed of propagation for this scalar balance law reads (1.2) which vanishes at the boundary r = 2M, so that no boundary condition is required in order to pose the Cauchy problem. Next, we consider the relativistic Euler-Schwarzschild model (as it is called in [15,16]): (1.3a) whose unknowns are the fluid density ρ = ρ(t, r ) ≥ 0 and the normalized velocity v = v(t, r ) ∈ (−1, 1). These functions are defined for all r > 2M and the limiting values v = ±1 being reached at the boundary r = 2M only, and (1.3d) Here, k ∈ (−1, 1) denotes the (constant) speed of sound, and this second model can be checked to be a strictly hyperbolic system. The eigenvalues of the corresponding Jacobian of the flux function read As usual, a state (ρ, v), by definition, is said to be sonic if one of the eigenvalues vanishes, i.e. if |v| = k, supersonic if both eigenvalues have the same sign, i.e. if |v| > k, or subsonic if the eigenvalues have different signs, i.e. if |v| < k. We will need to distinguish between these regimes in order to design a robust scheme for this model. Both eigenvalues μ ± vanish at the boundary r = 2M, so that no boundary condition is required in order to pose the Cauchy problem.
In order to be able of running reliable and accurate numerical simulations for these two models, in this paper we design shock-capturing, high-order, and well-balanced finite volume methods of first-and second-order of accuracy (and even third-order accurate for (1.1)). Specifically, we extend to the present problem the well-balanced methodology proposed recently by Castro and Parés [7] for nonlinear hyperbolic systems of balance laws. For earlier work on well-balanced schemes we also refer to [5,19,20] and, concerning the design of geometry-preserving schemes, we refer for instance to [1][2][3]6,[9][10][11]18,22] and the references therein.
The properties of the stationary solutions play a fundamental role in the design of well balanced schemes, as well as in the study of the long time behavior of solutions. We thus also built here upon earlier investigations by LeFloch and collaborators [14][15][16] on the theory and approximation of the relativistic Burgers-and Euler-Schwarzschild model (1.1) and (1.3). Remarkably, the stationary solutions to both models are available in explicit or implicit form.
An outline of the content of this paper is as follows.
In Sect. 2 we describe the methodology for this paper and indicate the challenges met with the two models above. The actual design of the schemes is the content of Sect. 3 (Burgers equation) and Sect. 5 (Euler equations). The proposed well-balanced algorithms, by construction, preserve all of the steady state solutions, which is an essential property since numerical methods without this property may lead to wrong conclusions. Furthermore, for each model we investigate the efficiency, accuracy, and robustness of the proposed algorithms first in Sect. 4 (Burgers equation) and Sect. 6 (Euler equations). Our numerical experiments below make comparisons between well-balanced and standard methods, and we demonstrate that the proposed schemes are numerically well-balanced and we emphasize the importance of this property in order to reach reliable results. Furthermore, we study the late time behavior of solutions to both models and discuss the role of the choice of the value of the initial data at the boundary. Finally, we also describe how steady shocks behave under small (or large) perturbations. In short, we demonstrate that the global dynamics can be accurately determined by the proposed algorithms and we reach some definite conclusions in Sects. 4.6 and 6.5, respectively. Further details concerning our methodology and conclusions are found in the corresponding sections for each model. with unknown V = V (t, r ) ∈ R N and N = 1 or 2. Systems of this form have non-trivial stationary solutions, which satisfy the ODE Our goal is to introduce a family of numerical methods that are well-balanced, i.e. that preserve the stationary solutions in a sense to be specified. We follow the strategy in [7] to which we refer for further details and arguments of proof. We consider semi-discrete finite volume numerical methods of the form 3) and the following notation is used: ] denote the computational cells, whose length r is taken to be a constant for the sake of simplicity in the presentation.
• V i (t) denotes the approximate average of the exact solution in the ith cell at the time t, that is, • P t i (r ) denotes the approximate solution in the ith cell, as given by a high-order reconstruction operator based on the cell averages Here, S i denotes the set of cell indices associated with the stencil of the ith cell.
• The flux terms are denoted by are the reconstructed states at the interfaces, i.e.
Here, F is a consistent numerical flux, i.e. a continuous function F : Furthermore, given a stationary solution V * of (2.1), we use the following terminology.
• The numerical method (2.3) is said to be well-balanced for V * if the vector of cell averages of V * is an equilibrium of the ODE system (2.3). • The reconstruction operator is said to be well-balanced for V * if, for every i, we have where P i is the approximation of V * obtained by applying the reconstruction operator to the vector of cell averages of V * .
It is easily checked that, if the reconstruction operator is well-balanced for a continuous stationary solution V * of (2.1) then the numerical method is also well-balanced for V * . The following strategy to design a well-balanced reconstruction operator P i on the basis of a standard operator Q i was introduced in ] we proceed as follows.
1. Seek, (whenever possible), a stationary solution V * i (x) defined in the stencil of cell I i (∪ j∈S i I j ) such that 1 r If such a solution does not exist, take V * i ≡ 0.
2. Apply the reconstruction operator to the cell values {W j } j∈S i given by It can be then easily shown that the reconstruction operator P i in (2.8) is well-balanced for every stationary solution provided that the reconstruction operator Q i is exact for the zero function. Moreover, if Q i is conservative then P i is conservative, in the sense that 1 r and P i has the same accuracy as Q i if the stationary solutions are sufficiently regular. Observe that, if there does not exists a stationary solution defined in the stencil and satisfying (2.6) then the standard reconstruction is used. Observe that this choice does not spoil the well-balanced character of the numerical method since, in this case, the cell values in the stencil cannot be the averages of a stationary solution (otherwise, there would be at least one solution V * i ) and therefore there exists no local equilibrium which would be required to preserve. On the other hand, if there exist more than one stationary solution defined on the stencil and satisfying (2.6), a criterion is needed in order to select one of them and indeed this sometimes happen for the Euler-Schwarzschild model of interest in the present paper. The relevant criterion in this regime depends on the particular problem and, in the case of the Euler-Schwarzschild model, we propose to distinguish between the (subsonic or supersonic) regimes of the flow, as we will explain later in Sect. 5.
If a quadrature formula (whose order of accuracy must be greater or equal to the one of the reconstruction operator) where α 0 , . . . , α q , r i,0 , . . . , r i,q represent the weights and the nodes of the formula, is used to compute the averages of the initial condition, namely V i,0 = q l=0 α l V 0 (r i,l ), the reconstruction procedure has to be modified to preserve the well-balanced property: Steps 1 and 2 have to be replaced by the following ones.
If this solution does not exist, take V * i ≡ 0. 2. Apply the reconstruction operator to the cell values {W j } j∈S i given by For first-or second-order methods, if the midpoint rule is selected to compute the initial averages, i.e. V i,0 = V 0 (r i ), then at the first step of the reconstruction procedure, the problem (2.10) reduces to finding the stationary solution satisfying The well-balanced property of the method can be lost if the quadrature formula is used to compute the integral appearing at the right-hand side of (2.3). In order to circumvent this difficulty, in [7] it is proposed to rewrite the methods as follows: where V t, * i is the function selected in Step 1 for the ith cell at time t. In this equivalent form, a quadrature formula can be applied to the integral without losing the well-balanced property, and this leads to a numerical method of the form: (2.13) First-order well-balanced methods are obtained by selecting the trivial constant piecewise reconstruction operator as the standard one, i.e. ]. (2.14) It can be easily checked that the numerical method then reduces to Further strategy for this paper. When applying the well-balanced reconstruction methodology, a main difficulty comes from the first step, in which we need to find a stationary solution (i.e. a solution of the ODE system (2.2)) satisfying (2.6), (2.10), or (2.11), and this depends upon the way that the relevant integrals are computed. This amounts, in general, to a non-trivial, nonlinear problem whose analysis and solution depend on the system of balance laws under consideration. For the Burgers-Schwarzschild model, the explicit form of the general solution of the ODE (2.2) is available (as we recall in Sect. 3) and the nonlinear problem to be solved in Step 1 has always a unique solution v * i . Nevertheless, it may not be possible to extend this solution to the whole stencil and, in such a case, the standard reconstruction operator must then be used.
For the Euler-Schwarzschild model, only first and second-order methods will be considered in the present paper, and the mid-point rule will be used in order to numerically compute the relevant integrals. As a consequence, only solutions of the ODE system (2.2) satisfying (2.11) (i.e. solutions of a standard Cauchy problem) should be found in Step 1. Now, the general solution of the ODE system is available in implicit form (as we will recall in Sect. 5) and, for this system, the Cauchy problem (2.2), (2.11) may have a single solution or no solutions as well as two or more solutions. If there exists no solution, the standard reconstruction is used. If there are two or more solutions, a criterion based on the (subsonic, supersonic) nature of the fluid flow will be used in order to select one of them, as mentioned above.
This methodology can be extended to systems of balance laws for which the solutions to the ODE system (2.2) are not available neither in explicit or implicit form. For such system, the nonlinear problems arising in Step 1 must be solved numerically. For instance, in [12] a control-based strategy combined with a standard ODE solver is used and solutions of (2.2) are computed that satisfy averaging conditions like (2.6) or (2.10).
Concerning the extension to multidimensional problems, the main challenge is again the problem in Step 1, which now consists in finding a solution to the nonlinear system of PDEs satisfied by the stationary solutions with prescribed average in the cell or from values at a point. Such a problem, clearly, is much more difficult to solve (either exactly or numerically) than an ODE system. Moreover, there may exist infinitely many stationary solutions satisfying the average or point value conditions. Nevertheless, the methodology can be still applied and allows one to design numerical methods that preserve a certain family of stationary solutions. In Step 1, the stationary solution belonging to the prescribed family that satisfies the imposed condition would be selected or the one that is closer in some sense to the cell values in the stencil. For instance, it is possible to apply this methodology and design numerical methods for the Euler-Schwarzschild model in three spatial dimension and impose spherical-symmetric stationary solutions are preserved.

Preliminaries
For the Burgers-Schwarzschild equation (1.1), the steady state solutions are of the form In Fig. 1 we plot the steady solutions for several values of K 2 . The domain of definition of these stationary solutions is It can be easily checked that, given a pair (K , r * ) such that r * ∈ D K , the discontinuous function defined in D K by is an entropy weak stationary solution of the Burgers-Schwarzschild model. (For further properties, see the study in [15,16].)

First-order Method
If the midpoint rule is used to compute the initial averages, at the first step of the reconstruction procedure one has to search for K 2 i such that There is always a unique solution K i given by In order to apply the numerical method (2.12), this stationary solution has to be computed at r i± 1 2 and this is only possible if r i+ 1 2 ∈ D K i , that is, provided If this condition is satisfied, then the numerical method (2.15) can be used. If this condition is not satisfied, then the standard trivial reconstruction is considered, i.e. Q i (r ) = v i . The numerical method writes then as follows: where ).
Summing up, the expression of the semi-discrete first-order method reads The forward Euler method is used for the time discretization. Furthermore, we emphasize that, if (3.7) is not satisfied, then v i cannot be the point value of a stationary solution defined in the computational domain, so that the use of the standard reconstruction does not destroy the well-balanced property of the method, since in this case there is no stationary solution to preserve.

Second-order Method
Let us suppose again that the midpoint rule is used to compute the cell averages and that the minmod reconstruction operator is considered: see [21]. The stationary solution v * i selected at the first stage of the reconstruction procedure is again (3.6) withK i given by (3.5). In order to compute the reconstructions, this stationary solution has to be computed at the points r i−1 , , r i+1 so that the following condition has to be satisfied r i+1 ∈ D K i , i.e.
If this condition is satisfied, the following step of the reconstruction procedure consists in computing the fluctuations: Then the reconstruction is defined as (3.14) Once the well-balanced reconstruction operator has been computed, the form (2.12) of the numerical method is considered and the midpoint rule is used again to approximate the integral. Observe however that, in this case: Therefore, the expression (2.13) reduces again to (2.15) with ).
If (3.11) is not satisfied, then the standard MUSCL reconstruction is applied, namely The expression of the numerical method is given again by (3.9)-(3.10) with the difference that the second-order reconstructions are used now to compute the numerical fluxes. The TVDRK2 method [13] is used in order to discretize the equations in time.
Observe that, according to the well-balanced reconstruction procedure described in the previous section, the fluctuations to be reconstructed should be in this case but in (3.12) the sign of v i has been replaced by that of v j . This modification allows one to preserve not only the continuous stationary solutions solution but also the discontinuous stationary solutions of the family (3.3).

Third-order Method
In order to design a third-order method the CWENO reconstruction of order 3 (see [8], [17]) will be considered and the two-point Gauss quadrature will be used to compute the initial averages and the integrals coming from the source term, namely v 0 Therefore, at the first step of the reconstruction procedure one has to look for K 2 i such that: If we define the function it is easily checked that g is a positive decreasing function defined in the interval [0, Therefore there are two possibilities: has to be computed in the points {r i−1,0 , r i−1,1 , r i+1,0 , r i+1,1 } in the reconstruction procedure. Therefore, these points have to be in the interval of definition of v * i , and this happens if r i+1,1 ∈ D K i . Therefore, the well-balanced reconstruction can be computed only if the following condition is satisfied: If this condition is satisfied, the fluctuations can be then computed: and the well-balanced reconstruction is given by The expression of the semi-discrete method is finally (3.9) where the numerical fluxes are computed at the reconstructed states and otherwise. (3.18) The TVDRK3 method of order 3 [13] will be used for the time discretization.

Preserving the Exact Averages of the Stationary Solutions
The three methods presented above can be modified to preserve the exact averages of the stationary solutions instead of its approximation computed with a quadrature formula. To do this, the problem to be solved at the first stage of the well-balanced reconstruction procedure is the following one: find K 2 i such that: If we define the function it can be easily checked that it is a decreasing function defined in [0, The explicit expression of g can be obtained: is a primitive function of 1 − x 1 − 2M r . Therefore g(K e,i ) can be explicitly computed. The well-balanced reconstruction can thus be computed if (3.20) is satisfied and the cells of the stencil S i are contained in the domain of definition D K i of the corresponding stationary solution. Otherwise, the standard reconstruction is applied. The expression of the numerical methods is the same as the ones above.

Preliminaries
In this section several numerical tests are applied to check the performance of the wellbalanced numerical methods introduced in the previous section. We consider the spatial interval [2M, L] with M = 1 and L = 4, a 256-point uniform mesh and the CFL number is set to 0.5. At r = 2M we impose F − 1 2 = 0 as boundary condition since 1 − 2M r = 0. At r = L we will use a transmissive boundary condition using ghost-cells if the initial condition is not a stationary solution; otherwise, the stationary solution is imposed in the ghost-cells. The following numerical flux will be used: is the self-similar solution of the Riemann problem for the standard Burgers equation with the initial condition In order to check the relevance of the well-balanced property, these methods will be compared with those based on the same numerical fluxes and the standard first-, second-, or third-order reconstructions.

Stationary Solutions
Positive stationary solution. We consider the initial condition v 0 (r ) = 3 4 + 1 2r (4.1)   Table 1 shows the error in L 1 norm between the initial condition and the numerical solution at time t = 50. Figure 2 compares the numerical solutions obtained with the well-balanced and the non-well-balanced methods: it can be seen how the latter are unable to capture the stationary solution. After a time that decreases with the order, the numerical solutions depart from the steady state. Negative stationary solution. Let us consider now as initial condition the negative stationary condition corresponding to K = 1 2 : Figure 3 shows the numerical solutions obtained with the different numerical methods.
Observe that the scale of the vertical axis is not the same as the one in Fig. 2: it has been changed so that the difference between the numerical solutions can be better seen. Table 2 shows the error in L 1 norm between the initial condition and the numerical solution at time t = 50. According to Fig. 3 and Table 2 we need more time to see the differences between the well-balanced and non-well-balanced schemes of order 1, 2 and 3 but the errors are again much smaller with the well-balanced schemes for this test. In this case we need more time to see these differences since this negative stationary solution is close to the constant state v(r ) = −1 where it seems that the non-well-balanced schemes converge.
that is a stationary entropy weak solution of the family (3.3). Table 3 shows the error in L 1 norm between the initial condition and the numerical solution at time t = 50. Figure 4 shows the differences between the numerical solutions obtained with well-balanced and non-wellbalanced methods: again the latter depart from the stationary solution at times that decrease with the order. The numerical results obtained for the equation with initial conditions (4.1), (4.2), and (4.3) clearly show the need of using well-balanced methods for this equation.

Moving Shocks Connecting two Steady Profiles
Right-moving shock. We consider now the initial condition The corresponding solution consists of a right-moving shock connecting two branches of stationary solutions. Figure 5 shows the numerical solutions obtained with the first-, second-, and third-order well-balanced methods and a reference solution computed with the first-order standard method using a mesh of 10000 cells. As it can be seen, the well-balanced methods capture correctly the shock with a resolution that increases with the order as expected. Left-moving shock. Similar conclusions can be drawn for the left-moving shock linking two branches of stationary solutions that generates from the initial condition: see Fig. 6. A reference solution computed with the first-order standard method has been computed again using a mesh of 10000 cells.

Perturbation of a Steady Shock Solution
Left-hand perturbation. In this test case we consider the initial condition: where v 0 is the steady shock solution given by (4.7) The first-, second-, and third-order well-balanced methods have been applied to this problem. In Fig. 7 it can be observed that, after the wave generated by the initial perturbation leaves the : with as displayed in Fig. 8. In this case we have used a 2000-point uniform mesh since the displacement of the shock is smaller in this case and more points in the mesh are needed in order to see that the steady shock is not recovered.
Left-hand perturbation with zero average. Now we consider an initial condition of the form (4.6) with a perturbation p L such that p L (r )dr = 0, in particular:  Left-hand and right-hand perturbations with zero average. In order to study the relation between the amplitude of the perturbation and the distance between the initial and the final stationary shocks, we consider the initial condition: where v 0 is the steady shock solution given by (4.3) and ( p L (r ) + p R (r ))dr = 0. In particular we take p L (r ) as in (4.7) and p R (r ) as in (4.9). In Fig. 11 it can be observed that, after the wave generated by the initial perturbation leaves the computational domain, the stationary solution (4.3) is not recovered: a different stationary solution of the family (3.3) with the shock placed at a different location is reached. This is a natural result since, as we saw Relation between the perturbation and the displacement of the shock. In order to study the relationship between the amplitude of the perturbation and the distance between the initial and the final shock locations, we consider the family of initial conditions: where v 0 is given again by (4.3) and with α > 0. The amplitude of the perturbation is measured by δ v (α, r ) dr and the distance between the shocks is measured by lim t→∞ |v(r , t) − v 0 (r )| dr. See Fig. 12. Table 4 and Fig. 13 show the relationship between those magnitudes that is clearly linear.

Long-time Behavior of the Solutions
In this section we consider different initial conditions and investigate the long-time behavior of the corresponding solutions using the first-order well-balanced scheme. A large number of tests have been performed with the first-order methods (that is the less costly one) consid-

the initial stationary solution is recovered. 2. If
In view of Figs. 14, 15, 16 and 17 we have reached the following.

Preliminaries
In the case of the Euler-Schwarzschild equations (1.3), the stationary solutions are implicitly given by the equations: where C 1 , C 2 are constants. The pair (v, ρ) of a stationary solution satisfies the following ODE system analyzed first in [15,16]:  Figure 18 shows the graph of v for some of them. When these functions are defined in (2M, ∞), they have a maximum or a minimum in that comes from solving dv dr = 0. In Fig. 18 the stationary solutions marked in red are those that take the value v = ±k at r = r c .
Given two constants C 1 and C 2 , in order to compute the stationary solution given by (5.1) in a point r = a, the following nonlinear system has to be solved: It is enough thus to solve, if it is possible, the nonlinear equation to compute v. Once this equation is solved, ρ is computed using the second equation of (5.5). It can be easily checked that g satisfies: Moreover, g is strictly monotone in [−1, −k), (−k, k), and (k, 1]. As a consequence we have the following conclusion. (For further properties, see the study in [15,16].) • If |K a | > g(k), the equation (5.6) has no solution, i.e. a stationary solution given by C 1 and C 2 cannot be defined at r = a. • If |K a | = g(k) then the equation (5.6) has only one solution (k if K a > 0, −k if K a < 0). Therefore, (5.5) has only one solution that is a sonic state. • Otherwise, (5.6) has two possible solutions. Therefore there are two states (ρ, v) that solve (5.5), one supersonic and one subsonic.
For the sake of clarity, together with the representation for the states, we will use V = [ρ, v] T . Here, V can be easily computed from V and, given V , V is also easily computed by (1.3d) that comes from solving a second-degree equation.

First-order Method
If the midpoint rule is used to compute the initial averages, given a family of cell values V i , in the first step of the well-balanced reconstruction procedure one has to find, if it is possible, a stationary solution Obviously such a stationary solution would correspond to the choice of constants: According to the discussion above, the corresponding stationary solution is defined in r i± 1 2 provided that: , one subsonic and one supersonic. Therefore, a criterion is needed to select one or the other. The following criterion will be used here: • if V i is not sonic, then the state whose regime (sub or supersonic) is the same as ).
• if V i is sonic, then the state whose regime is the same as V i+1 is selected for V * i (r i+ 1 2 ) and the state whose regime is the same as ).
Observe that this criterion aims to preserve the regime of the given cell values.
If condition (5.9) is satisfied, then the numerical method (2.15) is used. Otherwise the standard trivial reconstruction is considered.
The expression of the semi-discrete first-order method is then as follows: where The forward Euler method is used again for the time discretization.

Second-order Method
Let us use again the midpoint rule to compute cell averages and the minmod reconstruction operator. The stationary solution sought at the first stage of the well-balanced reconstruction procedure is again characterized by the constants (5.8). This time, this stationary solution has to be computed at the points r i−1 , r i− 1 2 , r i+ 1 2 , r i+1 so that the following condition has to be satisfied: where K i± 1 2 are given by (5.9) and If this condition is satisfied, the reconstruction is defined as follows: where the minmod function is applied component by component and Observe that the conserved variables V are used in the reconstruction procedure. On the other hand, if (5.12) is not satisfied, then the standard MUSCL reconstruction is applied: The expression of the numerical method is given again by (5.10)-(5.11) with the difference that the second-order reconstructions are used now to compute the numerical fluxes. The TVDRK2 method is used now to discretize the equations in time.

Preliminaries
In this section several tests are considered to check the performance of the first-and secondorder well-balanced numerical methods for Euler-Schwarzschild model introduced in the previous section. We consider the spatial interval In order to test the dependency of the results on the numerical method, two different first-order numerical fluxes are considered: the Lax-Friedrichs numerical flux and a HLL-like numerical flux in PVM form (see [4]): where λ 1 and λ 2 are the eigenvalues of some intermediate matrix J i+ 1 2 of the form where v m is some intermediate value between v n i and v n i+1 .
Given two states V L and V R , in order to choose an adequate intermediate value v m , we look for v such that the following Roe-type property is satisfied: i.e. the factor (1 − 2M/r ) is neglected for simplicity. Due to the form of the matrix, it is enough to find v such that This equality is equivalent to a second-order degree equation for v, namely αv 2 +βv+γ = 0, where 2 is always positive, there are always two real solutions: and it can be proven that: Finally, in the case α = 0 and V R = V L , we take v m = − γ β and in the case Once v m has been chosen, the expression of λ j , j = 1, 2 in (6.3) is as follows: Since for a 2-systems HLL and Roe methods are equivalent and the intermediate value chosen to compute the wave speeds satisfies a Roe-type property, this numerical flux will be called Roe-type numerical flux in what follows. The proposed numerical method will be compared with those based on the same numerical flux and the standard first-and secondorder reconstructions.    Table 6 shows the error in L 1 norm between the numerical solution at time t = 50. Figs. 21,22 show the difference between the numerical results given by well-balanced and non-wellbalanced methods using the Roe-type numerical flux. Again the numerical solutions obtained with non-well-balanced methods depart from the initial steady state. Discontinuous stationary entropy weak solution. We consider finally the initial condition where V * − (r ) is the supersonic stationary solution such that ρ * − (6) = 4, v * − (6) = 0.6 (6.9)  and V * + (r ) is the subsonic one such that . (6.10) V 0 is an entropy weak stationary solution of the system: see [15,16]. Table 7 shows the error in L 1 norm between the numerical solution at time t = 50 and Figs. 23, 24 show the comparison of the numerical results obtained with well-balanced and non-well-balanced methods at selected times. On the other hand, the numerical results of this section put on evidence, as for the Burgers-Schwarzschild system, the relevance of using well-balanced methods for the Euler-Schwarzschild model.

Perturbation of a Regular Stationary Solution
In this test we consider the initial condition V 0 (r ) = V * (r ) + δ(r ), (6.11) where V * is the supersonic stationary solution It can be observed in Fig. 25 that the stationary solution V * is recovered once the perturbation has left the domain. In this Figure, the numerical results obtained with the first-and secondorder well-balanced methods are compared with a reference solution computed using the first-order well-balanced method with a 5000-point mesh (the Roe-type numerical flux is used again). As expected, the second-order method is less diffusive.

Perturbation of a Steady Shock Solution
Left-hand perturbation. We consider the initial condition and V * (r ) is the stationary solution given by (6.8)-(6.10). In Figs. 26 and 27 the numerical results obtained with the first-and second-order well-balanced methods using the Lax-Friedrichs and the Roe-type numerical methods with different meshes are compared. As it happened for the Burgers-Schwarzschild model, the location of the stationary shock changes after the passage of the wave generated by the perturbation. Nevertheless in this case the displacement of the shock is slower. Different numerical methods have been applied to check the dependency of the motion on the scheme: although the evolution of the shock slightly depends on the number of points of the mesh, all the numerical solutions capture the same final location of the shock. In Fig. 28 the evolution of the shock given by the first-order WB method with different number of cells are compared. The location of the shock at every time step has been detected by using the condition Right-hand perturbation. Let us consider now two different initial conditions: on the one hand where V * (r ) is again the discontinuous stationary solution given by (6.8)-(6.10) and On the other hand, where δ R is given again by (6.16) and V * 2 (r ) is the steady shock of the form (6.8) satisfying ρ * − (6) = 5, v * − (6) = 0.6.   Observe that the definition of v is identical for both stationary solutions but ρ is different. After the passage of the perturbation, the shock starts moving leftward and, in both cases, the numerical solution converges to a smooth transonic stationary solution of the form: V * (r ) = V * − (r ), r ≤ r c , V * + (r ), otherwise, (6.20) where r c is given by (5.4); V * − (r ) and V * + (r ) are respectively a subsonic and a supersonic stationary solution satisfying v * ± (r c ) = −k: see Fig. 29. Nevertheless, the limits in time of   where V * is the steady shock solution given by (6.8)-(6.10) and with α > 0. In this case we will also use the Roe-type numerical flux and a 2000-point uniform mesh. Figures 31 and 32 show the numerical solution for different values of α and we observe that depending on the amplitude of the perturbation the numerical solutions converge in time to different steady shock solutions. The amplitude of the perturbation is measured with δ v (α, r ) dr and the distance between the shocks is measured by lim t→∞ |v(r , t) − v * (r )| dr, as we did for the Burgers-Schwarzschild model. Table 8 and Fig. 33 show the relationship between those magnitudes: the displacement of the shock seems to grow exponentially with the amplitude. Let us finally consider a family of initial conditions that generate leftward displacement of the initial steady shock: V 0 (r ) = V * (r ) + δ(β, r ), (6.23) where V * is again the steady shock solution given by (6.8)-(6.10) and with β < 0. In this case we will use the Roe-type numerical flux and a 2000-point uniform mesh. Figures 34 and 35 show the numerical solution for different values of β and we observe that it converges to the same stationary solution regardless of the perturbation.

Main Conclusions for the Euler-Schwarzschild Model
In view of Fig. 25 we arrive at the following observation. (1.3) is perturbed, the solution is restored once the wave generated by the perturbation goes away.  Table 8 we arrive at the following observation.
1. If the perturbation moves the steady shock to the right, then a different stationary solution of the form V (r ) = V * − (r ), r ≤ r 1 , V * + (r ), otherwise, with r 0 = r 1 , is obtained, and the distance between r 0 and r 1 seems to depend exponentially on the amplitude of the perturbation: see Table 8 and Fig. 33. 2. If the perturbation moves the steady shock to the left, then a steady shock solution of the form (6.20) is obtained.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.