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

Abstract

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.

1 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]):

$$\begin{aligned} v_t + F(v,r)_r = S(v,r), \quad t \ge 0, \quad r>2M, \quad \end{aligned}$$

(1.1a)

where $v=v(t,r) \in [-1,1]$ is the unknown function and the flux and source terms read

$$\begin{aligned} F(v,r) = \Big ( 1-\frac{2M}{r} \Big ) \frac{v^2-1}{2}, \quad S(v,r)= \frac{2M}{r^2} \, (v^{2}-1), \end{aligned}$$

(1.1b)

while the constant $M>0$ represents the mass of the black hole. Obviously, the speed of propagation for this scalar balance law reads

$$\begin{aligned} \partial _v F(v,r) = \Big ( 1-\frac{2M}{r} \Big ) v, \end{aligned}$$

(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]):

$$\begin{aligned} V_{t}+ F(V,r)_{r} = S(V,r), \quad t \ge 0, \quad r>2M, \quad \end{aligned}$$

(1.3a)

whose unknowns are the fluid density $\rho =\rho (t,r) \ge 0$ and the normalized velocity $v=v(t,r) \in (-1,1)$. These functions are defined for all $r>2M$ and the limiting values $v=\pm 1$ being reached at the boundary $r=2M$ only, and

$$\begin{aligned} V= & {} \left( \begin{array}{c} V^{0}\\ V^{1} \end{array}\right) = \left( \begin{array}{c} \displaystyle \frac{1+k^{2}v^{2}}{1-v^{2}}\rho \\ \displaystyle \frac{1+k^{2}}{1-v^{2}}\rho v\\ \end{array}\right) , \quad F(V,r)= \left( \begin{array}{c} \displaystyle \Big ( 1-\frac{2M}{r} \Big )\frac{1+k^{2}}{1-v^{2}}\rho v \\ \displaystyle \Big ( 1-\frac{2M}{r} \Big )\frac{v^{2}+k^{2}}{1-v^{2}}\rho \end{array}\right) , \end{aligned}$$

(1.3b)

$$\begin{aligned} S(V,r)= & {} \left( \begin{array}{c} \displaystyle -\frac{2}{r}\Big ( 1-\frac{2M}{r} \Big )\frac{1+k^{2}}{1-v^{2}}\rho v \\ \displaystyle \frac{-2r+5M}{r^{2}}\frac{v^{2}+k^{2}}{1-v^{2}}\rho -\frac{M}{r^{2}} \frac{1+k^{2}v^{2}}{1-v^{2}}\rho +2\frac{r-2M}{r^{2}}k^{2}\rho \end{array}\right) , \end{aligned}$$

(1.3c)

with

$$\begin{aligned} v=\frac{1+k^{2}-\sqrt{(1+k^{2})^{2}-4k^{2}\left( \frac{V^{1}}{V^{0}} \right) ^{2}}}{2k^{2}\frac{V^{1}}{V^{0}}}, \quad \rho =\frac{V^{1}(1-v^{2})}{v(1+k^{2})}. \end{aligned}$$

(1.3d)

Here, $k\in (-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

$$\begin{aligned} \mu _{\pm } = \Big ( 1-\frac{2M}{r} \Big )\frac{v \pm k}{1 \pm k^2 v}. \end{aligned}$$

(1.4)

As usual, a state $(\rho , 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 $\mu _{\pm }$ 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.

2 Well-Balanced Methodology and Strategy Proposed in this Paper

General methodology. Both problems of interest are of the form

$$\begin{aligned} V_t + F(V,r)_r = S(V,r), \quad r > 2M, \end{aligned}$$

(2.1)

with unknown $V =V(t,r)\in {\mathbb {R}}^N$ and $N = 1$ or 2. Systems of this form have non-trivial stationary solutions, which satisfy the ODE

$$\begin{aligned} F(V,r)_r = S(V,r). \end{aligned}$$

(2.2)

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

$$\begin{aligned} \frac{dV_{i}}{dt}= -\frac{1}{\Delta r} \Big ( F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}} - \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} S({\mathbb {P}}^t_{i}(r), r)\,dr \Big ), \end{aligned}$$

(2.3)

and the following notation is used:

• $I_i= [r_{i-\frac{1}{2}},r_{i+\frac{1}{2}}]$ denote the computational cells, whose length $\Delta 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,

  $$\begin{aligned} V_i(t) \cong \frac{1}{\Delta r} \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} V(r,t) \, dr. \end{aligned}$$

  (2.4)

• ${\mathbb {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 $\{ V_i(t) \}$, that is, ${\mathbb {P}}^t_{i}(r) = {\mathbb {P}}^t_i(r; \{V_j(t)\}_{j \in {\mathcal {S}}_i})$. Here, ${\mathcal {S}}_i$ denotes the set of cell indices associated with the stencil of the ith cell.

• The flux terms are denoted by $F_{i+\frac{1}{2}} = {\mathbb {F}}(V_{i+\frac{1}{2}}^{t, -}, V_{i+\frac{1}{2}}^{t, +}, r_{i+\frac{1}{2}})$, where $V_{i+\frac{1}{2}}^{t, \pm }$ are the reconstructed states at the interfaces, i.e.

  $$\begin{aligned} V_{i+\frac{1}{2}}^{t, -}={\mathbb {P}}_i^t(r_{i+\frac{1}{2}}), \quad V_{i+\frac{1}{2}}^{t, +}={\mathbb {P}}_{i+1}^t(r_{i+\frac{1}{2}}). \end{aligned}$$

  (2.5)

  Here, ${\mathbb {F}}$ is a consistent numerical flux, i.e. a continuous function ${\mathbb {F}}: {\mathbb {R}}^N \times {\mathbb {R}}^N \times (2M, + \infty ) \rightarrow {\mathbb {R}}^N$ satisfying ${\mathbb {F}}(V,V,r) = F(V,r)$ for all V, r.

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 ${\mathbb {P}}_i(r)=V^*(r)$ for all $r \in [r_{i-\frac{1}{2}}, r_{i+\frac{1}{2}}]$, where ${\mathbb {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 ${\mathbb {P}}_i$ on the basis of a standard operator ${\mathbb {Q}}_i$ was introduced in [5]:

Given a family of cell values $\{V_i\}$, in every cell $I_i=[r_{i-\frac{1}{2}}, r_{i+\frac{1}{2}}]$ we proceed as follows.

• 1. Seek, (whenever possible), a stationary solution $V_i^*(x)$ defined in the stencil of cell $I_i$ ($\cup _{j\in {\mathcal {S}} _{i}}I_j$) such that

  $$\begin{aligned} \frac{1}{\Delta r} \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} V_i^* (r) \, dr = V_i. \end{aligned}$$

  (2.6)

  If such a solution does not exist, take $V^*_i \equiv 0$.

• 2. Apply the reconstruction operator to the cell values $\{W_j\}_{j \in {\mathcal {S}}_i}$ given by

  $$\begin{aligned} W_j= V_j - \frac{1}{\Delta r} \int _{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}} V_i^* (r) \, dr, \quad j \in {\mathcal {S}}_i, \end{aligned}$$

  (2.7)

  in order to obtain ${\mathbb {Q}}_i(r)={\mathbb {Q}}_i(r;\{W_j\}_{j \in {\mathcal {S}}_i}).$

• 3. Define finally

  $$\begin{aligned} {\mathbb {P}}_i(r)=V_i^*(r)+{\mathbb {Q}}_i(r). \end{aligned}$$

  (2.8)

It can be then easily shown that the reconstruction operator ${\mathbb {P}}_i$ in (2.8) is well-balanced for every stationary solution provided that the reconstruction operator ${\mathbb {Q}}_i$ is exact for the zero function. Moreover, if ${\mathbb {Q}}_i$ is conservative then ${\mathbb {P}}_i$ is conservative, in the sense that

$$\begin{aligned} \frac{1}{\Delta r} \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}}{\mathbb {P}}_i (r) \, dr = V_i, \end{aligned}$$

(2.9)

and ${\mathbb {P}}_i$ has the same accuracy as ${\mathbb {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)

$$\begin{aligned} \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} f(x) \,dx \approx \Delta r\sum _{l=0}^q \alpha _l f(r_{i,l}) \end{aligned}$$

where $\alpha _0, \dots , \alpha _q$, $r_{i,0}, \dots , 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} = \sum _{l=0}^q \alpha _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.

• 1. Seek, if possible, the stationary solution $V_i^*(x)$ defined in the stencil of cell $I_i$ ($\cup _{j\in {\mathcal {S}} _{i}}I_j$) such that

  $$\begin{aligned} \sum _{l=0}^q \alpha _l V_i^* (r_{i,l}) = V_i. \end{aligned}$$

  (2.10)

  If this solution does not exist, take $V^*_i \equiv 0$.

