Numerical simulation of conservation laws with moving grid nodes

In the present article we describe a few simple and efficient finite volume type schemes on moving grids in one spatial dimension. The underlying finite volume scheme is conservative and it is accurate up to the second order in space. The main novelty consists in the motion of the grid. This new dynamic aspect can be used to resolve better the areas with high solution gradients or any other special features. No interpolation procedure is employed, thus an unnecessary solution smearing is avoided. Thus, our method enjoys excellent conservation properties. The resulting grid is completely redistributed according the choice of the so-called monitor function. Several more or less universal choices of the monitor function are provided. Finally, the performance of the proposed algorithm is illustrated on several examples stemming from the simple linear advection to the simulation of complex shallow water waves.


Introduction
The numerical simulation of conservation laws is one of the most dynamic and central parts of modern numerical analysis. The systems of conservation laws appear in many fields ranging from traffic modeling [46] to shallow [23] and compressible fluid mechanics [63]. The finite volume method was proposed in a pioneering work of S. Godunov (1959) [28,29] and developed later by Ph. Roe [49] and many other researchers (see [26] for the current state of the art).
Nowadays the complexity of the problems which arise in practice is such that we have to think about the optimal usage of available computational resources. In this way the community came with the idea of developing numerical methods on adaptive grids in order to have higher resolution only where it is needed [5]. Historically, hp-adaptive methods were first developed for the Finite Element-type (FEM) discretizations [6] (including the newest discontinuous Galerkin discretizations as well [34]). The main advantage of FEM is that this method is quite flexible in respect of local approximation spaces, thus allowing for relatively easy p-adaptivity. FEM was applied to hyperbolic problems under the guise of relaxation schemes [3]. Later h-adaptive methods have been adapted for finite volume methods as well [10]. However, the widely used approach nowadays consists in performing local grid refinement (or coarsening) with locally nested grids [9]. For instance, it has been successfully applied to tsunami propagation and run-up problems [27] as well as to more complex two-phase incompressible flows [47].
Adaptive mesh refinement techniques applied to the spatial variable after the computation of an approximate solution of the model equations at each time step usually consist of the following steps: • generation of a new spatial mesh according to prescribed adaptivity criteria, • reconstruction or interpolation of the numerical solution on the new mesh, • integration in time of the numerical solution on the new mesh. The main disadvantages of this methodology is that the reconstruction of the numerical solution on the new mesh increases the complexity of the method and also reduces the accuracy of the method by adding more dissipation or dispersion to the solution. Moreover, the variable number of nodes represents some difficulties while implementing the algorithm on a computer. The underlying data structures have to be flexible enough to insert or to suppress some elements.
The approach we propose in this study is different and does not require interpolation or reconstruction of the solution to the new mesh (a conservative interpolation is employed in [58]). Although the proposed method is closer to the redistribution method proposed in [2,58], is simpler and more elegant at the level of implementation details while it is based on the transformation of the model equations to some new equations where the solutions are independent of the new mesh. Namely, the method uses the same number of discretization points of the spatial mesh during the simulation and the adaptivity is achieved by moving the grid nodes to the places indicated by the so-called monitor function. The adapted grid is obtained as a solution of an elliptic (or parabolic) problem, which ensures the smoothness of the obtained grid. The main idea behind the equidistribution principle is to distribute the nodes over the computational domain so that the measure of a cell times the value of the monitoring function on it be approximatively a constant [19,20]. The heart of the matter in the equidistribution method is the construction of a suitable monitoring function. This method was described independently in the seminal paper [53], and more than twenty years later in [35]. Numerous subsequent developments were published in recent books [37,44]. The review of earlier works on the adaptive grid generation can be found in [30,43]. The 'Holy Grail' in choosing the monitoring function is to achieve ideally the reduction of error in orders of magnitude when the discretization step h → 0 , as illustrated in [38].
For simplicity, we focus on one-dimensional cases only and even in this simple case some open problems are outlined throughout the manuscript. According to our knowledge the pioneers in 1D case were [57] (for nonlinear shallow water equations) and [1] (for the gas dynamics). The generation of 2D curvilinear adaptive grids is a relatively classical topic [51]. A rigorous definition of a curvilinear grid was given in [54]. The generalization of the equidistribution method to two spatial dimensions was performed for the first time in [18]. There are a few studies which report recent modern implementations on moving redistributed two-dimensional grids (see e.g. [4,8,14,58]) with generalized monitor functions formalized in [37]. Some approaches for the inter-comparison of generated grids are described in [48].
The discretization of the model equations is based on explicit, fully discrete predictorcorrector Finite Volume schemes that are accurate up to the second order in space. Predictorcorrector schemes have been introduced by MacCormack in [45] and can be considered as a natural generalization of splitting schemes in higher spatial dimensions.
We have to mention that non-uniform grids can be used also to preserve some Lie symmetry group of the equation at the discrete level (see [15,Section 3.2.2] for a brief survey and Burgers-Hopf equation example). Originally this approach was proposed by Dorodnitsyn [21] and developed later in [36]. In our study the mesh motion is directed solely by the equidistribution principle as in [58] without taking into account the symmetry considerations for the moment.
The present article is organized as follows. In Section 2 we discuss in all details the implementation of the algorithm for the linear advection equation. Then, we generalize it to the nonlinear case in Section 3 and to systems of conservation laws in Section 4. Some open problems are outlined in Section 5. Finally, the main conclusions and perspectives of our study are given in Section 6.

