Implementation

The quantities we then need to calculate in order to obtain the Roe numerical flux at the boundary, x_iph (i.e. x_iplushalf) are as follows:

  1. We need the values of the conservative variable
    q[x]
    to the left and right of the boundary. Methods: recos_qL, recos_qR, with the help of dqL_1, dqR_1, etc.
  2. The idea is to estimate the value of q at the cell boundary xiph based on the averages in the adjacent cells. The most intuitive way to do this is to estimate the values of the slope on either side and extrapolate, so as follows

    $\displaystyle q^L_{i+\frac{1}{2}} = q_i + \frac{q_i - q_{i-1}}{x_i - x_{i-1}}
$

    $\displaystyle q^R_{i+\frac{1}{2}} = q_{i+1} - \frac{q_{i+1} - q_i}{x_{i+1} - x_i}
$

    However, this can lead to problems when there are discontinuities. Since the slopes around discontinuities are very steep, estimated slopes can sometimes overshoot true cell averages. This can lead to too much total variation, defined as

    $\displaystyle \sum_i \vert q_{i+1} - q_i\vert
$

    A high total variation indicates large amounts of oscillation, which is not desirable in a solution and indicates instability. Various "slope-limiter" reconstruction methods aim to limit the variations so the solution method is total variation diminishing, or TVD. The idea is to restrict the ratio between adjacent jumps in cell averages where such jumps seem to indicate oscillatory behavior, while (depending on the method) steepening this jumps in other places to sharpen resolution around discontinuities that are actually part of the solution. This way, unwanted variation is reduced without sacrificing too much accuracy. There are options for four different types of slope limiters this code (called minmod, superbee, MC, and van Leer). The methods dqL_1, dqR_1, dqL_2, dqR_2, etc. calculate limited jumps for use in the recos_qL and recos_qR methods. By default the minmod limiter (dq1) is used.
    private double dqL_1(double dq_imh, double dq_iph) {
    		return minmod(dq_iph, dq_imh);
    	}
    
    //minmod function
    private double minmod(double a, double b){
    	if (a*b < 0) return 0;
    	else if(Math.abs(a) < Math.abs(b)) return a;
    	else return b;
    	}
    
    In practical terms, this means that if the slops on either side of the boundary have different signs, which could indicate spurious oscillation, the value will be reconstructed using a slope of zero. Similarly, if the slopes on either side have the same sign, the less steep one is always chosen. This clearly reduces the total variation of the solution, but it may slightly smooth out discontinuities even when they are desired in the solution.
  3. We need the jumps between right and left values, omega_i. Method: calc_omega.
  4. We need the right and left fluxes based on (1). Method: calc_f
  5. We need the Jacobian matrix at x_iph, along with its eigenvalues and eigenvectors lambda and eta. Methods: calc_A (for the Jacobian), calc_eta, calc_lambda.

    In the case of Burger's equation, the conservative variable and the flux functions are both 1D, so the Jacobian and single eigenvalue lambda are just the scalar q, and the eigenvector is just 1. This makes the above methods trivial, but they would be non-trivial in a multidimensional situation.

Time-stepping method

This program uses a 2nd order Runge-Kutta time stepping method. This involves estimating the derivative twice at each time step, rather than once as in the first-order Euler method

$\displaystyle \textbf q^{n+1} = \textbf q^n + \textbf v(\textbf q^n, t_n) \triangle t + O(t^2)$$\displaystyle \mbox { where } \textbf v(\textbf q, t) = \frac{\partial \textbf q}{\partial t} \Big \vert _t
$

The Runge-Kutta method uses the slope at the nth time step to estimate q at $ t_{n+\frac{1}{2}}$ , and based on the value $ q^{n+\frac{1}{2}}$ the "midpoint" slope at the $ \frac{\partial \textbf q}{\partial t} \Big \vert _{t_{n+\frac{1}{2}}}$ level is obtained. This slope is used in place of $ \frac{\partial \textbf q}{\partial t} \Big \vert _{t_n}$ , and leads to a higher order of accuracy

$\displaystyle \tilde{\textbf q}^{n+\frac{1}{2}} = \textbf q^n +\textbf v(\textbf q^n, t_n) \frac{\triangle t}{2}
$

$\displaystyle \textbf q^{n+1} = \textbf q^n + \textbf v(\tilde{\textbf q}^{n + \frac{1}{2}}, t_{n+\frac{1}{2}}) \triangle t + O(t^3)
$

In this case,

$\displaystyle v(\textbf q, t) = \textbf F(
$

This is implemented as follows in the step() method.

//moving forward one time step
private double[] step(double[] q) {

	//half advanced time step to estimate v at the midpoint
	double[] FF = calc_flux(q, x, Nx, Ng);
	double[] q_nph = update_q(q, x, FF, delta_t/2, Nx, Ng);

	//full advanced time step
	FF = calc_flux(q_nph, x, Nx, Ng);
	double[] q_np1 = update_q(q_nph, x, FF, delta_t, Nx, Ng);

	return q_np1;

	}
where in each case updating uses the derivative v based on the fluxes as described above (note that FF[i] contains the flux at the boundary x_iph).
for (int i = Ng; i < Nx - Ng; i++) {
	qnp1[i] = qn[i] - delta_t/dx*(FF[i] - FF[i-1]);
}