• 2. Apply the reconstruction operator to the cell values $\{W_j\}_{j \in {\mathcal {S}}_i}$ given by

  $$\begin{aligned} W_j= V_j - \sum _{l=0}^q \alpha _l V_i^* (r_{j,l}), \quad j \in {\mathcal {S}}_i. \end{aligned}$$

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

$$\begin{aligned} V^*_i(r_i) = V_i. \end{aligned}$$

(2.11)

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:

$$\begin{aligned} \begin{aligned} \frac{dV_i}{dt}&=-\frac{1}{\Delta r} \left( F_{i+\frac{1}{2}}-F\left( V_i^{t,*}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}}\right) - F_{i-\frac{1}{2}} +F\left( V_i^{t,*}(r_{i-\frac{1}{2}}), r_{i-\frac{1}{2}}\right) \right) \\&\quad + \frac{1}{\Delta r} \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} \left( S({\mathbb {P}}_i^t (r), r) -S(V_i^{t,*}(r), r) \right) \, dr, \end{aligned} \end{aligned}$$

(2.12)

where $V_i^{t,*}$ 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:

$$\begin{aligned} \begin{aligned} \frac{dV_i}{dt}&=-\frac{1}{\Delta r} \left( F_{i+\frac{1}{2}}-F\left( V_i^{t,*}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}}\right) - F_{i-\frac{1}{2}} +F\left( V_i^{t,*}(r_{i-\frac{1}{2}}), r_{i-\frac{1}{2}}\right) \right) \\&\quad + \sum _{l =0}^q \alpha _l \left( S({\mathbb {P}}_i^t (r_{i,l} ), r_{i,l} ) -S(V_i^{t,*}(r_{i,l}), r_{i,l}) \right) . \end{aligned} \end{aligned}$$

(2.13)

First-order well-balanced methods are obtained by selecting the trivial constant piecewise reconstruction operator as the standard one, i.e.

$$\begin{aligned} {\mathbb {Q}}_i(r, V_i)=V_i, \quad r \in [r_{i-\frac{1}{2}}, r_{i+\frac{1}{2}}]. \end{aligned}$$

(2.14)

It can be easily checked that the numerical method then reduces to

$$\begin{aligned} \frac{dV_i}{dt} =-\frac{1}{\Delta r} \left( F_{i+\frac{1}{2}}-F\left( V_i^{t,*}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}}\right) - F_{i-\frac{1}{2}} +F\left( V_i^{t,*}(r_{i-\frac{1}{2}}), r_{i-\frac{1}{2}}\right) \right) , \end{aligned}$$

(2.15)

where $F_{i+\frac{1}{2}} = {\mathbb {F}}\left( V^*_i( r_{i+\frac{1}{2}}), V^*_{i+1}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}} \right) .$

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.

3 Burgers–Schwarzschild Model: Designing the Numerical Algorithm

3.1 Preliminaries

For the Burgers–Schwarzschild equation (1.1), the steady state solutions are of the form

$$\begin{aligned} v^{*}(r) = \pm \sqrt{1-K^{2}\Big ( 1-\frac{2M}{r} \Big )}, \ \ K>0. \end{aligned}$$

(3.1)

In Fig. 1 we plot the steady solutions for several values of $K^{2}$. The domain of definition of these stationary solutions is

$$\begin{aligned} D_{K} = {\left\{ \begin{array}{ll} [2M, \infty ), &{} { K^{2}\le 1,} \\ \displaystyle \left[ 2M, \frac{2MK^{2}}{K^{2}-1}\right] , &{} { K^{2}>1}. \end{array}\right. } \end{aligned}$$

(3.2)

It can be easily checked that, given a pair $(K, r^* )$ such that $r^* \in D_K$, the discontinuous function defined in $D_K$ by

$$\begin{aligned} w^*(r) = {\left\{ \begin{array}{ll} \sqrt{1-K^{2}\Big ( 1-\frac{2M}{r} \Big )}, &{} {r \le r^*,}\\ - \sqrt{1-K^{2}\Big ( 1-\frac{2M}{r} \Big )}, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(3.3)

is an entropy weak stationary solution of the Burgers–Schwarzschild model. (For further properties, see the study in [15, 16].)

Fig. 1

Steady solutions to the Burgers model

3.2 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_{i}^{2}$ such that

$$\begin{aligned} \sqrt{1-K_{i}^{2}\left( 1-\frac{2M}{r_{i}}\right) } = |v_{i}|. \end{aligned}$$

(3.4)

There is always a unique solution ${{\widetilde{K}}}_{i}$ given by

$$\begin{aligned} {{\widetilde{K}}}_{i}^{2} = \frac{1-v_{i}^{2}}{1-\displaystyle \frac{2M}{r_{i}}}, \end{aligned}$$

(3.5)

so that the stationary solution is

$$\begin{aligned} v_{i}^{*}(r) =\text {sgn}(v_i)\sqrt{1-{{\widetilde{K}}}_{i}^{2}\Big ( 1-\frac{2M}{r} \Big )}. \end{aligned}$$

(3.6)

In order to apply the numerical method (2.12), this stationary solution has to be computed at $r_{i\pm \frac{1}{2}}$ and this is only possible if $r_{i+\frac{1}{2}} \in D_{{{\widetilde{K}}}_i}$, that is, provided

$$\begin{aligned} {{\widetilde{K}}}_{i}^{2} \le 1 \text { or } \left( {{\widetilde{K}}}_{i}^{2}>1 \text { and } r_{i+\frac{1}{2}} \le \frac{2M{{\widetilde{K}}}_{i}^{2}}{{{\widetilde{K}}}_{i}^{2}-1}\right) . \end{aligned}$$

(3.7)

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.  ${\mathbb {Q}}_{i}(r) = v_{i}.$ The numerical method writes then as follows:

$$\begin{aligned} \frac{dv_{i}}{dt}= -\frac{1}{\Delta r}\left( F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}} - S(v_{i}, r_i)\right) , \end{aligned}$$

(3.8)

where $F_{i+\frac{1}{2}} = {\mathbb {F}}(v_{i}, v_{i+1}, r_{i+\frac{1}{2}}).$

Summing up, the expression of the semi-discrete first-order method reads

$$\begin{aligned} \frac{dv_{i}}{dt}= -\frac{1}{\Delta r}\left( F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}} - S_i \right) , \end{aligned}$$

(3.9)

where

$$\begin{aligned} S_ i = {\left\{ \begin{array}{ll} F(v_{i}^{*}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}}) - F(v_{i}^{*}(r_{i-\frac{1}{2}}), r_{i-\frac{1}{2}})), &{} \text { if }(3.7)\text { holds,}\\ S(v_i, r_i), &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

(3.10)

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.

3.3 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) with ${{\tilde{K}}}_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 -\frac{1}{2}}$, $r_{i+\frac{1}{2}}$, $r_{i+1}$ so that the following condition has to be satisfied $r_{i+1} \in D_{{\widetilde{K}}_i}$, i.e.

$$\begin{aligned} {{\widetilde{K}}}_{i}^{2} \le 1 \text { or } \left( {{\widetilde{K}}}_{i}^{2}>1 \text { and } r_{i+1} \le \frac{2M{{\widetilde{K}}}_{i}^{2}}{{{\widetilde{K}}}_{i}^{2}-1}\right) . \end{aligned}$$

(3.11)

If this condition is satisfied, the following step of the reconstruction procedure consists in computing the fluctuations:

$$\begin{aligned} \begin{aligned}&w_{i-1} = v_{i-1} - \text {sgn}(v_{i-1})\sqrt{1-{\widetilde{K}}_{i}^{2}\left( 1-\frac{2M}{r_{i-1}}\right) }, \\&w_{i} = v_{i} - \text {sgn}(v_{i})\sqrt{1-{\widetilde{K}}_{i}^{2}\left( 1-\frac{2M}{r_{i}}\right) } = 0, \quad \quad \\&w_{i+1} = v_{i+1} - \text {sgn}(v_{i+1})\sqrt{1-{\widetilde{K}}_{i}^{2}\left( 1-\frac{2M}{r_{i+1}} \right) }. \end{aligned} \end{aligned}$$

(3.12)

Then the reconstruction is defined as

$$\begin{aligned} {\mathbb {P}}_{i}(r) = v_{i}^{*}(r) + \text {minmod}\left( \frac{w_{i+1}-w_{i}}{\Delta r}, \frac{w_{i+1}-w_{i-1}}{2\Delta r},\frac{w_{i}-w_{i-1}}{\Delta r}\right) (r-r_{i}), \end{aligned}$$

(3.13)

where