Linear scalar equation
In order to present the moving grid method we shall start with the simplest scalar linear advection equation in this section and, then, we shall increase gradually the complexity by adding first the nonlinearities in the following section and moving to the systems by the end of our manuscript.
We pay special attention to the linear case for the following reasons: • The main properties of the scheme are the most transparent in the linear case • The generalization to the nonlinear and vectorial cases will be easier once the linear case is fully understood. So, it will allow us to go faster in the subsequent sections • An exhaustive error analysis is possible in the linear (and presumably only in the linear) case In order to construct an efficient finite volume scheme on a moving mesh, we have to choose first a robust and an accurate scheme on a fixed grid, which will be generalized later to incorporate the motion of mesh points. Such a scheme retained for our study is described in the next Section.

A predictor-corrector scheme on a fixed uniform mesh
Consider the following linear advection equation with a constant propagation velocity The subscripts in this study denote the partial derivatives, i.e. u t = ∂u ∂t , u x = ∂u ∂x , etc. We introduce a uniform discretization of the real space line R with nodes x j = jh, where h > 0, j ∈ Z is the spatial discretization step, and for simplicity we do not pay attention to the boundary conditions here. The interval C j = [x j , x j+1 ] will be referred as the cell C j . The time step is denoted by τ and the discrete solution is computed at t n = nτ , n ∈ Z + . The last equation (2.1) is discretized in space and time using an explicit fully discrete predictor-corrector scheme where τ is the time step. The intermediate quantities u * j±1 2 are evaluated in the middle of cells at x j±1 2 = x j + x j±1 2 = x j ± 1 2 h and at time instances t = t n + τ * j+1 2 . Moreover, the following notations have been introduced in (2.2): We underline the fact that the parameter θ has not been fixed for the moment and it may vary from one cell and one time layer to another one. The predictor-corrector scheme (2.2) can be recast as a one-step scheme by combining two equations together: The last fully discrete scheme is a canonical form of all two steps explicit schemes for equation (2.1), since many well-known schemes can be obtained from (2.3) by choosing carefully the parameter θ. Indeed, Lax-Wendroff scheme [41] (θ = 0): First order upwind [17] (θ = 1 C − 1 > 0): The number C ∶= a τ h < 1 is the famous Courant-Friedrichs-Lewy number [16], which controls the stability of the upwind scheme Lax-Friedrichs scheme [11] However, the schemes listed above are not suitable for accurate simulations, since the first order upwind and Lax-Friedrichs schemes are too diffusive in practice and Lax-Wendroff scheme is not monotonicity preserving. A good solution consists in choosing the parameter θ adaptively. In a previous work [52] the following strategy was proposed: x,j+1 2−s and u n x,j+1 2 ⋅ u n x,j+1 2−s ⩾ 0, θ 0 , u n x,j+1 2 ⋅ u n x,j+1 2−s < 0, . The first case corresponds to the classical Lax-Wendroff scheme, the third case is the classical first order upwind and the second case is detailed in the following Remark 1. However, the proposed combination of schemes produces a robust TVD-type scheme, which ensures a very good trade-off between the solution accuracy and monotonicity properties. Its performance on fixed and mobile meshes will be illustrated in numerical sections below.
Remark 1. We note that the canonical form (2.3) contains not only three points schemes on a symmetric stencil. For instance, if we assume for simplicity a > 0, and set then we shall obtain the following second order upwind scheme on an asymmetric stencil: Of course, the numerical properties such as the accuracy, stability and monotonicity depend crucially on the choice of the parameter θ and have to be studied carefully in each particular situation.

