Some Background on Finite Volume Methods

We are generally interested in solving PDE's of the form

$\displaystyle \frac{\partial \textbf u}{\partial t} + \frac{\partial \textbf f(\textbf u)}{\partial \textbf x} = 0

For the moment, let's focus our attention even further, on one of the simplest PDE's of that form, known as Burger's equation (inviscid form).

$\displaystyle \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0
Though this class of equations may appear relatively straightforward, it turns out that many such PDE's lead to the formation of discontinuities, or shocks, even from smooth initial data. In order to solve PDE's whose solutions contain shocks, we must turn to methods which do not assume the continuity of the solution, unlike the more intuitive finite difference methods. One way in which we can make sense of such solutions is to think of the average value of the solution over a given cell, rather than the value at specific grid points. Let the cell $C_{i}^{n+\frac{1}{2}}$ describe the spatial region between $x_{i-\frac{1}{2}}$ and $x_{i+\frac{1}{2}}$ and the temporal region between $t_i$ and $t_{i+1}$. The average value of the solution over this cell is

\begin{displaymath}\textbf u_i^n = \frac{1}{\triangle x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} u(x,t_n)\ \,dx\ \end{displaymath}

and let the numerical flux F be given by

\begin{displaymath}\textbf F^{n+\frac{1}{2}}_i =
\frac{1}{\triangle t} \int_{t_n}^{t_{n+1}} f(u(x_i,t))\ \,dt\ \end{displaymath}

Integrating the conservation form of Burger's eqn. and dividing by total volume gives,

\begin{displaymath}\frac{1}{\triangle x \triangle t}\int_{t_n}^{t_{n+1}} \int_{x...
...\frac{1}{2}}} \frac{\partial f(u)}{\partial x} \,dx\\ \,dt\ = 0\end{displaymath}

\begin{displaymath}\frac{1}{\triangle x \triangle t} \int_{x_{i-\frac{1}{2}}}^{x...
...}} [f(u_{i + \frac{1}{2}}) - f(u_{i - \frac{1}{2}})]\ \,dt\ = 0\end{displaymath}

\begin{displaymath}\frac{\textbf u_i^{n+1} - \textbf u_i^n}{\triangle t} +
...} - \textbf F^{n+\frac{1}{2}}_{i-\frac{1}{2}}}{\triangle x} = 0\end{displaymath}

If we can find a way to estimate the numerical fluxes $\textbf F^{n+\frac{1}{2}}_{i+\frac{1}{2}}$ and $\textbf F^{n+\frac{1}{2}}_{i-\frac{1}{2}}$ using information from the current time step, we will have a method for advancing the solution in time. Specifically

\begin{displaymath}\textbf u_i^{n+1} = \textbf u_i^n -
\frac{\triangle t}{\tri...
...}}_{i+\frac{1}{2}} - \textbf F^{n+\frac{1}{2}}_{i-\frac{1}{2}})\end{displaymath}

To see how this can be implemented in practice, go to the next section

Roe Solver

One technique for estimating these numerical fluxes, developed by Roe, involves linearizing the system by evaluating the Jacobian matrix $\textbf A = \frac{\partial \textbf f}{\partial \textbf u}$ at $(x_i, t_n)$ and approximating $\textbf f(\textbf u) = \textbf A \textbf u$. The original equation then becomes $\frac{\partial \textbf u}{\partial t} + \textbf A \frac{\partial \textbf u}{\partial x} = 0$. If the basis of eigenvectors $r_\alpha$ of A is chosen, the system is uncoupled and is reduced to a series of 1D equations of the form

\begin{displaymath}\frac {\partial u_{(\alpha)}}{\partial t} + \lambda_\alpha \frac {\partial u_{(\alpha)}}{\partial x} = 0\end{displaymath}

where $\lambda_\alpha$ is the eigenvalue corresponding to the eigenvector $r_\alpha$. Given initial data

\begin{displaymath}u_{(\alpha)}(x, 0) = \left\{
u_L & \quad ...
u_R & \quad x > x_{i+\frac{1}{2}}
\end{array} \right. \end{displaymath}

the value at $x_{i+\frac{1}{2}}$ at subsequent times is $q_L$ if $\lambda_\alpha > 0$ and $q_R$ if $\lambda_\alpha < 0$. The $\alpha$ component of the numerical flux at $x_{i+\frac{1}{2}}$ is then given by

\begin{displaymath}F_{i + \frac{1}{2}, \alpha} = \left\{
...f_\alpha (u_R) & \quad \lambda_\alpha < 0
\end{array} \right. \end{displaymath}

The two results can be combined in the following expression

\begin{displaymath}F_{i + \frac{1}{2}, \alpha} = \frac{1}{2}(f_\alpha (u_L) + f_\alpha (u_R) - \vert\lambda_\alpha\vert \omega_\alpha)\end{displaymath}

Summing over all $\alpha$, the final result is

\begin{displaymath}\textbf F_{i+\frac{1}{2}} = \frac{1}{2}(\textbf f_\alpha (\te...
... \sum \vert\lambda_\alpha\vert \omega_\alpha \textbf r_\alpha) \end{displaymath}