Conservation Properties of Iterative Methods for Implicit Discretizations of Conservation Laws

Conservation properties of iterative methods applied to implicit finite volume discretizations of nonlinear conservation laws are analyzed. It is shown that any consistent multistep or Runge-Kutta method is globally conservative. Further, it is shown that Newton’s method, Krylov subspace methods and pseudo-time iterations are globally conservative while the Jacobi and Gauss-Seidel methods are not in general. If pseudo-time iterations using an explicit Runge-Kutta method are applied to a locally conservative discretization, then the resulting scheme is also locally conservative. However, the corresponding numerical flux can be inconsistent with the conservation law. We prove an extension of the Lax-Wendroff theorem, which reveals that numerical solutions based on these methods converge to weak solutions of a modified conservation law where the flux function is multiplied by a particular constant. This constant depends on the choice of Runge-Kutta method but is independent of both the conservation law and the discretization. Consistency is maintained by ensuring that this constant equals unity and a strategy for achieving this is presented. Simulations show that this strategy improves the convergence rate of the pseudo-time iterations. Experiments with GMRES suggest that it also suffers from inconsistency but that this is automatically accounted for after some finite number of iterations. Similar experiments with coarse grid corrections based on agglomeration indicate no inconsistency.


Introduction
Conservation laws arise ubiquitously in the modelling of physical phenomena and their discretizations have been the subject of intense study; see e.g. [17,Chapter 1] and [18,Chapter 1]. They derive their name from the fact that they describe the conservation of some quantities of interest over time. In computational fluid dynamics (CFD), these quantities are usually mass, momentum and energy. For convenience we always refer to the conserved quantity as the "mass" in this paper.
Many numerical schemes have been designed to mimic mass conservation. Throughout, we refer to such schemes as globally conservative. As special cases, locally conservative schemes dictate that local variations in the mass propagate from one computational cell to neighboring ones without "skipping" any cells (typically referred to as "conservative schemes" in the literature). In particular, finite volume methods are designed in part upon this principle, although many other schemes possess the same property [13]. A major result on the topic is the Lax-Wendroff theorem, see e.g. [16], which provides sufficient conditions for numerical solutions of convergent and locally conservative discretizations to converge to weak solutions of the corresponding conservation law.
It is natural to ask whether such approximate solutions satisfy the conservation properties upon which the discretizations are based, and in particular, whether a Lax-Wendroff type result is available. It is typically claimed that this is not the case, e.g. in [15], where a version of a Newton-Jacobi iteration was noted to violate global conservation. In this paper we set out to answer this question in greater detail for some well-known (families of) iterative methods. We consider conservative implicit space-time discretizations whose solutions in each time step are approximated by a fixed number of iterations with different iterative methods. In practice, iterative methods are terminated using tolerances on the residual, set to match some predefined error tolerance on the entire scheme. Hence, the number of iterations is not known a priori. Nonetheless, in the pursuit of designing efficient iterative methods, it is of interest to understand their behaviour when terminated well before convergence. To this end, the assumption of a small number of iterations is beneficial.
After introducing relevant concepts and definitions in Sect. 2, we investigate the global conservation of iterative methods in Sect. 3. It is shown that Newton's method, Krylov subspace methods and pseudo-time iterations using Runge-Kutta (RK) methods are globally conservative if the initial guess has correct mass. However, the Jacobi and Gauss-Seidel methods in general are not. In Sect. 4 we focus on pseudo-time iterations using explicit RK methods. We show that such methods are locally conservative. By fixing the number of iterations and considering the limit of infinitesimal space-time increments, an extension of the Lax-Wendroff theorem is given. It turns out that the numerical solution in general converges to a weak solution of a modified conservation law for which the flux function is multiplied by a particular constant. An expression for this modification constant is given that only depends on the choice of Runge-Kutta method and the selected pseudo-time steps. A technique for ensuring that the constant equals one is presented, thereby ensuring consistency of the resulting scheme. Experiments indicate that the new technique improves the convergence rate of the pseudo-time iterations. Numerical tests corroborate these findings and further suggest that the results hold for systems of conservation laws in multiple dimensions. Further, experiments with GMRES suggest that they suffer from a similar modification as pseudo-time iterations, but that this is automatically accounted for after some finite number of iterations. Similar experiments with agglomeration coarse grid corrections indicate no presence of such modifications. Conclusions are drawn in Sect. 5.
Throughout the paper we separate quantities that belong to a spatial discretization from those that do not. Any vector living on a grid with m cells is denoted by a lower case bold letter, e.g. u ∈ R m . Similarly, any matrix operating on such a vector is represented by a bold upper case letter, e.g. A ∈ R m×m . In Sect. 4, we consider s-stage Runge-Kutta methods. Vector quantities spanning these stages are denoted by bold and underlined lower case letters, e.g. b ∈ R s . Matrices operating on the stages are represented by bold and underlined upper case letters e.g. A ∈ R s×s .