Predictor-corrector schemes on a moving grid
In this section we shall restrict our attention to bounded domains only since almost all numerical simulations in 1D are made on finite intervals. To be more specific we assume that the computational domain is I = [0, ℓ]. Consider also the reference domain Q = [0, 1] and a smooth bijective time-dependent mapping x(⋅, t) ∶ Q × R + ↦ I, which satisfies the following boundary conditions For our modest numerical purposes we do not even need to have the complete information about this mapping. Consider a uniform mesh of the reference domain Q h = {q j = jh} N j=0 , with a constant step h = 1 N . Strictly speaking we need to know only the image of nodes q j under the map x(q, t), since they constitute the nodes of the moving mesh: x(q j , t n ) = x n j .

The linear advection equation (2.1) can be rewritten on domain
whereā ∶= a − x t is the advection speed relative to the mesh nodes velocity and J ∶= x q > 0 is the Jacobian of the transformation x(q, t). Equation (2.5) does possess a conservative form as well Now we can write the predictor-corrector scheme for the linear advection equation on a moving mesh: Note that the first step is based on equation (2.5) and the second one on the conservative form (2.6). Here we introduced some new notations: As before, one can recast the scheme (2.7), (2.8) in an equivalent one-stage form: A stability study using the method of frozen coefficients [50] leads to the following restrictions on the scheme parameters where C n j+1 2 is the local CFL number [16] defined as In order to ensure the monotonicity of the numerical solution, it was shown in [52] that it is sufficient to use the following choice of the scheme parameter θ: where s = sign(ā n j+1 2 ) and .
Using (2.10), it can be shown that stability requirements (2.9) are fulfilled under the following single restriction on local CFL numbers: It is this value, which will be reported in Tables below. There is an additional relation, which is called the geometric conservation law. In 1D case it looks trivial since it translates a trivial fact that x qt ≡ x tq . However, it was demonstrated in [59] that it is absolutely crucial to respect the geometric conservation law at the discrete level as well in order to have a consistent fully discrete scheme which preserves all constant solutions exactly, for example. It can be shown that the scheme proposed above satisfies the discrete counterpart of the geometric conservation law in integer and half-integer nodes:

Construction and motion of the grid
In the previous sections we explained in some detail the predictor-corrector schemes on fixed and moving meshes. We saw that the mesh motion is parametrized by a bijective mapping x(q, t) ∶ Q×R + ↦ I. Now we explain how this map is effectively constructed, since so far it was quite general. Moreover, the predictor-corrector schemes presented above are valid for other constructions of the mesh motion.
The first step consists in choosing a positive valued function w(x, t) ∶ I ×R + ↦ R + , the socalled monitor function, which controls the distribution of the nodes. Without entering into the details for the moment, we suggest as one choice for a monitor function the following: where α > 0 is a real parameter. The monitor function (2.11) is a particular case of the class of monitor functions considered by Beckett, Mackenzie et al. [8]. Other choices of the monitor function are discussed in [12,14,55]. For instance, to give an example, another popular choice of the monitor function is In general, the choice of the monitor function and its parameters might be crucial for the accuracy of the scheme. The adaptive grid construction for a simple linear singularly perturbed Sturm-Liouville equation is deeply analysed in [38]. This question seems to be rather open for unsteady nonlinear problems. The current state-of-the-art in this field is described in [37]. It is noted that the monitor function varies in both space and time and takes large (positive) values in areas where the solution has important gradients. The monitor function w(x, t) specified in (2.11) can be readily discretized The proposed algorithm is slightly different on the very first step when the initial grid is generated. It is important to obtain a high-quality initial mesh, since the numerical error committed in the beginning cannot be corrected afterwards. So, the first step is explained in detail below and the algorithm for unsteady computations is presented in the following section.

Initial grid generation
It is of capital importance to produce a grid adapted to the initial condition as well, since no error committed initially can be corrected during dynamical simulations. At t = 0 we compute the monitor function w (x, 0) of the initial condition u (x, 0) and the mapping x (q, 0) is determined as the solution to the following nonlinear elliptic problem, which degenerates to a simple second order ordinary differential equation supplemented with Dirichlet boundary conditions (2.4). This elliptic problem can be readily discretized to produce the following finite difference approximation: with boundary conditions x 0 0 = 0, x 0 N = ℓ . The last nonlinear system of equations is solved iteratively (usually with some fixed point-type iterations). The iterations are initialized with a uniform grid as the first guess. One can notice also that its solution satisfies the equidistribution principle: where the constant C h can be computed as So, C h is the discrete version of the monitor function integral. The equidistribution principle gives meaning to the monitor function: in the areas where w 0 j+1 2 takes large values, the spacing between two neighbouring nodes x 0 j+1 and x 0 j has to be inversely proportionally small.

