Godunov’s theorem says that there are no linear schemes for the linear advection equation of accuracy two or higher. Thus, for oscillation-free and high-order of accuracy, the numerical scheme have to be non-linear method.

Let’s consider the initial-boundary value problem (IBVP) for the inhomogeneous \((s(q)\neq0)\) and non-linear hyperbolic equation as follows:

\[\partial_{t}q+\partial_{x}f(q)=s(q)\]

where initial and boundary conditions are given by

\[q(x, 0)=q^{(0)}(x),\] \[q(a, t) = b_{L}(t),\] \[q(b,t)=b_{R}(t).\]

Finite Volume Framework

Entire space domain \([a,b]\) is discretized into \(M\) cells having interval \(I_i=[x_{x-1/2},x_{x+1/2}]\), which is not a point defined as finite difference framework. The control volume is established as rectangular \(V=[x_{i-1/2}, x_{i+1/2}]\times[t_{n}, t_{n+1}]\).

Therefore, exact integration of the partial difference equation (PDE) on control volume is as follows:

\[\hat{q}^{n+1}_{i}=\hat{q}^{n}_{i}-\frac{\Delta t}{\Delta x}\left[\hat{f}_{i+1/2}-\hat{f}_{i-1/2}\right] + \Delta t \hat{s}_{i}\] \[\hat{q}^{n}_{i}=\frac{1}{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}}q(x,t_n)dx,\] \[\hat{f}_{i+1/2}=\frac{1}{\Delta t}\int_{t_n}^{t_{n+1}}f\left(q(x_{i+1/2},t)\right)dt,\] \[\hat{s}_{i} =\frac{1}{\Delta t}\frac{1}{\Delta x}\int^{t_{n+1}}_{t_{n}}\int^{x_{i+1/2}}_{x_{i-1/2}}s\left(q(x,t)\right)dx dt.\]

Here, \(\hat{}\) means the exact value. As mentioned, these value is not defined as point value, but interval. Thus integral average is used for representative value. This equation shows the finite volume framework is one of the conservative form.

By removing the hat \((\hat{})\), the approximate finite volume method is given by

\[{q}^{n+1}_{i}={q}^{n}_{i}-\frac{\Delta t}{\Delta x}\left[{f}_{i+1/2}-{f}_{i-1/2}\right] + \Delta t {s}_{i},\]

where the approximated, numerical flux \((f_{i+1/2})\) and approximated, numerical source \((s_i)\) have to be defined to fully determine a finite volume method. The accuracy of this numerical scheme depends crucially on the accuracy used to approximate the integrals.

Reconstruction

The integration of the value make the constant value, and thus remove some distribution of the original value in the interval. To recovering the removed information, the reconstruction can be used. The first-degree polynomial in cell \(i\) is defined as

\[p_i(x)=q^{n}_{i}+(x-x_i)\Delta_i,\]

where \(x\in[x_{i-1/2}, x_{i+1/2}]\) and \(\Delta_{i}\) is the slope. One of the simple slope is centred slope given by

\[\Delta_{i}=\frac{q^{n}_{i+1}-q^{n}_{i-1}}{2\Delta x}.\]

With this centred slope, the reconstruction becomes a linear reconstruction as it is fixed (same pattern) for all \(i\), making the solutions oscillatory. On the other hand, by choosing \(\Delta_i\) depending on the solution \((q)\) adaptively, the reconstruction becomes non-linear.

Procedure for Numerical Flux

Boundary values of the cell \([x_{i-1/2}, x_{x+1/2}]\) are as follows:

\[q^{L}_{i}=p_i\left(x_{x-1/2}\right)=q^{n}_{i}-\frac{1}{2}\Delta x\Delta_i,\] \[q^{R}_{i}=p_i\left(x_{x+1/2}\right)=q^{n}_{i}+\frac{1}{2}\Delta x\Delta_i.\]

For time evolution, taylor expansions are applied, i.e.,

\[\hat{q}^{L}_{i}=q^{L}_{i}+\frac{1}{2}\Delta t\partial_t q^L_i,\] \[\hat{q}^{R}_{i}=q^{R}_{i}+\frac{1}{2}\Delta t\partial_t q^R_i.\]

To invert time derivative into space derivative, Cauchy-Kowalevskaya procedure is used, i.e.,

\[\partial_t q=-\partial_x f(q) + s(q),\]

where the flux derivative is approximated as

\[\partial_x f(q)\approx\frac{f(q^R_i)-f(q^L_i)}{\Delta x}\equiv\Delta f_i,\]

and the source term \(s(q)\) can be obtained by simply substitution of boundary values \(q^L_i\) and \(q^R_i\).

Final evolved boundary values become

