Paper (PostScript).

Adding Implicit and Semi-Implicit Schemes to the Versatile Advection Code

R. Keppens, FOM-Institute for Plasma-Physics Rijnhuizen,
P.O. Box 1207, 3430 BE Nieuwegein, The Netherlands,
G. Tóth, Astronomical Institute, Utrecht University,
P.O. Box 80000, 3508 TA Utrecht, The Netherlands,
M. A. Botchev, Mathematical Institute, Utrecht University,
P.O. Box 80010, 3508 TA Utrecht, The Netherlands,
A. van der Ploeg, Centrum voor Wiskunde en Informatica,
P.O. Box 94079, 1090 GB Amsterdam, The Netherlands,

Subject Classification and Key Words:
35L65-Conservation Laws
65D99-Numerical approximation
65M99-PDEs, initial value and time-dependent
initial-boundary value problems


We formulate general guidelines to extend an explicit code with a great variety of implicit and semi-implicit time integration schemes. The discussion is based on their specific implementation in the Versatile Advection Code, which is a general purpose software package for solving systems of non-linear hyperbolic (and/or parabolic) partial differential equations, using standard high-resolution shock-capturing schemes. For all combinations of explicit high-resolution schemes with implicit and semi-implicit treatments, we show how second order spatial and temporal accuracy for the smooth part of the solutions can be maintained. We discuss strategies to obtain steady state and time accurate solutions implicitly.

The implicit and semi-implicit schemes require the solution of large linear systems containing the Jacobian matrix. The Jacobian matrix itself is calculated numerically to ensure the generality of our implementation. We discuss three options available in VAC in terms of applicability, storage requirements, and computational efficiency. One option is the easily implemented matrix-free approach, but the Jacobian matrix can also be calculated by using a general grid masking algorithm, or by an efficient implementation for a specific Lax-Friedrich type Total Variation Diminishing spatial discretization.

The choice of the linear solver depends on the dimensionality of the problem. In one dimension a direct block tridiagonal solver can be applied, while in more than one spatial dimension a Conjugate Gradient-type iterative solver is used. For advection dominated problems, preconditioning is needed to accelerate the convergence of the iterative schemes. We implemented the Modified Block Incomplete LU-preconditioner which performs very well. As an alternative approach, we describe the Minimal Residual Predictor-Corrector time stepping scheme, which can be combined with the matrix-free method.