Unsteady computations
Assume that the grid {x n j } N j=0 at t = t n is known and let us evolve it to the next time layer t n+1 . It is done by solving the following linear fully discrete problem: with the same boundary conditions x n+1 Here the real positive parameter β > 0 plays the rôle of the diffusion coefficient and it controls the smoothness of nodes trajectories. It can be easily seen that the scheme (2.13) is an implicit discretization of this nonlinear parabolic equation: However, it has to be stressed out that the fully discrete problem is linear since the monitor function w n j±1 2 is taken on the previous time layer. Thus, the overall complexity of the grid motion algorithm at every time step is equal to that of the tridiagonal matrix algorithm [25].

Smoothing step
In some problems where the solution contains a lot of oscillations (i.e. many local extrema), the monitor function w(x, t) has to be filtered out to ensure the smoothness of the mesh motion. Moreover, a smoothing step allows to enlarge the set of acceptable values of the parameter α in the definition (2.11) of the monitor function. Our numerical experiments show that the following implicit procedure (inspired by implicit discretizations of the heat equation) produces very robust results: where σ > 0 is an ad-hoc smoothing parameter to be set empirically. System (2.14) is completed with boundary conditions The smoothed discrete monitor function {w j } N j=0 is obtained by solving the linear system (2.14) with any favourite method. In the present study we used the simplest tridiagonal matrix algorithm [25]. In contrast, the authors of [58] used an explicit smoothing method.
We stress out that one should apply a smoothing operator to the monitor function and not to the numerical solution.

On the choice of the monitor function
In general, when one has a problem (i.e. a system of PDEs) to solve, there is a question of the best possible monitor function choice for the grid adaptation. One has to choose also the pertinent solution component as an input for the monitor function. It is extremely difficult to give some general recommendations which would work in all cases * . However, in choosing the monitor function one has to take into account the following aspects: Physical considerations: Any a priori knowledge on the solution can be used to underline important physical aspects (e.g. the presence of shock waves or the presence of multiple wave crests, etc. ) Geometrical considerations: Computational domain geometry (e.g. the presence of angles or fractal shorelines in wave run-up problems) and topology (i.e. the presence of islands or other obstacles) Mathematical considerations: Nature of equations (differential, partial differential, integrodifferential, singularly perturbed, etc. ), but also one has to pick up the most pertinent variable to compute the monitor function on it (the free surface elevation or the vorticity field in Navier-Stokes simulations) * It is probably as difficult as the general theory of PDEs.  Numerical considerations: Particularities of the employed discretization scheme may also hint the choice of the monitor function.
Once the monitor function is chosen, one has to determine also (nearly) optimal values of parameters which dependent not only on the problem, but also on the solution. Finally, the most important guideline nowadays is the experience in solving practical problems using moving grid techniques. It cannot be replaced by any guidelines.

Numerical results
In order to illustrate the performance of this algorithm we take a relatively serious testcase which consists in simulating the propagation of a discontinuous profile: The monitor function employed in this computation is defined in (2.11). The numerical parameters used in this simulations are given in Table 1. The simulation results at time t = T f is shown in Figure 1. From the comparison of profiles obtained on fixed and moving meshes, one can see that the moving mesh reduces significantly the "width" of the numerical shock wave profile.
The second test-case is a smooth bell-shaped function advected by the dynamics of (2.1). This test is used to measure the ability of the scheme to preserve the extrema (in other words we assess the strength of the numerical diffusion on a given mesh). The exact solution is given by Thus, the initial condition is simply given by u(x, 0) = e −25(x−x 0 ) 2 . Another particularity of the present test-case is that we use the monitor function specifically designed with emphasis on the fine resolution of the local extrema (those with high positive or high negative values): w(x, t) = 1 + α u(x, t) .
(2.15) It is noted that there is no differentiation operator in the definition of the last monitor function compared to the one defined in (2.11). The parameters used in the numerical simulation are provided in Table 2. The numerical results at t = T f are shown in Figure 2. For the sake of comparison we show also the numerical solution obtained with the same scheme, but on a uniform mesh. The use of a moving grid improves dramatically the preservation of the local maximum. One cannot even distinguish the exact solution from the numerical one up to the graphic resolution.