$$\begin{aligned} \text {minmod}(a,b,c) = {\left\{ \begin{array}{ll} \min \{a,b,c\}, &{} \text { } \ \ a,b,c>0, \\ \max \{a,b,c\}, &{} \text { } \ \ a,b,c <0, \\ 0, &{} \,\, \text { otherwise.} \end{array}\right. } \end{aligned}$$

(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:

$$\begin{aligned} \int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} \left( S({\mathbb {P}}_i^t (r), r) -S(V_i^{t,*}(r), r) \right) \, dr \cong \Delta r \left( S({\mathbb {P}}_i^t (r_i), r_i) -S(V_i^{t,*}(r_i), r_i) \right) = 0. \end{aligned}$$

Therefore, the expression (2.13) reduces again to (2.15) with $F_{i+\frac{1}{2}} = {\mathbb {F}}(v_{i+\frac{1}{2}}^{t, -}, v_{i+\frac{1}{2}}^{t, +}, r_{i+\frac{1}{2}}).$

If (3.11) is not satisfied, then the standard MUSCL reconstruction is applied, namely

$$\begin{aligned} {\mathbb {Q}}_{i}(r) = v_i + \text {minmod}\left( \frac{v_{i+1}-v_{i}}{\Delta r}, \frac{v_{i+1}-v_{i-1}}{2\Delta r},\frac{v_{i}-v_{i-1}}{\Delta r}\right) (r-r_{i}). \end{aligned}$$

(3.15)

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

$$\begin{aligned} w_j = v_j - v^*_i(r_j) = v_j - \text {sgn}(v_{i})\sqrt{1-{\widetilde{K}}_{i}^{2}\left( 1-\frac{2M}{r_{j}}\right) }, \quad j=i-1,i,i+1, \end{aligned}$$

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).

3.4 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_{i}^{0} = \frac{1}{2} (v_{0}(r_{i,0})+v_{0}(r_{i,1})),$ where

$$\begin{aligned} r_{i,0} = r_{i-\frac{1}{2}} + \frac{\Delta r}{2}\left( -\sqrt{\frac{1}{3}}+1\right) ,\ \ r_{i,1} = r_{i-\frac{1}{2}} + \frac{\Delta r}{2}\left( \sqrt{\frac{1}{3}}+1\right) . \end{aligned}$$

Therefore, at the first step of the reconstruction procedure one has to look for $K_{i}^{2}$ such that:

$$\begin{aligned} \frac{1}{2}\left( \sqrt{1-K_{i}^{2}\left( 1-\frac{2M}{r_{i,0}}\right) }+ \sqrt{1-K_{i}^{2}\left( 1-\frac{2M}{r_{i,1}}\right) }\right) =|v_{i}|. \end{aligned}$$

(3.16)

If we define the function

$$\begin{aligned} g(x) = \frac{1}{2}\left( \sqrt{1-x\left( 1-\frac{2M}{r_{i,0}}\right) }+ \sqrt{1-x\left( 1-\frac{2M}{r_{i,1}}\right) }\right) , \quad x \ge 0, \end{aligned}$$

it is easily checked that g is a positive decreasing function defined in the interval $[0, K_{i,c}^2]$ with $K_{i,c}^{2} = \big ( 1-\displaystyle \frac{2M}{r_{i,1}}\big )^{-1}.$ Therefore there are two possibilities:

• If $|v_{i}| \in [g(K_{i,c}^{2}), 1]$, there is a unique ${{\widetilde{K}}}^{2}_{i}$ satisfying (3.16).

• In other case, (3.16) has no solution.

If (3.16) is satisfied, the corresponding stationary solution $v_{i}^{*}(r) = \text {sgn}(v_{i})\sqrt{1-{{\widetilde{K}}}_{i}^{2}\Big ( 1-\frac{2M}{r} \Big )}$ 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} \in D_{{{\widetilde{K}}}_i}$. Therefore, the well-balanced reconstruction can be computed only if the following condition is satisfied:

$$\begin{aligned} |v_{i}| \in [g(K_{i,c}^{2}), 1] \text { and } \left( {{\widetilde{K}}}_{i}^{2} \le 1 \text { or } \left( {{\widetilde{K}}}_{i}^{2}>1 \text { and } r_{i+1,1} \le \frac{2M{{\widetilde{K}}}_{i}^{2}}{{{\widetilde{K}}}_{i}^{2}-1}\right) \right) . \end{aligned}$$

(3.17)

If this condition is satisfied, the fluctuations can be then computed:

$$\begin{aligned} w_j = v_j - \text {sgn}(v_{j})\frac{1}{2}\left( \sqrt{1-{{\widetilde{K}}}_{i}^{2} \left( 1-\frac{2M}{r_{j,0}}\right) }+ \sqrt{1-{{\widetilde{K}}}_{i}^{2}\left( 1-\frac{2M}{r_{j,1}}\right) }\right) ,\quad j=i-1,i,i+1, \end{aligned}$$

and the well-balanced reconstruction is given by ${\mathbb {P}}_{i}(r) = v_{i}^{*}(r) + {\mathbb {Q}}_{i}(r; w_{i-1}, w_i, w_{i+1}),$ where ${\mathbb {Q}}$ represents the CWENO approximation. On the other hand, if (3.17) is not satisfied, the standard CWENO reconstruction is applied, namely ${\mathbb {Q}}_{i}(r) = {\mathbb {Q}}_{i}(r; v_{i-1}, v_i, v_{i+1}).$ The expression of the semi-discrete method is finally (3.9) where the numerical fluxes are computed at the reconstructed states and

$$\begin{aligned} S_i = {\left\{ \begin{array}{ll} F(v_{i}^{*}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}}) - F(v_{i}^{*}(r_{i-\frac{1}{2}}), r_{i-\frac{1}{2}}) + \displaystyle \frac{\Delta r}{2}\\ \sum _{j=0,1} \left( S({\mathbb {P}}_{i}(r_{i,j}), r_{i,j})-S(v_{i}^{*}(r_{i,j}), r_{i,j}) \right) &{} \text { if }(3.17)\text { holds,} \\ \displaystyle \frac{\Delta r}{2}\sum _{j=0,1} S({\mathbb {Q}}_{i}(r_{i,j}), r_{i,j}), &{}\text { otherwise.} \end{array}\right. } \end{aligned}$$

(3.18)

The TVDRK3 method of order 3 [13] will be used for the time discretization.

3.5 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_i^{2}$ such that:

$$\begin{aligned} \frac{1}{\Delta r}\int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} \sqrt{1-K_i^{2}\Big ( 1-\frac{2M}{r} \Big )}dr = |v_{i}|. \end{aligned}$$

(3.19)

If we define the function

$$\begin{aligned} g(x) = \frac{1}{\Delta r}\int _{r_{i-\frac{1}{2}}}^{r_{i+\frac{1}{2}}} \sqrt{1-x\Big ( 1-\frac{2M}{r} \Big )}dr, \end{aligned}$$

it can be easily checked that it is a decreasing function defined in $[0, K_{e,i}^{2}]$ where $K_{e,i}^{2} = \left( 1-2M / r_{i+\frac{1}{2}} \right) ^{-1}$ and $g(0) = 1.$ Therefore, (3.19) has a unique solution ${{\widetilde{K}}}_i^2$ if

$$\begin{aligned} |v_i| \le g(K_{e,i}). \end{aligned}$$

(3.20)

The explicit expression of g can be obtained: $g(x) = \frac{1}{\Delta r}\left( f(x, r_{i+\frac{1}{2}})- f(x, r_{i-\frac{1}{2}})\right) ,$ where

