Simplicity is the ultimate sophistication.

― Clare Boothe Luce

urlWhen one is solving problems involving a flow of some sort, conservation principles are quite attractive since these principles follow nature’s “true” laws (true to the extent we know things are conserved!). With flows involving shocks and discontinuities, the conservation brings even greater benefits as the Lax-Wendroff theorem demonstrates ( In a nutshell you have guarantees about the solution through the use of conservation form that are far weaker without it. A particular set of variables is the obvious variables because they arise naturally in conservation form. For fluid flow these are density, momentum and total energy. The most seemingly straightforward thing to do is use these same variables to discretize the equations. This is generally a bad choice and should be avoided unless one does not care about the quality of results.

While straightforward and obvious, the choice of using conserved variables is almost always a poor one, and far better results can be achieved through the use of primitive variables for most of the discretization and approximation work. This is even true if one is using characteristic variables (which usually imply some sort of entirely one-dimensional character). The primitive variables have simple and intuitive meaning physically, and often equate directly to what can be observed in nature (conservedcsd240333fig7variables don’t). The beauty of primitive variables is that they trivially generalize to multiple dimensions in ways that characteristic variables do not. The other advantages are equally clear specifically the ability to extend the physics of the problem in a natural and simple manner. This sort of extension usually causes the characteristic approach to either collapse or at least become increasingly unwieldy. A key aspect to keep in mind at all times is that one returns to the conservation variables for the final approximation and update of the equations. Keeping the conservation form for the accounting of the complete solution is essential.

To keep the bulk of the discussion simple, I will focus on the Euler equations of fluid dynamics. These equations describe the conservations of mass, \rho_t + m_x = 0, momentum, m_t + (m^2/\rho + p)_x = 0 and total energy, E_t + \left[m/\rho(E + p) \right]_x = 0 in one dimension. Even in this very simple setting the primitive variables are immensely useful as demonstrated by HT Huynh, in another of his massively under-appreciated papers. In this paper he masterfully covers the whole of the techniques and utility of primitive variables. Arguably, the use of primitive variables went mainstream with the papers of Colella and Woodward. In spite of the broad appreciation of that paper, the use of primitive variables in work is still more a niche than common practice. The benefits become manifestly obvious whether one is analyzing the equations (which is equivalent to the more complex variable set!), or discretizing the solutions.

Study the past if you would define the future.

― Confucius

ClimateModelnestingThe use of the “primitive variables” came from a number of different directions. Perhaps the earliest use of the term “primitive” came from meteorology in terms of the work of Bjerknes (1921) whose primitive equations formed the basis of early work in computing weather in an effort led by Jules Charney (1955). Another field to use this concept is the solution of incompressible flows. The primitive variables are the velocities and pressure, which is distinguished from the vorticity-streamfunction approach (Roache 1972). In two dimensions the vorticity-streamfunction solution is more efficient, but lacks simple connection to measurable quantities. The same sort of notion separates the conserved variables from the primitive variables in compressible flow. The use of primitive variables as an effective approach computationally may have begun in the computational physics work at Livermore in the 1970’s (see e.g., Debar). The connection of the primitive variables to classical analysis of compressible flows and simple physical interpretation also plays a role.

What are the primitive variables? The basic conserved variables form compressible fluid flow are density, \rho, momentum, m=\rho u, and total energy, E = \rho e + \frac{1}{2} \rho u^2. Here the velocity is u and the internal energy is e. One also has the equation of state p=P(\rho,e) as the constitutive relation. Let’s take the Euler equations and rewrite them using the primitive variables, the conservations of mass, \rho_t + (\rho u)_x = 0, momentum, (\rho u)_t + (\rho u^2 + p)_x = 0 and total energy, \left[\rho (e + \frac{1}{2}u^2)\right]_t + \left[u(\left(rho (e + \frac{1}{2}u^2)+ p\right) \right]_x = 0. Except for the energy equation, the expressions are simpler to work with, but this is the veritable tip of the proverbial iceberg.

What are the equations for the primitive variables? The primitive variables can be expressed and evolved using simpler equations, which are primarily evolution equations dependent on differentiability, which must be present for any sort of accuracy to be in play anyway. The mass equation is the same although one might expand the derivative, \rho_t + u \rho_x + \rho u_x = 0. The momentum equation is replaced by an equation of motion, u_t + u u_x + \frac{1}{\rho} p_x = 0. The energy equation could be replaced with a pressure equation, p_t + u p_x + \gamma p u_x = 0 (\gamma is the generalized isentropic derivative \partial_\rho p|_S) or an internal energy equation, \rho e_t + \rho u e_x + p u_x = 0. One can use either energy representation to good measure, or better yet, use both and avoid having to evaluate the equation of state. Moreover if one wants you can evaluate the difference between the pressure from the evolution equation and the state relation as an error measure.

How does on convert to the primitive variables, and convert back to the conserved variables? If one is interested in analysis of the conservative equations, then one linearizes the equations about a point, U_t + \left(F(U)\right)_x = 0 \rightarrow U_t + \partial_U F(U) U_x = 0 where U is the vector of conserved varibles, and F(U) is the flux function. The matrix A_c = partial_U F(U) is the flux Jacobian. One does an eigenvalue decomposition, $ to analyze the equations. From this decomposition, A_c = R_c \Lambda L_c, one can get the eigenvalues, \Lambda, and the characteristic variables, L_c \Delta U. The analysis is difficult and non-intuitive with the conserved variables.

Here we get to the cool part of this whole thing, there is a much easier and more intuitive path through the primitive variables. One can get a matrix representation of the primitive variables which I’ll call V in vector form, V_t + A_p V_x = 0. One can get the terms in A_p easily from the differential forms, and recognizing that \gamma p = \rho c^2, with c being the speed of sound, the eigen-analysis is so simple that it can be done by hand (and it’s a total piece of cake for Mathematica). Using similar notation as the conserved form, A_p = R_p \Lambda L_p. The first thing to note is that \Lambda is exactly the same, i.e., the eigenvalues are identical. One then gets a result for the characteristics, L_p \Delta V that matches the textbooks, and that L_p \Delta V = L_c \Delta U. All the differences in the transformation are bound up in the right eigenvectors R_c and R_p, and the ease of physical insight provided by the analysis.

24-Figure17-1Now we can elucidate how to move between these two forms, and even use the primitive variables for the analysis of the conserved form directly. Using Huynh’s paper as a guide and repeating the main results one defines a matrix of partial derivatives of the conserved variables, U with respect to the primitive variables, V, M= \partial_V U. This matrix then can be inverted into M^{-1} and we then may define an identity, A_c = M A_p M^{-1}, which might allow the conserved eigen-analysis to be executed in terms of the more convenient primitive variables. The eigenvalues of A_c and A_p are the same. We can get the left and right eigenvectors through L_c = L_p M^{-1} and R_c = M R_p. All of this follows the simple application of the chain rule to the linearized versions of the governing equations.

The primitive variable idea can be extended in a variety of nifty and useful ways. One can augment the variable set in ways that can yield some extra efficiency to the solution by avoiding extra evaluations of the constitutive (or state) relations. This would most classically involve using both a pressure and energy equation in the system. Miller and Puckett provide a nice example of this technique in practice, building upon the work of Colella, Glaz and Ferguson where expensive equation of state evaluations are avoided. One must note that the system of equations being used to discretize the system is carrying redundant information that may have utility beyond efficiency.

One can go beyond this to add variables to the system of equations that are redundant, but carry information implicit in their approximation that may be useful in solving equations. One might add an equation for the specific volume of the fluid to compare with density. Similar things could be done with kinetic energy, vorticity, or entropy. In each case the redunency might be used to discover or estimate error or smoothness of the underlying solution and perhaps adapt the solution method on the basis of this information.

Using the primitive variables for discretization is almost as good as using characteristic variables in terms of solution fidelity. Generally if you can get away with 1-D ideas, the characteristic variables are unambiguously the best. The primitive variables are almost as good. The key is to use a local transformation to the primitive variables for the work of discretization even when your bookkeeping is all in conserved variables. Even if you are doing characteristic variables, the construction and use of them is enabled by primitive variables. The resulting expressions for the characteristics are simpler in primitive variables. Perhaps almost as important the expressions for the variables themselves are far more intuitively expressed in primitive variables.

A real source of power of the primitive variables comes when you extend past the simpler case of the Euler equations to things like magneto-hydrodynamics (MHD i.e., compressible magnetic fluids). Doing discretization of the MHD with conserved variables is a severe challenge and analysis of their mathematical characteristic structure can be a decent into utter madness. Doing the work in these more complex systems using the primitive variables is extremely advantageous. It is an approach that is far too often left out and the quality and fidelity of numerical methods suffers as a result.

Any intelligent fool can make things bigger, more complex, and more violent. It takes a touch of genius — and a lot of courage to move in the opposite direction.

― Ernst F. Schumacher

Lax, Peter, and Burton Wendroff. “Systems of conservation laws.”Communications on Pure and Applied mathematics 13, no. 2 (1960): 217-237.

Huynh, Hung T. “Accurate upwind methods for the Euler equations.” SIAM Journal on Numerical Analysis 32, no. 5 (1995): 1565-1619.

Colella, Phillip, and Paul R. Woodward. “The piecewise parabolic method (PPM) for gas-dynamical simulations.” Journal of computational physics 54, no. 1 (1984): 174-201.

Woodward, Paul, and Phillip Colella. “The numerical simulation of two-dimensional fluid flow with strong shocks.” Journal of computational physics54, no. 1 (1984): 115-173.

Van Leer, Bram. “Upwind and high-resolution methods for compressible flow: From donor cell to residual-distribution schemes.” Communications in Computational Physics 1, no. 192-206 (2006): 138.

Bjerknes, V. “The Meteorology of the Temperate Zone and the General Atmospheric CIRCULATION. 1.” Monthly Weather Review 49, no. 1 (1921): 1-3.

Charney, J. “The use of the primitive equations of motion in numerical prediction.” Tellus 7, no. 1 (1955): 22-26.

Roache, Patrick J. Computational fluid dynamics. Hermosa publishers, 1972.

DeBar, R. B. Method in two-D Eulerian hydrodynamics. No. UCID-19683. Lawrence Livermore National Lab., CA (USA), 1974.

Miller, Gregory Hale, and Elbridge Gerry Puckett. “A high-order Godunov method for multiple condensed phases.” Journal of Computational Physics128, no. 1 (1996): 134-164.

Colella, P., H. M. Glaz, and R. E. Ferguson. “Multifluid algorithms for Eulerian finite difference methods.” preprint (1996).