Nonlinear scalar equation
In this section we consider a general hyperbolic scalar nonlinear conservation law of the form:  As a chrestomathic example of such equations we can mention the celebrated Burgers-Hopf equation [13]. As above, we make a change of variables to the fixed reference domain Q: or in the conservative form: where v ∶= u ○ x and J ∶= x q as before. For nonlinear equations we need an extra transport equation for the flux function: where a(u) = f u (u). The last equation transforms to the following non-conservative form on the fixed domain Q: Now all the preparations are done and we can write down the predictor-corrector scheme directly on a moving grid (notice that in the nonlinear case three stages are required): where v n j+1 2 , f n j+1 2 were defined above. The last scheme has the second order approximation in space O(h 2 ) provided that θ = O(h). The only point which remains to be specified is the computation of the numerical Jacobian a n j+1 2 . In order to ensure good properties of the scheme, one has to ensure the compatibility condition f q = a(v)v q at the discrete level f n q,j+1 2 = (a v q ) n j+1 2 for finite differences analogue of the derivatives: This requirement can be satisfied with the following choice [33]: , v n j+1 = v n j . By noticing that the corrector step contains only the combination f * − x n t u * , we can write a compact two-stages form of the predictor-corrector scheme: whereā(v, q, t) = a(v) − x t . The parameter θ n j+1 2 is computed according to (2.10) which ensures the TVD property of the scheme. As a result we obtain a generalization of Harten scheme [33] to the case of moving meshes.

Numerical results
In order to illustrate the performance of the generalized Harten scheme on moving meshes we take the classical example of the inviscid Burgers-Hopf equation (f (u) = u 2 2 ). As the first illustration let us take the following initial condition on the computational domain x ∈ [0, ℓ]: u r x ⩾ x r . All numerical parameters are given in Table 3. The exact weak solution to this problem is given by the following formula: The gradient catastrophe takes place at t = t * and x = x * , which can be computed exactly as well After t = t * we obtain a travelling shock-wave profile: For the parameters chosen in our numerical simulation (see Table 3) the shock wave remains stationary at t = t * = 5 and x = x * = 15 since u l = −u r . The numerical solution and the grid motion are shown in Figure 3. It can be seen that the solution monotonicity is preserved on both fixed and moving meshes. However, on the moving grid we obtain a solution which cannot be distinguished from the exact one up to the graphical resolution. On the fixed grid the shock wave has a finite width equal to two control volumes. As another numerical illustration, we consider the same problem, but we change the value of u r ∶= 0 instead of −1. In this case the gradient catastrophe takes place at t = t * = 10 and x = x * = 20. The shock wave profile is not steady anymore since {u} = u l +ur 2 = 1 2 ≠ 0. So, it will move in the rightward direction with a constant speed equal to 1 2 . The results of numerical simulations are shown on Figure 4. It can be seen again that the adaptive solution is excellent, while the fixed grid produces a shock wave smeared out on four cells only.

Nonlinear system of equations
The predictor-corrector scheme for nonlinear systems on moving meshes will be presented on a practically important case of the Nonlinear Shallow Water Equations (NSWE) [56].      This choice reflects also the scientific interests of the authors of the present article. Moving grid techniques have been applied to other systems of conservation laws such as Euler equations in [55], the pressureless gas dynamics [40] and Magneto-Hydrodynamics [61]. So, consider one-dimensional nonlinear shallow water equations written in a vectorial form: where u is the vector of conservative variables, f (u) is the flux function and G (u, x) contains the source terms: Here u(x, t) is the depth-averaged velocity and H(x, t) = η(x, t) + h(x) is the total water depth defined as the sum of the bathymetry function h(x) and the free surface elevation η(x, t) over the mean water level. The sketch of the physical domain with all these definitions is shown in Figure 5. The system (4.1) is considered on a finite segment I = [0, ℓ] with appropriate boundary conditions. As before, we consider a smooth bijective timedependent transformation x = x(q, t) from the reference domain Q = [0, 1] into I. The NSWE on the reference domain are written in the form [7]: where v(q, t) = u ○ x = u x(q, t), t and the source term contains a derivative with respect to q : It is noted that the bathymetry becomes time-dependent under this change of variables, even if initially the bottom was fixed. There exists also a conservative form NSWE on the reference domain: Finally, the flux vector satisfies the following equation: where A is the advective flux's Jacobian matrix, which can be easily computed: Since the system of NSWE (4.1) is hyperbolic, the Jacobian A can be decomposed into the product of left L and right R eigenvectors: where the matrices L , R and D are defined as Here λ 1,2 = v ∓c (c = √ gH being the speed of gravity waves) are eigenvalues of the Jacobian A (v). They are real and distinct provided that H(q, t) > 0. Mathematically it implies the strict hyperbolicity property for the system (4.1).

