For scalar equations the front tracking algorithm is as follows:
The algorithm is the same for systems hyperbolic conservation laws, except that the flux function is not approximated by a piecewise linear function. Moreover, boundary conditions are easily included into the algorithm, and the algorithm then reads
- Make initial approximation: f piecewise linear, u0 piecewise constant.
- Solve initial Riemann problems and compute possible collisions between fronts.
- T = first.collision.time.
- WHILE (T < Tend)
- remove colliding fronts;
- solve corresponding Riemann problem;
- insert new fronts;
- compute possible new collisions;
- T = next.collision.time
- END
- Advance solution forward to Tend.
- Make initial approximation.
- Solve initial Riemann problems and compute possible collisions between fronts.
- T = first.collision.time.
- WHILE (T < Tend)
- if interior collision
- remove colliding fronts;
- solve corresponding Riemann problem;
- insert new fronts;
- else (boundary collision)
- resolve collision according to boundary type;
- compute possible new collisions;
- T = next.collision.time
- END
- Advance solution forward to Tend.
The discontinuities (or fronts) are typically represented by a data object in a linked list. The data object contains the left and right state, the start position (in x and t), the place of a possible collision with the rightmost neighbour, and pointers to the rightmost neighbours. The objects are then collected in a list (X-list), from left to right, corresponding to the spatial distribution of the fronts. Collisions must also be kept in a list (T-list). We have chosen to incorporate the T-list in the X-list, that is, each data object has an additional pointer, pointing to the next collision. For a pair of colliding fronts, only the leftmost is kept in the T-list. This completely describes the two list. The front tracking algorithm operates on this lists, and removes and adds new objects for each shock collision according to the above construction. In our experience implementing the front tracking algorithm can be simplified if both the X-list and the T-list are double-linked lists, and not single-linked as described above.
A boundary can be represented in the X-list and T-list by a front object with a special tag indicating that the object is a boundary. Thus the boundary object is assigned a position, a velocity (which is zero unless the boundary is moving), a left and a right state. Computing possible collisions with the boundary can then be carried out as for an ordinary front. An example of a the data structure is shown in Figure 1
Whether one uses pointers to implement the list, or `simulate' pointers by
the use of arrays is a matter of taste and programming language. The front
tracking algorithm has been implemented in C and C++ using pointers and in
FORTRAN and Matlab using arrays.
A C++
implementation by Jan Olav Langseth of the front tracking method is
available for download.
A Riemann Solver for Scalar Problems
Solving a Riemann problem is quite simple for a general non-convex flux
function. Consider the following Riemann problem
If uL>uR
the solution is determined from the upper concave envelope of the flux
function. Otherwise, the lower convex envelope is used. Assuming that the
picewise linear flux function is represented in the arrays
u0(x) =
uR, if x> 0
f
and u, and that start is the index of the left value
and stop the index of the right value. Then the following C code
determines the upper concave/lower convex envelope
/* assume that the piecewise linear flux function is given by
the two arrays u[] and f[]. Moreover, assume for simplicity
that u_L and u_R coincide with breakpoints in the flux. */
if (u_L > u_R)
{
start =
If the left state, for instance, does not coincide with a breakpoint in the piecewise
linear flux function, we temporarily insert an new breakpoint at the left
state before entering the above algorithm.
