High order unconditionally stable discretization for a scalar hyperbolic PDE

Question

In order to numerically solve the following differential equation: \begin{equation} \text{Fr}\{f\} := v(k)\dfrac{\partial f(z,k)}{\partial z} - F(z) \dfrac{\partial f(z,k)}{\partial k} = -\dfrac{f-f_0}{\tau} \end{equation} I have used the finite volume method, and have discretized the left hand side to: \begin{equation} \text{Fr}'\{f\} := v(k_j)\Delta k \Bigg[ f(z_{i+},k_j) - f(z_{i-},k_j) \Bigg] - F(z_i)\Delta z \Bigg[ f(z_{i},k_{j+}) - f(z_{i},k_{j-}) \Bigg] \end{equation} for every box.

The so-called flux averaging approximation, \begin{align} & f(z_{i+},k_j) = \dfrac{f(z_i,k_j)+f(z_{i+1},k_j)}{2} \\ & f(z_{i-},k_j) = \dfrac{f(z_i,k_j)+f(z_{i-1},k_j)}{2} \end{align} will lead to instability especially if the flux term is weak. This could be avoided by applying the upwind scheme. In this case: \begin{equation} f(z_{i+},k_j) = \begin{cases} f(z_{i},k_j) & v(k_j)>0 \\ f(z_{i+1},k_j) & v(k_j)<0 \end{cases} \end{equation} \begin{equation} f(z_{i-},k_j) = \begin{cases} f(z_{i-1},k_j) & v(k_j)>0 \\ f(z_{i},k_j) & v(k_j)<0 \end{cases} \end{equation} I have implemented the above upwinding method, and my results seem accurate for very small values of $F(z_i)$ throughout the system. However, when $F(z_i)$ gets large, the obtained results lose accuracy and deviate from the correct result.

The reason to this deviation, I guess, is that the analytical solution to $\text{Fr}\{f\}=0$ is an exponential function. Since the equations are linearly discretized, the discretization cannot follow the large exponential changes accurately.

Do you have any idea for increasing the accuracy of the above discretization while holding on to unconditional stability?

Answer 1

Your problem looks very much like "implicit upwind method for advection equation", let me comment it from this point of view.

Let us think that $f=f(z,k,t)$, $f(z,k,0)=f_0(z,k)$ and we search for $f(z,k,\tau)$ that fulfills for $t \in (0,\tau)$ the advection equation (I use a shorter notation for partial derivatives): \begin{equation} \partial_t f + v(k) \partial_z f - F(z) \partial_k f = 0 \end{equation}

Then your upwind scheme is so called implicit full upwind (the first order accurate) method, in which you try to solve for $f(z,k,\tau)$ the advection equation in one time step. For small values of $v$ and $F$ it can work, i.e. the results can have satisfactory accuracy, but for large values the accuracy must be very poor.

Although the implicit upwind scheme has no restriction on the time step (say your $\tau$) due to stability reasons, it has a natural restriction on it due to accuracy. In fact, the standard CFL condition is pretty good guess for such accuracy restriction, it means $$ \tau \le h \max \limits_{i,j} \left(|v(k_j)| + |F(z_i)| \right)^{-1} $$ where $h=\min \{\Delta k, \Delta z\}$.

If for any serious reasons your $\tau$ is given and fixed, you might look to your problem like the suggested advection equation and solve it with Lagrangian methods if you significantly violate the above condition for $\tau$. Such type of method can be implemented with arbitrary large $\tau$.