{'fraction_non_alphanumeric': 0.032908427461114824, 'fraction_numerical': 0.01346341841129727, 'mean_word_length': 4.761655141646219, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entropy Stable Finite Volume Approximations for Ideal Magnetohydrodynamics

This article serves as a summary outlining the mathematical entropy analysis of the ideal magnetohydrodynamic (MHD) equations. We select the ideal MHD equations as they are particularly useful for mathematically modeling a wide variety of magnetized fluids. In order to be self-contained we first motivate the physical properties of a magnetic fluid and how it should behave under the laws of thermodynamics. Next, we introduce a mathematical model built from hyperbolic partial differential equations (PDEs) that translate physical laws into mathematical equations. After an overview of the continuous analysis, we thoroughly describe the derivation of a numerical approximation of the ideal MHD system that remains consistent to the continuous thermodynamic principles. The derivation of the method and the theorems contained within serve as the bulk of the review article. We demonstrate that the derived numerical approximation retains the correct entropic properties of the continuous model and show its applicability to a variety of standard numerical test cases for MHD schemes. We close with our conclusions and a brief discussion on future work in the area of entropy consistent numerical methods and the modeling of plasmas.


Introduction
The focus of this summary article is to offer a comprehensive overview of mathematical entropy analysis applied to the equations of ideal magnetohydrodynamics (MHD). As such, we will motivate each step of the mathematical analysis with physical properties, and, where applicable, with everyday experience. In particular, we will discuss the three laws of thermodynamics and call their consequences for the mathematical model in question. Our analysis is first carried out at the continuous level, followed by a rigorous examination of the transition to the discrete level to demonstrate how grid-based numerical algorithms can be designed. The discrete form of the mathematical model will allow us to capture the physical structures and thermodynamic principles of a system of partial differential equations on a computer. By collecting the relevant physical and mathematical analysis into a single work we aim to demonstrate the development of thermodynamically consistent numerical approximations.
The world is full of complex (magneto-)hydrodynamic processes one wishes to understand. From the hurricanes that can devastate communities, the weather forecast for a given day, up to the formation and evolution of the Universe. Over the last two centuries, the conjoined scientific fields of physics and applied mathematics have provided a litany of tools to reproduce the world around us. The breadth of topics covered by the joint effort of these two disciplines is far reaching, from predicting the implication of quantum effects to cosmological simulations reproducing the distribution of galaxies in our Universe.
Fluid dynamics has a wide range of applications in science and engineering, including but not limited to, the prediction of noise and drag of modern aircraft, important for future mobility; the behavior of gas clouds, important for the understanding of the interstellar medium and star formation; the build up and propagation of tsunamis in oceans, important for prediction and warning systems; the impact of the solar wind on Earth, important for an uninterrupted operation of communication satellites, and predictions of instabilities in plasma, important for potential future energy sources like nuclear fusion power plant design.
In this work, we focus on the evolution of macroscopic fluid flows. For example, flows around an airfoil or the formation of new stars out of the diffuse gas in our Galaxy. One classifies a particular fluid into two important groups depending on the level of variation in the fluid's density: Incompressible fluids have a density that remains nearly constant as the fluid evolves, e.g., liquid water at a particular temperature; Compressible fluids have a density that varies as the fluid evolves, e.g., air flow around an automobile [36,61] or highly compressible gas forming the interstellar medium around us [43,44,104]. Due to the multitude of spatial and temporal scales involved with such fluid flows, laboratory experiments are very often impractical or even impossible to carry out. So we move away from a classical experimental setting and instead focus on mathematical models that can predict such complicated flow behaviors accurately. Interestingly, one has to invest much work into the mathematical model to ensure that the solutions generated actually remain physically relevant. Therefore, the applied mathematics we use and describe herein is built from a perspective that respects several fundamental physical principles by construction.
The remainder of this work is organized as follows: Sec. 2 outlines the background physical principles. Sec. 3 highlights how the physical model is represented in a mathematical context. Then, Sec. 3.1 describes the mathematical analysis of entropy on the continuous level. We present the equations of ideal MHD, their usefulness and an entropy analysis in Sec. 3.2. Sec. 4.1 describes one framework to design a numerical approximation that can approximate complicated flow configurations. The bulk of the work and main theorems are contained in Secs. 4.2-4.4, where we demonstrate how to derive a numerical method that is made to respect the thermodynamic properties of the governing equations. We present the utility of thermodynamically consistent algorithms with numerical results in Sec. 5. Finally, our conclusions and a brief outline of future research topics for the ideal MHD equations and mathematical entropy are given in Sec. 6.

A fairly brief introduction to thermodynamics
Thermodynamics is a branch of physics relating the heat, temperature, or entropy of a given physical system to energy and work. The definitions of most physical quantities used in the three laws of thermodynamics analysis are intuitive, but the entropy is a slightly esoteric quantity as we do not use this quantity in everyday life -we have no "feeling" for it. From a mathematical point of view, the entropy is defined as where θ is the specific inner energy of the system and T is the temperature of the system. For now, we will postpone any further details on the implications of entropy and thoroughly discuss it in Section 3.1.
The three laws of thermodynamics are some of the most important laws in all of physics, particularly because they separate possible physical processes, e.g., heat diffusion in a liquid, from the impossible, e.g., perpetual motion: First law of thermodynamics: Energy can neither be created nor destroyed. It can only change forms. For any physical process, the total energy of the closed system remains conserved.
Second law of thermodynamics: The entropy of a closed physical system not in equilibrium tends to increase over time, approaching its maximum value at equilibrium. It cannot shrink.
Third law of thermodynamics: It is impossible for any process, no matter how idealized, to reduce the temperature of a system to absolute zero in a finite number of steps.
All three laws of thermodynamics must be satisfied at the same time, as otherwise a fluid could exhibit strange, i.e., unintuitive (and obviously wrong) behavior. For example, a fluid governed solely by the first law of thermodynamics, i.e., conservation of total energy, could transfer all of its internal energy into kinetic energy. The result would be a very fast, very cold jet of air. This physical state is allowed under total energy conservation, yet has never been observed in nature. To explain this discrepancy requires the second law of thermodynamics. Under the second law, certain transfers between internal energy and kinetic energy are not possible (see also Fig. 1).
Although most thermodynamics calculations use only entropy differences and the zero point of the entropy scale is of no importance, we still mention the third law for purposes of completeness. Interestingly enough, the third law describes the minimum temperature (absolute zero temperature) as being unattainable in somewhat the same way as the maximums speed (speed of light). Special relativity states (and experiments have shown) that no matter how fast something is moving, it can always be accelerated, but it can never reach the speed of light. In a similar way, a system can always be cooled further down, but can never reach absolute zero.
We distinguish between reversible and irreversible thermodynamic processes. First, if work is applied to increase the fluid's kinetic energy at constant pressure and temperature, it can be completely recovered by decelerating the flow. This is called an isentropic 1 acceleration and/or deceleration. Second, if work is done to compress or expand the fluid isentropically, then the volume and therefore the temperature and pressure of the fluid change while its kinetic energy remains constant. Through such processes, the kinetic energy and internal energy (temperature) can be independently and reversibly varied. Other reversible state changes are a combination of these two processes. For example, the work required for an isentropic compression can be provided by an isentropic deceleration. In addition to these reversible state changes, there are irreversible state changes. For example, heat transfer across a finite temperature difference is irreversible as is viscous momentum transfer. Irreversibility implies that certain state changes (like the transferring of heat along a temperature gradient) cannot be reversed. The reversed processes of irreversible processes (e.g., like the transferring of heat against a temperature gradient) are unphysical and have never been observed in the laboratory at any macroscopic scales. Possible and impossible behavior of the fluid system judged by the second law of thermodynamics. Reversible processes conserve entropy over time ("isentropic"). In contrast, irreversible processes generate entropy over time ("entropic"). "Non-entropic" processes, i.e., processes that reduce the entropy are forbidden by the second law of thermodynamics in closed systems. They are never observed in nature on any macroscopic scale.
Hence, the three laws of thermodynamics provide an important role in selecting the tiny subset of physically feasible solutions from all imaginable state configurations. The first law of thermodynamics ensures that a closed system never has more or less total energy than it did at the beginning of its evolution. It brings into play the idea that energy cannot be destroyed, but can only be converted into different forms of itself (kinetic, thermal, magnetic, etc.) inside the system. The second law of thermodynamics provides a mechanism to distinguish further between possible state changes and the impossible ones. The fundamental concept of the entropy of a system is introduced to characterize such processes. For reversible processes the rate of change of the total system entropy with respect to time is zero, i.e., entropy is conserved for reversible state configurations. For irreversible processes, the entropy increases; fluid dynamics where the total system entropy shrinks in time (i.e. reversed processes of irreversible processes) are never observed. Finally, the third law of thermodynamics provides a limitation for the temperature -it can approach zero but never reach it, regardless of the cooling mechanism and the duration of the cooling.

Mission:
We now see the importance of the laws of thermodynamics and how entropy separates between possible and impossible realizations. So, the goal and focus of the remainder of this article are clear. We must first create a systemic approach to discuss the laws of thermodynamics in the language of mathematics on the continuous level. Then, we will mimic the continuous analysis on the discrete level such that we develop numerical algorithms that guarantee approximate solutions to the mathematical models obey the correct entropic behavior.

A fairly brief introduction to hyperbolic partial differential equations
We seek alternatives in theoretical predictions to overcome the limitations and tractability of laboratory experiments, which are often too costly, too time-consuming, too dangerous, or are even impossible to perform as they surpasses the capabilities of our technology. This is because the behavior of flows such as these are of utmost concern for modern industry, science, medicine and society.
It is possible for physicists to interpret the fundamental laws of physics, e.g., the first law of thermodynamics, into corresponding mathematical equations. Then another branch of the natural sciences, applied mathematics is concerned with developing methods to solve the given system of equations. The combined effort of both fields enables us to make predictive statements about how physical systems evolve in time. In particular, we are interested in mathematical models built from partial differential equations (PDEs). A PDE is an equation that involves the rate of change of unknown multivariable function(s) with respect to several continuous variables. To pose the PDE properly also requires information about the initial configuration of the physical system (initial condition) as well as the possible prescription of how the physics should behave at any boundary (boundary conditions), e.g., the surface of an airfoil [29].
Emmy Noether found an astounding relationship that symmetry in a given physical law implies that there must be a conserved quantity associated with that law [58,77]. We note that symmetry in this context refers to a transformation, e.g, rotational or translational, which does not alter the behavior of the physical system. A simple example is if we require that the results of an experiment performed today should remain the same if the experiment were performed tomorrow, i.e., symmetry in time, we arrive at the conservation of total energy. A consequence of Noether's theorem is that we can build mathematical conservation laws (which are governed by PDEs) and immediately relate them to a particular set of fundamental physical laws. This offers a powerful link between physics and mathematics as well as provides a multitude of practical computational tools for studying fluid mechanics.
Conservation laws that govern fluid dynamics are often members of the class of hyperbolic PDEs. A standard prototype of a hyperbolic PDE is the advection equation where a is a constant, t is the temporal variable, and x is the spatial variable. Without loss of generality we assume a > 0. Given an initial condition for the function we have an initial value problem where we want to determine the values of u(x, t) for positive values of time. We see that a possible solution to (3.1) is where it is known that (3.3) is actually the unique solution of (3.1) [29]. The formula (3.3) tells us two important things. First, the solution at any time t * is the original initial condition shifted to the right (as a > 0) by the amount at * . In other words, the solution at a point (x, t) depends only on the value of ξ = x − at. The lines in the (x, t) plane on which ξ is a constant are referred to as characteristics. The parameter a is the speed of propagation of the solution along each characteristic line. Thus, the solution of the advection equation (3.1) can be regarded as a wave that propagates with speed a without changing its shape. Second, it appears that the solution of (3.1) only makes sense if u is differentiable. However, the general solution (3.3) prescribes no differentiability constraint on the initial condition u 0 . In general, this means that discontinuous solutions are allowed for hyperbolic conservation laws [29]. An example of such discontinuous solutions are shock waves, which are a feature of non-linear hyperbolic conservation laws, that occur when two characteristic curves cross [29,91,99].
Many of the complicated, compressible fluid applications discussed in the previous section are modeled by the solution of non-linear conservation laws. However, it is precisely the possibility of discontinuous solutions of non-linear hyperbolic PDEs that introduces a challenge when computing their solution. To be explicit we now consider a general one-dimensional system of conservation laws of the form where q is the vector of conserved variables and f (q) is the non-linear vector flux. We know that regardless of the continuity of the initial data it is possible for the solution of a non-linear problem to develop discontinuities, e.g., shocks, in finite time [62,90]. Therefore, solutions of the system (3.4) are sought in the weak sense. A weak solution (sometimes called a generalized solution) of a PDE is a function for which the derivative may not exist everywhere, but nonetheless satisfies the equation in some sense [29]. Unfortunately, weak solutions of a PDE are, in general, not unique and must be supplemented with extra admissibility criteria in order to single out the physically relevant solution [21,48,59,64,75,96,97].
We previously discussed that the laws of thermodynamics provide guidelines in separating the possible state changes from the impossible. In the mathematical sense, solutions of reversible processes are smooth and differentiable. However, discontinuous solutions arise due to irreversible processes, such as within shock waves, where entropy must increase. It is straightforward for the mathematical model to recover the first law of thermodynamics provided one of the conserved variables in q is the total energy. It is trickier for the mathematical model to accurately capture the second law of thermodynamics because the mechanism that governs irreversible processes (entropy) is typically not explicitly built into the PDE [48,59,64,95,96,97]. Hence, we have to fulfill an auxiliary conservation law of entropy for reversible processes and must guarantee that entropy increases for irreversible processes.
3.1. Continuous mathematical entropy analysis for general systems. We outline some of the necessary additions to the mathematical analysis in order to incorporate the entropy inequality. The investigation of hyperbolic PDEs modeling compressible fluid dynamics that are mathematically entropy consistent has been the subject of research for approximately the past fifty years [45,46,48,59,60,62,64,75,93,95,96,97,98]. Here we indicate an important difference in the sign convention between the physical and mathematical communities. In physics entropy production is increasing with time, whereas in mathematics entropy is modeled with a decreasing function of time. This difference in notation is largely due to the fact that mathematicians want an upper bound on the entropy [73]. In contrast, physics -specifically statistical mechanicsrelate entropy to the number of microscopic configurations Ω that a thermodynamic system can have when it is in a state specified by some macroscopic variables. Since it seems natural to assume that each microscopic configuration is equally probable, the number of configurations can be denoted as where k B is the Boltzmann constant, a commonly found constant that can be interpreted as the bridge between macroscopic and microscopic physics [61].
Boltzmann's constant, and therefore also the entropy, have dimensions of energy divided by temperature as the number of microstates (and hence its logarithm) is a unitless quantity. From a different point of view, physicists -in contrast to mathematicians -interpret entropy as a measure of disorder within a macroscopic system. In fact, macroscopic systems tend to spontaneously evolve into the state with maximum entropy, commonly called thermodynamic equilibrium. As an example, imagine the temperature is zero everywhere in the system. Then every particle in the system is ordered (in a specific way) and the number of microstates Ω is minimal. However, if we increase the temperature, then the particles start to move because of their thermal energy, i.e., Brownian motion and the number of possible microstates increases. It is evident that any orderly arrangement would be very rare because all particles are in random thermal motion. Even if we would have the a special configuration at some time again, this would be for a brief instant and the particles would again move to some other configuration.
We consider the one-dimensional entropy analysis of the general hyperbolic conservation law (3.4). First, we define a scalar entropy function S(q) that satisfies [48,73,96,97] (i) The entropy function S is a strongly convex function of q.
(ii) The entropy function S is augmented with a corresponding entropy flux function F (q) in the x−direction such that These conditions are referred to as the convexity condition and the compatibility condition, respectively. The function S and its corresponding flux F form an entropy-entropy flux pair (S, F ). It is important to note that the entropy-entropy flux pair need not be unique, e.g., as was demonstrated by Harten [48] and Tadmor [96]. From the entropy function we define a new vector of entropy variables Next, we contract the hyperbolic conservation law (3.4) with the entropy variables v to obtain From the compatibility condition (3.6) and the definition of the entropy variables (3.7) we see that (3.8) transforms into the entropy conservation law for smooth solutions. If we account for non-smooth regimes as well then the entropy equality modifies to be the entropy inequality [97] (3.10) where we used the mathematical convention of a decreasing entropy. So, we see that with a particular choice of entropy function and satisfying a compatibility condition on the entropy flux we can build the second law of thermodynamics into the system in a mathematically consistent way.
There are other important quantities related to the mathematical entropy analysis necessary to discuss the PDE either in the space of conservative variables or in the space of entropy variables. Due to the strong convexity of the entropy function there exist symmetric positive definite (s.p.d.) Jacobian matrices we use to transform back and forth between the variable spaces. The Jacobian matrix to move from conservative space to entropy space is found by computing the Hessian of the entropy function (3.11) It is interesting to note that it was shown by Friedrichs and Lax [38] that multiplication of the conservation law "on the left" by the matrix H −1 symmetrizes the system (3.4). Because the matrix H −1 is s.p.d., we can immediately compute the other entropy Jacobian through inversion of (3.11) to obtain (3.12) H := ∂q ∂v .

3.2.
Example: Ideal magnetohydrodynamic equations. We are interested in capturing as much physics as possible with the mathematical model. Therefore, we consider a set of hyperbolic PDEs that offer applications to a wide range of physical phenomena. To do so we desire the ability to mathematically describe the evolution of plasmas (electrically ionized gases) such that we can model flow configurations local to our planet as well as to the Universe that surrounds us as a high electrical conductivity is ubiquitous in astrophysical objects. One example of interesting phenomena we wish to model is the star formation process.
To extend the description of fluid dynamics to gases that respond to electric and magnetic fields we couple the inviscid Euler equations [61] with a reduced set of Maxwell's equations [53] to obtain the ideal magnetohydrodynamic (MHD) equations [37]. The "ideal MHD" approximation involves some fundamental assumptions on the fluid we want to model with it: The fluid approximation: The essential approximation in the "ideal MHD" model is that variations in the macroscopic quantities we observe are slow compared to the time scale of microscopic processes in the fluid we want to model. We require that a sufficient number of gas particles is available to make meaningful definitions of macroscopic properties like density, velocity, and pressure. Instantaneous relation for the electric field: There is an instantaneous relation between current densities and the electric field (Ohm's law). The term "instantaneous" means that the involved processes have to take place on small scales that average out over temporal and spacial scales smaller than the ones of interest. For this to happen, we require perfect conductivity that can also be understood as vanishing electrical resistivity. Electric conductivity is given if the constituents of the plasma are at least partly ionized, which is sufficient in most astrophysical environments. Neutrality: The fluid is (overall) electrically neutral. This approximation is most often valid.
The coupling of gas dynamics and magnetic fields of a perfectly conducting fluid creates a three-dimensional hyperbolic system of conservation laws of the form where I is the 3 × 3 identity matrix. The density , momenta u, total energy E, and the magnetic fields B are the primary conserved quantities. The variables u and p are the fluid velocities and the thermal pressure, respectively. No evolution equations for the electric field are present since we know from Ohm's law that it vanishes in a perfect conductor. However, it should be noted that the electric field vanishes only in the co-moving or fluid reference frame we use that moves with the flow. In any other reference frame there is an electric field that would require attention. We close the system under the assumption of an ideal gas which relates the thermal pressure to the total energy by (3.14) p where the ratio of specific heats, γ, is the adiabatic constant.
We immediately see that the first law of thermodynamics, i.e., total energy conservation, is the fifth equation in (3.13). However, the second law of thermodynamics, i.e., the entropy inequality, is not explicitly built into the ideal MHD equations. As mentioned in Sec. 2, this can be problematic. It is well known (see e.g. [73]) that magnetized fluids governed solely by these conservation laws can, in fact, display very strange behavior. A prime example is an air conditioner that is able to transfer internal energy (heat) to kinetic energy (motion), resulting in a very fast, very cold jet of air. This strange air conditioner is allowed under (3.13), yet is never observed in practice. Thus, we see that special care must be taken to ensure that the mathematics respects the physics. We do this by explicitly adding the evolution of the thermodynamic entropy in our derivations. By this, we can incorporate the laws of thermodynamics explicitly into the numerical approximation. As this air conditioner requires state changes that are not valid in our scheme it is rendered impossible and the importance of our effort is immediately apparent.
An additional constraint for the ideal MHD equations to model the evolution of magnetized fluid dynamics we must ensure that the divergence on the magnetic field is zero, as is dictated by one of Maxwell's equations, specifically Gauss' law for magnetism, (3.15) ∇ · B = 0.
We refer to this property as the divergence-free condition (although it is sometimes called the involution [3,49] or solenoidal [53] condition in the literature). The geometrical meaning of (3.15) is that magnetic field lines have "no ends," i.e., regions of reduced field strength cannot be local since magnetic field lines are not allowed to meet any monopolar singularities, e.g., see Fig. 2, like seen with electric field lines and charges. Hence, in contrast to other macroscopic quantities like density, a change in magnetic field strength must be accommodated by changes in the field morphology on a larger scale causing additional difficulties in the numerical modeling.
On the continuous level (3.15) is assumed to always be satisfied. However, even if the divergencefree constraint is satisfied with the initial conditions, it is not necessarily true that it will remain satisfied through the discrete evolution of the equations [7]. So, we see that the divergence-free constraint provides an important indicator to decide if flows remain physically meaningful during their numerical approximation.
Just like entropy conservation, the satisfaction of the divergence-free condition (3.15) is not explicitly included in the PDE system of the ideal MHD equations (3.13) which introduces additional complexities in the (numerical) modeling of plasmas. This is problematic as mathematical analysis requires a single equation to model physical processes rather a set of uncoupled equations. We next demonstrate that the divergence-free condition can be built into the non-linear PDE system such that the equations can mathematically model all of the observed physics of the system simultaneously.
3.2.1. Continuous mathematical entropy and the divergence-free condition. Surprisingly, the issue of entropy conservation and satisfaction of the divergence-free condition are inextricably linked for the ideal MHD system. To demonstrate this we consider the entropy function of the ideal MHD system to be the physical entropy density [5] (3.16) where s( , p) = ln(p) − γ ln( ). The corresponding entropy fluxes are where u = (u, v, w) T . Therefore, we compute the entropy variables to be which contracts the ideal MHD equations into the entropy conservation law (3.9) or entropy inequality (3.10) depending on the smoothness of the solution q. The entropy Jacobian (3.12) is (3.19) Also, we compute the entropy potential defined as The usefulness of the entropy potential (3.21) will reveal itself later in Sec. 4.2 where the potential function φ plays an important role to mimic the continuous entropy analysis at the discrete level.
Just as in Sec. 3.1 we contract the PDE system into entropy space with the goal to recover entropy conservation. So, we premultiply (3.13) with the entropy variables and after some algebraic manipulation we arrive at the expression On the continuous level we know that the divergence-free constraint (3.15) is always satisfied and the last term in (3.22) vanishes, thus recovering entropy conservation. However, as previously discussed, a numerical approximation will not necessarily satisfy the divergence-free condition. Therefore, a discrete evolution of the ideal MHD equations cannot guarantee entropy conservation without addressing errors made in the divergence-free constraint.
In the course of an entropy analysis of the ideal MHD equations (3.13) Godunov [46] observed that the divergence-free condition can be incorporated into the ideal MHD equations as a source term proportional to the divergence of the magnetic field (which, on a continuous level, is a clever way of adding zero to the system). This additional source term not only allows the equations to be put in symmetric hyperbolic form [5,46], but it also restores Galilean invariance, i.e., Newton's laws of motion hold in all frames related to one another by a Galilean transformation. This source term, first investigated in the multi-dimensional numerical context by Powell et al. [81], augments the ideal MHD system to become (3.23) While it is now possible to recover entropy conservation from the form (3.23) there is a significant drawback. The new set of equations (3.23) loses conservation in all but the continuity equation. In numerics this is problematic as non-conservative approximations encounter difficulties obtaining the correct shock speeds or the physically correct weak solution, e.g., [100].
Fortunately, there is an alternative source term that incorporates the divergence-free condition into the ideal MHD system and maintains the conservation of the momenta and total energy equations. Janhunen [55] used a proper generalization of Maxwell's equations when magnetic monopoles are present and imposed electromagnetic duality invariance of the Lorentz force to derive an alternative MHD system From the structure of the Janhunen source term in (3.24) and the last three components of the entropy variables (3.18) we see that contracting the new MHD system (3.24) into entropy space will cancel the extraneous terms present in (3.22). Thus, the ideal MHD system with the Janhunen source term (3.24) can be used to recover entropy conservation on a discrete level, as demonstrated in [107]. Also, because the system remains conservative in all the hydrodynamic variables a numerical approximation will satisfy the Lax-Wendroff theorem [65] and retains the ability to capture the correct shock speeds.