Conservation and Conservative Discretizations
In this paper we consider numerical methods for conservation laws, meaning partial differential equations of the form possibly with additional boundary conditions. Here, the kth component u k of the vector u ∈ R q represents the concentration of a quantity, i.e. the total amount of that quantity in a domain is given by u k dx. The conservation law states that this amount is only changed by flow across the boundary of the domain. Assuming that the net flow across the boundary is zero, we thus have for all t ∈ (t 0 , t e ] that In this paper we focus on scalar conservation laws in 1D, i.e. the case when d = q = 1 in (1). Numerical experiments in Sect. 4 suggest that the results of the paper can be generalized to systems in multiple dimensions, however this task is deferred for future work. Further, we restrict our analysis to the Cauchy problem for (1), or the problem with periodic boundary conditions. Thus, the conservation law of interest becomes where either = (−∞, ∞) or = (a, b] and u(a) = u(b). Which case is considered will be clear from context. A very successful line of research in computational fluid dynamics is to construct numerical methods that respect (2) on a discrete level. Herein, we predominantly consider finite volume methods that utilize the implicit Euler method as temporal discretization. Let the computational grid be given by (x i , t n ) = (i x, n t) with n ≥ 0. This grid may be either infinite or finite and periodic in space, i.e. x i+m = x i for some positive integer m. We consider discretizations on the form Here, u n i approximates the solution u(x, t n ) of (3) in the computational cell i . The notation f i+ 1 2 (u) is used as a placeholder for a generic numerical flux of the formf (u n+1 i− p , . . . , u n+1 i+q ), where p, q ∈ N determine the bandwidth of the stencil. Multiplying (4) by x and summing over all i reveals that Comparing this with (2) shows that the finite volume scheme discretely mimics conservation of mass.

Local Conservation
It is the telescoping nature of the spatial terms in (4) that leads to the preceeding result. In the finite volume community, this property is typically referred to simply as conservation. However, to distinguish it from other mass conservative schemes, we will henceforth refer to this property as local conservation: Definition 1 A discretization of (1) that can be written in the form (4) is said to be locally conservative.
Local conservation plays a central role in the Lax-Wendroff theorem, which considers the Cauchy problem for (3). First, the notion of consistency must be specified: A numerical flux function is inconsistent if it fails to satisfy Definition 2. We will see examples of this in Sect. 4.
The Lax-Wendroff theorem applies to the Cauchy problem for (3) and considers locally conservative discretizations with consistent numerical flux. If the numerical solution of such a scheme converges to a function u in the limit of vanishing x and t, the theorem provides sufficient conditions for u to be a weak solution of the conservation law (1) [17,Chapter 12]. More precisely, consider a sequence of grids ( x , t ) such that x , t → 0 as → ∞. Let U (x, t) denote the piecewise constant function that takes the solution value u n i in (x i , x i+1 ] × (t n−1 , t n ] on the th grid. We make the following assumptions: The first of these assumptions asserts that the numerical solution is convergent in the L 1 -norm with limit u. The second one states that the total variation of the numerical solution remains bounded independently of the grid. The Lax-Wendroff theorem may now be stated as follows: Theorem 1 (Lax-Wendroff) Consider the locally conservative discretization (4), suppose that the numerical fluxf is consistent and that Assumption 1 is satisfied. Then, u(x, t) is a weak solution of (1).
The theorem and its assumptions require that the linear or nonlinear systems arising in (4) are solved exactly. However, in practice these systems will be solved approximately using iterative methods. The natural question that now arises is: Do iterative solvers maintain the convergence to weak solutions?

Global Conservation
To answer this question we must establish if iterative methods preserve local conservation. However, local conservation is a special case of the more general notion of global conservation. Let us therefore consider a finite and periodic grid. Discretizing (3) in space on a mesh with cells i of volume | i |, we arrive at an initial value problem with u(t),f ∈ R m .
Global conservation implies that (2) is fulfilled on a semidiscrete level by (6): Suppose next that a time integration method is used to solve (6). (6) is