$$\begin{aligned} f(x, r) = {\left\{ \begin{array}{ll} r\sqrt{1-x\Big ( 1-\frac{2M}{r} \Big )} +\displaystyle \frac{xM}{\sqrt{1-x}}\log \left( x(M-r) + r + \sqrt{1-x}r\sqrt{1-x\Big ( 1-\frac{2M}{r} \Big )}\right) , &{} {0 \le x < 1,}\\ 2r\sqrt{\displaystyle \frac{2M}{r}}, &{} { x = 1, }\\ r\sqrt{1-x}\Big ( 1-\frac{2M}{r} \Big ) -\displaystyle \frac{2xM}{\sqrt{x-1}}\tan ^{-1}\left( \frac{\sqrt{1-x\Big ( 1-\frac{2M}{r} \Big )}}{\sqrt{x-1}}\right) , &{} { x>1.} \end{array}\right. } \end{aligned}$$

is a primitive function of $\sqrt{1-x\Big ( 1-\frac{2M}{r} \Big )}$. 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 ${\mathcal {S}}_i$ are contained in the domain of definition $D_{{{\widetilde{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.

4 Burgers–Schwarzschild Model: A Numerical Study

4.1 Preliminaries

In this section several numerical tests are applied to check the performance of the well-balanced 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_{-\frac{1}{2}}=0$ as boundary condition since $\Big ( 1-\frac{2M}{r} \Big )= 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:

$$\begin{aligned} F_{i+\frac{1}{2}} = {\mathbb {F}}( v_{i}, v_{i+1},r_{i+\frac{1}{2}}) = \left( 1-\frac{2M}{r_{i+\frac{1}{2}}}\right) \frac{q^{2}(0;v_{i}, v_{i+1})-1}{2}, \end{aligned}$$

where $q(\cdot ; v_L, v_R)$ is the self-similar solution of the Riemann problem for the standard Burgers equation with the initial condition

$$\begin{aligned} v_{0}(r)=\left\{ \begin{array}{l} v_{L}, \quad \ r<{0},\\ v_{R}, \quad \ \ r>{0}. \end{array}\right. \end{aligned}$$

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.

4.2 Stationary Solutions

Positive stationary solution. We consider the initial condition

$$\begin{aligned} v_{0}(r) = \sqrt{\displaystyle \frac{3}{4} + \frac{1}{2r}} \end{aligned}$$

(4.1)

corresponding to the positive stationary solution with $K = \frac{1}{2}$. 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.

Table 1 Well-balanced versus non-well-balanced schemes: $L^{1}$ errors at $t=50$ for the Burgers model with the initial condition (4.1)

Fig. 2

Burgers–Schwarzschild model with the initial condition (4.1): first-, second- and third-order well-balanced and not-well-balanced methods at various times for variable v

Negative stationary solution. Let us consider now as initial condition the negative stationary condition corresponding to $K=\frac{1}{2}$:

$$\begin{aligned} v_{0}(r) = -\sqrt{\displaystyle \frac{3}{4} + \frac{1}{2r}}. \end{aligned}$$

(4.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.

Fig. 3

Burgers–Schwarzschild model with the initial condition (4.2): first-, second- and third-order well-balanced and non-well-balanced methods at selected times for variable v

Table 2 Well-balanced versus non-well-balanced schemes: $L^{1}$ errors at $t=50$ for the Burgers model with the initial condition (4.2)

Discontinuous stationary entropy weak solution. Let us consider finally the discontinuous initial condition

$$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} \sqrt{\displaystyle \frac{3}{4} + \frac{1}{2r}}, &{} { 2<r<3},\\ -\sqrt{\displaystyle \frac{3}{4} + \frac{1}{2r}}, &{} \text { otherwise}, \end{array}\right. } \end{aligned}$$

(4.3)

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-well-balanced 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.

Table 3 Well-balanced versus non-well-balanced schemes: $L^{1}$ errors at $t=50$ for the Burgers model with the initial condition (4.3)

Fig. 4

Burgers–Schwarzschild model with the initial condition (4.3): first-, second- and third-order well-balanced and non-well-balanced methods at selected times for variable v

4.3 Moving Shocks Connecting two Steady Profiles

Right-moving shock. We consider now the initial condition

$$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} \sqrt{\displaystyle \frac{1}{2} + \frac{1}{r}}, &{} { 2<r<2.5},\\ \sqrt{\displaystyle \frac{2}{r}}, &{} \text { otherwise}. \end{array}\right. } \end{aligned}$$

(4.4)

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.

Fig. 5

Burgers–Schwarzschild model with the initial condition (4.4): first-, second- and third-order well-balanced methods at selected times for variable v

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:

$$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} -\sqrt{\displaystyle \frac{2}{r}}, &{} { 2<r<2.5,}\\ -\sqrt{\displaystyle \frac{3}{4} + \frac{1}{2r}}, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(4.5)

see Fig. 6. A reference solution computed with the first-order standard method has been computed again using a mesh of 10000 cells.

Fig. 6

Burgers–Schwarzschild model with the initial condition (4.5): first-, second- and third-order well-balanced methods at selected times for variable v

4.4 Perturbation of a Steady Shock Solution

Left-hand perturbation. In this test case we consider the initial condition:

$$\begin{aligned} {{\tilde{v}}}_{0}(r) = v_0(r) + p_L(r), \end{aligned}$$

(4.6)

where $v_0$ is the steady shock solution given by (4.3) and

$$\begin{aligned} p_L(r) = {\left\{ \begin{array}{ll} \displaystyle -\frac{1}{5}e^{-200(r-2.5)^{2}}, &{} { 2.2<r<2.8,}\\ 0, &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

(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 computational domain, the stationary solution (4.3) is not recovered: a different stationary solution of the family (3.3) is obtained whose shock is placed at a different location. Observe that all the three methods capture the same stationary solution.

Fig. 7

Burgers–Schwarzschild model with the initial condition (4.6)–(4.3)–(4.7): first-, second- and third-order well-balanced methods at selected times for variable v

Right side perturbation. Similar conclusions can be drawn if a perturbation at the right side of the shock is superposed to the stationary solution $v_0$ given by (4.3):

$$\begin{aligned} {{\tilde{v}}}_0(r) = v_0(r) + p_R(r), \end{aligned}$$

(4.8)

with

$$\begin{aligned} p_R(r) = {\left\{ \begin{array}{ll} \displaystyle \frac{1}{5}e^{-200(r-3.5)^{2}}, &{} { 3.2<r<3.8},\\ 0, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(4.9)

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.

Fig. 8

Burgers–Schwarzschild model with the initial condition (4.8)–(4.3)–(4.9): first-, second- and third-order well-balanced methods at selected times for variable v

Left-hand perturbation with zero average. Now we consider an initial condition of the form (4.6) with a perturbation $p_L$ such that $\int p_L(r)dr = 0$, in particular:

$$\begin{aligned} p_L(r) = {\left\{ \begin{array}{ll} \displaystyle 0.1\text {cos}(-25.5 \pi + 10\pi x)e^{-200(x-2.8)^2}, &{} { 2.7<r<2.9,}\\ 0, &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

(4.10)

In Fig. 9 it can be observed that now, after the wave generated by the initial perturbation leaves the computational domain, the stationary solution (4.3) is recovered. Here we have used again a 2000-point uniform mesh to verify that the steady state does not move.

Fig. 9

Burgers–Schwarzschild model with the initial condition (4.6)–(4.3)–(4.10): first-, second- and third-order well-balanced methods at selected times and zoom of the initial and final stationary shocks (right-down) for variable v

Right side perturbation with zero average. Similar conclusions can be drawn if we consider an initial condition of the form (4.8) with $\int p_R(r)dr = 0$. In particular we take

$$\begin{aligned} p_R(r) = {\left\{ \begin{array}{ll} \displaystyle 0.1\text {cos}(-29.5 \pi + 10\pi x)e^{-200(x-3.2)^2}, &{} { 3.1<r<3.3,}\\ 0, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(4.11)

see Fig. 10. Here we have used again a 2000-point uniform mesh to verify that the steady state does not move.

Fig. 10

Burgers–Schwarzschild model with the initial condition (4.8)–(4.3)–(4.11): first-, second- and third-order well-balanced methods at selected times and zoom of the initial and final stationary shocks (right-down) for variable v

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:

$$\begin{aligned} {{\tilde{v}}}_{0}(r) = v_0(r) + p_L(r) + p_R(r), \end{aligned}$$

(4.12)

where $v_0$ is the steady shock solution given by (4.3) and $\int (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 before, the right perturbation creates a lower displacement than the left-hand perturbation. Here we have used again a 2000-point uniform mesh.

Fig. 11

Burgers–Schwarzschild model with the initial condition (4.8)–(4.3)–(4.11): comparison between the first-, second- and third-order well-balanced methods at selected times for variable v

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:

$$\begin{aligned} {{\tilde{v}}}_0 (r) = v_0(r) + \delta _v(\alpha , r), \end{aligned}$$

(4.13)

where $v_0$ is given again by (4.3) and

$$\begin{aligned} \delta _v(\alpha , r) = {\left\{ \begin{array}{ll} \alpha \cos (5\pi r-12\pi )e^{-200(r-2.8)^{2}}, &{} { 2.7<r<2.9,}\\ 0, &{} \text { otherwise}, \end{array}\right. } \end{aligned}$$

(4.14)

with $\alpha >0$. The amplitude of the perturbation is measured by $\int \delta _v(\alpha , r) \,dr$ and the distance between the shocks is measured by $\lim _{t \rightarrow \infty } \int |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.

Fig. 12

Burgers–Schwarzschild model with the initial condition (4.13)–(4.3)–(4.14): measures of the perturbation and the shock displacement for $\alpha = 1$

4.5 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) considering different initial conditions, different meshes, and different lengths of the computational domain: the observed behavior of the numerical solutions have been always one of the four ones shown here depending on the value at 2M (1 or lower) and at the right boundary (positive or negative).

• 1. Initial condition satisfying $v_{0}(2M)=1$ and $v_{0}(L)\ge 0$: let us consider the initial condition

  $$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} 1, &{} { 2<r<2.1},\\ \cos (30r)e^{\frac{-1}{(x-2.5)^{2}}}, &{} \text { otherwise}, \end{array}\right. } \end{aligned}$$

  (4.15)

  that takes value 1 in a neighborhood of $2M =2$ and a positive value at the right boundary of the computational domain $x = 4$. As it can be observed in Fig. 14 after a transient regime, the numerical solution takes the form of a right-moving shock linking the stationary solution $v \equiv 1$ with the negative stationary solution that takes value $-1$ at $x= 2M$ and value 0 at $x = 4$. Once this shock leaves the domain, the stationary solution $v\equiv 1$ is reached in the whole computational domain.

• 2. Initial condition satisfying $v_{0}(2M)=1$ and $v_{0}(L)< 0$: we consider now the initial condition

  $$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} 1, &{} { 2<r<2.1,}\\ \cos (20r)e^{\frac{-1}{(x-2.5)^{2}}}, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

  (4.16)

  that takes value 1 in a neigborhood of $2M =2$ and negative value at the right boundary of the computational domain $x = 4$. As it can be observed in Fig. 15 after a transient period, the numerical solution takes the form of a right-moving shock linking the stationary solution $v \equiv 1$ with the negative stationary solution that takes value $-1$ at $x= 2M$ and value $v_0(4)$ at $x = 4$. Once this shock leaves the domain, the stationary solution $v\equiv 1$ is reached in the whole computational domain.

• 3. Initial condition satisfying $v_{0}(2M)<1$ and $v_{0}(L)\ge 0$: we consider now the initial condition

  $$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} 0.8, &{} { 2<r<2.1,}\\ \cos (30r)e^{\frac{-1}{(x-2.5)^{2}}}, &{} \text { otherwise}. \end{array}\right. } \end{aligned}$$

  (4.17)

  In this case the numerical solution reaches in finite time the negative stationary solution $v^*$ such that $v^*(2) = -1$ and $v^*(4) = 0$: see Fig. 16.

• 4. Initial condition satisfying $v_{0}(2M)<1$ and $v_{0}(L)< 0$: we finally consider the initial condition

  $$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} 0.8, &{} { 2<r<2.1,}\\ \cos (20r)e^{\frac{-1}{(x-2.5)^{2}}}, &{} \text { otherwise}. \end{array}\right. } \end{aligned}$$

  (4.18)

  The numerical solution reaches in finite time the negative stationary solution $v^*$ such that $v^*(2) = -1$ and $v^*(4) = v_0(4)$: see Fig. 17.

Table 4 Burgers–Schwarzschild model with the initial condition (4.13)–(4.3)–(4.14): measures of the perturbation and the shock displacement for different values of $\alpha$

Fig. 13

Burgers–Schwarzschild model with the initial condition (4.13)–(4.3)–(4.14): values of $\lim _{t\rightarrow \infty } \int |v -v_{s}|$ as a function of $\int \delta _{v}$

Fig. 14

Burgers–Schwarzschild model with the initial condition (4.15): first-order well-balanced scheme at selected times for variable v

Fig. 15

Burgers–Schwarzschild model with the initial condition (4.16): first-order well-balanced scheme at selected times for variable v

Fig. 16

Burgers–Schwarzschild model with the initial condition (4.17): first-order well-balanced scheme at selected times for variable v

Fig. 17

Burgers–Schwarzschild model with the initial condition (4.18): numerical solution obtained with the first-order well-balanced scheme at selected times for variable v

4.6 Main Conclusions for the Burgers–Schwarzschild Model

From Figs. 7, 8, 9, 10, 11, 12 and 13 and Table 4 we can conclude the following:

Conclusion 1

If a perturbation $\delta _{v}$ is added to a steady shock solution of the form

$$\begin{aligned} v_{0}(r) = {\left\{ \begin{array}{ll} \sqrt{1-K_{0}^{2}(1-\frac{2M}{r})}, &{} { 2M<r<r_{0},}\\ -\sqrt{1-K_{0}^{2}(1-\frac{2M}{r})}, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

then the solution reaches at finite time another steady shock solution of the form:

$$\begin{aligned} v(r) = {\left\{ \begin{array}{ll} \sqrt{1-K_{0}^{2}(1-\frac{2M}{r})}, &{} { 2M<r<r_{1},}\\ -\sqrt{1-K_{0}^{2}(1-\frac{2M}{r})}, &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

• 1. If $\int _{2M}^{r_{0}} \delta _{v} = 0$ and $\int _{r_{0}}^{\infty } \delta _{v} = 0$, then $r_{1}=r_{0}$, i.e. the initial stationary solution is recovered.

• 2. If $\int _{2M}^{\infty } \delta _{v} = 0$ and $\int _{2M}^{r_{0}} \delta _{v} = -\int _{r_{0}}^{\infty } \delta _{v}$, then $r_{1}\ne r_{0}$ and a different stationary solution is obtained.

• 3. If $\int \delta _{v} \ne 0$, then $r_{1}\ne r_{0}$ and a different stationary solution is obtained. In this case the distance between $r_0$ and $r_1$ depends linearly on the amplitude of the perturbation: see Table 4 and Fig. 13.

In view of Figs. 14, 15, 16 and 17 we have reached the following.

Conclusion 2

• 1. For a bounded domain [2M, L]:

  • (a) If $v_{0}(r)= 1$ for $r\in [2M, 2M+\epsilon )$, with $\epsilon >0$, $v_{0}(L)\ge 0$ and $v_{0}\ne 1$, in finite time the solution has the form of a right-moving shock that links the stationary solution $v\equiv 1$ and the negative steady solution $v^*$ such that $v^*(2M) = -1$ and $v^{*}(L) = 0$, that is, $v^{*}_{0}(r)= -\sqrt{1-\frac{1}{1-\frac{2M}{L}}\Big ( 1-\frac{2M}{r} \Big )}.$

  • (b) If $v_{0}(r)= 1$ for $r\in [2M, 2M+\epsilon )$, with $\epsilon >0$ and $v_{0}(L)=a$, with $a<0$, then in finite time the solution has the form of a right-moving shock that links the stationary solution $v\equiv 1$ and the negative steady solution $v*$ such that $v^*(2M) = -1$ and $v^{*}_{0}(L) = a$, that is, $v^{*}_{0}(r)= -\sqrt{1-\frac{1-a^{2}}{1-\frac{2M}{L}}\Big ( 1-\frac{2M}{r} \Big )}.$

  • (c) If $v_{0}(2M)< 1$ and $v_{0}(L)\ge 0$, then in finite time the solution coincides with the negative steady solution such that $v^*(2M) = -1$ and $v^{*}_{0}(L) = 0$, that is, $v^{*}_{0}(r)= -\sqrt{1-\frac{1}{1-\frac{2M}{L}}\Big ( 1-\frac{2M}{r} \Big )}.$

  • (d) If $v_{0}(2M)< 1$ and $v_{0}(L)=a$, with $a<0$, then in finite time the solution coincides with the negative stationary solution $v^*$ such that $v^*(L) = a$, that is, $v^{*}_{0}(r)= -\sqrt{1-\frac{1-a^{2}}{1-\frac{2M}{L}}\Big ( 1-\frac{2M}{r} \Big )}.$

• 2. For the unbounded domain $[2M,\infty )$ the following conclusions can be drawn by passing to the limit when $L \rightarrow \infty$:

  • (a) If $v_{0}(r)= 1$ for $r\in [2M, 2M+\epsilon )$, with $\epsilon >0$, $\lim _{r\rightarrow \infty }v_{0}(r)\ge 0$ and $v_{0}\ne 1$, in finite time the solution has the form of a right-moving shock that links the stationary solution $v \equiv 1$ and the negative stationary solution $v^{*}_{0}(r)= -\sqrt{\frac{2M}{r}},$ corresponding to $K^2 = 1$.

  • (b) If $v_{0}(r)= 1$ for $r\in [2M, 2M+\epsilon )$, with $\epsilon >0$ and $\lim _{r\rightarrow \infty } v_{0}(r)=a$, with $a<0$, then in finite time $t_{0}$ the solution has the form of a right- moving shock that links the stationary solution $v \equiv 1$ and the negative stationary solution $v^*$ such that $v^*(2M) = -1$ and $\lim _{r\rightarrow \infty } v^{*}_{0}(r) = a$, that is, $v^{*}_{0}(r)= -\sqrt{1-(1-a^{2})\Big ( 1-\frac{2M}{r} \Big )}.$

  • (c) If $v_{0}(2M)< 1$ and $\lim _{r\rightarrow \infty } v_{0}(r)\ge 0$, then the solution converges as $t \rightarrow \infty$ to the negative stationary solution $v^*$ such that $v^*(2M) = -1$ and $\lim _{r\rightarrow \infty } v^{*}_{0}(r) = 0$, that is, $v^{*}_{0}(r)= -\sqrt{\frac{2M}{r}}.$

  • (d) If $v_{0}(2M)< 1$ and $\lim _{r\rightarrow \infty } v_{0}(r)=a$, with $a<0$, then the solution converges as $t \rightarrow \infty$ to the negative stationary solution $v^*$ such that $v^*(2M) = -1$ and $\lim _{r\rightarrow \infty } v^{*}_{0}(r) = a$, that is, $v^{*}_{0}(r)= -\sqrt{1-(1-a^{2})\Big ( 1-\frac{2M}{r} \Big )}.$

5 Euler–Schwarzschild Model: Designing the Numerical Algorithm

5.1 Preliminaries

In the case of the Euler–Schwarzschild equations (1.3), the stationary solutions are implicitly given by the equations:

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{\text {sgn}(v)(1-v^{2})|v|^{\frac{2k^{2}}{1-k^{2}}}r^{\frac{4k^{2}}{1-k^{2}}}}{\Big ( 1-\frac{2M}{r} \Big )} = C_{1}, \quad \displaystyle r(r-2M)\rho \frac{v}{1-v^{2}}=C_{2}, \end{array} \end{aligned}$$

(5.1)

where $C_1, C_2$ are constants. The pair $(v, \rho )$ of a stationary solution satisfies the following ODE system analyzed first in [15, 16]:

$$\begin{aligned} \frac{dv}{dr}= & {} v\frac{(1-v^{2})(1-k^{2})}{r(r-2M)}\left( \frac{2k^{2}}{1-k^{2}}(r-2M)-M\right) \Big / (v^{2}-k^{2}), \end{aligned}$$

(5.2)

$$\begin{aligned} \frac{d\rho }{dr}= & {} -\frac{2(r-M)}{r(r-2M)}\rho - \rho \frac{(1+v^{2})(1-k^{2})}{r(r-2M)}\left( \frac{2k^{2}}{1-k^{2}}(r-2M)-M\right) \Big / (v^{2}-k^{2}). \end{aligned}$$

(5.3)

Figure 18 shows the graph of v for some of them. When these functions are defined in $(2M,\infty )$, they have a maximum or a minimum in

$$\begin{aligned} r_{c} = \frac{M(1-k^{2})}{2k^{2}} +2M, \end{aligned}$$

(5.4)

that comes from solving $\frac{dv}{dr} =0$. In Fig. 18 the stationary solutions marked in red are those that take the value $v=\pm k$ at $r=r_{c}$.

Fig. 18

Euler–Schwarzschild model with $k = 0.3$: graph of the variable v for some stationary solutions

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:

$$\begin{aligned} \text {sgn}(v)(1-v^{2})|v|^{\frac{2k^{2}}{1-k^{2}}} = \displaystyle \frac{\left( 1-\frac{2M}{a}\right) }{a^{\frac{4k^{2}}{1-k^{2}}}}C_1, \quad \displaystyle \rho =\frac{1-v^{2}}{va(a-2M)}C_2. \end{aligned}$$

(5.5)

It is enough thus to solve, if it is possible, the nonlinear equation

$$\begin{aligned} g(v) = K_a, \end{aligned}$$

(5.6)

with

$$\begin{aligned} g(v) =\text {sgn}(v)(1-v^{2})|v|^{\frac{2k^{2}}{1-k^{2}}}, \quad v \in [-1,1], \quad K_a = \frac{\left( 1-\frac{2M}{a}\right) }{a^{\frac{4k^{2}}{1-k^{2}}}}C_1, \end{aligned}$$

(5.7)

to compute v. Once this equation is solved, $\rho$ is computed using the second equation of (5.5).

It can be easily checked that g satisfies:

$$\begin{aligned} -(1-k^{2})k^{\frac{2k^{2}}{1-k^{2}}} = g(-k) \le g(v) \le g(k) = (1-k^{2})k^{\frac{2k^{2}}{1-k^{2}}}, \quad v \in [-1,1]. \end{aligned}$$

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 $(\rho , v)$ that solve (5.5), one supersonic and one subsonic.

For the sake of clarity, together with the representation

$$\begin{aligned} V = [\rho (1+k^{2}v^{2})/(1-v^{2}), \rho v (1+k^{2})/(1-v^{2})]^T, \end{aligned}$$

for the states, we will use ${{\widetilde{V}}}= [\rho , v]^T.$ Here, V can be easily computed from ${{\widetilde{V}}}$ and, given V, ${{\widetilde{V}}}$ is also easily computed by (1.3d) that comes from solving a second-degree equation.

5.2 First-order Method

If the midpoint rule is used to compute the initial averages, given a family of cell values ${{\widetilde{V}}}_i$, in the first step of the well-balanced reconstruction procedure one has to find, if it is possible, a stationary solution ${{\widetilde{V}}}^*_i$ defined in the interval $[r_{i-\frac{1}{2}}$, $r_{i+\frac{1}{2}}]$ such that ${\widetilde{V}}_i^*(r_i) = {\widetilde{V}}_i = [\rho _i, v_i]^T.$ Obviously such a stationary solution would correspond to the choice of constants:

$$\begin{aligned} C_{i,1}= & {} \displaystyle \frac{\text {sgn}(v_{i})(1-v_{i}^{2})|v_{i}|^{\frac{2k^{2}}{1-k^{2}}}r_{i}^{\frac{4k^{2}}{1-k^{2}}}}{\left( 1-\frac{2M}{r_{i}}\right) }, \quad C_{i,2} = \displaystyle r_{i}(r_{i}-2M)\rho _{i}\frac{v_{i}}{1-v_{i}^{2}}. \end{aligned}$$

(5.8)

According to the discussion above, the corresponding stationary solution is defined in $r_{i \pm \frac{1}{2}}$ provided that:

$$\begin{aligned} |K_{i \pm \frac{1}{2}}| \le g(k), \quad K_{i \pm \frac{1}{2}} = \left( 1-\frac{2M}{r_{i \pm \frac{1}{2}}}\right) r_{i \pm \frac{1}{2}}^{-\frac{4k^{2}}{1-k^{2}}}C_{i,1}. \end{aligned}$$

(5.9)

When $|K_{i \pm \frac{1}{2}}| < g(k)$ there are two possible values for ${{\widetilde{V}}}^*_i(r_{i \pm \frac{1}{2}})$, one subsonic and one supersonic. Therefore, a criterion is needed to select one or the other. The following criterion will be used here:

• if ${{\widetilde{V}}}_i$ is not sonic, then the state whose regime (sub or supersonic) is the same as ${{\widetilde{V}}}_i$ is selected for ${\widetilde{V}}^*_i(r_{i \pm \frac{1}{2}}).$

• if ${\widetilde{V}}_i$ is sonic, then the state whose regime is the same as ${\widetilde{V}}_{i+1}$ is selected for ${\widetilde{V}}^*_i(r_{i + \frac{1}{2}})$ and the state whose regime is the same as ${\widetilde{V}}_{i-1}$ is selected for ${{\widetilde{V}}}^*_i(r_{i - \frac{1}{2}})$.

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:

$$\begin{aligned} \frac{dV_{i}}{dt}= -\frac{1}{\Delta r}\left( F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}} - S_i \right) , \end{aligned}$$

(5.10)

where

$$\begin{aligned} S_ i = {\left\{ \begin{array}{ll} F(V_{i}^{*}(r_{i+\frac{1}{2}}), r_{i+\frac{1}{2}}) - F(V_{i}^{*}(r_{i-\frac{1}{2}}), r_{i-\frac{1}{2}}), &{} \text { if }(5.9)\text { holds,} \\ S(V_i, r_i), &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

(5.11)

The forward Euler method is used again for the time discretization.

5.3 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 -\frac{1}{2}}$, $r_{i+\frac{1}{2}}$, $r_{i+1}$ so that the following condition has to be satisfied:

$$\begin{aligned} |K_{i+j}| \le g(k),\quad j = -1, -\frac{1}{2}, \frac{1}{2}, 1, \end{aligned}$$

(5.12)

where $K_{i\pm \frac{1}{2}}$ are given by (5.9) and

$$\begin{aligned} K_{i \pm 1} = \left( 1-\frac{2M}{r_{i \pm 1}}\right) r_{i \pm 1}^{-\frac{4k^{2}}{1-k^{2}}}C_{i,1}. \end{aligned}$$

(5.13)

If this condition is satisfied, the reconstruction is defined as follows:

$$\begin{aligned} {\mathbb {P}}_{i}(r) = V_{i}^{*}(r) + \text {minmod}\left( \frac{W_{i+1}-W_{i}}{\Delta r}, \frac{W_{i+1}-W_{i-1}}{2\Delta r},\frac{W_{i}-W_{i-1}}{\Delta r}\right) (r-r_{i}), \end{aligned}$$

where the minmod function is applied component by component and $W_j = V_j - V_i^*(r_j)$ for $j = i-1, i, i+1$. 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:

$$\begin{aligned} {\mathbb {Q}}_{i}(r) = V_i + \text {minmod}\left( \frac{V_{i+1}-V_{i}}{\Delta r}, \frac{V_{i+1}-V_{i-1}}{2\Delta r},\frac{V_{i}-V_{i-1}}{\Delta r}\right) (r-r_{i}). \end{aligned}$$

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.

6 Euler–Schwarzschild Model: A Numerical Study

6.1 Preliminaries

In this section several tests are considered to check the performance of the first- and second-order well-balanced numerical methods for Euler–Schwarzschild model introduced in the previous section. We consider the spatial interval [2M, L] with $M=1$ and $L=10$, a 500-point uniform mesh, $k=0.3$ and the CFL number is set to 0.5 again. At $r=2M$ we impose $F_{-\frac{1}{2}}=0$ as boundary condition since $\Big ( 1-\frac{2M}{r} \Big )= 0$. At $r=L$ we will use a transmissive boundary condition in the case we are not in a stationary solution or we will expand the stationary solution if we are in one.

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

$$\begin{aligned} F_{i+\frac{1}{2}} = \frac{1}{2}(F(V_{i})+F(V_{i+1}))-\frac{1}{2}\frac{\Delta t}{\Delta x}(V_{i+1}-V_{i}), \end{aligned}$$

(6.1)

and a HLL-like numerical flux in PVM form (see [4]):

$$\begin{aligned} F_{i+\frac{1}{2}} = \frac{1}{2}(F(V_{i})+F(V_{i+1}))-\frac{1}{2}\left( \alpha _{0}(V_{i+1}-V_{i}) +\alpha _{1}(F(V_{i+1})-F(V_{i}))\right) , \end{aligned}$$

(6.2)

with

$$\begin{aligned} \alpha _{0}=\frac{\overline{\lambda _{2}}|\overline{\lambda _{1}}| - \overline{\lambda _{1}}|\overline{\lambda _{2}}|}{\overline{\lambda _{2}}- \overline{\lambda _{1}}}, \quad \alpha _{1}=\frac{|\overline{\lambda _{2}}| - |\overline{\lambda _{1}}|}{\overline{\lambda _{2}}- \overline{\lambda _{1}}}, \end{aligned}$$

(6.3)

where $\overline{\lambda _{1}}$ and $\overline{\lambda _{2}}$ are the eigenvalues of some intermediate matrix $J_{i+\frac{1}{2}}$ of the form

$$\begin{aligned} J_{i+\frac{1}{2}} = \left( 1- \frac{2M}{r_{i+\frac{1}{2}}}\right) \left[ \begin{array}{cc} 0 &{} 1\\ \displaystyle \frac{k^2 - v_m^2}{1 - k^2 v_m^2} &{} \displaystyle \frac{2(1 - k^2)v_m}{1 - k^2 v_m^2} \end{array}\right] , \end{aligned}$$

(6.4)

where $v_m$ is some intermediate value between $v_i^n$ and $v_{i+1}^n$.

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:

$$\begin{aligned} \begin{aligned}&\left[ \begin{array}{cc} 0 &{} 1\\ \displaystyle \frac{k^2 - v^2}{1 - k^2 v^2} &{} \displaystyle \frac{2(1 - k^2)v}{1 - k^2 v^2} \end{array}\right] \cdot (V_R - V_L ) = {{\widehat{F}}}_R - {{\widehat{F}}}_L, \quad \\&\quad {{\widehat{F}}}_\alpha = \left( \begin{array}{c} \displaystyle \frac{1+k^{2}}{1-v_\alpha ^{2}}\rho _\alpha v_\alpha \\ \displaystyle \frac{v_\alpha ^{2}+k^{2}}{1-v_\alpha ^{2}}\rho _\alpha \end{array}\right) ,\quad \alpha = L, R, \end{aligned} \end{aligned}$$

(6.5)

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

$$\begin{aligned} \frac{k^2 - v^2}{1 - k^2 v^2}( V_{1, R} - V_{1,L}) + \frac{2(1 - k^2)v}{1 - k^2 v^2} ( V_{2,R} - V_{2,L}) = F_{2,R} - F_{2,L}. \end{aligned}$$

This equality is equivalent to a second-order degree equation for v, namely $\alpha v^2 + \beta v + \gamma = 0,$ where

$$\begin{aligned}&\alpha = \rho _R(1- v_L^2) - \rho _L(1- v_R^2), \\&\beta =- 2 \left( \rho _Rv_R(1- v_L^2) - \rho _Lv_L(1- v_R^2)\right) , \\&\gamma = \rho _Rv_R^2(1- v_L^2) - \rho _Lv_L^2(1- v_R^2). \end{aligned}$$

Since the discriminant $D = \rho _L \rho _R (1- v_L^2)(1- v_R^2) (v_R - v_L)^2$ is always positive, there are always two real solutions:

$$\begin{aligned} v_{\pm } = \frac{\rho _R v_R(1-v_L^2) - \rho _L v_L (1-v_R^2) \pm |v_R-v_L| \sqrt{\rho _L \rho _R (1- v_L^2)(1- v_R^2)}}{\rho _R(1-v_L^2)-\rho _L(1-v_R^2)} \end{aligned}$$

and it can be proven that:

• if $v_L<v_R$, then $v_{-} \in (v_L,v_R)$ and $v_{+} \notin (v_L,v_R)$, so that we will take $v_m=v_{-}$;

• if $v_L>v_R$ then $v_{-} \notin (v_R,v_L)$ and $v_{+} \notin (v_R,v_L)$, so that we will take $v_m=v_{+}$.

Finally, in the case $\alpha =0$ and $V_R\ne V_L$, we take $v_m = -\frac{\gamma }{\beta }$ and in the case $||V_R- V_L||_\infty < \epsilon$ we take $v_m=\frac{v_L+v_R}{2}$.

Once $v_m$ has been chosen, the expression of $\overline{\lambda _j}$, $j=1,2$ in (6.3) is as follows:

$$\begin{aligned}&\overline{\lambda _{1}} = \lambda _{1}\left( v_m\right) = \left( 1-\frac{2M}{r_{i+\frac{1}{2}}}\right) \frac{v_m-k}{1-k^{2}v_m}, \qquad \quad \overline{\lambda _{2}} = \lambda _{2}\left( v_m\right) = \left( 1-\frac{2M}{r_{i+\frac{1}{2}}}\right) \frac{v_m+k}{1+k^{2}v_m}. \end{aligned}$$

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 second-order reconstructions.

6.2 Stationary Solutions

Positive stationary solution. We take as initial condition the positive supersonic stationary solution satisfying

$$\begin{aligned} \rho ^{*}(10)=1, \quad v^{*}(10)=0.6. \end{aligned}$$

(6.6)

Table 5 shows the error in $L^1$ norm between the numerical solution at time $t = 50$ for the well-balanced and non-well-balanced methods using the Roe-type numerical flux. Figs. 19 and 20 compare the numerical solutions obtained with the well-balanced and the non-well-balanced methods: as it happened for the Burgers–Schwarzschild model, the numerical solutions obtained with non-well-balanced methods depart from the initial steady state.

Table 5 Well-balanced versus non-well-balanced schemes: $L^{1}$ errors at time $t=50$ for the Euler–Schwarzschild model with the initial condition (6.6)

Fig. 19

Euler–Schwarzschild model with the initial condition (6.6): first- and second-order well-balanced and non-well-balanced methods at selected times for the variable v

Fig. 20

Euler–Schwarzschild model with the initial condition (6.6): first- and second-order well-balanced and non-well-balanced methods at selected times for the variable $\rho$

Negative stationary solution. Let us consider now as initial condition the negative supersonic stationary solution $V^*$ that satisfies

$$\begin{aligned} \rho ^{*}(10)=1, \quad v^{*}(10)=-0.8. \end{aligned}$$

(6.7)

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-well-balanced methods using the Roe-type numerical flux. Again the numerical solutions obtained with non-well-balanced methods depart from the initial steady state.

Table 6 Well-balanced versus non-well-balanced schemes: $L^{1}$ errors at time $t=50$ for the Burgers–Schwarzschild model with the initial condition (6.7)

Fig. 21

Euler–Schwarzschild model with the initial condition the stationary solution satisfying (6.7): first- and second-order well-balanced and non-well-balanced methods at selected times for the variable v

Fig. 22

Euler–Schwarzschild model with the initial condition (6.7): first- and second-order well-balanced and non-well-balanced methods at selected times for the variable $\rho$

Discontinuous stationary entropy weak solution. We consider finally the initial condition

$$\begin{aligned} V_{0}(r) = {\left\{ \begin{array}{ll} V_{-}^{*}(r), &{} { r\le 6},\\ V_{+}^{*}(r), &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(6.8)

where $V_{-}^{*}(r)$ is the supersonic stationary solution such that

$$\begin{aligned} \rho _{-}^{*}(6)=4, \quad v_{-}^{*}(6)=0.6 \end{aligned}$$

(6.9)

and $V_{+}^{*}(r)$ is the subsonic one such that

$$\begin{aligned} \rho _{+}^{*}(6)=\frac{\rho _{-}^{*}(6)(v_{-}^{*}(6)^{2}-k^{4})}{k^{2} (1-v_{-}^{*}(6)^{2})}, \quad v_{+}^{*}(6)=\frac{k^{2}}{v_{-}^{*}(6)}. \end{aligned}$$

(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.

Table 7 Well-balanced versus non-well-balanced schemes: $L^{1}$ errors at time $t=50$ for the Burgers–Schwarzschild model with the initial condition (6.8)

Fig. 23

Euler–Schwarzschild model with the initial condition (6.8): first- and second-order well-balanced and non-well-balanced methods at selected times for the variable v

Fig. 24

Euler–Schwarzschild model with the initial condition (6.8): first- and second-order well-balanced and non-well-balanced methods at selected times for the variable $\rho$

6.3 Perturbation of a Regular Stationary Solution

In this test we consider the initial condition

$$\begin{aligned} {{\widetilde{V}}}_0(r) ={{\widetilde{V}}}^*(r) + \delta (r), \end{aligned}$$

(6.11)

where ${{\widetilde{V}}}^*$ is the supersonic stationary solution

$$\begin{aligned} \rho ^{*}(10)=1, \quad v^{*}(10)=0.9 \end{aligned}$$

(6.12)

and

$$\begin{aligned} \delta (r) = [\delta _{v}(r), \delta _{\rho }(r)]^T = {\left\{ \begin{array}{ll} [-0.01e^{-200(r-6)^{2}}, 0]^T, &{} { 5<r<7,}\\ {[}0, 0]^T, &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

(6.13)

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 second-order 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.

Fig. 25

Euler–Schwarzschild model with the initial condition (6.11): first- and second-order well-balanced at selected times for the variable v

6.4 Perturbation of a Steady Shock Solution

Left-hand perturbation. We consider the initial condition

$$\begin{aligned} {{\widetilde{V}}}_{0}(r) = {\widetilde{V}}^*(r) + \delta _L(r), \end{aligned}$$

(6.14)

where

$$\begin{aligned} \delta _L(r) = [\delta _{v,L}(r), \delta _{\rho ,L}(r)]^T = {\left\{ \begin{array}{ll} [0.2e^{-200(r-4)^{2}}, 0]^T, &{} { 3<r<5,}\\ {[}0, 0]^T, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(6.15)

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 $\displaystyle \frac{v_i-v_{i-1}}{v_{i+1}-v_{i}} \ge 0.8$.

Fig. 26

Euler–Schwarzschild model with the initial condition (6.14): comparison between the first-order well-balanced method with different meshes using the Roe-type and the Lax numerical fluxes at selected times for the variable v

Fig. 27

Euler–Schwarzschild model with the initial condition (6.14): first-order well-balanced method with different meshes using the Roe-type and the Lax numerical fluxes at selected times for the variable $\rho$

Fig. 28

Euler–Schwarzschild model taking as initial condition (6.14): evolution of the shock position with time obtained with the first-order well-balanced method using the Roe-type numerical flux with different meshes

Right-hand perturbation. Let us consider now two different initial conditions: on the one hand

$$\begin{aligned} {{\widetilde{V}}}_{0,1}(r) = {\widetilde{V}}^*(r) + \delta _R(r), \end{aligned}$$

(6.16)

where ${\widetilde{V}}^*(r)$ is again the discontinuous stationary solution given by (6.8)–(6.10) and

$$\begin{aligned} \delta _R(r) = [\delta _{v,R}(r), \delta _{\rho ,R}(r)]^T = {\left\{ \begin{array}{ll} [-0.05e^{-200(r-8)^{2}}, 0]^T, &{} { 7<r<9,}\\ {[}0, 0]^T, &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

(6.17)

On the other hand,

$$\begin{aligned} {{\widetilde{V}}}_{0,2}(r) = {\widetilde{V}}_2^*(r) + \delta _R(r), \end{aligned}$$

(6.18)

where $\delta _R$ is given again by (6.16) and ${\widetilde{V}}_2^*(r)$ is the steady shock of the form (6.8) satisfying

$$\begin{aligned} \rho _{-}^{*}(6)=5, \quad v_{-}^{*}(6)=0.6. \end{aligned}$$

(6.19)

Observe that the definition of v is identical for both stationary solutions but $\rho$ 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:

$$\begin{aligned} V^*(r) = {\left\{ \begin{array}{ll} V_{-}^{*}(r), &{} { r\le r_c},\\ V_{+}^{*}(r), &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(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_{\pm }^{*}(r_c)=-k$: see Fig. 29. Nevertheless, the limits in time of the approximations of $\rho$ are different: see Fig. 30. Observe that, in the Euler–Schwarzschild model (5.1), there are infinitely many stationary solutions with the same function v and different $\rho$.

Fig. 29

Euler–Schwarzschild model with the initial conditions (6.16) and (6.18): first-order well-balanced method with a 2000-point mesh using the Roe-type numerical flux at selected times for the variable v: the numerical solutions coincide

Fig. 30

Euler–Schwarzschild model with the initial conditions (6.16) and (6.18): first-order well-balanced method with a 2000-point mesh using the Roe-type numerical flux at selected times for the variable $\rho$

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:

$$\begin{aligned} {{\widetilde{V}}}_{0}(r) = {\widetilde{V}}^*(r) + \delta (\alpha , r), \end{aligned}$$

(6.21)

where ${{\widetilde{V}}}^*$ is the steady shock solution given by (6.8)–(6.10) and

$$\begin{aligned} \delta (\alpha , r) = [\delta _{v}(\alpha , r), \delta _{\rho }(\alpha , r)]^T = {\left\{ \begin{array}{ll} [\alpha e^{-200(r-4)^{2}}, 0]^T, &{} { 3<r<5,}\\ {[}0, 0]^T, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(6.22)

with $\alpha >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 $\alpha$ 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 $\int \delta _v(\alpha , r) \,dr$ and the distance between the shocks is measured by $\lim _{t \rightarrow \infty } \int |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.

Fig. 31

Euler–Schwarzschild model with the initial condition (6.21): first-order well-balanced method taking different values of $\alpha$ for variable v

Fig. 32

Euler–Schwarzschild model with the initial condition (6.21): first-order well-balanced method taking different values of $\alpha$ for variable $\rho$

Table 8 Euler–Schwarzschild model with the initial condition (6.21): measures of the perturbation and the shock displacement for different values of $\alpha$

Fig. 33

Euler–Schwarzschild model with the initial condition (6.21): values of $\lim _{t\rightarrow \infty } \int |v -v^{*}|$ as a function of $\int \delta _{v}$

Let us finally consider a family of initial conditions that generate leftward displacement of the initial steady shock:

$$\begin{aligned} {{\widetilde{V}}}_{0}(r) = {\widetilde{V}}^*(r) + \delta (\beta , r), \end{aligned}$$

(6.23)

where ${{\widetilde{V}}}^*$ is again the steady shock solution given by (6.8)–(6.10) and

$$\begin{aligned} \delta (\beta , r) = [\beta {v}(\beta , r), \delta _{\rho }(\beta , r)]^T = {\left\{ \begin{array}{ll} [\beta e^{-200(r-8)^{2}}, 0]^T, &{} { 7<r<8,}\\ {[}0, 0]^T, &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

(6.24)

with $\beta <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 $\beta$ and we observe that it converges to the same stationary solution regardless of the perturbation.

Fig. 34

Euler–Schwarzschild model with the initial condition (6.23): first-order well-balanced method taking different values of $\beta$ for variable v

Fig. 35

Euler–Schwarzschild model with the initial condition (6.23): first-order well-balanced method taking different values of $\beta$ for variable $\rho$

6.5 Main Conclusions for the Euler–Schwarzschild Model

In view of Fig. 25 we arrive at the following observation.

Conclusion 3

If a smooth stationary solution of the Euler system (1.3) is perturbed, the solution is restored once the wave generated by the perturbation goes away.

In view of Figures 26, 27, 28, 29, 30, 31, 32, 33, 34, 35 and Table 8 we arrive at the following observation.

Conclusion 4

Consider a perturbation $\delta = (\delta _v, \delta _\rho )$ added to a steady shock solution of the form

$$\begin{aligned} V_{0}(r) = {\left\{ \begin{array}{ll} V_{-}^{*}(r), &{} { r\le r_0},\\ V_{+}^{*}(r), &{} \text { otherwise.} \end{array}\right. } \end{aligned}$$

• 1. If the perturbation moves the steady shock to the right, then a different stationary solution of the form

  $$\begin{aligned} V(r) = {\left\{ \begin{array}{ll} V_{-}^{*}(r), &{} { r\le r_1},\\ V_{+}^{*}(r), &{} \text { otherwise,} \end{array}\right. } \end{aligned}$$

  with $r_0 \ne 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.