Finite volume discretization
The reference domain Q is discretized into N equal control volumes [q j , q j+1 ] of the width ∆q ≡ 1 N . The nodes x n j are obtained as images of uniform nodes under the mapping x(q, t): x n j = x(q j , t n ). In the predictor-corrector scheme on the first step we evaluate the flux vector in cell centers: The corrector step is based on the conservative form of the equation: H(x,t) where τ is the time step and the following discretizations are used The total water depth H and local depth gradient h q have to be computed in the following way: Despite the fact that H n+1 j±1 is present in the formulas above, the scheme remains explicit, since we compute first {H n+1 j } N −1 j=1 from the mass conservation and, then, this value is used in the correction of the momentum conservation equation.
The matrixD plays the rôle ofā (the wave speed relative to the grid nodes) in the scalar case:D Here λ n k, j+1 2 , k = 1, 2 are eigenvalues of the Jacobian matrix A n j+1 2 : the local gravity wave speed c n j+1 2 can be computed analytically: This choice of the Jacobian A discretization is dictated by the requirement to preserve the compatibility condition f q = A (v)v q at the discrete level as well. It is analogous to a similar algebraic condition on the Jacobian in Roe's scheme [49], but substantially different in details.
Finally, the scheme parameter θ k j+1 2 is computed according to the same strategy as above, with the only difference being that it is computed separately for each component k = 1, 2 : , s k ≡ s n k,j+1 2 ∶= sign λn k,j+1 2 . Hereλ n k,j+1 2 = (λ k − x t ) n j+1 2 are simply the diagonal elements of the matrixD. The coefficient θ k 0,j+1 2 is related to the local CFL number as Finally, we have to specify how to compute the indicators g k j+1 2 ∶= λ k (1 − C k )p k n j+1 2 , wherep k , k = 1, 2 are components of the following vector: The fully discrete scheme is stable under the following CFL-type restriction [16] on the time step τ : max k, j C n k, j+1 2 ⩽ 1.

Well-balanced property
For the scheme presented above one can show the following Theorem 1. Assume that we have a fixed, but possibly non-uniform, grid {x n j }. If on the n th time layer we have the lake at rest state, i.e.
then, the predictor-corrector scheme described above will preserve this state at the (n + 1) th time layer.
Proof. Since x t,j ≡ 0, ∀j and taking into account the definition (4.4) of the lake at rest state, we have: Consequently, the predictor step leadŝ Then, from the obvious identities h n+1 j ≡ h n j , J n+1 j ≡ J n j and the way how we compute G * j , it follows that η n+1 j ≡ 0, v n+1 j ≡ 0. Thus, the lake at rest state is preserved at the next time layer as well.
The well-balanced property is absolutely crucial to compute qualitatively correct numerical solutions to conservation (balance) laws [31]. The Theorem 1 shows that the discretization proposed above does possess this property regardless some complications coming from the grid motion.

Conservation property on moving grids
For conservation laws with source terms (the so-called balance laws), the conservation property of a numerical scheme is understood in the sense that the fully discrete solution verifies some thouroughly chosen finite difference analogues of continuous balance laws. For the NSWE we have the conservation of mass and the balance of horizontal (depthaveraged) momentum. Please, notice that on a flat bottom the balance of momentum becomes a conservation law as well. The conservation property is crucial as well to compute qualitatively correct numerical solutions [31].
On moving grids we have an additional difficulty: the computation of the speeds of mesh nodes has to be consistent with the calculation of the cell volumes. If it is done in the right way, the discrete geometric conservation law will be verified as well [59]. In 1D it can be written as We shall show below that other schemes can by put into the conservative form similar to (4.2). The pecularity here consists in the fact that the Jacobian matrix A is nonconstant, in contrast to the linear advection equation, for which any reasonable scheme is conservative (see also Remark 2 below about a nonlinear scalar equation). Moreover, the moving mesh is another complication, which makes this task rather non-trivial. For example, the scheme parameter θ k j±1 2 can be chosen as This choice corresponds to the first order upwind scheme (as it was discussed above in Section 2.1). So, the scheme (4.2) can be rewritten as One can see also that the following identities hold: By using these identities and the discrete geometric conservation law (4.5), one can derive also the following equivalent form of the scheme (4.2): The choice (4.3) of the scheme parameter θ k j±1 2 leads to a conservative scheme as well, but the underlying computations (to show it) are less elegant.
Remark 2. In order to illustrate even better the interplay between conservative and nonconservative forms, let us take again the inviscid Burgers-Hopf equation [42]: For the sake of simplicity let us discretize this equation on a uniform grid using the following non-conservative scheme: It is not difficult to see that this non-conservative scheme is equivalent to the following conservative upwind scheme: However, if one takes a different non-conservative scheme * u n+1 it is straightforwatd to see that one cannot transform it in a conservative form. Moreover, one can even show that scheme (4.6) is not converging. So, the conservation and convergence properties are intrinsically related.