Definition 4 A time integration method applied to the conservative semidiscretization
holds for each time step n = 0, 1, . . .
Note that the conservation property (5) derived for the locally conservative finite volume discretization (4) is a special case of Definition 4, where i = x for every i. Global conservation is thus necessary for local conservation, which in turn is necessary for the Lax-Wendroff Theorem.
In the CFD community, one hears from time to time the claim that implicit methods are not conservative. However, this is not true: All consistent linear multistep methods are globally conservative. To see this, consider a generic s-step linear multistep method s j=0 a j u n+ j = t s j=0 b jf (u n+ j ), where a j and b j are given method dependent coefficients. Suppose that m i=1 | i |u . . , s − 1. Multiplying (7) by | i | and summing over all cells gives Consistent multistep methods satisfy a s = 1 and a 0 + · · · + a s−1 = −1. Hence, global conservation is ensured. Similarly, all Runge-Kutta (RK) methods are globally conservative. Consider an s-stage RK method where a j,l and b j are given by the method. Sincef is globally conservative it follows from Definition 3 that the mass of k j is zero and therefore that Hence, RK methods are globally conservative. Exponential integrators have been proven to be conservative as well [9], essentially due to the following useful lemma: Lemma 1 For any y ∈ R m the Jacobianf (u) = ∂f ∂u ∈ R m×m of the globally conservative discretizationf(u) satisfies Proof See [27].
In the remainder we also need the following result: Lemma 2 Letf (u) ∈ R m×m be the Jacobian of a conservative discretizationf(u) and let α be any scalar. The solution x of a linear system of the form Proof Multiplying the ith element of (10) by | i | and summing over all i gives By Lemma 1 the final term on the right-hand side is zero.

Globally Conservative Iterative Solvers
We begin by asking: Do iterative solvers respect global conservation? Since global conservation is necessary for local conservation, which in turn is used in Theorem 1, a negative answer will immediately exclude the possibility to extend the Lax-Wendroff theorem to incorporate a particular iterative method. Discretizing (6) in time using the implicit Euler method results in a nonlinear equation system of the form u − αf(u) = u n .
We concern ourselves with a set of iterates u (k) with limit u n+1 as k → ∞.
In the following subsections we investigate the conservation properties of some of the most common iterative methods.

Newton's Method
Newton's method finds an approximate solution to (11) by solving the sequence of linear systems for k = 0, 1, . . . Sincef is conservative it has zero mass by Definition 3. If u (k) and u n have the same mass, the entire right-hand side of the first equation in (12) has zero mass. Hence, by Lemma 2, the increment u also has zero mass. Thus, the mass of u (k+1) is the same as that of u (k) , which by assumption is the same as that of u n . Hence, we have shown:

Proposition 1 Newton's method is globally conservative.
In large scale practical applications, each linear system in (12) is solved approximately using other iterative methods. In the following, we consider a generic linear system that may either represent a Newton iteration or a linear discretization of a conservation law.
In line with Lemma 1, it is assumed that the mass of Ay is zero for any vector y ∈ R m .

Stationary Linear Methods
Any iterative method for solving (13) that can be written in the form for some matrices M and N, is termed a stationary linear method. Well known examples include the Richardson, Jacobi and Gauss-Seidel iterations.

Richardson Iteration
For the Richardson iteration, M = I − θ(I − αA) and N −1 = θ I, where θ ∈ R is a fixed parameter. With some minor rearrangements, the iteration can thus be expressed as Multiplying the ith element by | i | and summing over i results in By assumption, the matrix-vector product Au (k) has zero mass, hence the final term vanishes. Further, if u (k) has the same mass as b, then the second term on the right-hand side vanishes. By induction we therefore conclude:

Proposition 2
The Richardson iteration is globally conservative.

Jacobi Iteration
Consider again the linear system (13). With the matrix decomposition I − αA = I − α(D + L+U), where D, L and U are the diagonal, lower and upper triangular parts of A respectively, the Jacobi iteration is obtained by setting M = α(I − αD) −1 (L + U) and Inserting this into (14), multiplying by I − αD, then adding and subtracting αDu (k) , gives after some rearrangement Multiplying the ith element by | i | and summing over i results in The final term vanishes by assumption. However, the second term on the right-hand side is in general non-zero and thus introduces a conservation error

Proposition 3 The Jacobi iteration is not globally conservative in general.
However, consider the special case when the diagonal D of A is a scalar multiple of the identity matrix, say D = aI. Then a i,i = a for i = 1, . . . , m. If u (k) has the same mass as b we find, after rearrangement of (15), that Thus, u (k+1) has correct mass, at least as long as αa = 1. Hence, the Jacobi iteration is conservative in this special case.

Gauss-Seidel Iteration
The Gauss-Seidel iteration is defined by M = α(I − α(D + L)) −1 U and N −1 = (I − α(D + L)) −1 . Inserting this into (14), multiplying by I − α(D + L), then adding and subtracting Following the same procedure as previously, we multiply the ith element by | i | and sum over i to obtain By assumption, the final term vanishes. However, the second term on the right-hand side is in general non-zero, and introduces a conservation error given by

Proposition 4 The Gauss-Seidel iteration is not globally conservative in general.
However, note that for problems where A is lower triangular the Gauss-Seidel iteration is indeed conservative. In fact, for such problems the method is by definition exact.

