Formulations of hyperbolic PDEs and Godunov's scheme
A wide variety of PDEs have the form,
,
where U is the set of variables we are interested in. This is known as the quasilinear formulation of PDEs. The governing equations of GRMHD can also be cast in a similar form (cf. Anile (1990)). Note that A(U) is an matrix that consists of functions of U where n is the length of U. We say that the PDE is hyperbolic if A is diagonalizable and has real entries (which are its eigenvalues,
). Of course this also means it has n eigenvectors corresponding to those eigenvalues and these are important as they physically mean something but we'll get to that in a moment.
Hyperbolicity is important because this ensures the problem is well-posed. This means that given some initial data (at t=0), there exists a unique solution for t>0 and the solution's behavior changes continuously with the initial conditions ie. small fluctuations in initial conditions lead to small fluctuations in the solution. The equations of GRMHD are hyperbolic (see Anile (1990)) but not strictly hyperbolic, ie. their eigenvalues are not distinct.
As noted above, the equations of GRMHD above are diagonalizable and can be written in a form,
,
where Li are the corresponding eigenvectors. One can then diagonalize A, , where L is the matrix of eigenvectors. Applying R-1 to the governing equations we obtain,
,
where w=R-1U is called the characteristic vector and the equation above is called a characteristic equation. By doing so we have decoupled the equations and effectively obtained 8 independent advection equations. These are wave equations and each equation has a characteristic wave velocity. The eigenvalues define characteristic curves , along which the waves propagate. Anile and Pennisi (1987) and Anile (1990) studied the eigenstructure of these equations in the 10-dimensional space (p, s, uµ, bµ) where 's' is the specific entropy. They find seven characteristic, distinct, wave velocities:
(entropy wave),
(Alfvén waves),
(fast magnetosonic waves) and
(slow magnetosonic waves). The first two have a closed form expression while the latter two involve solving a quartic expression and are computed numerically. The two fast magnetosonic velocities are used to compute HLL fluxes.
mhd_vchar in phys.c is responsible for computing these velocities.
In principle, this should solve problem as one can always compute U(x,t) from w(x,t) (which is obtained from solving the advection equation, given some initial data) by inverting the relation between w and U. However, for nonlinear hyperbolic systems, as it is the case in GRMHD, the characteristic curves intersect and the solution needn't be single-valued. We instead use the conservative formulation of hyperbolic PDEs.
By expressing the matrix A(U) as a gradient of a function F(U) w.r.t. U, hyperbolic PDEs can be recast into a conservative form,
,
where F is known as the flux vector and S is the source vector. In such a form, variables U are called conserved variables. The GRMHD equations when recast into a conservative form appear as,
This reveals the list of conserved variables in our system of equations which we discuss here along with what's known as primitive variables.
Godunov's scheme is a finite volume method for conservative equations where the fluxes are computed at zone faces by solving a local Riemann problem at the zone interface. A 'Riemann problem' consists of a set of hyperbolic PDEs expressed in a conservative form with some initial value data that is piecewise constant and has a single discontinuity at the spatial location that is of interest to us. Let us break down all the new terms introduced in this paragraph and physically/numerically understand what they mean and how they all fit together.
The finite volume method discretizes the conserved variables such that the grid/cell volume-averaged value is stored at the grid center,
,
where for simplicity we have considered a 1D problem (along x) with grid points indexed by i. n denotes the time step and ∆x = xi+1/2 - xi-1/2. The equation at the top of this section can then be expressed as,
where is the average numerical flux over ∆t = tn+1 - tn evaluated at the zone interfaces
and for Godunov's scheme this depends on the values of the conserved variables in the grid zones adjoining the zone face, for eg:
. The values of the conserved variables at the two adjoining zones are piecewise discontinuous, ie. the initial data has a jump from the left to the right face and we are trying to solve a set of hyperbolic PDEs. This is precisely what is known as a Riemann problem. Since the fluxes depend on the values of the conserved variables only in the adjoining zone, it is a local Riemann problem. This page discusses how one goes about using the zone-centered values to compute the flux at the interface.