Content area
Traditionally, to solve the hydrodynamic equations a Godunov method is used whose main stage is the solution of a Riemann problem to compute the fluxes of conservative variables through the interfaces of adjacent cells. Most numerical Riemann solvers are based on partial or full spectral decompositions of the Jacobian matrix with respect to the spatial derivatives. However, when using complex hyperbolic models and various types of equations of state, even partial spectral decompositions are quite difficult to find analytically. Such hyperbolic systems include the equations of special relativistic magnetohydrodynamics. In this paper, a numerical Riemann solver is constructed by means of a viscosity matrix based on Chebyshev polynomials. This scheme does not require any information about the spectral decomposition of the Jacobian matrix, while considering all types of waves in its design. To reduce the dissipation of the numerical solution, a piecewise parabolic reconstruction of the physical variables is used. The behavior of the numerical method is studied by using some classical test problems.
1. INTRODUCTION
Among the variety of physical processes in high-energy astrophysics, the magnetic field in relativistic gas flows plays a fundamental role. Such processes take place in supernovas with core collapse [1], gamma-ray bursts [2, 3], gravitational impacts [4], magnetars and pulsars [5], black holes [6], neutron stars [7], and jets [8]. Since these processes are highly complicated, corresponding problems can be described numerically only by using supercomputer modeling. Many codes have been proposed to solve such relativistic problems (see, for example, review [9]).
The codes for solving problems of special relativistic magnetohydrodynamics are based on a Godunov-type method and solution of a Riemann problem with partial or full spectral decompositions of the Jacobi matrix. In the general case, nonlinear Riemann problems [10], as well as full spectral problems (even for linearized equations) can only be exactly solved numerically [11], for example, by formulating symmetric hyperbolic systems [12]. In some cases, approximations for extreme waves can be used [13]. Methods based on solving full spectral problems for the Jacobi matrix include the Roe [14], Markin [15], and Glimm [16] methods. The HLLD [17], HLLC [18], HLL [19, 20], Rusanov [21, 22] and Lax–Friedrichs [23] methods are based on solving partial spectral problems. In addition to schemes for solving the Riemann problem which are based on spectral decompositions, methods with artificial viscosity terms have been actively developed [24]. An example of a code based on this approach is called ZEUS [25]. Paper [26] proposes a method based on approximation of the Jacobi matrix by using a Roe polynomial viscosity matrix. The method is extended to the Osher–Solomon method [27]. This approach has been used to develop a code for solving the equations of special relativistic magnetohydrodynamics [28], with radiation transfer in particular [29].
The present paper proposes an approach to constructing a numerical Riemann problem solver based on a polynomial viscosity matrix and a piecewise parabolic reconstruction of the physical variables. Section 2 presents the equations of special relativistic magnetohydrodynamics and a procedure of reconstruction of the physical variables from the conservative ones. In Section 3, a solution to the Riemann problem based on a viscosity matrix is described. In Section 4, the numerical method is verified by using some classical tests. Section 5 is devoted to a discussion of some relevant issues. Section 6 provides conclusions to the paper.
2. THE EQUATIONS OF SPECIAL RELATIVISTIC MAGNETOHYDRODYNAMICS
Let us introduce a vector of physical variables: , here is the density, the velocity components, the pressure, and are the magnetic field components. The model uses dimensionless speed of light and the Lorentz factor calculated by the formula
where . Next we consider the adiabatic index , special enthalpy , and sound speed as follows:
On the basis of the physical variables we introduce a vector of conservative variables , where is the relativistic density, , and denote the relativistic momentum, and is the total relativistic energy defined by the formulas
where and . We consider one-dimensional equations in direction . In this case, the flux vector of the conservative variables may be written as
In this case, we have a system of conservation equations that may be written in vector form as follows:
To solve the equations, we introduce a uniform calculation grid with step for which Godunov’s method is used in the form
To reduce dissipation, we use a piecewise parabolic reconstruction of the physical variables following [9]. The time step is determined by the Courant condition using the limitation of the wave propagation speed by the speed of light.
Since the reconstruction of the physical variables in problems of special relativistic hydrodynamics in the general case cannot be written in analytical form, we need to introduce a numerical procedure based on Newton’s method. For this we introduce the sought-for variable as and . The energy is used in the form
where . In this case Newton’s method for the variable may be written as
the function and its derivative are
where the variable , the pressure , and is calculated by the formula
Let us rewrite the variable in terms of :
Omitting some trivial cumbersome manipulation, we have
Since we consider an ideal gas model, , and, hence, . After the termination of the iterative process the velocity components are calculated as follows:
Let us recalculate the Lorentz factor and the variable by the formula
and the pressure and the density by the formulas
The vector components of the magnetic field are conservative variables.
3. SOLVING A RIEMANN PROBLEM BY USING A VISCOSITY MATRIX
To construct solutions to Riemann problems, we consider left (L) and right (R) cells. To denote the vectors of physical and conservative variables and their fluxes, we use the notation introduced earlier: , , and . Let denote a procedure of reconstructing the physical variables from the conservative functions defined as . The flux of conservative variables at the interface between adjacent cells is calculated as
where , , , and are the vectors of conservative variables and their left and right fluxes, respectively, and is a viscosity matrix. We search for the viscosity term by using Chebyshev polynomials for :
where are coefficients calculated by the formulas
The algorithm for calculating the viscosity term is as follows:
1. Make a piecewise parabolic reconstruction of the primitive variables and .
2. Calculate the vectors of conservative variables , and the vectors of fluxes of conservative variables and .
3. Determine the vector of conservative variables .
4. Calculate the vector of conservative variables .
5. Reconstruct the vector of primitive variables .
6. Determine the vector of conservative variables , where .
7. Reconstruct the vector of primitive variables .
8. Calculate the vectors of fluxes of conservative variables and .
9. Determine the vector of viscosity variables .
10. Reconstruct the vector of primitive variables .
11. Calculate the vector of flux of viscosity variables .
12. Calculate the vector of conservative variables .
13. Calculate the approximation for the vector of viscosity variables.
14. Cycle: for ,and repeat the next steps.
15. Determine the vector of conservative variables , where .
16. Reconstruct the vector of primitive variables .
17. Calculate the vector of flux of conservative variables .
18. Determine the vector of viscosity variables .
19. Reconstruct the vector of primitive variables .
20. Calculate the vector of flux of conservative variables .
Table 1.. Initial data for discontinuity breakdown problems
[See PDF for image]
21. Calculate the vector of conservative variables .
22. Calculate the approximation for the vector of viscosity variables.
23. Go to step 14.
The result is an approximation of the viscosity term.
4. VERIFICATION OF THE NUMERICAL METHOD
The numerical method has been verified by using a set of tests from [30] allowing an analytical solution described in [10]. In the figures, is the density, is the longitudinal velocity, and are the transverse velocities, is the longitudinal component, and are the transverse components of the magnetic field, and is the total pressure.
In all tests we used 800 cells and the Courant number . The initial discontinuity is at point . In all tests, discontinuous initial data are specified in the interval ( is the left state and is the right state), and the numerical solution is calculated at a given time . The initial data are given in Table 1.
[See PDF for image]
Fig. 1
Balsara test no. 1: exact solution (solid line), numerical solution (circles).
[See PDF for image]
Fig. 2
Balsara test no. 2: exact solution (solid line), numerical solution (circles).
[See PDF for image]
Fig. 3
Balsara test no. 3: exact solution (solid line), numerical solution (circles).
[See PDF for image]
Fig. 4
Balsara test no. 4: exact solution (solid line), numerical solution (circles).
[See PDF for image]
Fig. 5
Balsara test no. 5: exact solution (solid line), numerical solution (circles).
[See PDF for image]
Fig. 6
Alfven test: exact solution (solid line), numerical solution (circles).
Figure 1 shows the results of comparison of the numerical and analytical solutions for Balsara test no. 1. At the base of a contact discontinuity there is an artifact in the form of an oscillation obtained as a result of reflection of waves from a slow composite wave at . The slow wave appears in a zone where the magnetic field pressure dominates over the gas pressure, which causes a superposition of the waves. In the exact solution there is only a single wave in this zone [10]. This superposition of waves leads to the formation of a singularity, a density peak. The contact discontinuity, the shock wave, and the fast right magnetic wave are reproduced successfully.
Figure 2 shows the results of comparison of the numerical and analytical solutions for Balsara test no. 2. The aim of the test is to reproduce the weak shock wave in a static gas of constant density with a pressure discontinuity . In the zone of the contact discontinuity a slow magnetic wave is formed, which propagates with a Lorentz factor . Note that the method easily reproduces all components of the solution. The dissipation at the fast magnetic wave is over five cells (three cells for the density function), since the wave is “weak” and close to the slow wave, which “takes over” part of the numerical solution dissipation.
Figure 3 shows the results of comparison of the numerical and analytical solutions for Balsara test no. 3. The test is intended to reproduce the strong shock wave caused by an initial pressure jump of four orders in magnitude. The discontinuity in the magnetic field makes the problem much more difficult, since the fast and slow shock waves propagate close to each other creating a “thin” peak that is reproduced only over five cells with the grid we use. The numerical method successfully reproduces both the location and amplitude of the density peak, as well as the reverse flow along the -axis (this can be seen in the plot of the transverse speed ), which cannot be said about some of the existing methods.
Figure 4 shows the results of comparison of the numerical and analytical solutions for Balsara test no. 4. The purpose of the test is to verify whether the method can reproduce a collision of the gas at ultra-relativistic speeds with . The collision is difficult to reproduce because of high magnetic field pressure. In the figure, in addition to fast magnetic waves one can see two slow magnetic waves moving in different directions. Parasitic oscillations at the shock wave are characteristic of this test.
Figure 5 shows the results of comparison of the numerical and analytical solutions for Balsara test no. 5. The purpose of the test is to verify whether the method can correctly reproduce the Alfven waves formed by discontinuities in the transverse components of the magnetic field and velocity. The Lorentz factor is rather small, . In the figure one can see Alfven waves moving in different directions. On the whole, all components of the solution are reproduced satisfactorily.
Figure 6 shows the results of comparison of the numerical and analytical solutions for the Alfven test. The Alfven test aims to simultaneously reproduce a fast rarefaction wave, an Alfven wave, a left slow magnetic wave, a contact discontinuity, right slow and fast magnetic waves. The test is complicated by the fact that there is a thin shell in the transverse direction (along the -axis). Note that the numerical method successfully reproduces this singularity for and only over five cells.
The results of the computational experiments show that the numerical method successfully and with low dissipation reproduces various types of magnetic-hydrodynamic waves and their combinations.
5. DISCUSSION
Here one should note three important points from the results of the computational experiments:
1. In the present paper, we take a Chebyshev polynomial at according to [29], in which the choice of the degree of the polynomial is justified. Note that in our constructions we have already used such polynomials implicitly: we have the Lax-Friedrichs method at and the Rusanov method at . These methods perform quite well in test problems.
2. Although, on the whole, the verification basis has been worked out with the method rather well, the numerical solution may have some oscillations. Apparently, such oscillations can be reduced by decreasing the degree of the Chebyshev polynomial to or by using a Kolgan-type piecewise linear reconstruction.
3. An evident advantage of the method, which seems to outweigh all its disadvantages, is that there is no need to find solutions to spectral problems. This construction of the method can be easily extended, similarly to [29], to problems with radiation transfer.
These are debatable points based on the author’s own perspective.
6. CONCLUSIONS
An original numerical scheme for solving the Riemann problem for the equations of special relativistic magnetohydrodynamics has been developed. The scheme is based on a viscosity matrix and Chebyshev polynomials. To decrease the dissipation at discontinuities, a piecewise parabolic reconstruction of the physical variables is used. The robustness of the numerical scheme has been demonstrated in some classical test problems.
FUNDING
This work was supported by the Russian Science Foundation (project no. 23-11-00014); https://rscf.ru/project/23-11-00014/.
CONFLICT OF INTEREST
The author of this work declares that he has no conflicts of interest.
Translated from Sibirskii Zhurnal Vychislitel’noi Matematiki, 2025, Vol. 28, No. 1, pp. 75-87. https://doi.org/10.15372/SJNM20250106.
Publisher’s Note. Pleiades Publishing remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
AI tools may have been used in the translation or editing of this article.
REFERENCES
1 Burrows, Adam; Radice, David; Vartanyan, David. Three-dimensional supernova explosion simulations of 9-, 10-, 11-, 12-, and 13-M stars. Monthly Notices of the Royal Astronomical Society; 2019; 485, pp. 3153-3168. [DOI: https://dx.doi.org/10.1093/mnras/stz543] 07539702
2 Aloy, Miguel A.; Rezzolla, Luciano. A powerful hydrodynamic booster for relativistic jets. The Astrophysical Journal; 2006; 640, pp. L115-L118. [DOI: https://dx.doi.org/10.1086/503608] 1097.76073
3 Rivera-Paleo, F.J.; Guzman, F.S. Evolution of jets driven by relativistic radiation hydrodynamics as long and low-luminosity GRBs. Monthly Notices of the Royal Astronomical Society; 2018; 479, pp. 2796-2809. [DOI: https://dx.doi.org/10.1093/mnras/sty1603] 1482.93421
4 Krolik, Julian H.; Armitage, Philip J.; Jiang, Yanfei; Lodato, Giuseppe. Future Simulations of Tidal Disruption Events. Space Science Reviews; 2020; 216, pp. 1-15. [DOI: https://dx.doi.org/10.1007/s11214-020-00680-z] 0912.62117
5 Neutron Stars and Pulsars, W. Becker, Springer, 2009 ( Astrophysics and Space Science Library, vol. 357).
6 Zanotti, Olindo; Roedig, Constanze; Rezzolla, Luciano; Del Zanna, Luca. General relativistic radiation hydrodynamics of accretion flows–I. Bondi–Hoyle accretion. Monthly Notices of the Royal Astronomical Society; 2011; 417, pp. 2899-2915. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2011.19451.x] 1388.76002
7 Radice, David; Bernuzzi, Sebastiano; Perego, Albino; Haas, Roland. A new moment-based general-relativistic neutrino-radiation transport code: Methods and first applications to neutron star mergers. Monthly Notices of the Royal Astronomical Society; 2022; 512, pp. 1499-1521. [DOI: https://dx.doi.org/10.1093/mnras/stac589] 1465.83012
8 Ferrari, Attilio. Modeling extragalactic jets. Annual Review of Astronomy and Astrophysics; 1998; 36, pp. 539-598. [DOI: https://dx.doi.org/10.1146/annurev.astro.36.1.539] 1435.82048
9 Kulikov, Igor. A new code for the numerical simulation of relativistic flows on supercomputers by means of a low-dissipation scheme. Computer Physics Communications; 2020; 257, pp. 107532-107532.4133533 [DOI: https://dx.doi.org/10.1016/j.cpc.2020.107532] 1515.85001
10 Giacomazzo, Bruno; Rezzolla, Luciano. The exact solution of the Riemann problem in relativistic magnetohydrodynamics. Journal of Fluid Mechanics; 2006; 562, pp. 223-259.2263543 [DOI: https://dx.doi.org/10.1017/S0022112006001145] 1097.76073
11 Anton, Luis; Miralles, Juan A.; Marti, Jose M. et al. Relativistic magnetohydrodynamics: renormalized eigenvectors and full wave decomposition Riemann solver. The Astrophysical Journal Supplement Series; 2010; 188, pp. 1-31. [DOI: https://dx.doi.org/10.1088/0067-0049/188/1/1] 1525.85004
12 Freistuehler, Heinrich; Trakhinin, Yuri. Symmetrizations of RMHD equations and stability of relativistic current–vortex sheets. Classical and Quantum Gravity; 2013; 30, pp. 085012-085012.3044369 [DOI: https://dx.doi.org/10.1088/0264-9381/30/8/085012] 1267.83004
13 Leismann, T.; Anton, L.; Aloy, M.A. et al. Relativistic MHD simulations of extragalactic jets. Astronomy & Astrophysics; 2005; 436, pp. 503-526. [DOI: https://dx.doi.org/10.1051/0004-6361:20042520] 1176.46028
14 Koldoba, A.V.; Kuznetsov, O.A.; Ustyugova, G.V. An approximate Riemann solver for relativistic magnetohydrodynamics. Monthly Notices of the Royal Astronomical Society; 2002; 333, pp. 932-942. [DOI: https://dx.doi.org/10.1046/j.1365-8711.2002.05474.x] 1077.76049
15 Donat, Rosa; Marquina, Antonio. Capturing shock reflections: an improved flux formula. Journal of Computational Physics; 1996; 125, pp. 42-58.1381803 [DOI: https://dx.doi.org/10.1006/jcph.1996.0078] 0847.76049
16 Cannizzo, J.K.; Gehrels, N.; Vishniac, E.T. Glimm’s method for relativistic hydrodynamics. The Astrophysical Journal; 2008; 680, pp. 885-896. [DOI: https://dx.doi.org/10.1086/587164] 1456.53070
17 Mignone, A.; Ugliano, M.; Bodo, G. A five-wave Harten–Lax–van Leer Riemann solver for relativistic magnetohydrodynamics. Monthly Notices of the Royal Astronomical Society; 2009; 393, pp. 1141-1156. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2008.14221.x] 1425.76305
18 Honkkila, Ville; Janhunen, Pekka. HLLC solver for ideal relativistic MHD. Journal of Computational Physics; 2007; 223, pp. 643-656.2319227 [DOI: https://dx.doi.org/10.1016/j.jcp.2006.09.027] 1111.76036
19 Schneider, V.; Katscher, U.; Rischke, D.H. et al. New algorithms for ultra-relativistic numerical hydrodynamics. Journal of Computational Physics; 1993; 105, pp. 92-107.1210859 [DOI: https://dx.doi.org/10.1006/jcph.1993.1056] 0779.76062
20 Nakamura, Kouki; Miyoshi, Takahiro; Nonaka, Chiho; Takahashi, Hiroyuki R. Relativistic resistive magneto-hydrodynamics code for high-energy heavy-ion collisions. The European Physical Journal C; 2022; 83, pp. 229-229. [DOI: https://dx.doi.org/10.1140/epjc/s10052-023-11343-y]
21 Kulikov, Igor' Mikhailovich. Using piecewise parabolic reconstruction of physical variables in the Rusanov solver. I. The special relativistic hydrodynamics equations. Journal of Applied and Industrial Mathematics; 2023; 17, pp. 737-749.4730022 [DOI: https://dx.doi.org/10.1134/S1990478923040051] 1550.76180
22 Kulikov, Igor' Mikhailovich. Using Piecewise Parabolic Reconstruction of Physical Variables in Rusanov’s Solver. II. Special Relativistic Magnetohydrodynamics Equations. Journal of Applied and Industrial Mathematics; 2024; 18, pp. 81-92.4757212 [DOI: https://dx.doi.org/10.1134/S1990478924010083] 1550.76181
23 Kulikov, Igor' Mikhailovich; Karavaev, Dmitrii Alekseevich. Using a Low Dissipation Lax–Friedrichs Scheme for Numerical Modeling of Relativistic Flows. Numerical Analysis and Applications; 2023; 16, pp. 326-336.4684125 [DOI: https://dx.doi.org/10.1134/S1995423923040043] 1531.76064
24 Von Neumann, John; Richtmyer, Robert D. A method for the numerical calculation of hydrodynamic shocks. Journal of applied physics; 1950; 21, pp. 232-237.37613 [DOI: https://dx.doi.org/10.1063/1.1699639] 0037.12002
25 Clarke, David A. On the Reliability of ZEUS-3D. The Astrophysical Journal Supplement Series; 2010; 187, pp. 119-134. [DOI: https://dx.doi.org/10.1088/0067-0049/187/1/119] 0232.20012
26 Castro Dı́az, Manuel Jesús; Fernandez-Nieto, Enrique. A class of computationally fast first order finite volume solvers: PVM methods. SIAM Journal on Scientific Computing; 2012; 34, pp. A2173-A2196.2970401 [DOI: https://dx.doi.org/10.1137/100795280] 1253.65167
27 Castro, Manuel J.; Gallardo, José M.; Marquina, Antonio. Approximate Osher–Solomon schemes for hyperbolic systems. Applied Mathematics and Computation; 2016; 272, pp. 347-368.3421111 [DOI: https://dx.doi.org/10.1016/j.amc.2015.06.104] 1410.76325
28 Castro, Manuel J.; Gallardo, José M.; Marquina, Antonio. Jacobian-free approximate solvers for hyperbolic systems: Application to relativistic magnetohydrodynamics. Computer Physics Communications; 2017; 219, pp. 108-120.3677798 [DOI: https://dx.doi.org/10.1016/j.cpc.2017.05.013] 1411.65119
29 Lopez-Miralles, Jose; Marti, Jose Marı́a; Perucho, Manel. On the application of Jacobian-free Riemann solvers for relativistic radiation magnetohydrodynamics under M1 closure. Computer Physics Communications; 2023; 284, pp. 108630-108630. [DOI: https://dx.doi.org/10.1016/j.cpc.2022.108630] 1525.85004
30 Balsara, Dinshaw. Total Variation Diminishing Scheme for Relativistic Magnetohydrodynamics. The Astrophysical Journal Supplement Series; 2001; 132, pp. 83-101. [DOI: https://dx.doi.org/10.1086/318941] 1157.76369
© Pleiades Publishing, Ltd. 2025.