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:
- 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. - We need the jumps between right and left values, omega_i. Method: calc_omega.
- We need the right and left fluxes based on (1). Method: calc_f
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.
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
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
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; }
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
The Runge-Kutta method uses the slope at the nth time step to estimate q at , and based on the value the "midpoint" slope at the level is obtained. This slope is used in place of , and leads to a higher order of accuracy
In this case,
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; }
for (int i = Ng; i < Nx - Ng; i++) { qnp1[i] = qn[i] - delta_t/dx*(FF[i] - FF[i-1]); }
评论