Numerical results
In this section we present two test cases of complex wave propagation in shallow water equations. The first test case is a validation against an implicit analytic solution, derived below specifically for this purpose, while the second test case is related to real world applications.

Exact solution derivation
In order to validate the numerical scheme, we derive an exact solution to nonlinear shallow water equations which will serve as a reference solution. We already saw this particular solution in the literature [62], however, the derivation procedure is published for the first time to our knowledge.
First of all, we assume that the bottom is flat, i.e. h(x) = h 0 = const > 0 and we consider an initial value problem on the real line: The nonlinear shallow water equations are rewritten using Riemann invariants [60]: We seek for a particular solution in the form of an r-wave propagating in the leftward direction. In other words we assume s = s 0 = const. Assuming that at the infinity the water is at rest, we deduce that s 0 = 2c 0 = 2 √ gh 0 . Then, from formula u + 2c = s 0 we deduce that the speed u(x, t) is related to the total depth: (4.7) Consequently, taking the limit t → 0 we have a similar relation between the initial conditions: u 0 (x) = 2c 0 − 2 g(η 0 (x) + h 0 ). * The difference is in the estimation of propagation speed, cf. u n j vs. 1 2 (u n j + u n j−1 ).

If we introduce a new variable
p(x, t) ∶= 3r(x, t) + s 0 4 the governing equation for the r-invariant becomes simply the Burgers-Hopf equation: The solution can be obtained using the method of characteristics by the following implicit formula [50,63]: Consequently, the solution p(x, t) can be obtained at any space location x and time instance t by solving the following implicit equation Once, p(x, t) is found, the free surface elevation η(x, t) can be reconstructed by using this formula and the flow speed u(x, t) from (4.7).
Remark 3. Please, note that this method works only until the gradient catastrophe occurs at some time instance t = t * , which depends on the initial condition.
Remark 4. In order to guarantee the solvability of equation (4.8), it is enough to check that the following function

Validation
Using the method described above, we obtain the analytical solution for the following initial value problem: where x w is the wave crest position and a is the wave amplitude. For 0 < t < t * the profile η 0 (x) will propagate to the left with the speed c 0 preserving its wavelength λ and amplitude a. However, the gradient catastrophe will inevitably take place since the characteristics  compression wave followed by a rarefaction wave. With the parameters specified in Table 4 the gradient catastrophe takes place about t = t * ≈ 5.57. Numerical results are shown in Figure 6 right before the gradient catastrophe occurs. The adaptive grid was constructed using the following monitor function, which contains an additional term with respect to (2.11) w(x, t) = 1 + α 0 η(x, t) + α 1 η x (x, t) , (4.9) which is discretized as (4.10) Thin solid line in Figure 6 represents the exact solution computed by solving equation (4.8) as it was explained in the previous Section. The agreement between these two profiles is exemplary. On the right panel (b) one can see how the nodes follow the wave propagation and how their density increases around the steep gradient. On the other hand, the use of a fixed grid leads to a noticeable smearing of the numerical solution.

Solitary wave run-up
As the last test-case we consider a practically important problem of a solitary wave run-up on a sloping beach. The main difficulty here consists in the fact that the fluid domain may vary in time and its evolution is not known a priori. It has to be determined dynamically during the simulation. The wetting/drying algorithm was described in detail in our previous study [39]. Here we focus mainly on the mesh motion algorithm behaviour.
Consider a closed 1D channel with the bathymetry prescribed by the following function where θ is the sloping beach angle. The variable bathymetry region is located leftwards from x s = 2 cot θ ≈ 38.16. At the initial time the shoreline is located at x 00 = cot θ ≈ 19.08. The initial condition is a solitary wave prescribed by the following formulas (which come from the analytical solitary wave solution to the Serre equations, see [22] for example): Variables a and x 0,w = x s + 30 ≈ 68.16 are the initial solitary wave amplitude and position correspondingly. The right boundary is located at ℓ = x 0,w + 30 ≈ 98. 16. All numerical parameters used in this experiment are provided in Table 5. The considered situation is complex enough so that the analytical solution is not available. However, we can simulate it numerically with the predictor-corrector scheme on moving grids presented above. The wave run-up process simulated in silico is depicted in Figure 7(a), while the nodes motion is illustrated in the right panel (b). Please, notice that only every fifth node is depicted in Figure 7(b) in order to improve the visibility of the image. The rightmost node is fixed (i.e. its trajectory is a straight line), while th leftmost node's trajectory coincides with the shoreline motion. One can see that even on domains with variable spatial extent, the mesh motion algorithm is robust enough and it puts the resolution where it is needed. The monitor function used in this computation is It is important to use the wave elevation η(x, t) instead of the total water depth H(x, t) because of the solutions such as the lake at rest state. In these situations the total depth H varies in space, while the free surface elevation η(x, t) ≡ 0. Thus, the usage of the monitor function (4.11) leads to the well-balanced property shown in Theorem 1. Some further validation tests and numerical experiments with this moving grid technique were presented in our preceding work [39] (without giving many details about the algorithm in use). So, we invite the interested reader to consult the work [39] as well.