Krylov Subspace Methods
Given an initial guess u (0) for the solution of the linear system (13), the kth iteration of a Krylov subspace method belongs to the space We may express the kth iteration u (k) as for some coefficients c l , l = 0, . . . , k − 1 that depend on, and are chosen by the method. It follows that The final term vanishes by assumption. Note that the mass of r (0) is zero if the mass of u (0) equals that of b. Thus, the second term on the right-hand side also vanishes. We thus conclude: Proposition 5 All Krylov subspace methods are globally conservative when the initial guess has the same mass as the right-hand side.
Typically a preconditioner P −1 will pre-or post-multiply A in Krylov subspace methods. We will not delve into this subject here but merely observe that if y and P −1 y have the same mass for any choice of y, then the above analysis applies without modification.

Multigrid Methods
Multigrid methods combine two methods, namely a smoother and a coarse grid correction (CGC). The idea is to separate the residual into high and low frequency components. As a smoother, an iterative method is applied, designed to effectively damp the high frequency components.
The role of the CGC is to remove the low frequency components of the residual. The residual is mapped to a coarser grid using a restriction operator, R. Smoothing is then applied and a prolongation operator P is used to reconstruct the residual on the fine grid. Finally, the fine grid solution is updated by adding the correction to the previous iterate. The procedure can be applied to a hirearchy of grids, thereby resulting in a multigrid method. If we consider only two grid levels and assume that the system is solved exactly on the coarser one, then the coarse grid correction for the linear problem (13) has the form Here, the notation (·) and (·) −1 denote matrices operating on the fine and the coarse grid respectively. The easiest way to make a multigrid method globally conservative is to choose its components to be globally conservative. The smoother can be any of the globally conservative iterations discussed so far. By Lemma 2, the inverted matrices in (16) are mass conserving. Thus, it remains to choose globally conservative restriction and prolongation operators R and P. This is achieved using agglomeration. The restriction is performed by agglomerating a number of neighboring cells by constructing a volume weighted average of fine grid values. Conversely, the prolongation is done by injecting coarse grid values at multiple fine grid points. In fact, this is the standard choice of multigrid method in CFD since it results in faster convergence than other alternatives. For details, see [3].

Pseudo-Time Iterations
Pseudo-time iterations are obtained by adding a pseudo-time derivative term to the algebraic equations (11) (or (13)), Here, u 0 is initial data that must be provided, e.g. u 0 = u n . The idea is that the system eventually should reach a steady state as τ → ∞ where the pseudo-time derivative vanishes [14], resuting in a solution to the original nonlinear system (11). Any time integration method can in principle be applied to the initial value problem (17). Explicit RK methods are globally conservative by construction. Implicit methods result in new systems of equations whose solutions are once again approximated using iterative methods [5]. Whether the resulting approximation is conservative depends on the choice of method. Any one of the conservative methods discussed so far can be applied in principle. We will return to pseudo-time iterations in Sect. 4, where more details are provided.

Numerical Validation
To validate the results of the preceding section we consider two simple experiments. The first is the linear advection equation with periodic boundary conditions, We discretize using the finite volume scheme (4) with a central numerical fluxf This results in a linear system of the form (13) with α = − t/ x and A = Tridiag − 1 2 , 0, 1 2 . Here, Tridiag(·) refers to a tridiagonal matrix with periodic wraparound. Note that the diagonal of A is constant. We therefore expect the Jacobi method to be conservative for this particular discretization.
The second experiment is Burgers' equation with periodic boundary conditions, Again, we use the finite volume scheme (4), however with the upwind fluxf i+ 1 The result is a nonlinear system in the form (11) with α = − t and The Jacobian of this discretization is lower triangular except for a single element in the top right corner. Thus, we expect the Gauss-Seidel method to be nearly conservative if the grid is sufficiently fine.
Burgers' Eq. (19) is solved using Newton's method with initial guess u n . The resulting linear systems are either solved exactly or using the following iterative methods with 0 as initial guess: (H) Pseudo-time iterations using Heun's method with τ = 0.5.
As prolongation and restriction operators for the CGC, agglomeration gives The same methods are used directly on the advection problem (18) with u n as initial guess. First we consider very coarse discretizations with x = t = 0.5 for both problems. A single iteration is used with each method and the mass error and residual is computed after one time step. The results are shown in Table 1. Here, any error smaller than 10 −15 is denoted as zero. The Richardson iteration, GMRES, CGC and the pseudo-time iterations are conservative for both problems. Similarly, the exact Newton method (Exact) is conservative. As expected, the Jacobi method is conservative for the advection problem but not for Burgers' equation. Gauss-Seidel is not conservative for any of the two problems.
Next, we repeat the experiment with well resolved discretizations. For the advection problem we choose t = x = 0.006 and for Burgers' equation t = x = 0.003. In the latter case, two Newton iterations are used per time step. For each of the other methods, 5 iterations are used except for CGC where a single iteration is performed.
The total mass is computed in each time step and compared with the mass of the initial data. Figure 1 shows the results. The main difference to the previous test is that Gauss-Seidel has a non-detectible mass error for Burger's equation. This is in line with our expectations since the Jacobian only has a single nonzero element above the main diagonal.

