Initial conditions

We would like to create a "Burger object" which allows us the Burger equation on a domain from xmin to xmax, dividing this domain into Nx cells and updating with time step size dt. We will also have a choice of two types of "automatic" initial data; a Gaussian profile centered about the center of the domain, and a discontinuous initial profile with value qL on the left qR on the right of some point. Though it is possible to input other initial profiles, these two allow us to see how Burger's equation affects the evolution of discontinuities.

The constructor creates the Burger object given the number of cells Nx, time step size dt, and a boolean value "smooth" which determines whether the initial profile will be Gaussian or discontinuous when the init_fluid() method is invoked. The domain (xmin, xmax) is also automatically initialized to (0, 5), so to solve over a different domain it is only necessary to change xmin and xmax definitions here. Using xmin, xmax, and Nx, an array x[] of length Nx is created to hold the positions of the cell centers. The array q[] (also of length Nx) is also created to hold the average value of the function over each cell, but until init_fluid() has been called, every element of q[] is by default set to 0. (The number of ghost cells Ng is also automatically set to 2).

The init_fluid() method fills q[] with the relevant initial data, Gaussian if smooth = true and discontinuous otherwise. By default the maximum of the Gaussian and initial discontinuity are centered over the domain (with qL = 1.0 and qR = 0.1 in the latter case) but it is simple enough to modify this method if different initial profiles are desired.

Boundary conditions

It is clear that the Roe flux method cannot be used for the two cells adjacent to each boundary of the domain, since the calculation of the numerical flux in this scheme requires knowledge of values in two cells to the immediate left and right of the cell being updated. These cells are called "ghost cells", and are updated differently from cells in the interior of the domain depending on the desired boundary conditions. Here we use "outflow" conditions, which assume no waves travelling in from outside the domain will influence the solution within it. We can then choose to update the ghost cells by extrapolating from the interior solution, assuming continuity at the boundary. The most obvious way to do this is to copy the value from the nearest non-ghost cell (constant or zero-order extrapolation):

//constant extrapolation
for (int i = 0; i < Ng; i++) {
    qnp1[i] = qnp1[Ng];
    qnp1[Nx - i - 1] = qnp1[Nx - Ng - 1];

An equally intuitive choice would be to extrapolate linearly from the two closest non-ghost cells, or first-order extrapolation

//first order extrapolation
for (>int i = 0; i < Ng; i++) {
    qnp1[i] = qnp1[Ng]  - (Ng - i)*(qnp1[Ng + 1] - qnp1[Ng]);
    qnp1[Nx - i - 1] = qnp1[Nx - Ng - 1] + (Ng - i)*(qnp1[Nx - Ng - 1] - qnp1[Nx - Ng - 2]);

This seems as though it would always be more accurate than zero-order extrapolation, but in practice it may be lead to stability problems. In either case, we expect that this extrapolation will have a lower order than the scheme used in the interior, and therefore that the order of convergence at the boundaries will be less than around smooth interior sections. At later times, inaccuracy from the boundaries may "leak in" and corrupt the solution even at the interior. This is particularly evident in the case of the ultrarelativistic Gaussian evolution, though in this instance using first-order extrapolation appears to improve the situation somewhat.