KRYLOVâRIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS ...
c 2012 Society for Industrial and Applied Mathematics
SIAM J. SCI. COMPUT. Vol. 34, No. 4, pp. A2072–A2091
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS OF CONSERVATION LAWS∗ MANUEL TORRILHON† Abstract. This paper presents a Riemann solver for nonlinear hyperbolic systems of conservation laws based on a Krylov subspace approximation of the upwinding dissipation vector. In the general case, the solver does not require any detailed information of the eigensystem, except an estimate of the global maximal propagation speed. It uses successive flux function evaluations to obtain a numerical flux which is almost equivalent to that of a Godunov scheme with complete upwinding. The new Krylov–Riemann solver is particularly efficient when used for large systems with many nonlinear equations such that typically no explicit expression for the eigensystem is available. Also, no numerical procedures are necessary to compute the eigensystem. Numerical examples demonstrate the excellent performance of the solver with respect to other solvers. Key words. Riemann solvers, conservation laws, Krylov subspace methods, hyperbolic PDEs AMS subject classification. 65M08 DOI. 10.1137/110840832
1. Introduction. In the finite volume approach to numerical solutions of systems of hyperbolic conservation laws, the Riemann solver for cell interfaces is the building block of the numerical method; see, for example, the textbook [7]. After locally integrating the finite volume cells, a numerical flux function arises at the cell interfaces which, following Godunov [3], is typically evaluated by solving a local Riemann problem with initial states given by the adjacent cells. Over the last decades this has led to the development of many different Riemann solvers for various kinds of equations and applications; see, for example, the textbook [12]. The challenge is to develop a numerical flux function that satisfies all mathematical requirements, such as accuracy, high resolution, and monotonicity, and still allows a computationally efficient and easy implementation. The upwind scheme and its nonlinear extension, the approximate Riemann solver of Roe [9], rely on a characteristic decomposition based on the eigensystem of the flux Jacobian and typically yield the best resolution of the Riemann wave fan. Other flux functions such as Lax–Friedrichs or Harten–Lax–van Leer (HLL) can be viewed as approximations to the upwind scheme and reduce the use of characteristic information at the expense of less resolution of the wave fan. In many situations these simplified Riemann solvers give excellent results sufficient for the requirements at hand. However, the detailed upwinding method is preferable if waves need to be resolved specifically, for example, in the interaction with boundary conditions. Because upwind Riemann solvers are based on characteristic properties of the equations, they are typically developed specifically for a given system of equations, such as the shallow water system, Euler equations of gas dynamics, or magnetohydrodynamics (MHD); see [1], for example. Depending on the size and complexity of the system, the specific derivation may become increasingly complicated and eventually ∗ Submitted
to the journal’s Methods and Algorithms for Scientific Computing section July 14, 2011; accepted for publication (in revised form) April 17, 2012; published electronically July 19, 2012. This work was supported by the EURYI award of the European Science Foundation (ESF). http://www.siam.org/journals/sisc/34-4/84083.html † Center for Computational Science, RWTH Aachen University, 52062 Aachen, Germany (mt@ mathcces.rwth-aachen.de). A2072
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2073
impossible. For example, consider moment equations in extended fluid dynamics as in [17], where for a system with 13 variables no analytical expression for the eigenvalues is available. In these cases an upwind scheme is typically implemented by using generic numerical libraries to compute the eigensystem and decomposition of the given equations. In this paper we develop an upwind Riemann solver for nonlinear flux functions that does not require the implementation or direct computation of characteristic information except a possibly global estimate of the largest propagation speed. It relies on a specific approximation of the eigensystem by means of a Krylov subspace sequence and uses only iterative flux function evaluations. As such, the Riemann solver can be implemented as an almost black-box procedure applicable to general hyperbolic systems of nonlinear conservation laws. While it will not be competitive with existing upwind or Roe solvers for small systems, the new solver provides an efficient method for upwinding of large and complicated systems by avoiding the direct numerical computation of an eigensystem. The aim of the paper is to derive a Riemann solver for large and complex systems that is efficient and easy to implement without knowledge of the eigensystem and that reproduces the result of a complete upwind Godunov flux with best accuracy. We will not address the problem of how to increase the resolution of particular waves beyond what upwinding would give. Additionally, it will turn out that full equivalence to the Godunov upwinding scheme cannot be achieved as the Krylov–Riemann solver will not be able to resolve stationary waves, that is, waves with speed zero. The paper is organized as follows: In the next section we give a brief overview of Riemann solvers and present them in a combining framework by identifying their dissipation matrices. The Krylov–Riemann solver is developed in section 3 by constructing an appropriate polynomial approximation to the upwind dissipation function. The matrix evaluation of this function is then written in a Jacobian-free way. In section 4 some numerical results for large systems of conservation laws demonstrate the performance of the new solver. The paper ends with a conclusion, and an appendix shows the code of an implementation of the Krylov–Riemann solver used in the numerical examples. 2. Riemann solver. We consider an initial value problem for a system of conservation laws in one space dimension: (1) (2)
∂t U + ∂x F (U ) = 0
in R × R+ ,
U (x, 0) = U0 (x)
for the unknown variable vector U : R × R+ → RN . The flux function F : RN → RN and the initial conditions U0 are given. The system is assumed to be hyperbolic such that the Jacobian DF (U ) has only real-valued eigenvalues for U in the solution space. 2.1. Finite volume method. The equations are discretized according to a finite volume method [7] based on equidistant cells located at xi = i Δx, i ∈ Z, with width Δx, and the cell averages are given by xi+1/2 1 U (x, tn )dx (3) Uin = Δx xi−1/2 at time level tn . The update after one time step of size Δt for the cell average reads as Δt (num) (num) (4) Uin+1 = Uin + Fi−1/2 − Fi+1/2 Δx
A2074
MANUEL TORRILHON
and contains a numerical flux function (num) (5) Fi+1/2 = F˜ UL,i+1/2 , UR,i+1/2 , which takes as input variables the left and right limiting values of the variable vector U at the interface i + 1/2. The numerical flux function is the essential ingredient of the finite volume method, and we will discuss standard choices below. The numerical flux function is also referred to as the Riemann solver of the method. In this paper we will assume the general form for the numerical flux function given by (6)
1 1 F˜ (UL , UR ) = (F (UL ) + F (UR )) + D (UL , UR ) (UL − UR ) 2 2
with some dissipation matrix D depending on the left and right states and properties of the flux function of the system. We will denote the Jacobian of the flux function by A (U ) = DF (U ) and define its maximal and minimal eigenvalues or characteristic speed by (7)
λmax (U ) = max {λ : λ eigenvalue of A (U )} ,
(8)
λmin (U ) = min {λ : λ eigenvalue of A (U )} .
Correspondingly, the spectral radius at U , that is, the maximal absolute characteristic ¯ (U ) = max{|λmax | , |λmin |}. speed, is given by λ 2.2. Lax–Friedrichs solver. The simplest numerical flux function results from the choice DLF (UL , UR ) =
(9)
Δx I, Δt
with I the identity matrix, and yields the Lax–Friedrichs solver. It is a very robust but also very diffusive solver. The diffusivity can be reduced by choosing ¯ (UL ) , λ ¯ (UR ) I (10) DLLF (UL , UR ) = max λ based on the fastest speeds to the left and right of the interface. In this form, the Lax–Friedrichs flux is also known as the local Lax–Friedrichs or Rusanov flux; see [10]. Several refinements of the choice of speeds are possible. 2.3. HLL solver. By solving a simplified Riemann problem at the interface containing only two waves, Harten, Lax, and van Leer [6] arrived at a numerical flux function which can be written as (11) 1 α0 α1 (F (UL ) − F (UR )) − (UL − UR ) , F˜ HLL (UL , UR ) = (F (UL ) + F (UR )) + 2 2 2 where the coefficients are given by (12)
α0 =
|λL | − |λR | , λL − λR
α1 =
|λL | λR − |λR | λL λL − λR
with λL = λmin (UL ) the minimal left-hand speed, and λR = λmax (UR ) the maximal right-hand speed of the system. If we assume the existence of a so-called Roe matrix A˜ [9] satisfying (13)
A˜ (UL , UR ) (UL − UR ) = F (UL ) − F (UR ) ,
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2075
we can substitute it for the flux difference in the above formula. This allows us to identify the dissipation matrix of the HLL solver as (14)
DHLL (UL , UR ) =
|λL | λR − |λR | λL |λL | − |λR | ˜ I− A (UL , UR ) , λL − λR λL − λR
where we inserted the definitions of α0,1 . 2.4. Upwind Godunov solver. In the upwind scheme, one-sided spatial derivatives are chosen according to the wave propagation of the system. For a linear system with F (U ) = A U with some constant A ∈ RN ×N , the upwind scheme reads as (15)
Uin+1 = Uin +
Δt + A (Ui−1 − Ui ) + A− (Ui − Ui+1 ) Δx
with A± as the matrices containing only the positive/negative eigenvalues. This scheme solves the complete Riemann problem on the interface respecting all N resulting waves. We have A+ + A− = A and A+ − A− = |A| such that the numerical flux function of the upwind Riemann solver is (16)
1 1 F˜ up (UL , UR ) = A (UL + UR ) + |A| (UL − UR ) 2 2
and the dissipation matrix Dup (UL , UR ) = |A| is the modulus or absolute value of the flux matrix A, containing only positive eigenvalues. The finite volume approach with the solution of the Riemann problem on the interfaces is known as the Godunov method [3]. For a nonlinear system with general flux function, a similar expression can be derived when solving the approximative Riemann problem with local Roe matrix (13). The numerical flux function reads as (17)
1 1 F˜ Roe (UL , UR ) = (F (UL ) + F (UR )) + A˜ (UL , UR ) (UL − UR ) 2 2
and is generally referred to as the Roe solver. Its dissipation matrix is (18) DRoe (UL , UR ) = A˜ (UL , UR ) . Note that to evaluate this function, typically a characteristic decomposition of the Roe matrix is performed, which requires the knowledge of all eigenvalues and eigenvectors. These are either hard-coded for smaller systems or obtained via an external numerical method. 2.5. Lax–Wendroff method. When the exact solution is Taylor-expanded, for the linear system with F (U ) = A U , we arrive at the Lax–Wendroff method (19)
Uin+1 = Uin +
Δt Δt2 2 A (Ui−1 − Ui+1 ) − A (Ui−1 − 2Ui + Ui+1 ) , 2Δx 2Δx2
which exhibits second order accuracy in space and time. The numerical flux function can be identified as (20)
1 1 Δt 2 A (UL − UR ) . F˜ LW (UL , UR ) = A (UL + UR ) + 2 2 Δx
A2076
MANUEL TORRILHON
In the nonlinear setting this flux is typically obtained by
1 Δt LW ˜ (UL + UR ) − (F (UR ) − F (UL )) , (UL , UR ) = F (21) F 2 2Δx and the dissipation matrix given by (22)
DLW (UL , UR ) =
Δt ˜ 2 A (UL , UR ) Δx
contains the square of the Roe matrix. Note that the Lax–Wendroff flux function yields nonmonotone solutions especially for discontinuities. 2.6. FORCE flux. A possible monotone variant of the Lax–Wendroff flux is the so-called FORCE flux; see the textbook [12] for details. This flux is obtained from an average of the Lax–Wendroff method and the original Lax–Friedrichs method, viz. (23)
1 1 F˜ FORCE (UL , UR ) = F˜ LW (UL , UR ) + F˜ LF (UL , UR ) . 2 2
Consequently, the dissipation matrix is given by
Δt Δx I , (24) DFORCE (UL , UR ) = A2 + 2Δx Δt where the matrix A must be considered as a Roe matrix in the nonlinear setting. 2.7. Comparison. In general, we view the dissipation matrix as a function of the Jacobian or a Roe matrix either explicitly or through the characteristic speeds. The dissipation function D (A) is then written as function D (λ) of the eigenvalues of the Jacobian A, which is a scalar function. To compare the different functions we define the dimensionless characteristic speed ν = λΔt/Δx (Courant number). In dimensionless formulation we finally obtain the function
Δx Δt D ν (25) d (ν) = Δx Δt for the different schemes discussed above. The plot of Figure 1 gives a visualization. The figure represents the scalar setting, but equivalently holds for each wave in a sys¯ tem. Under the CFL condition ν¯ = λΔt/Δx ≤ 1 for the scheme (4), the Courant number ν based on any eigenvalue of the Jacobian ranges between the bounds −1 ≤ ν ≤ 1. Furthermore, it holds that d (ν) ≥ 0 for linear stability. Lax–Friedrichs, Rusanov, and HLL give linear functions for d (ν). HLL is shown for an exemplary situation in which the speed of the leftgoing wave is given by the Courant number νL = λL Δt/Δx, while the rightgoing wave moves with νR = λR Δt/Δx, both shown on the abscissa in the plot. The HLL flux can be obtained as a linear interpolation of the modulus function (upwind scheme) based on the values νR,L . This interpolation becomes exact if the values νR,L have the same sign, that is, all waves in the Riemann problem move in the same direction. The Lax–Wendroff method gives the least dissipation, but it results in oscillations for discontinuous solutions. The figure also shows how the average of Lax–Wendroff and Lax–Friedrichs gives rise to the FORCE flux. For scalar equations the theory of TVD-schemes [2] gives the condition (26)
d (ν) ≥ |ν|
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2077
1.0
d (º)
Lax-Friedrichs
Rusanov
HLL
FORCE upwind Lax-Wendroff 0.0
-1.0
ºL
0.0
ºR
1.0
Courant number º
Fig. 1. Comparison of dissipation functions of different numerical fluxes, namely Lax– Friedrichs, local Lax–Fiedrichs, HLL, upwind, and Lax–Wendroff. Nonmonotone schemes lie in the gray area.
for a monotone method, that is, nonoscillatory solutions. This means that the dissipation function with minimal dissipation and monotone behavior is that of the upwind or Roe scheme. The nonmonotonicity region is shaded in gray in Figure 1. The dissipation functions of Lax–Friedrichs, local Lax–Friedrichs, HLL, and FORCE can be viewed as successively better monotone approximations of the upwind function d(ν) = |ν|. In a similar spirit the HLL scheme can be tuned for specific systems like the Euler equations or MHD to become closer to upwind by introducing additional waves. These solvers are known as HLLC or HLLM; see [13] and [19], for example. 3. A Krylov–Riemann solver. In this section we derive a Riemann solver that avoids the computation of an eigensystem and yet reproduces almost all waves in a way equivalent to the upwind Godunov method. It therefore gives almost optimal dissipation for a monotone scheme, but requires only flux function evaluations and in the general case only an estimate for the globally fastest speed in the domain. The fundamental idea is to consider the Krylov subspace sequence of the Jacobian or Roe matrix A based on the vector ΔU = UL − UR and construct successive approximations for the dissipation vector |A| ΔU ; that is, we look for (27)
D (A) ΔU = |A| ΔU ∈ span(ΔU, AΔU, A2 ΔU, A3 ΔU, . . .)
with certain side constraints for stability and monotonicity. In a second step the resulting scheme will be formulated in a Jacobian-free way. 3.1. Polynomial approximation. The Krylov subspace approximation of the dissipation matrix represents a polynomial approximation to the modulus function, namely (28)
D (A) ΔU =
n
ak Ak ΔU = pn (A) ΔU
k=0
with coefficients ak . We construct an appropriate approximation for d (ν) in the scalar setting, where we use abs (x) = |x|.
A2078 1.0
MANUEL TORRILHON
dissipation function
pn (x) ¡ abs(x)
0.2
pn 0.1
n=4,...12
0.0
n=2,...12 0.0 -1.0
x
0.0
1.0
-0.1 -1.0
0.0
x
1.0
Fig. 2. Left: result for the dissipation function after minimization of the L1 -error. Right: error with respect to the modulus function.
Correspondingly, the absolute value abs (x) should be approximated by a polynomial with unknown coefficients ak , k = 0, 1, . . . , n, according to (29)
abs (x) ≈ pn (x) =
n
ak xk ,
k=0
where the degree n of the polynomial is considered fixed. In dimensionless variables we restrict x ∈ [−1, 1] and fix the extreme values pn (−1) = pn (−1) = 1. We also assume that the polynomial is symmetric with respect to the ordinate, i.e., n is even. These conditions can be built into some of the coefficients ak . The remaining coefficients of the polynomial pn are then the result of the minimization (30)
pn (x) − abs (x) → min
in an appropriate norm. For monotonicity of the resulting numerical flux function we require (31)
pn (x) ≥ abs (x)
such that the dissipation function does not drop below the upwind method. This constraint complicates the minimization problem; however, it cannot be ignored, since (30) alone typically gives an oscillatory result for pn . As an example we consider the minimization problem (30) with L1 -norm which already gives fewer oscillations than the L2 -norm. Figure 2 gives the result of this L1 -minimization for polynomial degrees n = 2, 4, . . . , 12 and the respective errors for the cases n = 4, 6, . . . , 12. Even though the approximation becomes better with increasing n, the deviation pn (x) − abs (x) shows negative values in violation of (31), which was not explicitly enforced. To enforce the constraint we consider the minimization 1 with |f (x)| = κ (f (x)) dx, (32) |pn (x) − abs (x)| → min −1
where κ (y) is given by (33)
κ (y) =
y, − 1ε y,
y ≥ 0, y < 0,
for some ε > 0. For finite ε this represents a nondifferentiable, but convex, minimization, which can be solved by standard computer algebra software. The result
A2079
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
dissipation function
1.0
pn (x) ¡ abs(x)
0.2
pn n=4,...12
0.1
0.0
n=2,...12 0.0 -1.0
0.0
x
1.0
-0.1 -1.0
x
0.0
1.0
Fig. 3. Left: result for dissipation function after minimization of the error with positivity constraint. Right: error with respect to the modulus function.
with ε = 10−2 is shown in Figure 3 for n = 2, 4, . . . , 12 together with the errors for the cases n = 4, 6, . . . , 12. Due to the constrained minimization the deviation pn (x) − abs (x) is essentially positive within an error of ε. The graph of p2 shows the exact result of the minimization, which is a symmetric parabola with unit slope at x = 1, that is, p2 (x) = 12 1 + x2 . Due to the nondifferentiability of abs (x) in the origin, convergence of pn is very slow close to x = 0. The same is true for the unconstrained L1 -minimization in Figure 2, but the constraint additionally increases the errors around x = 0. However, for larger n the constraint loses influence, and the results in Figure 3 become similar to the unconstrained ones above. The resulting coefficients of the optimal polynomials of Figure 3 are shown in Table 1. Only coefficients ak for k even are given; odd coefficients vanish. Due to the exact interpolation at x = ±1 the coefficients sum up to unity. The coefficients also give the typical behavior of larger coefficients with oscillating sign. Hence, in the implementation of these polynomials a different representation will be used. 3.2. Jacobian-free implementation. The numerical flux function of the Krylov–Riemann solver should give an accurate flux equivalent to the upwind flux after sufficiently many iterations, but the iterations should be designed such that at any stage a valid stable flux of a monotone method is computed. Following this requirement, we consider a recursive representation of the optimal polynomial pn in the form (m)
p2m (x) = α0
(34)
+ x2
m−1 k=0
(m)
αk+1 p2k (x)
Table 1 k Coefficients of the optimal polynomial dissipation functions pn (x) = n k=0 ak x . Odd coefficients vanish.
n 2 4 6 8 10 12
a0 0.5 0.201079 0.158864 0.119830 0.102425 0.086850
a2 0.5 1.33254 1.72916 2.32356 2.73512 3.24176
a4
a6
a8
a10
a12
-0.533619 -1.483070 -3.975520 -6.759890 -11.61570
0.595048 4.206920 11.46820 30.38110
-1.67478 -9.65750 -44.2558
3.11166 32.7875
-9.62571
A2080
MANUEL TORRILHON (m)
Coefficients αk
m 0 1 2 3 4 5 6
(m)
Table 2 of the recursive formula for the polynomial dissipation function (34).
(m)
α0
(m)
α1
1.0 0.5 0.201079 0.158864 0.119831 0.102425 0.086186
(m)
α2
0.5 1.86616 1.95052 1.80998 1.61711 1.32195
-1.06724 0.005735 1.94606 3.28221 4.89137
(m)
α3
-1.115120 -0.061604 0.950356 2.100740
α4
-2.81427 -3.09398 -2.45167
(m)
α5
-1.85812 -1.69004
(m)
α6
-3.25854
(m)
with coefficients αk , m ≥ 0, k = 0, 1, . . . , m − 1. The sum is empty for m = 0. The (m) coefficients αk of the recursion are given in Table 2. They have been found after recursively comparing expression (29) with (34) for increasing m. The new coefficients (m) αk remain in the order of unity and do not show strong sign oscillations. The above formula computes the value of p2m (x) based on the values of p2k (x) for k = 0, 1, . . . , m − 1. In the numerical flux function pn (x) is used as dissipation function d (ν), which is evaluated on the Jacobian or Roe matrix. Hence, x needs to be replaced by the matrix A scaled by some maximal speed that is not larger than the spectral radius of A. At this point we use a generic cmax for that purpose and discuss its role in the next section. The result of cmax p2m (A/cmax ) is then multiplied by the state difference ΔU = UL − UR , and we write V (m) = p2m (A/cmax ) ΔU . The recursive procedure reads for m ≥ 0 as (35)
(m)
V (m) = α0
(UL − UR ) +
1 c2max
(m) 2 m−1 ¯ A U αk+1 V (k) , k=0
where again the sum is empty for m = 0. The final numerical flux function of the Krylov–Riemann solver is given by 1 cmax (m) F˜ KR(m) (UL , UR ) = (F (UL ) + F (UR )) + V 2 2 with approximation order m. For any m > 0 the result of F˜ KR(m) gives a stable and monotone flux which for large m converges to the upwind or Godunov flux. Each iteration step requires a matrix-vector multiplication of the square of the Jacobian or Roe matrix A with a linear combination of earlier vectors V (k) . Very often the Jacobian is difficult to obtain or implement analytically. Since the Krylov–Riemann flux requires only the result of matrix-vector multiplications ¯ V = involving the Jacobian, it is easy to use the finite difference formulation DF U ¯ + εV ) − F (U ¯ ))/ε to replace the Jacobian A in the above formula. With limε→0 (F (U this, the recursion for the dissipation vector V (m) is given by (36)
(37)
V (0) = UL − UR ,
V˜ (m) =
m−1 k=0
(38)
V (m) =
(m) α0 V (0)
(m)
αk+1 V (k) ,
1 ¯ ¯ + εV˜ (m) ) − F U ¯ −F U ¯ , F U + F (U + 2 ε cmax
¯ . The finite differwhere ε needs to be chosen small in comparison to the norm of U ence expression for the directional derivative of the Jacobian is used here in a nested
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2081
2 ¯ V˜ . If the flux evaluation way in order to reproduce the square expression DF U (m) ¯ F U is stored separately, the computation of V requires 2m additional evaluations of F . 3.3. The role of cmax . The maximum speed that scales the polynomial p2m to become a dissipation function must be chosen such that all eigenvalues of the Jacobian at the local interface satisfy |λj | ≤ cmax . A value of cmax which is very close to the local maximal characteristic speed reduces the dissipation of the scheme. The choice of cmax depends on how much information about the eigenvalues is available. Any explicit finite volume scheme based on grid cells of size Δx must respect the ¯ (U ) ≤ Δx/Δt for the spectral radius global CFL condition to be convergent, namely λ of all cell values U . Thus, the simplest choice for cmax is cmax = Δx/Δt. With this choice no characteristic information enters the Krylov–Riemann solver at all, and it becomes fully black-box. For m = 0 the flux reduces to the classical Lax–Friedrichs flux, and for m = 1 the FORCE flux is recovered. Note, however, that in order to satisfy the global CFL condition some global ¯max for the spectral radius, that is, maximal characteristic speed from all estimate λ cells, must be available. The dissipation of the Krylov–Riemann solver can be reduced ¯ max is used globally for all cells. Additional reduction is obtained if a cmax = λ ¯ if some λ (U ) in each cell is available and is used to define cmax = local estimate ¯ ¯ max λ (UL ) , λ (UR ) for each interface evaluation of the Krylov–Riemann solver. In this case m = 0 gives the local Lax–Friedrichs or Rusanov flux. It follows that some characteristic information about the system is always needed. But very often the local estimates are allowed to be rough at the expense of some additional dissipation. These estimates often can also be based on a priori physical intuition; see section 4.3. 3.4. Discussion. The aim of the Krylov–Riemann solver is to replace an upwind Godunov solver with a procedure that requires only flux evaluations and no details about the Jacobian or its eigensystem. As seen above, up to 12 flux evaluations are necessary to achieve sufficient accuracy, and more may be needed in extreme situations. It becomes clear that for a small system such an approach may not be very efficient especially if analytical expressions are available for the eigensystem. However, for large systems and/or when explicit expressions for the eigensystem are not available the Krylov–Riemann solver turns into an efficient tool as the computation of an eigensystem and the decomposition of the difference of left and right states UL − UR becomes expensive. Even if analytical expressions are available the Krylov–Riemann solver might be more efficient if the time needed for implementation is considered. The use of only the flux function of the system and an estimate for the maximal local speed is identical to local Lax–Friedrichs. Hence, the Krylov–Riemann solver can be viewed as an extension of the local Lax–Friedrichs solver to incorporate full upwinding by using more flux evaluations. It can be implemented as a black-box solver for arbitrary hyperbolic systems as shown with the code in the appendix. For any order m in (36) the flux function gives a stable numerical flux which converges to the optimal upwind flux for large m. Depending on the actual distribution of the local characteristic speed of waves in the diagram of Figure 3, even a small value of m might be sufficient. For instance, if all characteristic speeds are clustered close to the maximal speed (positive and negative away from zero) even m = 4 might give a flux which is very close to the upwind scheme. The Krylov–Riemann solver in principle allows adaptivity in the order m, such that the iteration is stopped if
A2082
MANUEL TORRILHON
the dissipation vector V (m) does not change significantly in the next step. This approach will be exploited in future works, while this paper will use a constant m in the computations below. The nonlinear solver of order m as presented in (36)–(38) needs not only a total of 2(m + 1) flux evaluations but also memory storage of m + 1 arrays of the dimension of the system. The storage is due to the requirement that the solver give a valid flux at any iteration stage and the fact that the polynomials in Table 1 do not exhibit a three-term recursion for their construction. However, if the order is fixed for any usage of the solver, a low storage implementation is possible by implementation of a Horner rule for the particular polynomial of Table 1. Note, however, that the polynomial coefficients are less stable than in Table 2. Also, the additional memory storage is used only locally on an interface and is negligible in comparison to a full grid function. In contrast to a local Lax–Friedrichs solver, the HLL numerical flux requires the estimate of a maximal and a minimal characteristic speed. If both speeds have the same sign, the HLL flux produces exact upwinding by using the proper left or right flux evaluation. In the diagram of Figure 1 this means the left or right half of the modulus function is exactly approximated by a linear function. If this detailed information about the local speeds is available, such a procedure could also be incorporated in the Krylov–Riemann solver for reduced dissipation. Note, however, that this full upwinding is relatively rare in a general simulation and also that the improvement would act only on very particular waves. Even though the optimal dissipation function has the shape of an absolute value and corresponds to the upwind scheme, that function also has disadvantages. In the nonlinear situation upwinding is typically performed on a local linearization given by ˜ which solves the Riemann problem with only discontinuities. It is a Roe matrix A, well known that this yields an entropy violating scheme with incorrect rarefaction waves, for example, for gas dynamics in regions where the flow is close to the sonic point such that one characteristic speed is close to zero; see [7]. Several appropriate entropy fixes have been proposed, for instance, in [5], which essentially add dissipation to the slow waves. Ultimately, the dissipation functions of these corrected schemes are very similar to that of the Krylov–Riemann solver in Figure 3 for large m. In that sense, the weak accuracy of the Krylov–Riemann dissipation function for small values of characteristic speeds can be viewed as a built-in entropy fix of the solver.
3.5. Comparison to MUSTA flux. In [11] and [14] a multistage approach has been proposed to solve the Riemann problem at an interface to higher resolution. The idea is to choose a simple base flux and apply it in a repeated fashion assuming that both adjacent cells contain a subgrid on which the Riemann fan evolves. In this way the waves of the system can be resolved without additional characteristic information. A scheme using k steps beyond the base scheme is called MUSTAk . In principle the Riemann wave fan is completely resolved for the theoretical k → ∞; in practice only moderate values of k are used. In the following we compare in detail the relation between the MUSTA schemes and the Krylov–Riemann solver above. As base scheme the FORCE flux has been advocated in [11] and [14]. Starting (0) (0) with the values U1 = UL and U0 = UR , the MUSTA iteration can be written as
(39)
(k)
U1/2 =
Δt 1 (k) (k) (k) (k) − − F U0 , U0 + U1 F U1 2 Δx
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
(40)
(k)
A2083
Δx 1 (k) (k) (k) (k) (k) + 2F U1/2 + F U1 − , F U0 U1 − U0 4 4Δt Δt (k) (k) (k) , F1/2 − F U0 = U0 − Δx Δt (k) (k) (k) − F1/2 F U1 = U1 − Δx
F1/2 = (k+1)
(41)
U0
(42)
U1
(k+1)
for k = 0, 1, 2, . . . , and the corresponding flux is given by (43)
(k) F˜ MUSTAk (UL , UR ) = F1/2
for order k. Note that order k = 0 reproduces the FORCE flux. After some algebra it is easy to find the dissipation matrices for first order as a fourth order polynomial,
2 4
Δt 1 Δt Δx 1 I+ A − A (44) DMUSTA1 (UL , UR ) = , Δt 4 Δx 4 Δx and for second order as a sixth order polynomial, (45) Δx DMUSTA2 (UL , UR ) = Δt
1 11 I+ 8 8
2 4 6
Δt 5 Δt 1 Δt A − A + A . Δx 8 Δx 8 Δx
It makes sense to compare the corresponding scalar dissipation functions with those of the Krylov–Riemann solvers KR2 and KR3 , which are of the same polynomial degree. Note that the MUSTA fluxes are based on minimal characteristic information; i.e., they essentially use cmax = Δx/Δt. Hence, this setting is also used in the Krylov– Riemann case. Figure 4 displays the four polynomials together with the optimal dissipation function represented by the modulus function. The right half of the plot shows the fourth order polynomials of MUSTA1 and KR2 . In comparison to MUSTA, the Krylov–Riemann result clearly shows smaller dissipation towards zero, but at the expense of some increased dissipation for faster waves. Interestingly, the MUSTA flux drops slightly below the monotonicity boundary given by the upwind scheme. This phenomenon is more pronounced on the left-hand side of the plot, which shows the sixth order polynomials of MUSTA2 and KR3 . While the KR-polynomial follows the upwind method very accurately, the MUSTA flux enters the nonmonotone region and yields lower dissipation for very slow waves. The MUSTA fluxes are far away from the strong oscillations of the nonmonotone Lax–Wendroff flux, but the nonmonotone behavior can be observed in numerical solutions. We consider the scalar wave equation ut + ux = 0 and use as initial condition a discontinuity with jump from 0 to 1 on a sufficiently large interval. The discontinuity is evolved in time on a grid with 200 grid points according to the Lax–Wendroff method as well as the MUSTA2 and KR3 schemes. To demonstrate the effect we choose a Courant number of ν = 0.43 which shows the maximal deviation from the monotonicity region for MUSTA2 . As a measure for oscillations in the solution, Figure 5 shows the maximum value of the solution u for the initial 50 time steps. For Lax–Wendroff an overshoot of 20% quickly occurs and remains in the solution. The Krylov–Riemann solver does not introduce any oscillations. The profile is diffused with an amount equivalent to an upwind method (not shown in the figure). The MUSTA solver shows nonmonotone behavior with oscillations of up to 5% of the initial discontinuity. Interestingly, however, the oscillations are damped for longer times, and the profile becomes essentially monotone.
A2084
MANUEL TORRILHON
dissipation functions
1.0
d (º) KR3
MUSTA1 MUSTA2 KR2 0.0 -1.0
0.0
Courant number º
1.0
Fig. 4. Comparison of the dissipation functions of the Krylov–Riemann solver with the MUSTA flux. The right-hand side shows KR2 and M U ST A1 ; the left-hand side shows KR3 and M U ST A2 . In both cases the polynomial degrees are the same, 4 and 6, respectively.
maximum value ku (¢ ; tk )k1 1.20
Lax-Wendroff
1.15
1.10
MUSTA22 MUSTA 1.05
KR3 1.00 0
10
20
30
40
50 time step k
Fig. 5. Maximum value at a discontinuity with an initial jump from 0 to 1 as a function of the number of time steps. The result is essentially independent of the spatial resolution.
Altogether, the performance of the MUSTA and KR schemes is similar. The biggest differences lie in the philosophy of the derivation where the Krylov–Riemann solver uses an approach based on the abstract mathematical approximation of the modulus function and MUSTA follows the insight about the Riemann problem. Note, however, that due to the optimized character of the Krylov–Riemann solver the KR3 flux uses only 8 flux evaluations, while MUSTA2 requires 9. In general, a MUSTA flux of order k + 1 uses 3k flux evaluations, which must be compared to 2(k + 1) evaluations of the KRk flux yielding a dissipation polynomial of the same order. 4. Numerical examples. We will investigate the numerical performance of the Krylov–Riemann solver on a large linear system and the nonlinear systems of magnetohydrodynamics and extended thermodynamics. To focus on the Riemann solver and following the aims of this paper, the tests are conducted with first order methods and compared to the results obtained with an upwind Godunov solver. An extension
A2085
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
to higher order is straightforward and typically further increases the resolution of the waves. The Krylov–Riemann solver uses a nonadaptive, straightforward implementation with a fixed order m. The upwind method or Roe solver used in the following is based on an implementation for general equations. It uses finite differences to calculate the Jacobian from the flux function and a numerical library to obtain the eigensystem. Finally, a generic linear solver from the numerical library is employed to compute the dissipation vector. For many popular systems much more efficient implementations of the Godunov upwind schemes are available, typically relying on analytical expressions. However, in the case of extended thermodynamics (section 4.3) no analytical expressions are available. 4.1. Large linear system. We consider the linear system (46)
where A = (aij ) ∈ RN ×N
∂t U + A∂x U = 0,
with aij = tan (2ij) /50 and N = 20. Due to the symmetry of A, the system is hyperbolic, and the 20 characteristic speeds are nonuniformly distributed between λmin = −1.012 and λmax = 0.954. The initial conditions are those of a Riemann problem, viz. (47)
Ui (x < 0, t = 0) = 3,
Ui (x > 0, t = 0) = 1,
i = 1, 2, . . . , 20.
The time step is chosen such that for the maximal Courant number ν¯ = 0.97 holds and the simulation runs up to t = 0.9 on a computational domain with x ∈ [−1, 1]. An exact solution is available after diagonalization of the Jacobian and solving numerically for the jumps of the waves. Figure 6 shows the result for components U3 and U7 obtained with the Krylov– Riemann solver on 400 cells compared to the upwind method and local Lax–Friedrichs using the same number of grid points. We use cmax = 1.1 for the Krylov–Riemann solver, and the left-hand plot uses m = 6, while the right-hand side uses m = 2. The Krylov–Riemann result for m = 6 (left) follows the upwind result almost exactly except for the wave close to zero. This is explained by the discrepancy between the dissipation function of both methods close to zero as shown in Figure 3. The Krylov– Riemann solver adds more dissipation and produces a less steep wave at the origin
Ui
Ui
Upwind / Krylov
3.0
Upwind
3.0
exact
exact
2.0
2.0
1.0
1.0
local LF
0.0 -1.0
Krylov
local LF
0.0 0.0
x
1.0
-1.0
0.0
x
1.0
Fig. 6. Solution of a Riemann problem for a large linear system with 20 waves. The plot compares the result of the upwind method to the results of the Krylov–Riemann solver and local Lax–Friedrichs based on 400 cells. The left-hand side uses m = 6, and the right-hand side uses m = 2 for the Krylov–Riemann solver.
A2086
MANUEL TORRILHON
in Figure 6. In the right plot (m = 2) the dissipation is increased for slower waves, and the deviation spreads to waves more distant from the origin. This corresponds to the discrepancy between the dissipation functions p4 and p12 in Figure 3. The result of the local Lax–Friedrichs is given for additional comparison. It approximates the fast waves with high accuracy. Again this is where the dissipation functions of the methods coincide. 4.2. Magnetohydrodynamics. The equations of ideal magnetohydrodynamics describe the flow of highly conducting plasma and read in one space dimension as
(48)
∂t ρ ∂t ρvx ∂t ρvt ∂t Bt ∂t E
+ + + + +
∂x (ρv x) ∂x ρvx2 + p + 12 B2t ∂x (ρvx vt − Bx Bt ) ∂x (v x Bt − Bx vt ) ∂x E + p + 12 B2t vx − Bx Bt · vt
= = = = =
0, 0, 0, 0, 0
for the density ρ, normal velocity vn , tangential velocity vt = (vy , vz ), tangential magnetic field Bt = (By , Bz ), and the total energy E given in terms of the pressure p by (49)
E=
1 1 1 1 p + ρvx2 + ρvt2 + B2t , γ−1 2 2 2
where γ is the adiabatic constant. The normal magnetic field Bx is constant in onedimensional processes and appears as a parameter in the equations. Thus the system contains 8 equations for 8 unknowns U = (ρ, vx , vt , Bt , p). There exist various explicit Riemann solvers that are based on the characteristic properties of the MHD system, that is, explicit eigenvalues and eigenvectors; see [1]. However, the expressions exhibit a significant complexity and are probably at the limit of what is doable analytically. Even though 8 equations are not considered a very large system in the setting of this paper, we use the MHD system to demonstrate the performance of the new Krylov– Riemann solver. We consider the Riemann problem given by the initial conditions U (1) = const for x < 0 and U (0) = const for x > 0, with ρ1 , vx(1) , vy(1) , vz(1) , By(1) , Bz(1) , p1 = (3, −1.2, 0, 0, 1, 0, 3) , (50) ρ1 , vx(1) , vy(1) , vz(1) , By(1) , Bz(1) , p1 = (1, −1.2, 0, 0, cos α, sin α, 1) , (51) where α = 1.5 (in radians). The normal magnetic field is set to Bx = 1.5. The solution of this Riemann problem triggers all 7 MHD waves as shown in [16]. Note that the normal velocity is constant and nonvanishing in the initial conditions. The value is chosen such that one of the slow magnetic waves moves with a very small velocity. The computational domain is x ∈ [−1.0, 0.5] with end time t = 0.4 and maximal Courant number ν¯ = 0.95. The Riemann problem is solved with an upwind solver as in (17). However, for simplicity the Roe matrix A˜ is chosen as the Jacobian evaluated at the arithmetic average of the left and right states. Similarly, the solution using the Krylov–Riemann ¯ in (38). The value of cmax solver (36) uses the arithmetic mean for the state U is computed from the local values of the maximal characteristic speed for which an explicit expression is available, and the order is set m = 6. The result for the magnetic field component By is shown in Figure 7 for both Riemann solvers together with an
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2087
By 1.0 0.8
exact
0.6 0.4
Upwind / Krylov
0.2 0.0 -1.0
-0.5
0.0
x
0.5
Fig. 7. Solution of a Riemann problem for magnetohydrodynamics.
exact solution [15]. The solution shows 7 waves in total of which 5 are visible in the By -component. To the right there are a fast shock, a rotational linear Alfven wave, and a slow shock. The contact discontinuity is not visible, but a leftgoing slow rarefaction and a leftgoing rotational wave can be observed. A leftgoing fast rarefaction wave does not affect By and is not visible. The solution obtained by the Krylov–Riemann solver is virtually identical to that of the full upwind scheme except for the slow shock. This shock moves slowly to the right due to the imposed normal velocity in the initial conditions. Similar to the slow waves in the large linear system above, this wave shows more smearing due to the dissipation function of the Krylov– Riemann solver; see Figure 3. However, the difference is clearly smaller than in the linear situation above. This is due to nonlinear steepening of the shock waves. 4.3. Extended thermodynamics in three dimensions. Moment equations as used in [17] are derived in kinetic gas theory in order to describe gas processes with thermal nonequilibrium with higher physical accuracy than in classical models. This higher accuracy is achieved by extending the set of variables, e.g., considering not only density, velocity, and pressure, but, additionally, the pressure tensor and the heat flux vector. These extended thermodynamic models typically contain a hyperbolic core which is supplemented by relaxational and diffusive expressions. The hyperbolic core of the 13-moment equations of Grad [4] in full three-dimensional formulation is given by ∂ ∂ ( vk ) = 0, ρ+ ∂t ∂xk ∂ ∂ (ρvi ) + (53) (ρvi vk + pik ) = 0, ∂t ∂xk
∂ ∂ 6 (pij + ρvi vj ) + (54) ρvi vj vk + 3p(ij vk) + δ(ij qk) = 0, ∂t ∂xk 5 2 2 14 1 (55) ∂t Qi + ∂j E vi vj + 2vk pk(i vj) + 5 qk vk δij + 5 q(i vj) + 2 pij vk + 5θpδij = 0 (52)
and uses the variables density ρ, velocity vi , pressure tensor pij , and heat flux qi — a total of 13 fields. The equations can be viewed as a coarse approximation to a collisionless monatomic ideal gas. In the equations the total energy flux Qi = E vi +
A2088
MANUEL TORRILHON
pik vk +qi with total energy E = 12 ρvj2 + 23 p is used as well as the tensorial abbreviations (56)
a(i bj) =
1 (ai bj + aj bi ) , 2
δ(ij qk) =
1 (δij qk + δjk qi + δik qj ) . 3
Due to its complexity and nonlinearity, only very limited analytical investigations are possible for this system; see, e.g., [8]. In particular, explicit general expressions for the eigenvalues and eigenvectors of the Jacobian are not available. Still, precise upwinding is sometimes necessary, for example, when dealing with nontrivial boundary models as described in [18]. We will consider the system (52)–(55) for a one-dimensional Riemann problem in the direction of the vector n = ( 23 , 13 , 23 )T and with initial data U (1) = const for x < 0, and U (0) = const for x > 0, given by ⎛ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎞ 0 3 0 0 0 (1) (1) (1) ρ1 , vi , pij , qi = ⎝3, ⎝ 0.1 ⎠ , ⎝ 0 3 0 ⎠ , ⎝ 0 ⎠⎠ , (57) 0 0 0 3 0 ⎛ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎞ 0 1 0 0 0 (0) (0) (0) = ⎝1, ⎝ 0 ⎠ , ⎝ 0 1 0 ⎠ , ⎝ 0 ⎠⎠ . (58) ρ0 , vi , pij , qi 0.1 0 0 1 0 The system contains a total of 9 waves, including a 5-fold contact wave, and the initial conditions trigger 8 of these waves. The computational domain is x ∈ [−2, 2.4] discretized with 200 cells, and the end time is t = 0.9. We use a maximal Courant number of ν¯ = 0.8. In the present case the precise evaluation of the maximal char¯ (U ) = v + α cs , acteristic speed is very difficult. A rough estimate is of the form λ where cs is the speed of sound and α some safety factor; see [17]. We use α = 2, and ¯ (U ). the value of cmax in the Krylov–Riemann solver is based on the local values of λ As in the MHD case we choose m = 6. Figure 8 compares the solutions of the Krylov–Riemann scheme with one obtained from the generic implementation of the Godunov method using the arithmetic mean as Roe values. The usage of numerical libraries to compute the eigensystem for the Roe flux is mandatory for these equations, and consequently the Krylov–Riemann solver showed a computational time which was roughly 15 times faster than the generic upwind scheme. However, the number can be viewed only as a tendency, as no special effort on code optimization has been made. 0.4
0.4
vz
qz
0.2
0.3
vy
0.0
0.2
pxz -0.2
0.1 -0.4 0.0 -2.0
-1.0
0.0
1.0
2.0
x
-2.0
-1.0
0.0
1.0
2.0
x
Fig. 8. Riemann problem for the system of collisionless extended thermodynamics with high resolution reference solution as thin line. The figure shows the fields of velocity components vy and vz (left) and heat flux component qz together with shear stress pxz (right). The results of both the Roe scheme and the Krylov–Riemann method are shown but are indistinguishable to the eye.
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2089
0.004 0.002 0.0 -0.002 -0.004 -2.0
qz(KR) ¡ qz(Upwind) -1.0
0.0
1.0
2.0
x
Fig. 9. Difference between the solutions based on the Roe upwind method and the Krylov– Riemann solver for the heat flux component qz .
The plots in Figure 8 show only selected fields and a reference solution as thin lines obtained on 6000 cells with a Rusanov flux. The difference between the solutions of both methods is not visible to the eye. To get an impression of the deviation, Figure 9 shows the difference of both solutions for the heat flux component qz plotted for each space point. Relative to the global maximum of the heat flux, the plot shows a maximum deviation of 1%. The other fields show similar behavior. 5. Conclusion. This paper presented an approximative upwind Riemann solver for nonlinear hyperbolic systems of conservation laws based on a Krylov subspace approximation of the upwinding dissipation vector. In this way, the solver does not require any explicit characteristic information, such as eigenvalues or eigenvectors, of the equation, except for an estimate of the maximal propagation speed. It uses successive flux function evaluations to obtain a numerical flux which is almost equivalent to that of an approximate Roe scheme or complete Godunov upwinding. Extremely slow waves can be captured in the same way as in Godunov upwinding only if very high orders m of the new solver are used. Stationary waves will always be approximated with slightly higher diffusivity than in the case of complete upwinding. The new Krylov–Riemann solver has been designed to be efficient when used for large complex systems with many nonlinear equations such that, typically, no explicit expression for the eigensystem is available. Also, no numerical procedures are necessary to compute the eigensystem. Numerical examples demonstrated the excellent performance of the solver with respect to complete upwinding. The increased dissipation for very slow waves can be interpreted as a built-in entropy fix. Appendix. Example code of the solver. Below we give an example of a code for the Krylov–Riemann solver as described in (36) and (38) for a variable order m. For readability the code does not contain the explicit data of Table 2. #define N // number of equations #define NKR // available levels of Krylov-Riemann steps // input: left state UL[N], right state UR[N], maximal speed cMax // function fcn(U[],F[]) should give flux F[N] for state U[N] // output: numerical flux F[N] //
A2090
MANUEL TORRILHON
void FnumKrylovRiemann( double *UL, double *UR, double *F, void (*fcn)(double*,double*), double cMax ) { int n, m, k, idx, M = 6; // M is the order of the solver double alpha[] = { /* Table 2 */ }; double V[NKR][N], U[N], F0[N], F1[N], eps = 1.0e-7; for( n=0; n
SIAM J. SCI. COMPUT. Vol. 34, No. 4, pp. A2072–A2091
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS OF CONSERVATION LAWS∗ MANUEL TORRILHON† Abstract. This paper presents a Riemann solver for nonlinear hyperbolic systems of conservation laws based on a Krylov subspace approximation of the upwinding dissipation vector. In the general case, the solver does not require any detailed information of the eigensystem, except an estimate of the global maximal propagation speed. It uses successive flux function evaluations to obtain a numerical flux which is almost equivalent to that of a Godunov scheme with complete upwinding. The new Krylov–Riemann solver is particularly efficient when used for large systems with many nonlinear equations such that typically no explicit expression for the eigensystem is available. Also, no numerical procedures are necessary to compute the eigensystem. Numerical examples demonstrate the excellent performance of the solver with respect to other solvers. Key words. Riemann solvers, conservation laws, Krylov subspace methods, hyperbolic PDEs AMS subject classification. 65M08 DOI. 10.1137/110840832
1. Introduction. In the finite volume approach to numerical solutions of systems of hyperbolic conservation laws, the Riemann solver for cell interfaces is the building block of the numerical method; see, for example, the textbook [7]. After locally integrating the finite volume cells, a numerical flux function arises at the cell interfaces which, following Godunov [3], is typically evaluated by solving a local Riemann problem with initial states given by the adjacent cells. Over the last decades this has led to the development of many different Riemann solvers for various kinds of equations and applications; see, for example, the textbook [12]. The challenge is to develop a numerical flux function that satisfies all mathematical requirements, such as accuracy, high resolution, and monotonicity, and still allows a computationally efficient and easy implementation. The upwind scheme and its nonlinear extension, the approximate Riemann solver of Roe [9], rely on a characteristic decomposition based on the eigensystem of the flux Jacobian and typically yield the best resolution of the Riemann wave fan. Other flux functions such as Lax–Friedrichs or Harten–Lax–van Leer (HLL) can be viewed as approximations to the upwind scheme and reduce the use of characteristic information at the expense of less resolution of the wave fan. In many situations these simplified Riemann solvers give excellent results sufficient for the requirements at hand. However, the detailed upwinding method is preferable if waves need to be resolved specifically, for example, in the interaction with boundary conditions. Because upwind Riemann solvers are based on characteristic properties of the equations, they are typically developed specifically for a given system of equations, such as the shallow water system, Euler equations of gas dynamics, or magnetohydrodynamics (MHD); see [1], for example. Depending on the size and complexity of the system, the specific derivation may become increasingly complicated and eventually ∗ Submitted
to the journal’s Methods and Algorithms for Scientific Computing section July 14, 2011; accepted for publication (in revised form) April 17, 2012; published electronically July 19, 2012. This work was supported by the EURYI award of the European Science Foundation (ESF). http://www.siam.org/journals/sisc/34-4/84083.html † Center for Computational Science, RWTH Aachen University, 52062 Aachen, Germany (mt@ mathcces.rwth-aachen.de). A2072
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2073
impossible. For example, consider moment equations in extended fluid dynamics as in [17], where for a system with 13 variables no analytical expression for the eigenvalues is available. In these cases an upwind scheme is typically implemented by using generic numerical libraries to compute the eigensystem and decomposition of the given equations. In this paper we develop an upwind Riemann solver for nonlinear flux functions that does not require the implementation or direct computation of characteristic information except a possibly global estimate of the largest propagation speed. It relies on a specific approximation of the eigensystem by means of a Krylov subspace sequence and uses only iterative flux function evaluations. As such, the Riemann solver can be implemented as an almost black-box procedure applicable to general hyperbolic systems of nonlinear conservation laws. While it will not be competitive with existing upwind or Roe solvers for small systems, the new solver provides an efficient method for upwinding of large and complicated systems by avoiding the direct numerical computation of an eigensystem. The aim of the paper is to derive a Riemann solver for large and complex systems that is efficient and easy to implement without knowledge of the eigensystem and that reproduces the result of a complete upwind Godunov flux with best accuracy. We will not address the problem of how to increase the resolution of particular waves beyond what upwinding would give. Additionally, it will turn out that full equivalence to the Godunov upwinding scheme cannot be achieved as the Krylov–Riemann solver will not be able to resolve stationary waves, that is, waves with speed zero. The paper is organized as follows: In the next section we give a brief overview of Riemann solvers and present them in a combining framework by identifying their dissipation matrices. The Krylov–Riemann solver is developed in section 3 by constructing an appropriate polynomial approximation to the upwind dissipation function. The matrix evaluation of this function is then written in a Jacobian-free way. In section 4 some numerical results for large systems of conservation laws demonstrate the performance of the new solver. The paper ends with a conclusion, and an appendix shows the code of an implementation of the Krylov–Riemann solver used in the numerical examples. 2. Riemann solver. We consider an initial value problem for a system of conservation laws in one space dimension: (1) (2)
∂t U + ∂x F (U ) = 0
in R × R+ ,
U (x, 0) = U0 (x)
for the unknown variable vector U : R × R+ → RN . The flux function F : RN → RN and the initial conditions U0 are given. The system is assumed to be hyperbolic such that the Jacobian DF (U ) has only real-valued eigenvalues for U in the solution space. 2.1. Finite volume method. The equations are discretized according to a finite volume method [7] based on equidistant cells located at xi = i Δx, i ∈ Z, with width Δx, and the cell averages are given by xi+1/2 1 U (x, tn )dx (3) Uin = Δx xi−1/2 at time level tn . The update after one time step of size Δt for the cell average reads as Δt (num) (num) (4) Uin+1 = Uin + Fi−1/2 − Fi+1/2 Δx
A2074
MANUEL TORRILHON
and contains a numerical flux function (num) (5) Fi+1/2 = F˜ UL,i+1/2 , UR,i+1/2 , which takes as input variables the left and right limiting values of the variable vector U at the interface i + 1/2. The numerical flux function is the essential ingredient of the finite volume method, and we will discuss standard choices below. The numerical flux function is also referred to as the Riemann solver of the method. In this paper we will assume the general form for the numerical flux function given by (6)
1 1 F˜ (UL , UR ) = (F (UL ) + F (UR )) + D (UL , UR ) (UL − UR ) 2 2
with some dissipation matrix D depending on the left and right states and properties of the flux function of the system. We will denote the Jacobian of the flux function by A (U ) = DF (U ) and define its maximal and minimal eigenvalues or characteristic speed by (7)
λmax (U ) = max {λ : λ eigenvalue of A (U )} ,
(8)
λmin (U ) = min {λ : λ eigenvalue of A (U )} .
Correspondingly, the spectral radius at U , that is, the maximal absolute characteristic ¯ (U ) = max{|λmax | , |λmin |}. speed, is given by λ 2.2. Lax–Friedrichs solver. The simplest numerical flux function results from the choice DLF (UL , UR ) =
(9)
Δx I, Δt
with I the identity matrix, and yields the Lax–Friedrichs solver. It is a very robust but also very diffusive solver. The diffusivity can be reduced by choosing ¯ (UL ) , λ ¯ (UR ) I (10) DLLF (UL , UR ) = max λ based on the fastest speeds to the left and right of the interface. In this form, the Lax–Friedrichs flux is also known as the local Lax–Friedrichs or Rusanov flux; see [10]. Several refinements of the choice of speeds are possible. 2.3. HLL solver. By solving a simplified Riemann problem at the interface containing only two waves, Harten, Lax, and van Leer [6] arrived at a numerical flux function which can be written as (11) 1 α0 α1 (F (UL ) − F (UR )) − (UL − UR ) , F˜ HLL (UL , UR ) = (F (UL ) + F (UR )) + 2 2 2 where the coefficients are given by (12)
α0 =
|λL | − |λR | , λL − λR
α1 =
|λL | λR − |λR | λL λL − λR
with λL = λmin (UL ) the minimal left-hand speed, and λR = λmax (UR ) the maximal right-hand speed of the system. If we assume the existence of a so-called Roe matrix A˜ [9] satisfying (13)
A˜ (UL , UR ) (UL − UR ) = F (UL ) − F (UR ) ,
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2075
we can substitute it for the flux difference in the above formula. This allows us to identify the dissipation matrix of the HLL solver as (14)
DHLL (UL , UR ) =
|λL | λR − |λR | λL |λL | − |λR | ˜ I− A (UL , UR ) , λL − λR λL − λR
where we inserted the definitions of α0,1 . 2.4. Upwind Godunov solver. In the upwind scheme, one-sided spatial derivatives are chosen according to the wave propagation of the system. For a linear system with F (U ) = A U with some constant A ∈ RN ×N , the upwind scheme reads as (15)
Uin+1 = Uin +
Δt + A (Ui−1 − Ui ) + A− (Ui − Ui+1 ) Δx
with A± as the matrices containing only the positive/negative eigenvalues. This scheme solves the complete Riemann problem on the interface respecting all N resulting waves. We have A+ + A− = A and A+ − A− = |A| such that the numerical flux function of the upwind Riemann solver is (16)
1 1 F˜ up (UL , UR ) = A (UL + UR ) + |A| (UL − UR ) 2 2
and the dissipation matrix Dup (UL , UR ) = |A| is the modulus or absolute value of the flux matrix A, containing only positive eigenvalues. The finite volume approach with the solution of the Riemann problem on the interfaces is known as the Godunov method [3]. For a nonlinear system with general flux function, a similar expression can be derived when solving the approximative Riemann problem with local Roe matrix (13). The numerical flux function reads as (17)
1 1 F˜ Roe (UL , UR ) = (F (UL ) + F (UR )) + A˜ (UL , UR ) (UL − UR ) 2 2
and is generally referred to as the Roe solver. Its dissipation matrix is (18) DRoe (UL , UR ) = A˜ (UL , UR ) . Note that to evaluate this function, typically a characteristic decomposition of the Roe matrix is performed, which requires the knowledge of all eigenvalues and eigenvectors. These are either hard-coded for smaller systems or obtained via an external numerical method. 2.5. Lax–Wendroff method. When the exact solution is Taylor-expanded, for the linear system with F (U ) = A U , we arrive at the Lax–Wendroff method (19)
Uin+1 = Uin +
Δt Δt2 2 A (Ui−1 − Ui+1 ) − A (Ui−1 − 2Ui + Ui+1 ) , 2Δx 2Δx2
which exhibits second order accuracy in space and time. The numerical flux function can be identified as (20)
1 1 Δt 2 A (UL − UR ) . F˜ LW (UL , UR ) = A (UL + UR ) + 2 2 Δx
A2076
MANUEL TORRILHON
In the nonlinear setting this flux is typically obtained by
1 Δt LW ˜ (UL + UR ) − (F (UR ) − F (UL )) , (UL , UR ) = F (21) F 2 2Δx and the dissipation matrix given by (22)
DLW (UL , UR ) =
Δt ˜ 2 A (UL , UR ) Δx
contains the square of the Roe matrix. Note that the Lax–Wendroff flux function yields nonmonotone solutions especially for discontinuities. 2.6. FORCE flux. A possible monotone variant of the Lax–Wendroff flux is the so-called FORCE flux; see the textbook [12] for details. This flux is obtained from an average of the Lax–Wendroff method and the original Lax–Friedrichs method, viz. (23)
1 1 F˜ FORCE (UL , UR ) = F˜ LW (UL , UR ) + F˜ LF (UL , UR ) . 2 2
Consequently, the dissipation matrix is given by
Δt Δx I , (24) DFORCE (UL , UR ) = A2 + 2Δx Δt where the matrix A must be considered as a Roe matrix in the nonlinear setting. 2.7. Comparison. In general, we view the dissipation matrix as a function of the Jacobian or a Roe matrix either explicitly or through the characteristic speeds. The dissipation function D (A) is then written as function D (λ) of the eigenvalues of the Jacobian A, which is a scalar function. To compare the different functions we define the dimensionless characteristic speed ν = λΔt/Δx (Courant number). In dimensionless formulation we finally obtain the function
Δx Δt D ν (25) d (ν) = Δx Δt for the different schemes discussed above. The plot of Figure 1 gives a visualization. The figure represents the scalar setting, but equivalently holds for each wave in a sys¯ tem. Under the CFL condition ν¯ = λΔt/Δx ≤ 1 for the scheme (4), the Courant number ν based on any eigenvalue of the Jacobian ranges between the bounds −1 ≤ ν ≤ 1. Furthermore, it holds that d (ν) ≥ 0 for linear stability. Lax–Friedrichs, Rusanov, and HLL give linear functions for d (ν). HLL is shown for an exemplary situation in which the speed of the leftgoing wave is given by the Courant number νL = λL Δt/Δx, while the rightgoing wave moves with νR = λR Δt/Δx, both shown on the abscissa in the plot. The HLL flux can be obtained as a linear interpolation of the modulus function (upwind scheme) based on the values νR,L . This interpolation becomes exact if the values νR,L have the same sign, that is, all waves in the Riemann problem move in the same direction. The Lax–Wendroff method gives the least dissipation, but it results in oscillations for discontinuous solutions. The figure also shows how the average of Lax–Wendroff and Lax–Friedrichs gives rise to the FORCE flux. For scalar equations the theory of TVD-schemes [2] gives the condition (26)
d (ν) ≥ |ν|
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2077
1.0
d (º)
Lax-Friedrichs
Rusanov
HLL
FORCE upwind Lax-Wendroff 0.0
-1.0
ºL
0.0
ºR
1.0
Courant number º
Fig. 1. Comparison of dissipation functions of different numerical fluxes, namely Lax– Friedrichs, local Lax–Fiedrichs, HLL, upwind, and Lax–Wendroff. Nonmonotone schemes lie in the gray area.
for a monotone method, that is, nonoscillatory solutions. This means that the dissipation function with minimal dissipation and monotone behavior is that of the upwind or Roe scheme. The nonmonotonicity region is shaded in gray in Figure 1. The dissipation functions of Lax–Friedrichs, local Lax–Friedrichs, HLL, and FORCE can be viewed as successively better monotone approximations of the upwind function d(ν) = |ν|. In a similar spirit the HLL scheme can be tuned for specific systems like the Euler equations or MHD to become closer to upwind by introducing additional waves. These solvers are known as HLLC or HLLM; see [13] and [19], for example. 3. A Krylov–Riemann solver. In this section we derive a Riemann solver that avoids the computation of an eigensystem and yet reproduces almost all waves in a way equivalent to the upwind Godunov method. It therefore gives almost optimal dissipation for a monotone scheme, but requires only flux function evaluations and in the general case only an estimate for the globally fastest speed in the domain. The fundamental idea is to consider the Krylov subspace sequence of the Jacobian or Roe matrix A based on the vector ΔU = UL − UR and construct successive approximations for the dissipation vector |A| ΔU ; that is, we look for (27)
D (A) ΔU = |A| ΔU ∈ span(ΔU, AΔU, A2 ΔU, A3 ΔU, . . .)
with certain side constraints for stability and monotonicity. In a second step the resulting scheme will be formulated in a Jacobian-free way. 3.1. Polynomial approximation. The Krylov subspace approximation of the dissipation matrix represents a polynomial approximation to the modulus function, namely (28)
D (A) ΔU =
n
ak Ak ΔU = pn (A) ΔU
k=0
with coefficients ak . We construct an appropriate approximation for d (ν) in the scalar setting, where we use abs (x) = |x|.
A2078 1.0
MANUEL TORRILHON
dissipation function
pn (x) ¡ abs(x)
0.2
pn 0.1
n=4,...12
0.0
n=2,...12 0.0 -1.0
x
0.0
1.0
-0.1 -1.0
0.0
x
1.0
Fig. 2. Left: result for the dissipation function after minimization of the L1 -error. Right: error with respect to the modulus function.
Correspondingly, the absolute value abs (x) should be approximated by a polynomial with unknown coefficients ak , k = 0, 1, . . . , n, according to (29)
abs (x) ≈ pn (x) =
n
ak xk ,
k=0
where the degree n of the polynomial is considered fixed. In dimensionless variables we restrict x ∈ [−1, 1] and fix the extreme values pn (−1) = pn (−1) = 1. We also assume that the polynomial is symmetric with respect to the ordinate, i.e., n is even. These conditions can be built into some of the coefficients ak . The remaining coefficients of the polynomial pn are then the result of the minimization (30)
pn (x) − abs (x) → min
in an appropriate norm. For monotonicity of the resulting numerical flux function we require (31)
pn (x) ≥ abs (x)
such that the dissipation function does not drop below the upwind method. This constraint complicates the minimization problem; however, it cannot be ignored, since (30) alone typically gives an oscillatory result for pn . As an example we consider the minimization problem (30) with L1 -norm which already gives fewer oscillations than the L2 -norm. Figure 2 gives the result of this L1 -minimization for polynomial degrees n = 2, 4, . . . , 12 and the respective errors for the cases n = 4, 6, . . . , 12. Even though the approximation becomes better with increasing n, the deviation pn (x) − abs (x) shows negative values in violation of (31), which was not explicitly enforced. To enforce the constraint we consider the minimization 1 with |f (x)| = κ (f (x)) dx, (32) |pn (x) − abs (x)| → min −1
where κ (y) is given by (33)
κ (y) =
y, − 1ε y,
y ≥ 0, y < 0,
for some ε > 0. For finite ε this represents a nondifferentiable, but convex, minimization, which can be solved by standard computer algebra software. The result
A2079
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
dissipation function
1.0
pn (x) ¡ abs(x)
0.2
pn n=4,...12
0.1
0.0
n=2,...12 0.0 -1.0
0.0
x
1.0
-0.1 -1.0
x
0.0
1.0
Fig. 3. Left: result for dissipation function after minimization of the error with positivity constraint. Right: error with respect to the modulus function.
with ε = 10−2 is shown in Figure 3 for n = 2, 4, . . . , 12 together with the errors for the cases n = 4, 6, . . . , 12. Due to the constrained minimization the deviation pn (x) − abs (x) is essentially positive within an error of ε. The graph of p2 shows the exact result of the minimization, which is a symmetric parabola with unit slope at x = 1, that is, p2 (x) = 12 1 + x2 . Due to the nondifferentiability of abs (x) in the origin, convergence of pn is very slow close to x = 0. The same is true for the unconstrained L1 -minimization in Figure 2, but the constraint additionally increases the errors around x = 0. However, for larger n the constraint loses influence, and the results in Figure 3 become similar to the unconstrained ones above. The resulting coefficients of the optimal polynomials of Figure 3 are shown in Table 1. Only coefficients ak for k even are given; odd coefficients vanish. Due to the exact interpolation at x = ±1 the coefficients sum up to unity. The coefficients also give the typical behavior of larger coefficients with oscillating sign. Hence, in the implementation of these polynomials a different representation will be used. 3.2. Jacobian-free implementation. The numerical flux function of the Krylov–Riemann solver should give an accurate flux equivalent to the upwind flux after sufficiently many iterations, but the iterations should be designed such that at any stage a valid stable flux of a monotone method is computed. Following this requirement, we consider a recursive representation of the optimal polynomial pn in the form (m)
p2m (x) = α0
(34)
+ x2
m−1 k=0
(m)
αk+1 p2k (x)
Table 1 k Coefficients of the optimal polynomial dissipation functions pn (x) = n k=0 ak x . Odd coefficients vanish.
n 2 4 6 8 10 12
a0 0.5 0.201079 0.158864 0.119830 0.102425 0.086850
a2 0.5 1.33254 1.72916 2.32356 2.73512 3.24176
a4
a6
a8
a10
a12
-0.533619 -1.483070 -3.975520 -6.759890 -11.61570
0.595048 4.206920 11.46820 30.38110
-1.67478 -9.65750 -44.2558
3.11166 32.7875
-9.62571
A2080
MANUEL TORRILHON (m)
Coefficients αk
m 0 1 2 3 4 5 6
(m)
Table 2 of the recursive formula for the polynomial dissipation function (34).
(m)
α0
(m)
α1
1.0 0.5 0.201079 0.158864 0.119831 0.102425 0.086186
(m)
α2
0.5 1.86616 1.95052 1.80998 1.61711 1.32195
-1.06724 0.005735 1.94606 3.28221 4.89137
(m)
α3
-1.115120 -0.061604 0.950356 2.100740
α4
-2.81427 -3.09398 -2.45167
(m)
α5
-1.85812 -1.69004
(m)
α6
-3.25854
(m)
with coefficients αk , m ≥ 0, k = 0, 1, . . . , m − 1. The sum is empty for m = 0. The (m) coefficients αk of the recursion are given in Table 2. They have been found after recursively comparing expression (29) with (34) for increasing m. The new coefficients (m) αk remain in the order of unity and do not show strong sign oscillations. The above formula computes the value of p2m (x) based on the values of p2k (x) for k = 0, 1, . . . , m − 1. In the numerical flux function pn (x) is used as dissipation function d (ν), which is evaluated on the Jacobian or Roe matrix. Hence, x needs to be replaced by the matrix A scaled by some maximal speed that is not larger than the spectral radius of A. At this point we use a generic cmax for that purpose and discuss its role in the next section. The result of cmax p2m (A/cmax ) is then multiplied by the state difference ΔU = UL − UR , and we write V (m) = p2m (A/cmax ) ΔU . The recursive procedure reads for m ≥ 0 as (35)
(m)
V (m) = α0
(UL − UR ) +
1 c2max
(m) 2 m−1 ¯ A U αk+1 V (k) , k=0
where again the sum is empty for m = 0. The final numerical flux function of the Krylov–Riemann solver is given by 1 cmax (m) F˜ KR(m) (UL , UR ) = (F (UL ) + F (UR )) + V 2 2 with approximation order m. For any m > 0 the result of F˜ KR(m) gives a stable and monotone flux which for large m converges to the upwind or Godunov flux. Each iteration step requires a matrix-vector multiplication of the square of the Jacobian or Roe matrix A with a linear combination of earlier vectors V (k) . Very often the Jacobian is difficult to obtain or implement analytically. Since the Krylov–Riemann flux requires only the result of matrix-vector multiplications ¯ V = involving the Jacobian, it is easy to use the finite difference formulation DF U ¯ + εV ) − F (U ¯ ))/ε to replace the Jacobian A in the above formula. With limε→0 (F (U this, the recursion for the dissipation vector V (m) is given by (36)
(37)
V (0) = UL − UR ,
V˜ (m) =
m−1 k=0
(38)
V (m) =
(m) α0 V (0)
(m)
αk+1 V (k) ,
1 ¯ ¯ + εV˜ (m) ) − F U ¯ −F U ¯ , F U + F (U + 2 ε cmax
¯ . The finite differwhere ε needs to be chosen small in comparison to the norm of U ence expression for the directional derivative of the Jacobian is used here in a nested
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2081
2 ¯ V˜ . If the flux evaluation way in order to reproduce the square expression DF U (m) ¯ F U is stored separately, the computation of V requires 2m additional evaluations of F . 3.3. The role of cmax . The maximum speed that scales the polynomial p2m to become a dissipation function must be chosen such that all eigenvalues of the Jacobian at the local interface satisfy |λj | ≤ cmax . A value of cmax which is very close to the local maximal characteristic speed reduces the dissipation of the scheme. The choice of cmax depends on how much information about the eigenvalues is available. Any explicit finite volume scheme based on grid cells of size Δx must respect the ¯ (U ) ≤ Δx/Δt for the spectral radius global CFL condition to be convergent, namely λ of all cell values U . Thus, the simplest choice for cmax is cmax = Δx/Δt. With this choice no characteristic information enters the Krylov–Riemann solver at all, and it becomes fully black-box. For m = 0 the flux reduces to the classical Lax–Friedrichs flux, and for m = 1 the FORCE flux is recovered. Note, however, that in order to satisfy the global CFL condition some global ¯max for the spectral radius, that is, maximal characteristic speed from all estimate λ cells, must be available. The dissipation of the Krylov–Riemann solver can be reduced ¯ max is used globally for all cells. Additional reduction is obtained if a cmax = λ ¯ if some λ (U ) in each cell is available and is used to define cmax = local estimate ¯ ¯ max λ (UL ) , λ (UR ) for each interface evaluation of the Krylov–Riemann solver. In this case m = 0 gives the local Lax–Friedrichs or Rusanov flux. It follows that some characteristic information about the system is always needed. But very often the local estimates are allowed to be rough at the expense of some additional dissipation. These estimates often can also be based on a priori physical intuition; see section 4.3. 3.4. Discussion. The aim of the Krylov–Riemann solver is to replace an upwind Godunov solver with a procedure that requires only flux evaluations and no details about the Jacobian or its eigensystem. As seen above, up to 12 flux evaluations are necessary to achieve sufficient accuracy, and more may be needed in extreme situations. It becomes clear that for a small system such an approach may not be very efficient especially if analytical expressions are available for the eigensystem. However, for large systems and/or when explicit expressions for the eigensystem are not available the Krylov–Riemann solver turns into an efficient tool as the computation of an eigensystem and the decomposition of the difference of left and right states UL − UR becomes expensive. Even if analytical expressions are available the Krylov–Riemann solver might be more efficient if the time needed for implementation is considered. The use of only the flux function of the system and an estimate for the maximal local speed is identical to local Lax–Friedrichs. Hence, the Krylov–Riemann solver can be viewed as an extension of the local Lax–Friedrichs solver to incorporate full upwinding by using more flux evaluations. It can be implemented as a black-box solver for arbitrary hyperbolic systems as shown with the code in the appendix. For any order m in (36) the flux function gives a stable numerical flux which converges to the optimal upwind flux for large m. Depending on the actual distribution of the local characteristic speed of waves in the diagram of Figure 3, even a small value of m might be sufficient. For instance, if all characteristic speeds are clustered close to the maximal speed (positive and negative away from zero) even m = 4 might give a flux which is very close to the upwind scheme. The Krylov–Riemann solver in principle allows adaptivity in the order m, such that the iteration is stopped if
A2082
MANUEL TORRILHON
the dissipation vector V (m) does not change significantly in the next step. This approach will be exploited in future works, while this paper will use a constant m in the computations below. The nonlinear solver of order m as presented in (36)–(38) needs not only a total of 2(m + 1) flux evaluations but also memory storage of m + 1 arrays of the dimension of the system. The storage is due to the requirement that the solver give a valid flux at any iteration stage and the fact that the polynomials in Table 1 do not exhibit a three-term recursion for their construction. However, if the order is fixed for any usage of the solver, a low storage implementation is possible by implementation of a Horner rule for the particular polynomial of Table 1. Note, however, that the polynomial coefficients are less stable than in Table 2. Also, the additional memory storage is used only locally on an interface and is negligible in comparison to a full grid function. In contrast to a local Lax–Friedrichs solver, the HLL numerical flux requires the estimate of a maximal and a minimal characteristic speed. If both speeds have the same sign, the HLL flux produces exact upwinding by using the proper left or right flux evaluation. In the diagram of Figure 1 this means the left or right half of the modulus function is exactly approximated by a linear function. If this detailed information about the local speeds is available, such a procedure could also be incorporated in the Krylov–Riemann solver for reduced dissipation. Note, however, that this full upwinding is relatively rare in a general simulation and also that the improvement would act only on very particular waves. Even though the optimal dissipation function has the shape of an absolute value and corresponds to the upwind scheme, that function also has disadvantages. In the nonlinear situation upwinding is typically performed on a local linearization given by ˜ which solves the Riemann problem with only discontinuities. It is a Roe matrix A, well known that this yields an entropy violating scheme with incorrect rarefaction waves, for example, for gas dynamics in regions where the flow is close to the sonic point such that one characteristic speed is close to zero; see [7]. Several appropriate entropy fixes have been proposed, for instance, in [5], which essentially add dissipation to the slow waves. Ultimately, the dissipation functions of these corrected schemes are very similar to that of the Krylov–Riemann solver in Figure 3 for large m. In that sense, the weak accuracy of the Krylov–Riemann dissipation function for small values of characteristic speeds can be viewed as a built-in entropy fix of the solver.
3.5. Comparison to MUSTA flux. In [11] and [14] a multistage approach has been proposed to solve the Riemann problem at an interface to higher resolution. The idea is to choose a simple base flux and apply it in a repeated fashion assuming that both adjacent cells contain a subgrid on which the Riemann fan evolves. In this way the waves of the system can be resolved without additional characteristic information. A scheme using k steps beyond the base scheme is called MUSTAk . In principle the Riemann wave fan is completely resolved for the theoretical k → ∞; in practice only moderate values of k are used. In the following we compare in detail the relation between the MUSTA schemes and the Krylov–Riemann solver above. As base scheme the FORCE flux has been advocated in [11] and [14]. Starting (0) (0) with the values U1 = UL and U0 = UR , the MUSTA iteration can be written as
(39)
(k)
U1/2 =
Δt 1 (k) (k) (k) (k) − − F U0 , U0 + U1 F U1 2 Δx
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
(40)
(k)
A2083
Δx 1 (k) (k) (k) (k) (k) + 2F U1/2 + F U1 − , F U0 U1 − U0 4 4Δt Δt (k) (k) (k) , F1/2 − F U0 = U0 − Δx Δt (k) (k) (k) − F1/2 F U1 = U1 − Δx
F1/2 = (k+1)
(41)
U0
(42)
U1
(k+1)
for k = 0, 1, 2, . . . , and the corresponding flux is given by (43)
(k) F˜ MUSTAk (UL , UR ) = F1/2
for order k. Note that order k = 0 reproduces the FORCE flux. After some algebra it is easy to find the dissipation matrices for first order as a fourth order polynomial,
2 4
Δt 1 Δt Δx 1 I+ A − A (44) DMUSTA1 (UL , UR ) = , Δt 4 Δx 4 Δx and for second order as a sixth order polynomial, (45) Δx DMUSTA2 (UL , UR ) = Δt
1 11 I+ 8 8
2 4 6
Δt 5 Δt 1 Δt A − A + A . Δx 8 Δx 8 Δx
It makes sense to compare the corresponding scalar dissipation functions with those of the Krylov–Riemann solvers KR2 and KR3 , which are of the same polynomial degree. Note that the MUSTA fluxes are based on minimal characteristic information; i.e., they essentially use cmax = Δx/Δt. Hence, this setting is also used in the Krylov– Riemann case. Figure 4 displays the four polynomials together with the optimal dissipation function represented by the modulus function. The right half of the plot shows the fourth order polynomials of MUSTA1 and KR2 . In comparison to MUSTA, the Krylov–Riemann result clearly shows smaller dissipation towards zero, but at the expense of some increased dissipation for faster waves. Interestingly, the MUSTA flux drops slightly below the monotonicity boundary given by the upwind scheme. This phenomenon is more pronounced on the left-hand side of the plot, which shows the sixth order polynomials of MUSTA2 and KR3 . While the KR-polynomial follows the upwind method very accurately, the MUSTA flux enters the nonmonotone region and yields lower dissipation for very slow waves. The MUSTA fluxes are far away from the strong oscillations of the nonmonotone Lax–Wendroff flux, but the nonmonotone behavior can be observed in numerical solutions. We consider the scalar wave equation ut + ux = 0 and use as initial condition a discontinuity with jump from 0 to 1 on a sufficiently large interval. The discontinuity is evolved in time on a grid with 200 grid points according to the Lax–Wendroff method as well as the MUSTA2 and KR3 schemes. To demonstrate the effect we choose a Courant number of ν = 0.43 which shows the maximal deviation from the monotonicity region for MUSTA2 . As a measure for oscillations in the solution, Figure 5 shows the maximum value of the solution u for the initial 50 time steps. For Lax–Wendroff an overshoot of 20% quickly occurs and remains in the solution. The Krylov–Riemann solver does not introduce any oscillations. The profile is diffused with an amount equivalent to an upwind method (not shown in the figure). The MUSTA solver shows nonmonotone behavior with oscillations of up to 5% of the initial discontinuity. Interestingly, however, the oscillations are damped for longer times, and the profile becomes essentially monotone.
A2084
MANUEL TORRILHON
dissipation functions
1.0
d (º) KR3
MUSTA1 MUSTA2 KR2 0.0 -1.0
0.0
Courant number º
1.0
Fig. 4. Comparison of the dissipation functions of the Krylov–Riemann solver with the MUSTA flux. The right-hand side shows KR2 and M U ST A1 ; the left-hand side shows KR3 and M U ST A2 . In both cases the polynomial degrees are the same, 4 and 6, respectively.
maximum value ku (¢ ; tk )k1 1.20
Lax-Wendroff
1.15
1.10
MUSTA22 MUSTA 1.05
KR3 1.00 0
10
20
30
40
50 time step k
Fig. 5. Maximum value at a discontinuity with an initial jump from 0 to 1 as a function of the number of time steps. The result is essentially independent of the spatial resolution.
Altogether, the performance of the MUSTA and KR schemes is similar. The biggest differences lie in the philosophy of the derivation where the Krylov–Riemann solver uses an approach based on the abstract mathematical approximation of the modulus function and MUSTA follows the insight about the Riemann problem. Note, however, that due to the optimized character of the Krylov–Riemann solver the KR3 flux uses only 8 flux evaluations, while MUSTA2 requires 9. In general, a MUSTA flux of order k + 1 uses 3k flux evaluations, which must be compared to 2(k + 1) evaluations of the KRk flux yielding a dissipation polynomial of the same order. 4. Numerical examples. We will investigate the numerical performance of the Krylov–Riemann solver on a large linear system and the nonlinear systems of magnetohydrodynamics and extended thermodynamics. To focus on the Riemann solver and following the aims of this paper, the tests are conducted with first order methods and compared to the results obtained with an upwind Godunov solver. An extension
A2085
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
to higher order is straightforward and typically further increases the resolution of the waves. The Krylov–Riemann solver uses a nonadaptive, straightforward implementation with a fixed order m. The upwind method or Roe solver used in the following is based on an implementation for general equations. It uses finite differences to calculate the Jacobian from the flux function and a numerical library to obtain the eigensystem. Finally, a generic linear solver from the numerical library is employed to compute the dissipation vector. For many popular systems much more efficient implementations of the Godunov upwind schemes are available, typically relying on analytical expressions. However, in the case of extended thermodynamics (section 4.3) no analytical expressions are available. 4.1. Large linear system. We consider the linear system (46)
where A = (aij ) ∈ RN ×N
∂t U + A∂x U = 0,
with aij = tan (2ij) /50 and N = 20. Due to the symmetry of A, the system is hyperbolic, and the 20 characteristic speeds are nonuniformly distributed between λmin = −1.012 and λmax = 0.954. The initial conditions are those of a Riemann problem, viz. (47)
Ui (x < 0, t = 0) = 3,
Ui (x > 0, t = 0) = 1,
i = 1, 2, . . . , 20.
The time step is chosen such that for the maximal Courant number ν¯ = 0.97 holds and the simulation runs up to t = 0.9 on a computational domain with x ∈ [−1, 1]. An exact solution is available after diagonalization of the Jacobian and solving numerically for the jumps of the waves. Figure 6 shows the result for components U3 and U7 obtained with the Krylov– Riemann solver on 400 cells compared to the upwind method and local Lax–Friedrichs using the same number of grid points. We use cmax = 1.1 for the Krylov–Riemann solver, and the left-hand plot uses m = 6, while the right-hand side uses m = 2. The Krylov–Riemann result for m = 6 (left) follows the upwind result almost exactly except for the wave close to zero. This is explained by the discrepancy between the dissipation function of both methods close to zero as shown in Figure 3. The Krylov– Riemann solver adds more dissipation and produces a less steep wave at the origin
Ui
Ui
Upwind / Krylov
3.0
Upwind
3.0
exact
exact
2.0
2.0
1.0
1.0
local LF
0.0 -1.0
Krylov
local LF
0.0 0.0
x
1.0
-1.0
0.0
x
1.0
Fig. 6. Solution of a Riemann problem for a large linear system with 20 waves. The plot compares the result of the upwind method to the results of the Krylov–Riemann solver and local Lax–Friedrichs based on 400 cells. The left-hand side uses m = 6, and the right-hand side uses m = 2 for the Krylov–Riemann solver.
A2086
MANUEL TORRILHON
in Figure 6. In the right plot (m = 2) the dissipation is increased for slower waves, and the deviation spreads to waves more distant from the origin. This corresponds to the discrepancy between the dissipation functions p4 and p12 in Figure 3. The result of the local Lax–Friedrichs is given for additional comparison. It approximates the fast waves with high accuracy. Again this is where the dissipation functions of the methods coincide. 4.2. Magnetohydrodynamics. The equations of ideal magnetohydrodynamics describe the flow of highly conducting plasma and read in one space dimension as
(48)
∂t ρ ∂t ρvx ∂t ρvt ∂t Bt ∂t E
+ + + + +
∂x (ρv x) ∂x ρvx2 + p + 12 B2t ∂x (ρvx vt − Bx Bt ) ∂x (v x Bt − Bx vt ) ∂x E + p + 12 B2t vx − Bx Bt · vt
= = = = =
0, 0, 0, 0, 0
for the density ρ, normal velocity vn , tangential velocity vt = (vy , vz ), tangential magnetic field Bt = (By , Bz ), and the total energy E given in terms of the pressure p by (49)
E=
1 1 1 1 p + ρvx2 + ρvt2 + B2t , γ−1 2 2 2
where γ is the adiabatic constant. The normal magnetic field Bx is constant in onedimensional processes and appears as a parameter in the equations. Thus the system contains 8 equations for 8 unknowns U = (ρ, vx , vt , Bt , p). There exist various explicit Riemann solvers that are based on the characteristic properties of the MHD system, that is, explicit eigenvalues and eigenvectors; see [1]. However, the expressions exhibit a significant complexity and are probably at the limit of what is doable analytically. Even though 8 equations are not considered a very large system in the setting of this paper, we use the MHD system to demonstrate the performance of the new Krylov– Riemann solver. We consider the Riemann problem given by the initial conditions U (1) = const for x < 0 and U (0) = const for x > 0, with ρ1 , vx(1) , vy(1) , vz(1) , By(1) , Bz(1) , p1 = (3, −1.2, 0, 0, 1, 0, 3) , (50) ρ1 , vx(1) , vy(1) , vz(1) , By(1) , Bz(1) , p1 = (1, −1.2, 0, 0, cos α, sin α, 1) , (51) where α = 1.5 (in radians). The normal magnetic field is set to Bx = 1.5. The solution of this Riemann problem triggers all 7 MHD waves as shown in [16]. Note that the normal velocity is constant and nonvanishing in the initial conditions. The value is chosen such that one of the slow magnetic waves moves with a very small velocity. The computational domain is x ∈ [−1.0, 0.5] with end time t = 0.4 and maximal Courant number ν¯ = 0.95. The Riemann problem is solved with an upwind solver as in (17). However, for simplicity the Roe matrix A˜ is chosen as the Jacobian evaluated at the arithmetic average of the left and right states. Similarly, the solution using the Krylov–Riemann ¯ in (38). The value of cmax solver (36) uses the arithmetic mean for the state U is computed from the local values of the maximal characteristic speed for which an explicit expression is available, and the order is set m = 6. The result for the magnetic field component By is shown in Figure 7 for both Riemann solvers together with an
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2087
By 1.0 0.8
exact
0.6 0.4
Upwind / Krylov
0.2 0.0 -1.0
-0.5
0.0
x
0.5
Fig. 7. Solution of a Riemann problem for magnetohydrodynamics.
exact solution [15]. The solution shows 7 waves in total of which 5 are visible in the By -component. To the right there are a fast shock, a rotational linear Alfven wave, and a slow shock. The contact discontinuity is not visible, but a leftgoing slow rarefaction and a leftgoing rotational wave can be observed. A leftgoing fast rarefaction wave does not affect By and is not visible. The solution obtained by the Krylov–Riemann solver is virtually identical to that of the full upwind scheme except for the slow shock. This shock moves slowly to the right due to the imposed normal velocity in the initial conditions. Similar to the slow waves in the large linear system above, this wave shows more smearing due to the dissipation function of the Krylov– Riemann solver; see Figure 3. However, the difference is clearly smaller than in the linear situation above. This is due to nonlinear steepening of the shock waves. 4.3. Extended thermodynamics in three dimensions. Moment equations as used in [17] are derived in kinetic gas theory in order to describe gas processes with thermal nonequilibrium with higher physical accuracy than in classical models. This higher accuracy is achieved by extending the set of variables, e.g., considering not only density, velocity, and pressure, but, additionally, the pressure tensor and the heat flux vector. These extended thermodynamic models typically contain a hyperbolic core which is supplemented by relaxational and diffusive expressions. The hyperbolic core of the 13-moment equations of Grad [4] in full three-dimensional formulation is given by ∂ ∂ ( vk ) = 0, ρ+ ∂t ∂xk ∂ ∂ (ρvi ) + (53) (ρvi vk + pik ) = 0, ∂t ∂xk
∂ ∂ 6 (pij + ρvi vj ) + (54) ρvi vj vk + 3p(ij vk) + δ(ij qk) = 0, ∂t ∂xk 5 2 2 14 1 (55) ∂t Qi + ∂j E vi vj + 2vk pk(i vj) + 5 qk vk δij + 5 q(i vj) + 2 pij vk + 5θpδij = 0 (52)
and uses the variables density ρ, velocity vi , pressure tensor pij , and heat flux qi — a total of 13 fields. The equations can be viewed as a coarse approximation to a collisionless monatomic ideal gas. In the equations the total energy flux Qi = E vi +
A2088
MANUEL TORRILHON
pik vk +qi with total energy E = 12 ρvj2 + 23 p is used as well as the tensorial abbreviations (56)
a(i bj) =
1 (ai bj + aj bi ) , 2
δ(ij qk) =
1 (δij qk + δjk qi + δik qj ) . 3
Due to its complexity and nonlinearity, only very limited analytical investigations are possible for this system; see, e.g., [8]. In particular, explicit general expressions for the eigenvalues and eigenvectors of the Jacobian are not available. Still, precise upwinding is sometimes necessary, for example, when dealing with nontrivial boundary models as described in [18]. We will consider the system (52)–(55) for a one-dimensional Riemann problem in the direction of the vector n = ( 23 , 13 , 23 )T and with initial data U (1) = const for x < 0, and U (0) = const for x > 0, given by ⎛ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎞ 0 3 0 0 0 (1) (1) (1) ρ1 , vi , pij , qi = ⎝3, ⎝ 0.1 ⎠ , ⎝ 0 3 0 ⎠ , ⎝ 0 ⎠⎠ , (57) 0 0 0 3 0 ⎛ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎞ 0 1 0 0 0 (0) (0) (0) = ⎝1, ⎝ 0 ⎠ , ⎝ 0 1 0 ⎠ , ⎝ 0 ⎠⎠ . (58) ρ0 , vi , pij , qi 0.1 0 0 1 0 The system contains a total of 9 waves, including a 5-fold contact wave, and the initial conditions trigger 8 of these waves. The computational domain is x ∈ [−2, 2.4] discretized with 200 cells, and the end time is t = 0.9. We use a maximal Courant number of ν¯ = 0.8. In the present case the precise evaluation of the maximal char¯ (U ) = v + α cs , acteristic speed is very difficult. A rough estimate is of the form λ where cs is the speed of sound and α some safety factor; see [17]. We use α = 2, and ¯ (U ). the value of cmax in the Krylov–Riemann solver is based on the local values of λ As in the MHD case we choose m = 6. Figure 8 compares the solutions of the Krylov–Riemann scheme with one obtained from the generic implementation of the Godunov method using the arithmetic mean as Roe values. The usage of numerical libraries to compute the eigensystem for the Roe flux is mandatory for these equations, and consequently the Krylov–Riemann solver showed a computational time which was roughly 15 times faster than the generic upwind scheme. However, the number can be viewed only as a tendency, as no special effort on code optimization has been made. 0.4
0.4
vz
qz
0.2
0.3
vy
0.0
0.2
pxz -0.2
0.1 -0.4 0.0 -2.0
-1.0
0.0
1.0
2.0
x
-2.0
-1.0
0.0
1.0
2.0
x
Fig. 8. Riemann problem for the system of collisionless extended thermodynamics with high resolution reference solution as thin line. The figure shows the fields of velocity components vy and vz (left) and heat flux component qz together with shear stress pxz (right). The results of both the Roe scheme and the Krylov–Riemann method are shown but are indistinguishable to the eye.
KRYLOV–RIEMANN SOLVER FOR LARGE HYPERBOLIC SYSTEMS
A2089
0.004 0.002 0.0 -0.002 -0.004 -2.0
qz(KR) ¡ qz(Upwind) -1.0
0.0
1.0
2.0
x
Fig. 9. Difference between the solutions based on the Roe upwind method and the Krylov– Riemann solver for the heat flux component qz .
The plots in Figure 8 show only selected fields and a reference solution as thin lines obtained on 6000 cells with a Rusanov flux. The difference between the solutions of both methods is not visible to the eye. To get an impression of the deviation, Figure 9 shows the difference of both solutions for the heat flux component qz plotted for each space point. Relative to the global maximum of the heat flux, the plot shows a maximum deviation of 1%. The other fields show similar behavior. 5. Conclusion. This paper presented an approximative upwind Riemann solver for nonlinear hyperbolic systems of conservation laws based on a Krylov subspace approximation of the upwinding dissipation vector. In this way, the solver does not require any explicit characteristic information, such as eigenvalues or eigenvectors, of the equation, except for an estimate of the maximal propagation speed. It uses successive flux function evaluations to obtain a numerical flux which is almost equivalent to that of an approximate Roe scheme or complete Godunov upwinding. Extremely slow waves can be captured in the same way as in Godunov upwinding only if very high orders m of the new solver are used. Stationary waves will always be approximated with slightly higher diffusivity than in the case of complete upwinding. The new Krylov–Riemann solver has been designed to be efficient when used for large complex systems with many nonlinear equations such that, typically, no explicit expression for the eigensystem is available. Also, no numerical procedures are necessary to compute the eigensystem. Numerical examples demonstrated the excellent performance of the solver with respect to complete upwinding. The increased dissipation for very slow waves can be interpreted as a built-in entropy fix. Appendix. Example code of the solver. Below we give an example of a code for the Krylov–Riemann solver as described in (36) and (38) for a variable order m. For readability the code does not contain the explicit data of Table 2. #define N // number of equations #define NKR // available levels of Krylov-Riemann steps // input: left state UL[N], right state UR[N], maximal speed cMax // function fcn(U[],F[]) should give flux F[N] for state U[N] // output: numerical flux F[N] //
A2090
MANUEL TORRILHON
void FnumKrylovRiemann( double *UL, double *UR, double *F, void (*fcn)(double*,double*), double cMax ) { int n, m, k, idx, M = 6; // M is the order of the solver double alpha[] = { /* Table 2 */ }; double V[NKR][N], U[N], F0[N], F1[N], eps = 1.0e-7; for( n=0; n