Local Conservation of Pseudo-Time Iterations
Having covered global conservation in the previous section, we now switch focus to local conservation. The Jacobi and Gauss-Seidel methods cannot be locally conservative since they are not globally conservative in general. For the other methods considered so far, the question remains open.
In what follows, we restrict our attention to the Cauchy problem for (3). We focus on pseudo-time iterations, show that these methods are locally conservative and prove an extension of the Lax-Wendroff theorem. Noting that the Richardson iteration is equivalent to pseudo-time iterations with explicit Euler, this method is also covered.

Pseudo-Time Iterations in Locally Conservative Form
We once again consider the finite volume method (4) where the implicit Euler method is used as time discretization. In order to apply pseudo-time iterations to this scheme, we introduce a pseudo-time derivative, where the nonlinear function g i is given by Several different methods are available for iterating in pseudo-time [5,28]. Herein, we use an explicit s-stage Runge-Kutta method (ERK). Let (A, b, c) denote the coefficient matrix and vectors of the ERK method. We denote the kth pseudo-time iterate by u where the stage vectors U (k) As previously, p and q determine the bandwidth of the finite volume stencil.
In the remainder we always use a fixed number N of iterations. A step in physical time is taken by setting u n+1 i is chosen as initial iterate. Recall that the stability function φ(z) of an RK method (A, b, c) is given by (see e.g. where I is the s × s identity matrix. The stability region of the RK method is defined as the subset of the complex plane for which |φ(z)| < 1. The proof of the following lemma is rather lengthy and is therefore deferred to Appendix A.
where the fluxĥ
Proof Firstly, the locally conservative form (27) follows from setting u n+1 (25). Secondly we recall that the numerical fluxf i+ 1 2 depends on p + q + 1 parameters, e.g. (w i− p , . . . , w i+q ). Inspecting the fluxĥ (N ) i+ 1 2 we see that it is additionally dependent on u n i so that we may writeĥ To establish the consistency we must therefore show thatĥ N (u, . . . , u; u) = c f (u). To do this, we first note from (21) that g i (u, . . . , u; u) = 0 due to the consistency off i+ 1 2 . Using the fact that A is lower triangular for any ERK method, it follows from (22) and (23) where the consistency off has been used in the final equality. Inserting this into the flux functionĥ (N ) i+ 1 2 in (26) and using (24) giveŝ where we have utilized the fact that the sum in the third equality is telescoping. Finally, note thatĥ

Suppose that the numerical fluxf in (4) is consistent with f and that Assumption 1 is satisfied. Then, u(x, t) is a weak solution of the conservation law
Proof We follow the proof of the Lax-Wendroff theorem given in [17,Chapter 12] with changes and details added where necessary. Throughout the proof we let U (x, t) denote the piecewise constant function that takes the solution value u n i in (x i , x i+1 ] × (t n−1 , t n ] on the th grid. Similarly, for j = 1, . . . , s we let U (k) j (x, t) be the piecewise constant function that takes the value U (k) The discretization (27) can equivalently be expressed as Let ϕ ∈ C 1 0 be a compactly supported test function. Multiply (30) by ϕ(x i , t n ) and sum over all n and i. Due to the compact support of ϕ, these sums can be extended arbitrarily beyond the bounds of x and t, hence we obtain At this point we will use summation by parts on the sums in (31). Recall that for two sequences (a i ) and (b i ), the summation by parts formula can be expressed as Applied to the n-sum in the first term and to the i-sum in the second term in (31), this results in Here we have used the fact that ϕ has compact support in order to eliminate all boundary terms except the one at n = 0. We now let → ∞ and investigate the convergence of the terms arising in (32). The first and third terms are identical to those in the proof of the Lax-Wendroff theorem [17,Chapter 12]. Thus, it can immediately be concluded that the first term converges to To this expression, add and subtract The bracketed part of this expression evaluates to c(μ 0 , . . . , μ N −1 ), as in the proof of Theorem 2. Thus, the first half of this expression equals Apart from the factor c(μ 0 , . . . , μ N −1 ), which is independent of , this term appears identically in the proof of the Lax-Wendroff theorem [17,Chapter 12] and converges to It remains to show that vanishes in the limit → ∞. Since U and U U (x, t)) .
Here,f (k) i+ 1 2 (x, t) is placeholder notation for the vector whose jth element iŝ To establish that this term indeed vanishes, it suffices to show that tends to zero for almost every x. By the Cauchy-Schwarz and triangle inequalities, (34) is bounded by where the Euclidean norm in R s is used. The only term in (35) that depends on is f (k) , t)) and it therefore suffices to show that this term vanishes for almost every x.
Sincef i+ 1 2 is Lipschitz continuous in each argument there is a constant L such that the final expression is bounded by The second of these terms appear in the proof of the Lax-Wendroff theorem [17,Chapter 12] and vanishes in the limit for almost every x due to the bounded total variation of U . Further, from (23) and (22) and the fact that u A few remarks about Theorem 3 are in place: First, we demand from a useful iterative method that it converges to the correct solution as N → ∞. Thus, the pseudo-time steps should be chosen in a way that ensures that c → 1, or equivalently, The simplest way to do this is to choose each pseudo-time step so that |φ(−μ l )| < 1, i.e. to stay within the stability region of the RK method. Secondly, observe that if φ(−μ l ) = 0 for any l, then c = 1 irrespective of how many further iterations that are carried out. For some RK methods such a root exists while for others it does not. For instance, the explicit Euler method has stability function φ(−μ) = 1 − μ, hence μ = 1 is a root of φ. On the other hand, Heun's method has stability function φ(−μ) = 1 − μ + μ 2 /2 > 0 for all μ ∈ R and therefore does not have any real roots. A strategy is thus to choose a RK method with a root, begin the pseudo-time iterations with a step that corresponds to this root, then resort to a conventional method for choosing the remaining pseudo-time steps. The initial iteration will not change the limit as N → ∞ if the remaining pseudo-time steps are selected within the stability region.