Open problems
Only in very simple problems (linear and steady) one can study theoretically the choice of the optimal monitoring function. Recently we performed such an analysis for a singularly perturbed Sturm-Liouville problem [38]. In particular, we showed that a judicious choice of the monitoring function allows to increase the convergence order from the 2 nd to the 4 th order (the so-called supraconvergence phenomenon), even if the scheme is just 2 nd order accurate on uniform meshes. So, the first open problem consists in providing a similar analysis for, at least, 1D unsteady or nonlinear problem. The method employed in [38] is suitable for smooth solutions. For discontinuous ones a different technique has to be proposed. Even a linear transport equation is not completely understood from the perspective of constructing a suitable moving grid. It is even possible that for every discretization scheme there is its own "optimal" monitoring function. So, one choice clearly does not fit all possible schemes. Moreover, even if we choose the monitoring function as in (4.11), it was shown in [38] that the error may depend non-monotonically on the coefficient α. Thus, even if the monitoring function form is fixed, one has to set optimal values of coefficients in this function. In the present work our goal was to present the moving grid methodology, which does not involve any interpolation procedure (thus, avoiding the unnecessary smearing of the numerical solution). We showed some reasonable choices of the monitoring function, which are currently used in practice, but these choices have to be refined further in the forthcoming works. The current state-of-the-art in this field is described in [37].

Discussion
The main conclusions and perspectives of this study are outlined below.

Conclusions
In this manuscript we presented a detailed description of a family of predictor-corrector schemes up to the second order in space, which were generalized to the case of non-uniform moving meshes, while having all these properties simultaneously: • Conservation in space • Second order accurate (in smooth regions) • TVD property for scalar problems • Well-balanced character • Preservation of the discrete conservation law in cell centers and cell interfaces • Adaptivity in space and in time * • No interpolation between the meshes, thus, we do not add any unnecessary numerical diffusion • Implementation is trivial, despite all good properties listed above This is the main novelty and contribution of our study to the growing field of adaptive numerical methods. This method was illustrated on a number of relatively serious tests coming from linear, nonlinear and vectorial cases. Every time the gain of using a moving grid (with respect to the fixed one) was clearly visible. Consequently, we can only encourage the community to use these (or similar) adaptive strategies to improve the quality and efficiency of their numerical simulations.

Perspectives
The adaptive moving mesh technology presented in this study relies on the choice of the monitoring function along with several free real parameters. Numerical experiments show that the method performance may depend crucially on the choice of these real numbers. Consequently, it would be highly desirable to obtain some theoretical estimations of the optimal choice of these parameters at least in the case of the simplest linear advection equation (2.1). It constitutes one of the main challenges raised by our study. In the absence of such theoretical indications, one might envisage extensive numerical studies in order to issue some practical recommendations on how to choose reasonable values of parameters. Most probably the case of nonlinear systems of conservation laws will remain * The time step is chosen the largest possible while meeting the stability requirement. essentially empiric in the foreseeable future due to the important theoretical difficulties in its analysis.
From more practical point of view, the predictor-corrector schemes on moving grids presented in this study have to be generalized at least to two-dimensional case as it was done for the Godunov scheme in [4]. Moreover, the spatial adaptivity has to be coupled with time-adaptive strategies such as the PI step size control [32] in order to achieve the full adaptivity.
Finally, there is another goal behind this study. Ultimately, we would like to generalize moving grid methods to some dispersive wave equations [22], which are going to be addressed with the operator splitting approach. In this way, the hyperbolic part would be addressed with methods outlined above. The particularity of dispersive wave equations is that they possess localised solitary wave solutions. In other words, all the dynamics is concentrated in small portions of space, which move perpetually. Consequently, the grid redistribution techniques seem to be the natural choice to simulate complex solitonic gas dynamics [24].