\[\hat{q}^{L}_{i}=q^{L}_{i}-\frac{1}{2}\Delta t\Delta f_i+\frac{1}{2}\Delta ts(q^L_i),\] \[\hat{q}^{R}_{i}=q^{R}_{i}-\frac{1}{2}\Delta t\Delta f_i+\frac{1}{2}\Delta ts(q^R_i).\]

The conventional Riemann problem can be solved for solution at the local interface \(q_{i+1/2}(0)\) with Centred flux (Riemann problem-free method) or Godunov flux (Riemann solution) using two boundary values as

if \(x<0\),

\[q(x,0)=h(x)=\hat{q}^R_i\]

else if \(x>0\),

\[q(x,0)=h(x)=\hat{q}^L_{i+1},\]

It is worth noting that these two fluxes are first-order methods when without reconstruction; however, they become second-order methods with reconstruction.

Finally, the numerical flux can be obtained by

\[f_{i+1/2}=\frac{1}{\Delta t}\int^{t_{n+1}}_{t_{n}}f\left(q\left(x_{i+1/2},\frac{1}{2}\Delta t\right)\right)dt =f\left(q_{i+1/2}(0)\right).\]

Centred Flux (FORCE)

One of the centred flux is FORCE, which expressed by

\[f_{i+1/2}=\frac{1}{\Delta t}\int^{t_{n+1}}_{t_n}f^\mathrm{force}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)dt=f^\mathrm{force}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)\] \[f^\mathrm{force}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)=\frac{1}{2}\left[f^\mathrm{LF}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)+f^\mathrm{LW}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)\right]\] \[f^\mathrm{LF}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)=\frac{1}{2}\left(f\left(\hat{q}^R_i\right)+f\left(\hat{q}^L_{i+1}\right)\right)-\frac{1}{2}\frac{\Delta x}{\Delta t}\left(\hat{q}^L_{i+1}-\hat{q}^R_i\right)\] \[q^\mathrm{LW}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)=\frac{1}{2}\left(\hat{q}^R_i+\hat{q}^L_{i+1}\right)-\frac{1}{2}\frac{\Delta t}{\Delta x}\left(f\left(\hat{q}^L_{i+1}\right)-f\left(\hat{q}^R_{i}\right) \right)\] \[f^\mathrm{LW}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)=f\left(q^\mathrm{LW}_{i+1/2}\left(\hat{q}^R_i,\hat{q}^L_{i+1}\right)\right)\]

Godunov Flux (Upwind Flux)

If the conventional Riemann problem is applied for linear advection equation, the solution and numerical flux are respectively given by

if \(\lambda > 0\),

\[q(x_{i+1/2}, t)=\hat{q}^R_i\] \[f_{i+1/2}=\lambda \hat{q}^R_i\]

else if \(\lambda < 0\),

\[q(x_{i+1/2}, t)=\hat{q}^L_{i+1}\] \[f_{i+1/2}=\lambda \hat{q}^L_{i+1}\]

Procedure for Numerical Source

To approximate the volume integral of source term, mid-point rule in space and time is applied with Taylor expansion as follows:

\[q\left(x_i, t_n+\frac{1}{2}\Delta t\right)=q^n_i+\frac{1}{2}\Delta t\partial_tq\left(x_i, t_n\right)\]

To change the time derivative into space derivative, the Cauchy Kowaleswkaya is used as

\[q\left(x_i,t_n+\frac{1}{2}\Delta t\right)=q^n_i-\frac{1}{2}\Delta t\Delta f_i+\frac{1}{2}\Delta ts\left(q^n_i\right)\]

Finally, the numerical source is

\[s_i=s\left[q^n_i-\frac{1}{2}\Delta t\Delta f_i+\frac{1}{2}\Delta ts\left(q^n_i\right)\right]\]

It is worth noting that the equation for numerical flux includes the source term and the equation for numerical source also includes the flux term. Therefore, flux and source terms make a interaction each other, so have to addressed together, not separately.

Non-Linear schemes

TVD slope

With left and right polynomial slopes,

\[\Delta_{i-1/2}=\frac{q^n_i - q^n_{i-1}}{\Delta x}\] \[\Delta_{i+1/2}=\frac{q^n_{i+1} - q^n_{i}}{\Delta x}\]

The TVD slope (minmod slope) is given by

\[\Delta _i =\mathrm{minmod}\left(\Delta_{i-1/2},\Delta_{i+1/2}\right)\]

where \(\mathrm{minmod}\) function means that if signs of left and right slopes are same, the slope becomes smaller absolute value between the left and right slopes; however, if the signs are different, the value is zero.

ENO slope

The ENO means Essentially Non Oscillatory. This is the same as minmod slope except for not considering the different signs. In other words, ENO considers only smaller absolute value.