Numerical Results
Next we validate Theorem 3 by numerically solving a series of linear and nonlinear conservation laws.

Linear Advection
The first setting is the linear advection Eq. (18). The computational domain is x ∈ (−1, 1], t ∈ (0, 0.25] and periodic boundary conditions are used. The upwind fluxf i+ 1 2 = u n+1 i is used for the spatial discretization. The resulting finite volume method becomes Throughout the experiments the temporal and spatial increments are chosen to be equal: t = x. We validate Theorem 3 by studying the convergence of the numerical scheme to the solution of the original advection problem (18) as well as to the solution of the modified version u t + cu x = 0. As pseudo-time iteration we use the explicit Euler method, Heun's method and the third order strong stability preserving RK method SSPRK3 [26]. Theorem 3 predicts that these methods respectively will modify the propagation speed by the factor Here, we fix N = 4 and consider two different sequences of pseudo-time steps, one with constant and one with variable step sizes. The first is given by μ l = 1/20 and the second by μ l = 2 −l for l = 0, . . . , 3. The corresponding modification constants c(μ 0 , . . . , μ 3 ) are shown in Table 2. Note that with the second sequence, c = 1 for the explicit Euler method. This is due to the fact that μ 0 is a root of the stability polynomial in this case. The advection problem is solved on a sequence of grids with grid spacing x = 2/(40 × 2 j ) for j = 1, . . . , 12. The L 2 error is calculated with respect to the exact solution of the original and modified equations. The results are shown in Fig. 2. When constant pseudo-time steps are used (Fig. 2a), the lines overlap. None of the methods converge to the solution of the original problem since the numerical solution is moving with the incorrect speed. Instead, all   (Fig. 2b). Solutions of the advection Eq. (18) propagate from left to right at unit speed. In this setting it follows from Theorem 3 that the numerical solution will propagate with the speed c(μ 0 , . . . , μ N −1 ) in the limit x, t → 0. Thus, if c = 1 in each physical time step, the numerical solution will drift out of phase. However, it should be noted that Theorem 3 is an asymptotic result and does not necessarily imply that c(μ 0 , . . . , μ N −1 ) alone captures the propagation speed error on coarse grids. In fact, we may generally expect a contribution to the speed from dispersion errors built into the discretization; see e.g. [19-22, 29, 30] for details and remedies.
In the next experiment we verify that c depends on N as predicted by Theorem 3 by measuring the propagation speedc of the numerical solution while varying N . To measurē c we track the x-coordinate of the maximum of the propagating pulse through time. To get accurate measurements we extend the computational domain to x ∈ (−1/5, 1/5], t ∈ (0, 6].
For this experiment we fix x = 0.003 and set μ l = μ to be fixed for each l = 0, . . . , N − 1. The measured propagation speed error 1 −c for different μ are shown in Fig. 3a for Heun's method and in Fig. 3b for SSPRK3. The solid lines show the theoretically predicted error, The measured and theoretical speed errors agree well for μ = 0.05 and μ = 0.2. For μ = 0.5 a discrepancy between theory and measurement is seen for small errors. This suggests that the dispersion error intrinsic to the finite volume scheme is starting to dominate the propagation speed error.