Discrete entropy analysis for the ideal MHD equations
We have discussed the continuous components of the mathematical entropy analysis for the ideal MHD equations. With the necessary physical and mathematical background complete, we are prepared to embark on the mission set out at the end of Sec. 2. That is, we are equipped to describe a discrete numerical approximation of a system of hyperbolic PDEs that remains entropy consistent. We divide the discussion of the discrete entropy analysis into four parts: Sec. 4.1 provides background details on finite volume methods, a particularly useful numerical approximation for systems of hyperbolic conservation laws. Next, Sec. 4.2 presents a baseline entropy conservative scheme for the ideal MHD equations. We have seen that physics dictates that entropy must dissipate in the presence of discontinuities. Therefore, we augment the baseline entropy conservative scheme with numerical dissipation in Sec. 4.3 to create an entropy stable scheme. Finally, Sec. 4.4 provides detailed information on how to evaluate the numerical dissipation term to ensure that the approximation avoids unphysical states, e.g., negative values of the density.
For clarity, we provide a brief summary of the overall approximation in Sec. 4.5.

4.1.
Finite volume numerical discretization. Based on the current state-of-the-art of scientific knowledge, we investigate the behavior of fluids by using the conservation laws for mass, momentum, total energy, magnetic field, and entropy. However, it turns out that only a minority of very simplified hyperbolic PDEs have a known explicit form of an analytic solution. Although the knowledge of such simplified solutions are extremely valuable in helping understand how fluids behave, they can only rarely be used for practical applications in physics and engineering. Most often, the simplified equations are based on a combination of approximations and dimensional analysis; empirical input is almost always required [30]. This approach is very successful when the system can be described by one, or at most, two parameters. However, many flows require several parameters for their specification, rendering the required experiments extremely difficult or even impossible as, often, the necessary quantities are simply not measurable with present techniques or can only be measured with insufficient accuracy.
Hence, apart from rare exceptions, it is mandatory to be able to solve the underlying equations in a different way to be able to investigate the behavior of flows accurately. A solution came with the invention of sufficiently powerful computers. Although many key numerical techniques had already been developed more than a century ago, they were of little use before the invention of the computer. The performance of computers has increased spectacularly since the 1950s. In the early days of computing, computers were able to work on a few hundred instructions per second, whereas modern computers are designed to perform in the 10 12 operations per second range, with no apparent sign of this trend to slow down. At the same time, the ability to store large amounts of data has increased in a similarly rapid fashion.
It requires little to no imagination to see that computers might be a very valuable tool in trying to understand complex fluid flows. Out of this recognition the field of computational fluid dynamics (CFD) emerged. This field solves the governing equations numerically in order to approximate the solution of a particular problem. The range of applications extends from small scale simulations running on personal computers within seconds to large scale simulations running for months on the largest supercomputers that exist to date. This, of course, depends on the complexity and the desired accuracy of the flows to be simulated. For example, if we are interested in the design of a new car, we would like to simulate the flow around a moving car in a wind tunnel, we would have to build a car model and blow air in a wind tunnel at it. However, to obtain accurate results, the floor of the wind tunnel would need to move with the same velocity as the air flow to avoid boundary effects in the flow pattern. At the same time, however, we would want to have the car stand perfectly still, which is of course technically difficult. However, using a computer simulation we neither have to build a car model nor is setting the boundary condition for the moving floor boundary a difficulty.
It may have become apparent that it is not practical to try to cover the whole field of CFD in a single work. Also, the field is evolving so rapidly, that we risk the discussion becoming out of date in a short while. Hence, we focus here on the specific methods designed to solve the system of equations of interest numerically.
To numerically approximate, we have to discretize the mathematical model, i.e., the set of partial differential equations (3.13) themselves. Besides advanced methods like spectral methods, e.g., [12,13,89], the most important methods are finite difference (FD), e.g., [69,92], finite element (FE), e.g., [8,57], and finite volume (FV), e.g., [67,68]. The introduction of numerical approximations of fluid flow introduces the notion of resolution. Once we move to the discrete level, the approximation of a particular fluid is only known at a particular set of points. Therefore, it is intuitive to note that we need "enough" discrete points to capture the flow features we wish to study. The issue of resolution for a numerical method is discussed later in this section. However, we note that each of the numerical methods introduced yield the same result if the numerical resolution is sufficiently high.
Due to their simplicity of presentation and particular relevance to the solution of systems of hyperbolic conservation laws [67,68] we focus mainly on the development of a FV scheme in this work. The basic idea can be broken into four steps: 1. Divide the domain, the simplest being a square box, into a number of Cartesian cells.
2. Approximate the solution q by the average of each quantity at the center of each cell.
3. Approximate the fluxes and any source term at each of the interfaces in between the cells. 4. The spatial discretization creates a systems of ordinary differential equations (ODEs) to evolve in time. We evolve the approximate solution in time using an explicit time integration technique, e.g., Adams-Bashforth [82], Runge-Kutta [10], or strong stability preserving Runge-Kutta [47]. We note it is also possible to use implicit ODE intergration techniques [82], but they will not be discussed further.
We provide a sketch of the first three steps of the finite volume approximation in Fig. 3. With these four steps it is possible to take a hyperbolic PDE, including a source term, from the continuous level where n denotes the discrete level in time, q n i is the average solution in a grid cell,f is a numerical flux function where i ± 1 2 indicates the value at cell edges, and s i is a consistent discretization of the source term. Also, ∆t is the time step for the explicit time integration method and ∆x is the size of the grid cell. We present the fully discrete form of the FV method (4.2) to identify the final form of the numerical scheme. Next, we provide specific details of how each of the terms in (4.2) arise.
A first question that arises in the FV discretization comes from the fact that the value of the flux is multiply defined at each of the interfaces coupling the individual cells. At a given interface do we evaluate the flux from the centroid value on the left or on the right? To resolve this issue and create a unique flux at the interface, while maintaining discrete conservation of the variables q, we introduce the notion of a numerical flux function denoted asf [68,85,99]. We delay the explicit form of the numerical flux to Sec. 4.2. Presently, we provide a more in-depth background on the construction of FV schemes for hyperbolic conservation laws.
The first step in designing a numerical model for the flows of interest is to divide the continuous space around us into small portions at discrete locations in space, i.e., we discretize space. This discretization of space can be achieved by placing a suitable numerical grid onto our continuous geometry (see Fig. 4 for an illustration in two spatial dimensions). The spatial approximation that we select is of special importance as it will affect the quality of the result significantly. If the discretization is too coarse, we risk oversimplifying the approximated flow of interest and lose most of the details, obtaining a solution which may not be sufficient enough to gain the required conclusions. This amplifies, that, although the field of CFD and its applications are thrilling, one always has to bear in mind that numerical results are always only approximate solutions. If we, Divide the domain into a set of elements labeled with integers.
Approximate the solution with cell-centered averages.
Approximate the fluxes at cell interfaces labeled with half integer values. however, discretize with a resolution that is unreasonably high, we increase the computational complexity, and hence the costs of our simulation exponentially without gaining a noticeably higher degree of accuracy in the result.
We briefly compare the FD and FV methods to further motivate the choice of the FV method for the solution of hyperbolic conservation laws. FD type methods are the oldest of the mentioned methods and are accounted to the mathematician Leonhard Euler, who introduced this method in the mid-eighteenth century [30]. The starting point is the conservation equations in differential form as given by (3.4). Derivatives in the original equations are replaced by difference quotients, resulting in a set of algebraic equations for each nodal point on the grid (cf.  the quantities (like density, velocity, etc.) at the current and a number of neighboring cells are used to approximate the derivatives. In principle, this method can be applied to any type of grid. On regular Cartesian grids, this method is particularly easy to understand and implement.
i+2, j+1 i j (c) FV: Cell centers, polygons Figure 5. Illustration of the discretization using the FD and FV method in 2D. i and j are the indices counting in x and y direction, respectively.
The FV approach is perhaps the most popular method used to approximate hyperbolic conservation law as it is particularly easy to understand. It uses an integral form of the conservation equations as a starting point. Here, the solution is subdivided into a finite number of contiguous cell volumes. The conservation equations can then directly be applied to each volume. The discrete points are located at the centroid of each volume (cf. Fig. 5B). Interpolation is used to extrapolate the value of the volume centered values to the cell faces where the numerical fluxes are computed. Surface and volume integrals can easily be approximated using appropriate quadrature methods. The grid need not be aligned to a coordinate system. Hence, the FV method is suitable for any type of grid, as well as for discretizations of a domain that contain complex geometries. The cells can be constructed using a polyhedron of any dimension (in 2D polygons like triangles, rectangles, pentagons, etc., in 3D polyhedrons like tetrahedra, prisms, hexahedra, etc.). The only requirements are that the sum of the individual volumes must cover the whole computational domain and that the overlap is limited to the boundaries of the cells. On any given volume the flux entering through a particular boundary is identical to the flux leaving the adjacent volume. Thus, the method is conservative by construction, provided surface integrals are the same for the volumes sharing the same boundary [30,68]. An interesting detail is that all terms that need to be approximated have a direct physical meaning, which is a major reason for its widespread use and popularity among physicists and engineers.
We seek a numerical discretization of the ideal MHD equations. Additionally, because the divergence-free condition is important to guarantee entropy conservation we include a source term. So, let us look at the general conservation law with a source term, is the vector valued flux function depending on the physical variables of the solution, s is the source term, and x is the vector of the independent spatial variables. This differential conservation law can be understood as the local formulation of the integral conservation law equation which may be obtained if we take the volume integral over the total volume of a particular cell i, V i , and apply the divergence theorem to the second term. This gives, where S is the surface of the volume, and n is the (unit) surface normal in outward direction (see Fig. 6). Figure 6. Illustration for the integral conservation law.
On integrating the first term in (4.4) to get the volume average of q, q, this yields which is a general result equivalent to (4.3). It is immediately clear that the FV scheme, which is derived from this equation, is conservative by construction as cell averages change through the edge fluxes. By integrating over time, we obtain the evolution equation for the cell averaged For the derivation of the FV scheme, we now assume that the boundary of grid cell i with volume V i is given by k i edges K i,j with one certain neighboring grid cell. Further, we make the assumption that the source term contributes at each of the cells edges. We can then rewrite (4.6) into the form For the sake of convenience, we assume a one-dimensional geometry with equidistant discretization of space and time with uniform step sizes ∆x and ∆t. The grid intervals shall be given by Note that this can be done without loss of generality.
From (4.7) we immediately obtain is the source term approximation found by applying the trapezoid rule to the right hand side integral in (4.5). One particular difficulty that arises is that the numerical flux has to be computed at the interfaces, whereas the FV scheme is concerned with cell-centered quantities only. Hence, we have to reconstruct the values at the cell faces to compute the numerical fluxes.
The most simple approximation for the numerical flux is when we assume that the face values are identical to the center values, i.e.
(4.10)f n i+ 1 2 =f (q n i , q n i+1 ). Since the approximate solution is constant in each grid interval, this approximation is commonly called piecewise constant reconstruction: Unfortunately, it turns out that the piecewise constant reconstruction gives poor accuracy in smooth regions of the flow [68]. Moreover, shocks tend to be heavily smeared and poorly resolved by the numerical grid. A more sophisticated and significantly more accurate reconstruction technique can be derived if we take neighboring cells into account for the approximation of the interface quantities. A particularly simple higher-order method is the piecewise linear reconstruction [102]: The interface values can now be computed using Piecewise linear (4.13) with (4.15) q n The interface values of cell i are denoted with q i,L and q i,R at x = x i−1/2 , and x = x i+1/2 , restrictively.
The slope may be computed by using a generalized difference quotient with a free parameter α ∈ [0, 1]: For α = 0.5, this approximation is second-order accurate in space: However, it has been shown that this reconstruction is not unconditionally stable for any constant value of α [76]. Therefore, one has to take special measures to ensure that the monotonicity of the cell averaged data is preserved. Using the ansatz given by (4.15) we see, e.g, in Fig. 7B, that the monotonicity at the x i+3/2 interface is violated. A commonly used slope limiter that avoids spurious oscillations and violations of the original monotonicity is the minmod limiter [86]. We define where the minmod function is defined as In practical terms, this means that if the slopes on either side of the boundary have different signs, which could either indicate spurious oscillation or a (local) extremum, the value will be reconstructed using a slope of zero in order to avoid unphysical overshoots of the interpolation. Similarly, if the slopes on either side have the same sign, the slope that is less steep is always chosen. The general idea of the minmod algorithm is to restrict the ratio between adjacent jumps in cell averages where such jumps could indicate oscillatory behavior. Clearly, such a scheme will not develop spurious oscillations near a jump. Note that this, however, may smooth out extrema even when they are desired in the solution.
We note that the particular choice of the reconstruction is a very important for the accuracy and robustness of the resulting numerical FV scheme. In this work we restrict ourselves to the use of the simple minmod reconstruction (4.16), however many other reconstructions are available in the literature, e.g., [4,15,86,88,94,102].
To complete the construction of the FV scheme we must find a suitable numerical flux function f . In addition, we seek to construct a numerical flux function that conserves (or preserves the correct sign) of the entropy predicted by the continuous equations. The construction of such a numerical flux is the focus of the next section. As a final note on the discretization, it is important to guarantee that the numerical flux is a consistent approximation of the physical flux. This requirement can be understood as a desire that the approximation will be exact if the value of q is identical on the left and right of the interface, i.e., We now have a numerical scheme for the system of hyperbolic conservation laws with a finite volume discretization in space and some type of explicit time integrator. The full discretization is completed once we prescribe a suitable, consistent numerical flux functionf . We apply the additional constraint that the resulting numerical flux must recover the correct entropic properties of the continuous problem. Much work has been done over the past 20 years to create entropy consistent numerical approximations for hyperbolic conservation laws (both scalar and systems) [5,6,14,16,17,18,31,32,33,34,35,40,42,50,51,52,56,66,74,83,84,96,97,98,103,105,107,108].

4.2.
Entropy conservative numerical flux for the ideal MHD equations. In the classical framework of Tadmor the discrete entropy conserving numerical flux is defined as an integral in phase space [96,97], i.e., For non-linear problems, like the ideal MHD equations, the expression of the physical flux components in terms of the entropy variables are quite complicated. Thus, evaluating the integral (4.19) exactly is often impossible. Plus, evaluation of (4.19) with a sufficiently high accuracy quadrature rule to createf i+ 1 2 is computationally intensive. To avoid so much computational overhead, we will describe an affordable entropy conservative numerical flux function. Affordable entropy conservative flux functions are available for the shallow water equations [33,105], the compressible Euler equations [16,52], and the ideal MHD equations [17,107].
To derive an entropy conserving flux for the ideal MHD equations we will mimic the continuous entropy analysis from Sec. 3.1 on the discrete level through the use of the discrete entropy analysis tools developed by Tadmor [96,97]. We noted in Sec. 3 that hyperbolic PDEs which model fluid mechanics are rotationally invariant. Consider an arbitrary hyperbolic system of conservation laws where we explicitly indicate that the fluxes f , g, and h are functions of the conserved quantities q. Given a rotation matrix S, if it is true that that for any non-zero vector n = [n x , n y , n z ] T then the conservation law (4.20) is said to satisfy the rotational invariance property. For example, the ideal MHD equations possess the rotational invariance property, with the rotation matrix defined as where n, , and m are three orthogonal unit vectors. As the rotation matrix (4.22) is orthogonal it is trivial to compute the inverse matrix to be (4.23) That is, the explicit form of fluxes g and h do not need to be known. We recover their action from f and three orthogonal unit vectors. Thus, from the rotational invariance property we can find a flux vector projected in any direction completely from the flux vector f .
Rotational invariance offers a significant advantage in the context of the FV method. That is, it is sufficient to consider each of the spatial directions to be split, i.e., we can immediately apply algorithms developed in the one-dimensional context to the multi-dimensional context, e.g., [25]. To apply the one-dimensional numerical flux into the y or z−directions, any direction dependent quantities, i.e., the velocity and magnetic field components, are rotated such that the x−direction numerical flux remains valid. Quantities not listed in Table 1, e.g., density, do not dependent on the directionality and consequently, do not require rotation. Table 1. Rotation of multi-dimensional quantities to compute the numerical flux using the one-dimensional formulation in x−direction.
This rotational invariance greatly simplifies the discrete entropy analysis of the ideal MHD equations as it is sufficient to design the numerical flux for the one-dimensional case. Note that we have to rotate three times for any update. First, we rotate a copy of the physical quantities to compute the numerical fluxes so that we can apply the fluxes as if they were to be computed in x−direction. Applying the rotation given in Table 1 twice results in the original reference frame of the computational domain. Hence, we apply the same rotation both on the computed numerical fluxes and the source term to rotate them back into the reference frame of the computational simulation. Although this procedure sounds computationally intense, it actually is not, since rotating means only interchanging four quantities. Moreover, we gain the noteworthy advantage that we can apply the same flux and source term computation algorithms in all three directions. Not only does this save us from repeating the (algorithmically quite intense) derivation for the other directions, it also allows us to thoroughly test our implementation in comparable simple and standardized one-dimensional test cases. Once we are comfortable with the obtained performance, we can immediately apply this code to two and three-dimensional systems with a comparably negligible chance of introducing mistakes in the code that would come with the implementation of new y and z−direction fluxes.
To begin let's assume that we have two adjacent states (L, R) with cell areas (∆x L , ∆x R ). We discretize the ideal MHD equations semi-discretely and examine the approximation at the i + 1 2 interface. We suppress the interface index unless it is necessary for clarity. Also note the factor of one half from the source term discretization.
We can interpret the update (4.26) as a finite volume scheme where we have left and right cell-averaged values separated by a common flux interface.
We premultiply the expressions (4.26) by the entropy variables to convert the system into entropy space. From the definition of the entropy variables (3.18) we know that S t = v T q t , hence a semi-discrete entropy update is If we denote the jump in a state as · = (·) R − (·) L and the average of a state as {{·}} = 1 2 ((·) R + (·) L ), then the total update will be We want the discrete entropy update to satisfy the discrete entropy conservation law. To achieve this, we require [33,97,107] We combine the known entropy potential φ in (3.21) and the linearity of the jump operator to rewrite the entropy conservation condition (4.29) as , and denote the constraint (4.30) as the discrete entropy conserving condition. This is a single condition on the numerical flux vectorf , so there are many potential solutions for the entropy conserving flux. Recall, however, that we have the additional requirement that the numerical flux must be consistent (4.18). In the following, we develop the expression for an entropy conserving fluxf for the ideal MHD equations. We will see that the discretization of the source term s i+ 1 2 plays an important role to ensure that discrete entropy conservation (4.30) is maintained.
With the definition of the entropy variables and the formulation of the discrete entropy conserving condition (4.30) we are ready to derive a computationally affordable entropy conserving numerical flux. We have previously defined the arithmetic mean. To derive an entropy conserving flux we also require the logarithmic mean A numerically stable procedure to compute the logarithmic mean is described by Ismail and Roe [52,App. B].
Just as Chandrashekar [16] did for the compressible Euler equations, we develop a baseline numerical flux function that is both kinetic energy preserving [54] and entropy conserving (KEPEC) for the ideal MHD equations. We explicitly derive the flux in the x−direction, but generalization to higher dimensions is straightforward through symmetry arguments [25]. To do so we introduce notation for the inverse of the temperature therefore the entropy variables (3.18) are rewritten to be

Theorem 1 (Kinetic Energy Preserving and Entropy Conserving (KEPEC) Numerical Flux).
If we define the logarithmic mean (·) ln (4.31), the arithmetic mean {{·}}, and discretize the source term in the finite volume method to contribute to each element as (4.34) then we determine a discrete entropy and kinetic energy conservative flux to be where (4.36) Proof. To derive an affordable entropy conservative flux for the one-dimensional ideal MHD equations we first expand the discrete entropy conserving condition (4.30) componentwise to find (4.37) To determine the unknown components off we must expand each jump term in (4.37) into linear jump components. This will provide us with a system of eight equations from which we determine f . To do so we utilize two identities of the linear jump operator We first note that the jump in the physical entropy s written in terms of logarithmic mean (4.31) of and β is .
Thus the left hand side of the expression (4.37) expands to become (4.40) Next we expand the right hand side of (4.37) into a combination of linear jumps: Finally, we expand the source term contribution on the right hand side of (4.30) where we leave the individual components of the source term general as a consistent approximation will reveal itself in the later analysis. Note that we used ∂B1 ∂x ≈ B1,R−B1,L ∆x = B1 ∆x in the discrete case.
Every term in the discrete entropy conservation condition (4.37) is now rewritten into linear jump components. Though algebraically laborious this provides us with a set of eight equations for which we can determine the yet unknown components in the entropy conserving numerical flux. Next we gather the like terms of each jump component. Once we have grouped all the like terms for each linear jump it will become clear how to discretize the Janhunen source term in order to guarantee consistency. Gathering terms from (4.40), (4.41), and (4.42) we determine the system of eight equations to be (4.43) With the collection of equations (4.43) we find a rather alarming result. We know that the sixth component of the physical flux for the ideal MHD system is zero, i.e., f 6 = 0. However, we have found in our entropy conservation condition that the sixth component of the numerical flux must satisfy In general, we cannot guarantee that the right hand side of (4.44) will vanish. In one spatial dimension the argument could be made that B 1 = 0 (as it is a constant [37]) and there is, in fact, no issue. However, this assumption is too restrictive to discuss higher dimensional extensions of the entropy conservative flux formulae. The assumption that B 1 = constant revealed extra terms which otherwise would have been hidden from the analysis.
To remove this inconsistency introduced by the B 1 equation we discretize the source term to cancel the problematic terms (4.44). We compare the structure of the extra terms in (4.44) and the degrees of freedom s 6 , s 7 , and s 8 to determine a consistent discretization to cancel the extraneous terms in the B 1 equation in (4.44): The source term at each interface of the cell has an identical structure. We collect the total source term discretization in cell i for clarity: It is straightforward to check the consistency of the source term discretization (4.46).
We substitute the source term discretization (4.45) into the entropy constraint (4.43) and find that the source term components exactly cancel the extraneous terms in the B 1 equation (4.44). Thus, we recover a consistent term for the sixth numerical flux component and it is now true that (4.48)f 6 = 0.
Finally, we are able solve the remaining seven equations (4.43) and determine the components of the numerical flux to be (4.49) The newly derived numerical flux (4.49) conserves the discrete entropy by construction. Next, we verify that the numerical fluxf KEPEC is consistent to the physical flux. It will make the demonstration of consistency more straightforward if we simplify the fifth component of the physical flux to be (4.50) Now, if we assume that the left and right states are identical in the numerical flux (4.49), we find that Thus, we have shown that the numerical flux given by (4.49) is consistent and entropy conservative.

4.3.
Entropy stable numerical flux for the ideal MHD equations. The KEPEC numerical flux presented so far conserves the entropy in the semi-discrete approximation. However, the solution of the ideal MHD equations (or any hyperbolic PDE for that matter) may develop discontinuities, e.g., shock waves or contact discontinuities, in finite time even for smooth initial data. We know in the presence of discontinuities that the conservation law for the entropy function (3.9) must be replaced by the entropy inequality (3.10), e.g., [96]. Thus, we seek to add numerical dissipation to the KEPEC flux so that the entropy is guaranteed to be dissipated (or conserved for smooth, well-resolved solutions), thereby ensuring that a discrete version of the entropy inequality holds. A typical way to add dissipation in an FV approximation is via the definition of the numerical flux function.
To create a kinetic energy preserving and entropy stable (KEPES) numerical flux function we use the entropy conserving flux (4.35) as a baseline and subtract a general form of numerical dissipation, e.g., where D is a dissipation matrix. The form of dissipation (4.52) is motivated by the fact that we want less dissipation in regions of the flow that are smooth (or well-resolved) and more dissipation in regions containing discontinuities. Of utmost concern for entropy stability of the approximation is to formulate the dissipation term (4.52) such that it is guaranteed to cause a negative contribution in the discrete entropy equation (4.28).
We select the dissipation matrix to be the absolute value of the flux Jacobian for the ideal MHD 8-wave formulation [46,80,81] in the x−direction: where A is the flux Jacobian for the homogeneous ideal MHD equations and P is the Powell source term written in matrix form, i.e., The fact that we force the positivity of D = |A + P| does not guarantee that the dissipative term (4.52) is negative [5]. Thus, in the remaining sections we will motivate a reformulation of the dissipation term (4.52) in order to restore negativity.
It is important to distinguish that we use the Janhunen source term (3.24) to derive an entropy conservative numerical flux function in Sec. 4.2. However, to design an entropy stable approximation we require that the eigendecomposition of the flux Jacobian matrix can be related to the entropy Jacobian (3.19). This particular scaling, first examined by Merriam [73] and explored more thoroughly by Barth [5], requires the PDE system to be symmetrizable. Previous analysis of the ideal MHD equations [5,46] has demonstrated that the Powell source term is necessary to restore a symmetric ideal MHD system. We reiterate that the altered flux Jacobian is used only to derive the dissipation term. Just as Lax-Friedrichs [67] differs from Roe [85] in the construction of a dissipation term, we use the Powell source term only to build our dissipation term. Thus, no inconsistency with the previous entropy conserving flux derivations is introduced. 1. We compute the eigendecomposition for the symmetrizable MHD system written in the primitive variables.
2. We use previous results from Roe and Balsara [87] and rescale the eigenvectors to remove degeneracies.
3. We recover the, now stabilized, eigendecomposition for the matrix D.
For brevity we provide minimal details of the derivation and present the final result for the eigendecomposition of D. However, complete details can be found in [107].
To discuss the eigenstructure of the matrix D it is easiest to work with primitive variables, which we denote ω, and convert back to conservative variables, denoted by q, when necessary. We first write the ideal MHD system modified by the Powell source term in terms of the conservative variables (4.55) ∂q ∂t + D ∂q ∂x = 0.
We are free to move between primitive and conservative variables in the system (4.55) with the matrix and from conservative to primitive variables with M −1 . Then we can rewrite the system (4.55) in terms of the vector of primitive variables ω To describe the eigenstructure of the symmetrizable ideal MHD system flux Jacobian in conservative variables D we first investigate the eigendecompostion of B, the flux Jacobian in primitive variables: From (4.58) we see that we can convert the resulting eigendecomposition to conservative variables with straightforward matrix multiplication and the identity So, once we compute the eigendecomposition we can recover the eigendecomposition of the flux Jacobian matrix in conservative variables as Forgoing a large amount of algebra, we present the eigendecomposition of the matrix D. The 8-wave formulation supports eight traveling wave solutions • two fast magnetoacoustic waves (±f ), • two slow magnetoacoustic waves (±s), • two Alfvén waves (±a), • an entropy wave (E), • a divergence wave (D).
with eigenvalues where c f , c s are the fast and slow magnetoacoustic wave speeds and c a is the Alfvén wave speed. The double eigenvalue u represent the entropy wave and divergence wave. The divergence wave is a direct result of including the Powell source. This is because the Powell source term turns the divergence wave into an advected scalar which is directly reflected in the eigenstructure. The values for the characteristic wave speeds may be written as In (4.64) the plus sign is for the fast speed c f and minus sign is the slow speed c s . It is known that the right eigenvectors of D may exhibit several forms of degeneracy that are carefully described by Roe and Balsara [87]. We follow the same rescaling procedure of Roe and Balsara for the fast/slow magnetoacoustic eigenvectors. The algebra is simplified greatly if we introduce the parameters The parameters (4.66) have several useful properties: The parameters α f,s measure how closely the fast/slow waves approximate the behavior of acoustic waves. In the rescaling process we utilize several identities that arise from the quartic equation for the magnetoacoustic wave speeds ±c f,s Applying the identities (4.67) and (4.69) it is possible to rewrite the eigenvectors for the fast/slow waves in a more stable form in terms of the parameters α f,s [87].
The matrix of right eigenvectors is given by (4.70) R = [r +f |r +a |r +s |r E |r D |r −s |r −a |r −f ] , with the eigenvectorsr, and corresponding eigenvalues λ [5,25,87] Entropy and Divergence Waves: where we introduced several convenience variables (4.74)

Entropy scaled right eigenvectors.
We have a symmetrizable matrix D with a complete set of eigenvalues and right eigenvectors. We next utilize a previous result from Barth [5] which provides a systematic approach to restructure a general eigenvalue problem to a symmetric eigenvalue problem. To do so we rescale the right eigenvectors of an eigendecomposition with respect to a right symmetrizer matrix in the following way: Lemma 1 (Eigenvector Scaling).
Let A ∈ R n×n be an arbitrary diagonalizable matrix and S the set of all right symmetrizers: Further, let R ∈ R n×n denote the right eigenvector matrix which diagonalizes A, i.e., A = RΛR −1 , with r distinct eigenvalues. Then for each B ∈ S there exists a symmetric block diagonal matrix T that block scales columns of R, R = RT, such that which implies Proof. The proof of the eigenvector scaling lemma is given in [5].

Theorem 2 (Kinectic Energy Preserving and Entropy Stable (KEPES) Numerical Flux).
If we apply the diagonal scaling matrix to the matrix of right eigenvectors R (4.70), then we obtain the identity [5,73] (4.79) that relates the right eigenvectors of D to the entropy Jacobian matrix (3.19). For convenience, we introduce the diagonal scaling matrix Z = T 2 in (4.79). We then have the guaranteed entropy stable flux interface contribution Proof. We define the dissipation term in the numerical flux (4.52) to be where the eigendecomposition of D is given by (4.63) and (4.70). We define entropy stability to mean the approximation guarantees that the entropy within the system is a decreasing function, satisfying the following inequality From the previously computed discrete entropy update (4.28) and the condition (4.29) we find the total entropy update (now including the dissipation term (4.81)) to be (4.83) due to the design condition (4.29) on the kinetic energy preserving and entropy conserving flux, To ensure entropy stability, we must guarantee that the right hand side in (4.83) is non-positive. Unfortunately, it was shown by Barth [5] that the term may become positive in the presence of very strong shocks. However, we know from entropy symmetrization theory, e.g [5,73], that the entropy Jacobian H, given by (3.19), is a right symmetrizer for the flux Jacobian that incorporates the Powell source term D. Therefore, with the proper scaling matrix T we acquire the identity The rescaling of the right eigenvectors of D to satisfy the identity (4.85) is sufficient to guarantee the negativity of (4.84). We see from (4.84) and (4.85) So with the appropriate diagonal scaling matrix Z we have shown that (4.86) is guaranteed negative because the product is a quadratic form scaled by a negative. We use the right eigenvectors from (4.70), the constraint (4.85), and are able to determine the diagonal scaling matrix to be from which the the diagonal scaling matrix T (4.78) follows directly.

4.4.
Discrete evaluation of the dissipation term. We now know the structure of the dissipation term at the interface between two cells. However, in the discrete setting, there is still the question of where to evaluate the dissipation operator. We know the specific averaging for the baseline entropy conserving flux (4.35). However, there is an open question of how to evaluate the dissipation term in (4.80) discretely at some mean state. Much care is taken in the baseline fluxf KEPEC by using very specific averages to guarantee discrete entropy conservation, e.g., [17,97,107]. We will see that an equal amount of care must be taken for the dissipation term to guarantee that the numerical flux remains applicable to a wide variety of flow configurations [106].
Specifically, this section provides detailed derivations on the discrete forms of the matrices R, Z, and Λ needed to evaluate the entropy stable numerical flux function (4.80).

4.4.1.
Evaluation of the entropy Jacobian H. Somewhat counterintuitively, the first step to evaluate the dissipation term in the entropy stable flux (4.80) is to determine a discrete evaluation of the entropy Jacobian matrix H. This is due to the critical first step of (4.86) where we took To avoid unphysical dissipation we want to build the average state of H in such a way that equality holds in (4.88) whenever possible [24].

Lemma 2 (Failure of equality).
It is not possible to construct a consistent discrete matrix H such that holds for all eight components of the ideal MHD system.
Proof. We seek a proof by contradiction. Thus, we assume that averages in the entries of H exist in such a way that the relation q = H v holds discretely. We derive the entries of the matrix H through a step-by-step process computing the solution of 64 linear equations: The procedure is to multiply each row of H with the expanded jump in the entropy variables. By examining each equation individually, we are able to determine all unknown entries of the discrete matrix.
The derivation of the first row of H is straightforward and therefore serves as an excellent example for the derivation technique. First, we use properties of the linear jump operator (4.38) to expand the jump in the conservative (4.91) as well as the entropy variables According to (4.90), the entries of the first row of H can be obtained by solving From this equation, we directly obtain the entries of the first row of the discretized entropy Jacobian, The same procedure is repeated over each row to identify the remaining entries.
Unfortunately, we find that such a forthright solution of (4.90) leads to an asymmetric discrete matrix H. This asymmetry indicates that the scheme will not be entropy stable nor will it satisfy the necessary eigenvector scaling realtionship (4.79). Therefore, it is not possible to derive a discrete symmetric matrix such that the equality q = H v holds exactly for all components of q.

Theorem 3 (Discrete entropy Jacobian H).
If we take the discrete dissipation matrix to be of the form (4.97) then we obtain exact equality in (4.89) for each component except the total energy, i.e., Proof. Due to result in Lemma 2 we know that equality cannot hold for each component in (4.89). However, if special care is taken during the expansion of the total energy term, a matrix H that obeys the required property can be found. It guarantees equality in all but the jump in the total energy term where the equality reduces to an asymptotic one. The modified jump in total energy reads (4.99) We replace E with E in (4.91). Now, we apply the previously described tactic from the proof of Lemma 2 to solve (4.90). Forgoing the computational details we arrive at the discrete matrix (4.100) Clearly, the discrete entropy Jacobian matrices (4.97) is symmetric. Furthermore, it has been shown using Sylvester's criterion that the discrete matrix H is symmetric positive definite (s.p.d.) [24].

4.4.2.
Discrete right eigenvector, R, and scaling, Z, matrices. With the knowledge of the discrete symmetric matrix H the next goal is to determine the discrete evaluation of the matrix of right eigenvectors R and the diagonal scaling matrix Z such that the identity (4.101) H = RZ R T , will hold discretely.
Theorem 4 (Discrete right eigenvector R and scaling Z matrices).
If we take the discrete right eigenvector matrix (4.104) and discrete scaling matrix (4.108), then we have the discrete entropy scaled eigenvector relation Due to their complexity we present the specific form of R and Z in three parts. First, we give the specific averages of the convenience variables (4.74) (4.103) Next, we give the discrete right eigenvector matrix Alfvén Waves: Magnetoacoustic Waves: and finally the discrete diagonal scaling matrix (4.108) Proof. We have already computed the discrete entropy Jacobian H given in (4.97). The strategy to obtain the discrete evaluation of the right eigenvector and scaling matrices is similar to that used to create H. That is, we seek average states in the right eigenvectors (4.70) and scaling matrix (4.87) such that we have the discrete entropy scaled eigenvector relationship Once this relationship is satisfied discretely it will define unique averaging procedures for three of the terms in the dissipation operator (4.80) while retaining the almost equal property (4.98).
It is straightforward to relate the entries of the matrix H and determine the 64 individual components of the matrices R and Z. We will explicitly demonstrate two computations to outline the general technique and justify the somewhat unconventional averaging strategies employed in the final form. After this brief outline of the derivation, we will arrive at the complete, discrete evaluation for the right eigenvector and scaling matrices.
We begin by computing the first entry of the first row of the system, which must satisfy (4.110) The α f,s variables satisfy many useful identities [87] that we must recover discretely. Namely, Thus, there is still some freedom in the underlying averaging of theα f,s and the wave speedsĉ f,s as long as the identities (4.111) hold. We use the first identity in (4.111) and chooseˆ = ln in (4.110) to guarantee that H 1,1 = ( R Z R T ) 1,1 . Therefore, we have determined the first row of the discrete eigenvector matrix, denoted by R 1,: , as well as five entries of the scaling matrix Z to be (4.112) R 1,: = α f ln 0α s ln 1 0α s ln 0α f ln , The second example is the second entry of the second row of the system given by (4.113) We now select the particular average for the sound speed in the second identity of (4.111) to be and also takeˆ = ln andû = {{u}} to guarantee equality in (4.113). So, we have determined some important terms necessary for the second row of the eigenvector matrix R 2,: (4.115) where the averages of the two wave speedsĉ f,s are still arbitrary. We apply this same process to the remaining unknown values from the condition (4.102) and, after a considerable amount of algebraic manipulation, determine the unique averaging procedure for the discrete eigenvector and scaling matrices, given in (4.103)-(4.108). We note that the derivations as well as the constraint (4.102) have been verified using the symbolic algebra program Maxima [72].    Therefore, the discrete evaluation of the diagonal matrix of eigenvalues is Proof. We begin with the augmented ideal MHD equations written in primitive variables Note that the flux Jacobian D in conservative variables and the flux Jacobian in primitive variables B are similar. That is, they have the same eigenvalues, but different eigenvectors. Therefore, we can use the much simpler matrix B (4.59) to determine the discrete evaluation of the eigenvalues.
We make the discrete ansatz, where we discretize the update of the primitive variables in the spatial dimension (4.57) to be (4.122) where we used the trapezoidal rule to approximate the integral on the right hand side of (4.122) using the left and right states. We immediately see that the discretized version of the coefficient matrix is the continuous coefficient matrix, B, arithmetically averaged in each entry, We compute the discrete eigenvalues of B to be (4.124) After many algebraic manipulations and introducing the new average wave speeds (4.118), we find the greatly simplified form of the discrete eigenvalues given in (4.117). The discrete diagonal matrixΛ follows immediately.

4.5.
Summary of the discrete entropy analysis. The derivations and proofs in this section were quite technical. Therefore, we found it appropriate to collect the main results in a short summary for easier reference. From the derivations performed we have constructed a fully discrete and unique evaluation for the KEPES numerical flux function (4.80) given by where the KEPEC flux is given by (4.35), the discrete right eigenvector matrix R is found in (4.104), the discrete scaling matrix Z is (4.108), the discrete eigenvalues matrix Λ (4.120), and the jump in the entropy variables v (4.92).
So, we can now describe the complete FV update for each cell in the numerical approximation to be of the form (in one spatial dimension) wheref KEPES is given by (4.125) and the source term discretization s i (4.34).
The numerical discretization and proofs provided throughout this section creates an approximation for the ideal MHD equations that is conservative in the hydrodynamic variables and discretely obeys the correct entropic behavior predicted by the continuous hyperbolic PDEs, a fact we demonstrate in Sec. 5.2.
We close our derivation by looking back at the three laws of thermodynamics defined at the beginning of Sec. 2 and scrutinize how they are incorporated into the numerical approximation we derived: First law of thermodynamics: Energy can neither be created nor destroyed.
We fulfill this law by having the total energy of the system be a conserved quantity that can only change due to fluxes. The corresponding entry in the source term is zero. Hence, mathematically, there is no way that energy could be either created or destroyed in a closed box simulation.
Second law of thermodynamics: Entropy is either conserved or increases over time. It cannot shrink.
We have made sure that our model guarantees the correct sign in the entropy evolution. Entropy can only be conserved (KEPEC scheme) or "generated" (KEPES scheme, i.e. entropy conserving scheme + dissipation term). However, entropy can never be "destroyed," as the dissipation scheme is of quadratic form as shown in (4.86).
Third law of thermodynamics: It is impossible for any process to cool a system to absolute zero in a finite number of steps.
From the ideal gas law, we have T ∝ p. If we assume = const., we can re-write (3.16) to isolate the dependence of the entropy on the temperature, where we used the physical sign convention of an increasing entropy function. It can easily be seen that for vanishing temperature, the entropy diverges: As our scheme ensures that the entropy evolution is well-behaved, it is clear that T = 0 cannot be reached in any finite amount of time.
Thus, the Mission we gave at the end of Sec. 2 is complete!

Numerical results
We apply the described entropy stable scheme to several numerical examples. Section 5.1 provides a brief discussion of the implementation and its place in a state-of-the-art numerical framework. Next, Sec. 5.2 verifies the primary results of this work, e.g., that the baseline FV scheme is entropy conservative (from Sec. 4.2) and improved robustness of the discrete evaluation for the dissipation operator (from Sec. 4.4). Finally, we apply the entropy stable numerical solver to explore its utility for several standard test problems for the ideal MHD equations in Sec. 5.3.

5.1.
A fairly brief discussion of software. For completeness, we provide a few details regarding the implementation of the KEPES FV numerical approximation. This is partially motivated due to the fact that this aspect of CFD can feel overwhelming because there are so many choices. Thus, we note again, that the discussion is limited to a very small sector of the computational mathematics community.
We implement the three-dimensional FV solver in the FORTRAN programming language. As mentioned earlier, we are most interested in the numerical scheme to reach a broad user base and maintain the ability to model a large range of physical phenomena. Therefore, we choose to implement the entropy stable framework into the multi-scale multi-physics simulation code FLASH [26,39]. A chief motivation of this choice is because FLASH is publicly available (http://flash.uchicago.edu) and has a wide international user base. Complete details on the implementation of the entropy stable ideal MHD solver is provided by the authors in [25].
We know that complex flow phenomena, particularly in multiple spatial dimensions, exhibit a large range of spatial and temporal scales. Increasing the accuracy in time might be straightforward, e.g., selecting a higher order explicit time integration scheme. The increase in spatial accuracy has more features. We outlined how to increase spatial accuracy through the use of sophisticated reconstruction techniques in Sec. 4.1. However, accuracy can also be increased by adding more cells into the computational grid. Obviously, it is best to do this grid refinement locally and automatically near the flow complexities which may develop, e.g., vortical structures, and de-refine the grid where the flow is less complex, to reduce the overall computational cost. The simulation code we are using, FLASH, is an adaptive mesh refinement (AMR) code that uses the PARAMESH library [78]. Here the grid is organized in a block structured, oct-tree adaptive grid, supporting local refinement and de-refinement as can be seen, e.g., in Fig. 12 [25].
A final aspect of the numerical approximation is how to determine the time step for an explicit time integration scheme such that the method is stable and convergent [69]. The selection of the time step for hyperbolic systems is equation dependent as it involves the fastest wave speed. In order for the results of the numerical approximation to remain relevant we select the time step using the Courant, Friedrichs, and Lewy (CFL) condition [19,20,63]. One version of the general CFL condition for the three spatial dimension approximation on Cartesian meshes is given by [25] where λ d max is the largest wave speed at the current time traveling in the d = {x, y, z} direction and CFL ∈ (0, 1] is an adjustable coefficient. If λ max is known exactly, then the choice CFL = 1.0 may be adequate for stability [99, pg. 222]. However, λ max is usually computed in some approximate way. For example, the discrete eigenvalues of the ideal MHD system discussed in Thm. 5. Thus, a more conservative choice for the CFL coefficient is typically used in practice. For example, in the numerical results discussed in Secs. 5.

5.2.1.
Verification of spatial convergence. The smooth Alfvén wave test [100] is used to compare the accuracy of MHD schemes for smooth flows. The initial circularly polarized Alfvén wave propagates across a periodic domain [25]. The Alfvén wave speed in x-direction is |v A | = B 1 / √ = 1 and thus, the wave is expected to return to its initial state at each time t ∈ N. The initial conditions listed in Table 2 ensure that the magnetic pressure is constant.   [25].
For sufficiently smooth fields, i.e., in cases where discontinuous features are absent, the reconstruction technique used is designed to achieve more than first-order accuracy. To test the accuracy of our scheme, we run several simulations with varying resolutions and compute the L 1 and L 2 errors for the quantity B ⊥ = B 2 . We fix the time step to a very small value (∆t = 1 × 10 −5 ) in order to diminish the influence of the time integration scheme on this test. The obtained errors are plotted as a function of the number of grid points in logarithmic scale in Fig. 8 and are listed in Table 3.   As can be seen, the standard finite volume scheme is at most first-order accurate in space.
With second order reconstruction, close to second-order accuracy is achieved already at a coarse resolution. Given a sufficient numerical resolution (N 64), the minmod limited interface values converge to higher accuracy. We do not expect pure second-order accuracy to be achieved using minmod reconstruction in this test, as the slope limiter will always interpret the extrema of the setup Alfvén wave as discontinuities and drop to first-order accuracy at these points. However, this effect becomes less significant the higher the numerical resolution. Note that we used 2 nd order reconstructed interface values in this test as we knew that the solution will always stay smooth. However, we remind the reader that it cannot be used in the general case as discontinuous solutions may develop even from smooth initial data. We note that SSP RK time integration schemes are weighted convex combinations of explicit Euler steps.
On the basis of all these considerations, we measure the temporal accuracy and independently confirm the order of each time integration schemelisted above. In order to do so, we disable the entropy stabilization, i.e., the resulting scheme is entropy conserving. As our test of choice, we run the two-dimensional version of the Brio and Wu magnetohydrodynamical shock tube problem [9] with varying fixed time step lengths ∆t. This test includes discontinuities, a magnetic field and is performed in multiple dimensions and starts from discontinuous initial conditions. Hence, it utilizes the full set of features of the entropy aware scheme we derived in this work. We keep the previously used periodic boundary conditions to eliminate any possible influence from the boundaries of the domain and to ensure that we observe a closed system. We construct the two-dimensional initial conditions by rotating the one-dimensional conditions (see Table 4) at a 45 • angle. The fluid is initially at rest on either side of the interface.
Note that the entropy conserving scheme cannot describe systems with discontinuities as it cannot add the physically needed dissipation. We limit the end time step to t end = 0.001 when the oscillations have not grown too large as to cause numerical instabilities. Our sole intention is to show that even under high stresses our scheme is still capable of conversing the thermodynamic entropy correctly. We point out that the obtained solution is physically not correct by intention.
We ran a number of simulations, for each time integration scheme, using logarithmically equally spaced time steps and plot the measured entropy conservation error. We see that for the first Domain size x min = 0, x max = 1 Initial shock position x shock = 0.5 Boundary conditions all: periodic Simulation end time t max = 0.1 Adiabatic index γ = 2.0 (artificial)  Figure 9. Entropy conservation test, explicit Euler time integration ( Fig. 9) and second order (Fig. 10) temporal approximations, we obtain the expected order of convergence in time. For the third order time integrator (Fig. 11) we see that for very fine temporal resolution (i.e. very small time steps) the approximation tapers around 10 −14 . This is caused by the finite precision of the computer and hence a natural limit for the accuracy of the entropy conservation in our scheme.

Wind Tunnel with a
Step (2D). The wind tunnel that contains a step was first described in [27], who used it to compare several hydrodynamical schemes. Woodward and Colella [109] later reused it to compare more advanced methods including their piece-wise parabolic method (PPM) scheme. It exercises the scheme's ability to handle strong unsteady shock interactions in multiple dimensions. Furthermore, it can be used to verify that a code solves problems with non-trivial boundaries correctly. Note, that we include this test to show that our MHD scheme performs well in the purely hydrodynamical limit (i.e. vanishing magnetic fields). The flow is unsteady and exhibits multiple shock reflections and interactions between different discontinuities.
We use this test case to demonstrate the increased robustness gained from the somewhat unusual average states that we use to evaluate the discrete dissipation operator. If we use a naive evaluation of the discrete operator, for example the simple arithmetic mean for each component, . Between x = 0.6 and x = 3.0 we insert a step with height ∆y = 0.2. We initialize a uniform Mach 3 flow using = 1.4, p = 1.0, and u = 3.0, with a ratio of specific heats of γ = 1.4. The front of the step is treated as a reflecting boundary, as are the upper and lower y-boundaries. The left-hand x-boundary is an inflow boundary constantly feeding in gas with density, pressure, and velocity like given by the initial conditions. The right-hand side uses an outflow (zero-gradient) boundary. Since the exit velocity is always supersonic the exit boundary conditions have no effect on the flow. Immediately, a shock forms in front of the step and curves around the corner. The general shape and position of the shocks are accurately represented. The shocks are thin, however, numerical instabilities are not observed. The corner of the step evolves into the center of a rarefaction fan and into a singular point of the flow. In contrast to previous works (e.g. [109]), where special boundary conditions near the corner of the step had to be applied, no special boundary condition is applied in our work. The shock curves around the corner of the step and expand to the right (downstream). It grows in size until it strikes the upper reflecting boundary just after about t = 0.55 and is reflected back to the bottom where it again reflected at about t = 1.1. Throughout the computation (until t ≈ 12.0), the flow is unsteady, exhibits multiple shock reflections and interactions between different types of discontinuities.

5.3.
Standard numerical test problems. The properties of the numerical approximation are now verified and the EC/ES schemes are implemented correctly into FLASH. We next demonstrate the utility, robustness, and accuracy of the new solver by computing the solution to several well-known MHD test problems.

5.3.1.
Orzsag-Tang vortex (2D). The Orszag-Tang vortex problem [79] is a two-dimensional, spatially periodic problem well suited for studies of MHD turbulence. Thus, it has become a classical test for numerical MHD schemes. It includes dissipation of kinetic and magnetic energy, magnetic reconnection, the formation of high-density jets, dynamic alignment and the emergence and manifestation of small-scale structures [25]. The Orszag-Tang MHD vortex problem starts from non-random, smooth initial data as shown in Table 5. As the flow evolves it gradually becomes increasingly complex, forming intermediate shocks. Thus, this problem demonstrates the transition from initially smooth data to compressible, supersonic MHD turbulence. The initial data is chosen such that the root mean square values of the velocity and the magnetic fields as well as the initial Mach number are all one. all: periodic Adaptive refinement on density, magnetic field Simulation end time t max = 0.5 Adiabatic index γ = 5/3 Table 5. Initial conditions and runtime parameters: Orszag-Tang MHD vortex [25]. Fig. 13 shows the density of the plasma at t = 0.5 given the initial conditions listed in Table 5. We superimpose the magnetic field using the technique of line integral convolution (LIC) as this method, in contrast to an ordinary streamline plot, preserves all details of the vector field geometry [11]. As the solution evolves in time, the initial vortex splits into two vortices. Sharp gradients accumulate and the vortex pattern becomes increasingly complex due to highly non-linear interactions between multiple intermediate shock waves traveling at different speeds. The result compares very well with results given in the literature, e.g. [1,22,70].

MHD rotor test (2D).
The MHD rotor problem [2] describes a rapidly spinning dense cylinder embedded in a magnetized, homogeneous medium at rest. Due to centrifugal forces, the dense cylinder is in non-equilibrium. As the rotor spins with the given initial rotating velocity, the initially uniform magnetic field will wind up the rotor. The wrapping of the rotor by the magnetic field leads to strong torsional Alfvén waves launched into the ambient fluid. The initial conditions are listed in Table 6.
At t = 0.15 the rotor is, up to a certain radial distance, still in uniform rotation. Beyond this radius, the rotor has exchanged momentum with its environment and decelerated. In Fig. 14 we present the magnetic field superimposed on the fluid density, . It is clearly seen that the magnetic field basically maintains its initial shape outside of the region of influence of the Alfvén waves. Inside, the magnetic field is refracted by the MHD discontinuities. In the center, however, we can clearly see that the magnetic field remains unchanged and only gets rotated together with the flow (the cylinder has rotated almost 180 • ). The shown snapshot is taken shortly before the torsional Alfvén waves leave the computational domain and demonstrates that our solver is able to resolve torsional Alfvén waves.

MHD explosion wave (2D/3D
). The two-dimensional version of the MHD blast wave problem was studied by [2]. We use an extended three-dimensional version to demonstrate the robustness of our scheme in simulations involving regimes with low thermal pressures and high kinetic as well as magnetic energies in three dimensions. This test problem leads to the onset of strong MHD discontinuities, relevant to astrophysical phenomena where magnetic fields can have strong dynamical effects. It describes an initially circular pressure pulse. We choose here a relative magnitude of 10 4 for comparison with [2]. The initial conditions used are listed in Table 7.  Table 7. Initial conditions and runtime parameters: MHD blast wave test [25].
The MHD explosion wave is highly anisotropic and the explosion bubble is strongly distorted as the displacement of the gas in the transverse y and z−direction is inhibited due to the magnetic field (set in the x−direction). The out-going blast wave shows no grid alignment effects or other numerical defects as seen for various schemes, e.g., in [25].

Concluding remarks and outlook
In this work, we collected background information and demonstrated the inextricable link between physics and applied and computational mathematics. In particular, we examined the fundamental laws of thermodynamics and built partial differential equations (PDEs) to mathematically model the evolution and propagation of wave phenomena. Of utmost importance was the concept of entropy conservation (or stability) depending on the smoothness of the solution. This is because it is the entropy (introduced by the second law of thermodynamics) that separates  possible state configurations of systems from impossible states. However, entropy conservation is almost always an auxiliary conservation law of a given PDE system.
On the continuous level, the link between the primary conserved quantities of a PDE and the conservation of entropy has been well understood since the late 1970s. However, it wasn't until the landmark work of Tadmor in the 1980s that the design of numerical algorithms that could discretely remain entropy conservative was known. It was within this framework, built from a finite volume type discretization, that we described the construction of a baseline entropy conservative numerical flux function for the ideal MHD equations.
A very subtle part of the continuous as well as discrete entropy analysis is the treatment of the divergence-free constraint on the magnetic fields, ∇ · B = 0. Entropy conservation was shown in the continuous case by Godunov. To do so Godunov augmented the ideal MHD system to include a source term proportional to the divergence-free condition (essentially adding zero in a clever way). However, on the discrete level, this addition of zero strategy no longer works. This is due to inherent numerical inaccuracies in the divergence-free constraint for a discretization. Even if the initial flow is divergence-free, as the numerical approximation of the flow evolves inaccuracies can build and drive instabilities. We showed that a particular discretization of the Janhunen source term (similar to the Godunov or Powell source term but remains conservative in the hydrodynamic variables) was enough to guarantee a consistent, entropy conservative numerical flux for ideal MHD.
The PDEs considered are of a hyperbolic nature. This is problematic as discontinuities may develop, e.g., shock waves, in the solution regardless of the initial solution smoothness. Physically, it is known that entropy must be dissipated in the presence of shocks. Therefore, we presented the details of a particular dissipation term to be added to the baseline entropy conservative flux to create an entropy stable approximation. That is, the discrete entropy in the numerical approximation will always possess the correct sign. Furthermore, we demonstrated at which mean state the dissipation term should be evaluated to avoid unphysical transfers in the conserved quantities, e.g., mass.
The final section of this work applied the entropy stable algorithm to five two-and threedimensional numerical examples. The high Mach number flow in a wind tunnel was a pure Euler calculation, i.e., no magnetic fields. Through the battery of numerical examples considered we demonstrated that the numerical algorithm is able to model the fluid dynamics of a wide range of flows, e.g., low or high Mach number.
Two main avenues guide our future research: 1. There are many approaches to address the errors in the divergence-free constraint. The three most popular are projection methods [7,71], constrained transport [2,28], and hyperbolic divergence cleaning [23,101]. An extensive comparison of projection and constrained transport type methods is provided by Tóth [100]. The most computationally efficient of these choices is hyperbolic divergence cleaning. The idea is to introduce an additional variable proportional to the divergence error. This additional variable is governed by a simple advection equation and coupled into the flux of the magnetic field equations such that errors in the divergence are moved away from where they are generated in the approximation. However, at this time, it is unclear if hyperbolic divergence cleaning can lay "on top" of an entropy stable algorithm and remain entropy stable.
2. High-order spatial discretizations for the approximation of the ideal MHD equations that retain entropy stability are also of interest. This is because high-order spatial methods are able to capture complex flow features with fewer degrees of freedom, e.g., [41]. The results presented in this work used a simple linear reconstruction with a minmod limiter to increase spatial accuracy to second order. There are higher order reconstructions available in the literature, e.g., [4,86,94,102] but it might be that they fail to preserve the entropy stability we have worked so hard to obtain for the low order approximation. This issue has been previously addressed for general systems of conservation laws and there are arbitrarily high-order entropy stable schemes available [32,34,35,66]. However, all the high-order extensions referenced assume that there was no source term. We demonstrated the importance of the source term discretization for the ideal MHD equations and its crucial role to guarantee discrete entropy conservation. Therefore, an interesting question is to determine how the source term can be discretized at high-order while maintaining entropy consistency.