An implicit high order finite volume scheme with a posteriori limiting for advection networks

In this paper, we describe an implicit high order finite volume scheme for networks of advection equations. Possible oscillations are eliminated using an a posteriori limiting approach. A continuous transition between the different orders is achieved by a convex combination of the numerical fluxes. The scheme shows accurate results even for large CFL numbers.


Introduction
In the modeling of district heating systems, networks of advection equations arise [1]. These models offer despite their simplicity some challenges in their numerical treatment. Often the considered time horizon is large compared to the time steps enforced by explicit methods. As first order implicit methods suffer from large numerical diffusion high order implicit schemes are to be preferred. On a network the flow along one edge is modeled by the equation where the boundary input g is provided by coupling conditions at the nodes. At these points the information can not be transported upstream and thus a corresponding scheme should be upwind based only. This is also beneficial for designing implicit methods iteratively. In high order methods spurious oscillations have to be suppressed by suitable limiters. Here we follow the a posteriori MOOD-type approach. A more detailed analysis of the constructed scheme is in preparation [2].

Finite Volume Scheme
For advection problems the method of characteristics provides the exact solution. If the time step ∆t is chosen large compared to the spacial resolution ∆x, the characteristics reach back over several branches in a network. Therefore we choose a classical finite volume approach, implicit in time to avoid time step restrictions. In order to construct a 4th order scheme at least the values of four cells are required to compute the numerical flux F i+ 1 2 . The most compact upwind stencil uses the cell means u n i , u n+1 i , u n i−1 , u n+1 i−1 , where i is the index of the cell and n the index of the time step. These values are transported to the interface at x i+ 1 2 and serve as interpolating conditions for a cubic polynomial p. This is used to evaluate the numerical flux as F i+ 1 2 = t n+1 t n ap(τ )dτ . The resulting scheme is given by whith the CFL number c = a ∆t ∆x . A different formulation of the same scheme can be found in the CTCS scheme in [3]. Similarly, a scheme of third order can be constructed using the values u n i , u n+1 i and u n+1 i−1 . When applied to Dirichlet problems (1), we can show numerically that both schemes are stable for c ≥ 1, but unstable in the explicit regime c < 1. Note that due to the given upwind direction theses methods can be implemented with a marching technique and no linear systems have to be solved. This makes the numerical costs comparable to those of explicit schemes.

A Posteriori Limiting
High order schemes tend to show spurious oscillations in the presence of discontinuities. To suppress such oscillations appropriate limiter have to be chosen. In [4] implicit WENO schemes are proposed, but these result in large nonlinear systems, which become hard to solve for large time steps. As the scheme (2) is solved marching in the upwind direction, a posteriori limiting as proposed for explicit methods in the MOOD framework [5] can be applied. First the new value is computed with the high order scheme. This solution has to pass several checks, e.g. for physical bounds and local oscillation. If one of these properties is not satisfied the proposed value is recomputed with a scheme of lower order. We use a cascade of schemes with the orders 4 → 3 → 1, where the first order scheme is the implicit upwind method. As the actual domain of dependence exceeds the neighboring cells, this procedure does not perform as good as for explicit schemes. Especially next to hard transitions of high to first order updates some spurious behavior can occur. Such artifacts are avoided using a continuous transition between the different orders. At a neighboring cell interfaces, the numerical fluxes of two schemes with different orders are mixed as a convex combination to further improve the result. Nevertheless with a pure sequential process not all of the desired properties can be enforced, even with a first order update. This results from possible order decisions, that might cause mild overshoots more downstream in the domain. Such a behavior could be cured using several solves or backtracking. For our practical applications however this was not necessary as the errors were always below our preset thresholds.

Numerical Results
As a test case we consider Shu's linear test. Problem (1) with a = 1 and g ≡ 0 is solved on [0, 2]. As the method is just designed for Dirichlet problems we impose the initial condition of  (green). Compared to the exact solution (black) the first order scheme looses any contour even for small CFL numbers. The higher order schemes show oscillations next to discontinuities or kinks. These oscillations can be removed by the proposed limiting strategy applied in the hybrid scheme (red). All four shapes are captured accurately and even for a CFL number of 10 they still can all be distinguished. With the hybrid scheme the height of the peaks is approximated better compared to the third order method and the overshoots of the fourth order method are eliminated.

Conclusion
We extended the approach of a posteriori limiting for explicit schemes [5] to implicit finite volume schemes for the advection equation. By using upwind stencils the implicit scheme can be solved with a marching technique such that the computational costs are comparable to those of explicit methods. By using a convex combination of numerical fluxes artifacts of the limiting procedure are avoided. The scheme shows accurate results even for large CFL numbers as oscillations are suppressed by the limiter and the numerical diffusion is reduced due to the high order accuracy.