Burgers' Equation
Next we consider a triangular shock wave propagating under the 1D Burgers' equation with periodic boundary conditions: The exact solution to this problem is given by Modifying the conservation law to u t + (cu 2 /2) x = 0 while retaining the same initial condition modifies the exact solution to For this problem we therefore expect that pseudo-time iterations to affect both the speed and the amplitude of the shock front.
As for the advection problem, we investigate the L 2 convergence of the numerical scheme to the original and modified conservation laws. We once again use an upwind numerical flux and implicit Euler in time; We run the simulation to time t = 0.1 with t = x using the sequence of grids x = 1/(25 × 2 j ) for j = 1, . . . , 15. The explicit Euler method is used as pseudo-time iterator with N = 12 and μ l = 1/4 for l = 0, . . . , 11. The corresponding modification constant is c ≈ 0.9683. Figure 4a shows that the numerical solution converges to the solution of the modified conservation law as expected. Fixing x = 0.004 and setting μ l = μ = 1/4, we run the simulation to time t = 1 with different choices of N . Figure 5 shows the numerical solutions together with the initial data (dashed line) and the exact solution (dotted line). The locations of the tips of the shock waves as predicted by Theorem 3 are indicated as crosses in the figure. There is good agreement between theory and experiment, although the shocks appear slightly smeared due to the built-in dissipation in the numerical scheme.  Next, we repeat the experiments but change the initial data to a step function, Since this data is not periodic, we impose the boundary condition u(0, t) = 1. Theorem 3 does not treat boundary conditions and it is therefore interesting to see if it still provides useful predictions of the behaviour of the numerical solution in this setting. The exact solution is the initial step function travelling to the right. The shock speed is given by the Rankine-Hugoniot condition as where u l,r denote left and right states of the discontinuity respectively. The shock speed of the modified conservation law is instead c/2. We use the exact same grids and pseudo-time iterations as for the triangular shock and measure the L 2 error with respect to the exact and modified conservation laws. The results are shown in Fig. 4b. Convergence is once again seen towards the modified equation, thus verifying that Theorem 3 predicts the propagation speed error correctly despite the added boundary condition.
Next we extend the time domain to t ∈ (0, 1], fix x = t = 4 τ = 1/100 and vary the number of iterations, N . The computed solutions using N = 1, N = 3 and N = 12 are shown in Fig. 6 together with the predicted shock locations (dashed lines). Once again, there is good agreement between prediction and experiment, although numerical dissipation smears the shock fronts somewhat.
As mentioned previously, the shock speed error can be eliminated entirely for the explicit Euler method by noting that μ = 1 is a root of the stability polynomial 1 − μ. To highlight the effect of this, we introduce two strategies for choosing the pseudo-time steps: The first strategy is the same that we have used in the experiments so far. The second one ensures that c = 1 by taking a large initial pseudo-time step. Note that both strategies correspond to integration in pseudo-time to the same point; τ = 3 t.
The relative residuals of the pseudo-time iterates in the first physical time step are shown in Fig. 7 for the two strategies. The large initial pseudo-time step in Strategy 2 results in a considerably greater residual reduction than the corresponding iteration using Strategy 1. Interestingly, subsequent iterations yield faster convergence of the residual using Strategy 2 as seen by the steeper gradient. This suggests that the incorrect shock speed makes the dominant and slowest converging contribution to the residual for this problem when using Strategy 1. We also conclude that the point to which we march in pseudo-time, here τ = 3 t, have less of an impact on the convergence than the choice of pseudo-time steps used to reach this point.

The Euler Equations
As a second nonlinear problem we consider the 2D compressible Euler equations, ⎡ where γ = 1.4 and e = E − (u 2 + v 2 )/2 is the internal energy density. The domain is taken to be periodic in both spatial coordinates. The setting is the isentropic vortex problem [25] with initial conditions where r = 1 − x 2 − y 2 . Here, = 5 is the circulation and M ∞ = 0.5 is the Mach number. As the solution evolves in time, the initial vortex propagates in the horizontal direction with unit speed. As in previous experiments, we use implicit Euler in time, yielding a finite volume scheme of the form Along the x-coordinate we use a fourth order centered flux Both strategies integrate in pseudo-time to the point τ = 1.8 t in each time step. However, Theorem 3 predicts that Strategy 1 will give a speed modification c(μ 0 , . . . , μ 8 ) ≈ 0.866. On the other hand, Strategy 2 will give c(μ 0 , . . . , μ 4 ) = 1, i.e. the correct propagation speed, due to the large initial pseudo-time step. Figure 8 shows the numerical solutions at time T = 10 for the case where x = y = 0.2, t = 0.05. The predicted and observed vortex locations agree very well. Figure 9a shows the convergence of the numerical solutions, measured at time T = 1. Convergence to the correct solution is observed for Strategy 2 (S2) but not for Strategy 1 (S1). This is expected due to the incorrect location of the vortex in the latter case. However, convergence to the solution of the modified conservation law is seen (S1 mod). Thus, Theorem 3 accurately predicts the behavior of the numerical solution despite the violated assumptions.
In practical applications, the pseudo-time iterations are terminated when the residual has decreased beneath some tolerance. It is interesting to see in what way the convergence is affected by the choice of strategy. Returning to the setting in Fig. 8, the residual in each pseudo-time iteration for all 200 physical time steps are shown for the two strategies in Fig.  9b. The 200 lines overlap nearly perfectly, suggesting that the residual behaves similarly in each physical time step. The initial iteration in Strategy 2 evidently has a large impact on the reduction of the residual that supercedes those of the other iterations put together. The remaining iterations appear to reduce the residual by comparable amounts for the two strategies, as seen by the similar gradients. In contrast to the shock problem considered previously, this suggests that the convergence rate of the residual is dictated by other factors than the propagation speed for this particular problem. Nonetheless, a correct propagation speed is visibly very beneficial, here with a drop in relative residual of more than an order of magnitude.

