Invariant conservative finite-difference schemes for the one-dimensional shallow water magnetohydrodynamics equations in Lagrangian coordinates

Invariant finite-difference schemes for the one-dimensional shallow water equations in the presence of a magnetic field for various bottom topographies are constructed. Based on the results of the group classification recently carried out by the authors, finite-difference analogues of the conservation laws of the original differential model are obtained. Some typical problems are considered numerically, for which a comparison is made between the cases of a magnetic field presence and when it is absent (the standard shallow water model). The invariance of difference schemes in Lagrangian coordinates and the energy preservation on the obtained numerical solutions are also discussed.


Introduction
The consideration of the shallow water equations in the presence of a magnetic field (SMHD) is a relatively new area of magnetohydrodynamics.One of the first models to describe SMHD was presented in [1] (the Gilman model), and since then, many papers by various authors have been devoted to studying this model, including its stability [2], numerical simulation [3][4][5][6][7][8], and conservation laws [9,10] (see, for example, a more detailed review in [9]).
]ocnmp[ E I Kaptsov and V A Dorodnitsyn Initially, the SMHD equations were used to model the behavior of the solar tachocline, which is a thin layer at the base of the Sun's convection layer.Toroidal magnetic fields are usually considered in this case, which can be assumed to act approximately in the tachocline plane, so that the model [1] describes two-dimensional flows.In some papers (e.g., [9,11]), one-dimensional SMHD equations are considered.In [9], one-dimensional SMHD equations with arbitrary bottom topography were studied in Lagrangian coordinates.In this case, most of the equations of the SMHD system are integrated, and only one equation remains unintegrated.For this equation in [9], a group classification was carried out according to the function describing the topography of the bottom.Conservation laws and invariant solutions in Lagrangian and Eulerian coordinates were obtained.It was found that in Lagrangian coordinates, the one-dimensional SMHD equations differ from the standard one-dimensional shallow water equations [12,13] only by a linear term.Thus, based on the known results for the shallow water equations in Lagrangian coordinates, the conservation laws for SMHD was constructed, which can be done both for differential equations and for their discretizations.
It is known [14,15] that the presence of conservation laws for equations is closely related to the symmetries of the equations.The inheritance of symmetries and conservation laws is also important in the discretization of equations and numerical modeling [16][17][18][19][20].The study of symmetries and conservation laws of discrete models, including finite-difference equations, is the subject of a significant number of works, such as [16][17][18][19][20][21][22].We would like to emphasize significant contributions of Professor Decio Levi to the study of discrete dynamical system symmetries [18] and nonlinear differential difference equations as Bäcklund transformations, elucidating their connection to continuous evolution equations [23].D. Levi has also stated conditions for the existence of higher symmetries [24] and contact transformations [25,26] for discrete equations.Comprehensive research in the area of symmetries and integrability of discrete equations has been reflected in a recent book [27] by D. Levi, P. Winternitz and I. R. Yamilov.
As part of the development of methods of group analysis of discrete equations, in [28,29], finite-difference analogues of the Lagrangian and Hamiltonian formalisms were proposed, which simplify the construction of invariant (symmetry-preserving) schemes and the derivation of their conservation laws.If no invariant Lagrangian or Hamiltonian exists, the alternative method was introduced in [30].Practically, to construct invariant schemes for partial differential equations it is often more convenient to utilize the so-called direct method [31][32][33], which was recently used to obtain invariant conservative schemes for various equations of hydrodynamic type [31,[34][35][36][37][38][39].
In the recent papers [31,34] the authors constructed finite-difference schemes for the one-dimensional shallow water equations in Lagrangian coordinates, which admit symmetries and possess finite-difference analogues of the conservation laws of the original model.On their basis, new schemes were also constructed for various extended models, such as the Green-Naghdi equations [36] and the modified shallow water equations [35].In the present paper, based on the results of [31,34] and the paper [32] devoted to the construction of schemes for wave equations, the authors construct and implement new finite-difference schemes for the one-dimensional SMHD equations in Lagrangian coordinates and in mass Lagrangian coordinates.
The paper is organized as follows.In Section 2, the SMHD equations in Eulerian and Lagrangian coordinates are considered and their dimensionless forms are given.Sections 3 and 4 are devoted to the construction of schemes for the SMHD equations in Lagrangian and mass Lagrangian coordinates and their conservation laws are listed for various bottom topographies according to the group classification results.In Section 5 invariance of finitedifference schemes in Lagrangian coordinates and invariance the constructed schemes in particular are discussed.Section 6 is devoted to the numerical implementation of the schemes.The results are summarized in Conclusion.

The one-dimensional SMHD equations in Lagrangian coordinates
The SMHD model proposed in [1] in Eulerian coordinates has the form where u = (u, v, 0) and H = (H x , H y , 0) are the two-dimensional velocity and magnetic filed vectors, k = (0, 0, 1) denotes the unit vector in the vertical direction, ∇ ′ • is the horizontal divergence operator, k • ∇× is the vertical component of the curl operator, and the constant g ̸ = 0 characterizes the gravitational acceleration.Here h = h 0 + η is the depth of the fluid, where η characterizes a deviation of the free surface from the undisturbed level h 0 > 0. It is also considered |η| < h 0 .
By analogy with the standard shallow water equations, one introduces the function b(x, y) characterizing topography of the bottom [9,40].Then, one writes system (1) with uneven bottom in coordinate form as ]ocnmp[ E I Kaptsov and V A Dorodnitsyn Assuming that h, b, u and H only depend on the single space variable x, the latter system is brought to the form Using (3a) and (3f) one can rewrite (3d) as Thus, by means of (3f), where a is constant that characterizes magnitude of the magnetic force.
Further we consider the model in mass Lagrangian coordinates that is where the Lagrangian coordinates are introduced similarly to the gas dynamics equations [41] through the relations in such a way that the following relation for the differentials dt, ds and dx holds [42] ds = h dx − hu dt, and, The d'Alembert solution of the acoustic equations (4d), (4e) has the form where f 1 and f 2 are arbitrary functions of their arguments.Thus, as they are integrated, we can exclude them from consideration when performing numerical modeling.Notice that in the variables φ, t and s the system is reduced to the single equation [9] which is more suitable for constructing numerical schemes [31,[34][35][36].
The dimensionless form of system (4) is the following ]ocnmp[ where α and g 1 are dimensionless constants characterizing intensity of magnetic and gravity fields.Then, the solution ( 6) becomes The dimensionless form of ( 7) is Further we consider the dimensionless equations, and symbol˜is omitted for brevity.We also assume g 1 = 2 by means of equivalence transformations [9].
The local conservation laws of ( 9) have been obtained in [9].They are listed below depending on the bottom topography according to the results of the group classifications with respect to the function b ′ .
• Case b ′ is arbitrary.Conservation laws of momentum and energy are and the conservation law of mass is just the relation φ ts = φ st ; • Case b ′ = 0 (a horizontal bottom).Center-of-mass motion law and the following alternative form of the conservation law of momentum.
• The case of inclined bottom topography (b ′ = C) is reduced to the case of horizontal bottom topography by means of the transformation [43,44] Further, to avoid confusion with the standard finite-difference notation, we denote h = ρ and φ = x.The standard notation [45] is used where f and f denote the values of f = f (t, s) at the point shifted up and down along the time axis.Similarly f + and f − (or f + and f − ) denote right and left shifts along the space axis.The total differentiations are defined as follows This notation should not cause confusion with standard partial derivatives, since the rest of the discussion is dedicated to finite differences.
The authors would like to emphasize here that the main goal of the present publication is the construction of symmetry-preserving finite-difference schemes that possess the largest possible number of conservation laws.The study of such issues as convergence and stability of the constructed schemes, their well-posedness and regularity of solutions goes beyond the scope of our work.The interested reader can find a detailed discussion of nonlinear stability of schemes for hyperbolic systems and related issues, e.g., in [46].
Considering the form of ( 9), one notices that in case α = 0 (i.e., magnetic field is absent) one gets the one-dimensional shallow water equations in Lagrangian coordinates (see, e.g., [13]).Finite-difference schemes for the shallow water equations have been constructed by the authors in their recent publications [31,34].Here, equation ( 9) differs from the shallow water equations only by a linear term α 2 φ ss , which can be easily approximated.The construction can be accomplished by combining the schemes previously obtained in [31] and the scheme for the linear wave equation derived in [32].Thus, one writes the following finite-difference scheme on an orthogonal uniform mesh ]ocnmp[ E I Kaptsov and V A Dorodnitsyn where h and τ are steps along the space and time axes, and B is some approximation for b ′ .Scheme (11) has the second order of approximation in h and τ .
The conservation laws for the scheme are derived by straightforward extension of the previous results [31,32].Recall that any (local) conservation law of ( 11) can be represented in the divergent form where T and S are conserved quantities of the conservation law, and Λ is called a conservation law multiplier.
The finite-difference counterparts of the conservation laws listed in Section 3 possessed by scheme (11) and the corresponding conservation law multipliers are listed below.One can verify by direct calculations that these conservation laws are satisfied on the solutions of scheme (11) (and therefore can also be represented in form (12)).To do this, it suffices, for example, to solve the equations of the scheme with respect to x t , h + , and τ and substitute the result and its finite-difference shifts into the conservation laws.
• In case b(x) = const ( B = 0), the conservation laws are the following.
• The case of inclined bottom ( B = C) is reduced to the case of horizontal bottom topography by means of the following finite-difference analogue of (10) (see [43]).
• For the parabolic bottom b(x) = x 2 2 (see also [34]), there are two additional conservation laws • For the parabolic bottom b(x) = − x 2 2 (see also [34]), there are two additional conservation laws , one can preserve the conservation law of energy by modifying the first equation of the scheme as follows (see also [31]).
Then, the conservation law of energy becomes 4 Schemes for the SMHD equations in mass Lagrangian coordinates Instead of the three-layer scheme in Lagrangian coordinates one can consider a two-layer scheme in mass Lagrangian coordinates.This can be achieved by introducing a specific approximation of the 'state equation' p = ρ 2 and for the transformation (5) as it was previously done in [31,[34][35][36].The resulting scheme is ]ocnmp[ where the quantity Q is given by and the equation approximates the equation p = ρ 2 .Scheme (15) along with equation ( 17) have the first order of approximation in h and τ .
In mass Lagrangian coordinates, the conservation laws are the following.
• In case b(x) = const ( B = 0): 1. mass • As mentioned above, the case of an inclined bottom is reduced to the case of a horizontal bottom.
• For the parabolic bottom b(x) = x 2 2 : • For the parabolic bottom b(x) = − x 2 2 : • In case the bottom is arbitrary (b = b(x)), the conservation law of energy is As in the case of scheme (11), the validity of the listed conservation laws on solutions of ( 15) can be verified by direct calculations.

Discussion on the invariance of the constructed schemes
In the previous sections we focused on the conservation laws of the constructed finitedifference schemes.Here we mention that the constructed schemes are not only conservative (i.e., possessing conservation laws), but also invariant, that is, they preserve the symmetries of the original differential model.Recall that there is a close relationship between conservation laws and symmetries of an equation, which can be expressed in the form of Noether's theorem [14,15].In the finite-difference case for ordinary difference equations, there is a finite-difference analogue of Noether's theorem [20,28], which makes it possible to derive conservation laws from known symmetries admitted by the equation.In the case of finite-difference schemes for partial differential equations, conservation laws are usually established using the direct method [31,32] or algebraic transformations [37-39, 45, 47].In the latter case, knowledge of the symmetries and conservation laws of the original differential model may suggest the form of finite-difference conservation laws [37][38][39].As far as the authors know, there is no rigorous proof of the converse statement, i.e., from the existence of a conservation law for a scheme, in general, it does not follow that the scheme admits a symmetry corresponding to the conservation law.At the same time, a large number of examples [31,32,[34][35][36][37][38] indirectly indicate the existence of such a connection.In particular, as can be verified, scheme (11) for various bottom topographies admits the same sets of symmetries (Lie algebras) as equation ( 9).Here we refer to [9] where the Lie algebras admitted by ( 9) are stated as a result of the group classification procedure.
The second significant factor that should be taken into account when constructing invariant finite-difference schemes is the preservation of the orthogonality and uniformity of the finite-difference mesh by group transformations.The corresponding mesh invariance criteria are given in [20,48].In contrast to Eulerian coordinates, in Lagrangian coordinates in the one-dimensional case, these criteria are usually satisfied (this was the case for the Lie algebras admitted by the schemes for shallow water equations [31], modified shallow water equations [35], Green-Naghdi equations [36], and this is also true for the obtained SMHD schemes (11) and ( 15)).Thus, the advantage of using Lagrangian coordinates for onedimensional equations is the possibility of constructing symmetry-preserving schemes on simple orthogonal uniform meshes.We note that in higher-dimensional cases the situation becomes more complicated and some symmetries, such as relabeling symmetries, may not satisfy the criteria as it was stated in [49].As simplicial meshes are often used in practical applications for numerical calculations in two and three spatial dimensions, the orthogonality requirement in this case turns out to be quite severe.It can be weakened with the loss of some symmetries (for example, relabeling symmetries seem not to be ]ocnmp[ E I Kaptsov and V A Dorodnitsyn generally admitted by simplicial meshes).Although such meshes may naturally arise in the construction of schemes admitting certain Lie algebras, this is primarily refers to the finite element or finite volume method, the discussion of which is beyond the scope of our study.

Numerical implementation of the constructed schemes
In the present section, some typical problems are solved numerically using schemes of the form (15) in mass Lagrangian coordinates.The numerical implementation is carried out using pseudo-viscosity, which often leads to a loss of accuracy at discontinuities, but which is quite enough to evaluate the qualitative picture of solutions.Scheme ( 15) is linearized using the Newton method (see, e.g., [45]), after which it is transformed to a form suitable for solving through an iterative algorithm and the tridiagonal matrix method [50].Since the linearized version of the scheme may include derivatives of approximations of the function describing the bottom topography, we first discuss the possible forms of these approximations.

Remark on approximations for the function describing the bottom topography
As it was mentioned above, B is some approximation of the derivative b ′ of the function b that describes the bottom topography.For the linearization purposes, we need to know the set of variables on which approximations B may depend.According to the group classification [9] 3 Thus, in general we can restrict our consideration to the function B = B(x, x, u), or B = B(x, x, û).

Linearization of the schemes
We utilize Newton's linearization method, for which we first augment scheme (15) with the equation [45] ω − Ω(ρ, u s ) = 0, where The grid function ω introduces linear and quadratic pseudo-viscosities characterized by the corresponding coefficients ν and µ.When implementing scheme (15), for the equation of motion, instead of ( 16) we also consider the quantity Θ = Q − ω.
Linearizing the scheme, one derives the equations x i − (k) q denotes the value of the quantity q on the k-th iteration, and Notice that in case µ 2 + ν 2 ̸ = 0 conservation law ( 18) cannot be preserved.Therefore, when verifying the conservation law of energy, one should take the results of calculations obtained without pseudo-viscosity or with very small values of ν and µ.

Introducing variables
x , • • • , one derives the three-point equation where A i , B i , D i , and F i depend on the values obtained on the previous iterations.The latter equation can be solved using tridiagonal matrix algorithm [45,50].Then, provided .
Further we consider only problems in which velocity of the flow at the boundaries is held constant.Such problems correspond to the simplest left (L) and right (R) boundary conditions of the form [45] δu L = δu R = 0.

Numerical calculations
First, we consider the dam break problem over horizontal, parabolic and logarithmic bottoms.At the initial moment of time, the liquid is divided in the middle of the computational domain 0 ⩽ s ⩽ S into two regions with heights ρ L and ρ R < ρ L , and the liquid height profile is slightly smoothed near the initial discontinuity.We chose the following calculation parameters.τ = 0.05h, S = 4.0, ρ L = 1.0, ρ R = 0.5, α 2 = 1.6, The viscosity coefficients ν 0 = 1 ÷ 2, µ 0 = 3 ÷ 4, and the step h = 0.04 ÷ 0.1 are vary depending on the problem and the bottom topography under consideration.
By means of equations ( 8b) and (8e), one can approximately write This means that the change in the velocity of the flow depends linearly on the slope of the profile of the longitudinal component H x of the magnetic field.Further, we are only interested in the shape of this profile, therefore, to visualize the profile, we introduce the quantity where the coefficient κ is selected for each problem so that B changes within close to the area of change in the liquid height ρ.
The obtained solutions are depicted in Figure 1, Figure 3 and  ]   The second problem is related to the calculation of the collapse of a liquid column over the inclined bottom (Figure 5) given by the function b In Figure 4, the initial column profile is shown as a dotted gray line.One sees that in the presence of a magnetic field the column collapses faster (gray solid line) than in the case of standard shallow water model (black solid line).Notice that the same result is obtained by calculating for a horizontal bottom followed by applying transformation (13).
In all the experiments, an acceleration of fluid motion in the areas of growth of the gradient of H x is observed (see also trajectories of motion depicted in Figure 2), as predicted by formula (19).In case b ′ = 0, the shock wave velocity D = (α/h 0 ) 2 + g 1 h 0 is obtained using the Rankine-Hugoniot conditions (see, e.g., [5] for details).Here h 0 is the height of the liquid column in front of the wave.Accordingly, the characteristics of the flow are given in Figure 2.
Notice that the difference analogue (15) of ( 5  where, according to the conservation law (18), and the index notation is used: f n k = f (t + nτ, s + kh).H(n) gives the total energy in the computational domain, and its value should tend to constant in the continuous limit [29].
The result for the time interval 0 ⩽ t ⩽ 10 is given in Figure 6.The figure demonstrates that the total energy practically does not change with time.
]ocnmp[ E I Kaptsov and V A Dorodnitsyn

Conclusion
In the present paper, finite-difference schemes have been constructed for one-dimensional shallow water equations in the presence of a magnetic field (the Gilman model) in Lagrangian coordinates for various bottom topographies.According to the recent group classification [9], the cases of an arbitrary, horizontal, inclined, parabolic, and logarithmic bottom are distinguished, while the case of an inclined bottom is reduced to the case of a horizontal bottom by a simple point change of variables.The constructed schemes possess finite-difference analogues of the conservation laws of the original differential model for all the listed cases of bottom topography, and are also invariant, i.e. preserve the symmetries of the original model.The schemes are constructed on uniform orthogonal meshes.In Lagrangian coordinates, they are given on three time layers, and in mass Lagrangian coordinates, they can be given on two time layers by means of a specially selected finite-difference equation of state.
In mass Lagrangian coordinates, the schemes are implemented numerically.Typical one-dimensional problems of dam break and liquid column collapse above an inclined bottom are considered.To demonstrate the effect of a magnetic field on the motion of fluid particles, the Gilman model is compared to the standard shallow water model.Numerical experiments show that the magnetic field accelerates the movement of compression waves along the various bottom topographies, and the destruction of the liquid column above the inclined bottom.For one of the problems calculated with no pseudo-viscosity, the control of the conservation law of energy was carried out, from which it is seen that the invariant conservative scheme preserves the total energy almost with no loss.
, there the following possibilities are of interest: b = C, b = Cx, b = ± x 2 2 , b = ln x, and b = b(x).As in (14), we consider approximations B = [b(x) + b(x)] t x t + xt .The approximations for the listed topographies are the following.1) In case b = C one gets B = 0. 2) In case b = Cx one gets B = C.

Figure 4
for the time t = 0.92.The parabolic bottom profile is given by the function b(x) = 0.5k p (x − 0.5S) 2 , k p = 0.05, and the logarithmic bottom profile is described by b

Figure 1 .
Figure 1.Dam break problem for the horizontal bottom.'SW' (black) is the solution for the standard shallow water equations (α = 0), 'Gilman' is the solution for the case α 2 > 0 (SMHD).Initial profile is given by the dotted line, and the magnetic field gradient B is denoted by the dashed line.

Figure 2 .
Figure 2. Dam break problem for the horizontal bottom.Trajectories x(t, s) for the case α > 0 (gray) and α = 0 (black).The characteristics of the flow are outlined with thick dashed lines.

Figure 3 .
Figure 3. Dam break problem for the parabolic bottom.

Figure 4 .
Figure 4. Dam break problem for the logarithmic bottom.

Figure 5 .
Figure 5. Collapse of the liquid column over the inclined bottom.