Experiments with GMRES and Agglomeration CGC
The analysis and experiments with pseudo-time iterations reveals that iterative methods can cause (predictable) alterations to the physics described by a conservation law, even if conservation is preserved. Although we have not presented a similar analysis for multigrid and Krylov subspace methods, it is interesting to experimentally investigate whether they display similar modifications to the propagation speed. To this end, we return to the linear advection problem in Sect. 4.2.1. Throughout we set t = 4 x and solve the linear system in each time step either with an agglomeration based coarse grid correction or with GMRES. The initial guess is taken to be the solution in the previous time step. The prolongation and restriction operators used are given in (20). Figure 10a shows the L 2 error with respect to the exact solution. The solutions to the linear systems have been approximated either with one coarse grid correction (CGC) or with a fixed number of GMRES iterations per time step, N . The errors appear to converge both for the coarse grid correction and for GMRES with N = 10. However, with N = 1 and N = 3, convergence is lost as the grid is refined. Figure 10b shows the numerical solutions computed with N ∈ {1, 3, 10} in the subdomain x ∈ [0, 1] at time t = 0.25, here with x = 1/5120. The exact and numerical solution computed with N = 10 overlap, with the pulse centred at x = 0.5. However, the solutions computed with N = 1 and N = 3 are seen at incorrect locations, suggesting reduced propagation speeds.
The result suggest that, for this problem, agglomeration coarse grid corrections retain consistency with the conservation law. GMRES on the other hand, behaves similarly to pseudo-time iterations in the sense that it slows down the speed of the solution. However, this problem appears to be solved automatically after some finite number of iterations.

Summary and Conclusions
In this paper we have studied conservation properties of a selection of iterative methods applied to 1D scalar conservation laws. The fact that conservation is a design principle behind many numerical schemes motivates such a study. We have established that Newton's method, the Richardson iteration, Krylov subspace methods, coarse grid corrections using agglomeration, as well as ERK pseudo-time iterations preserve the global conservation of a given scheme if the initial guess has correct mass. However, the Jacobi and Gauss-Seidel methods do not, unless the linear system being solved possesses particular properties. The stronger requirement of local conservation has been investigated for ERK pseudotime iterations. We have shown that local conservation is preserved for finite volume schemes that employ the implicit Euler method in time. However, the resulting modified numerical flux may be inconsistent with the governing conservation law. An extension of the Lax-Wendroff theorem shows that this inconsistency leads to convergence to a weak solution of a conservation law modified by a particular constant. We give an exact expression for the modification constant, which depends only on the stability function of the ERK method and the pseudo-time steps. Depending on the problem solved, this modification can alter both the propagation speed and amplitude of the numerical solution if the constant differs from unity. We present a strategy for ensuring that the constant equals one and show numerically that the strategy results in faster convergence.
Numerical tests corroborate the theoretical findings and suggest that the results hold even for systems of conservation laws in multiple dimensions. Further, experiments with GMRES suggest that they suffer from a similar modification as pseudo-time iterations, but that this is automatically accounted for after a finite number of iterations. Similar experiments with coarse grid corrections based on agglomeration indicate no presence of such modifications.
Author Contributions Philipp Birken and Viktor Linders both contributed to the conceptualization of this article and to the mathematical analysis. Viktor Linders was the main contributor to the numerical experiments and production of the graphics. Philipp Birken and Viktor Linders jointly contributed in writing the article.
Funding Open access funding provided by Lund University.

Conflict of interest
The authors declare that they have no conflict of interest.

Code availability N/A
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Proof of Lemma 3
The purpose of this appendix is to provide a detailed proof of Lemma 3.
Proof Note first that (22)-(23) can be equivalently written in the form where 1 = (1, . . . , 1) ∈ R s . Here, we use the notation g i (U k ) = g i (U holds. We will show that it also holds for N + 1. From (21), (25) and (41) it follows that Note that (I + μ N A) −1 exists since A is lower triangular.
Evaluating g i U (N ) using (21) gives In the last equality we have used the fact that Inserting the above expression for g i U (N ) into u (N +1) i as given in (41) leads to Rearranging, using the induction hypothesis and the expression (26) forĥ (N ) results in