On the efficiency and robustness of implicit methods in computational astrophysics

Most of the numerical methods used in astrophysical fluid dynamics rely on explicit time-stepping schemes, whereas the higher robustness of implicit m ethods which constitute the core of modern computational fluid dynamics is rarely explored. In this paper, we survey some modem implicit solvers which are specially adapted to multi-dimensional problems and discuss their pote ntial and range of application in comparison to common explicit methods. Sp ecial emphasis is put on the aspect of efficiency and robustness. Our refer ence set of equations are those corresponding to radiative magnetohydrodyna mics (MHD) modeling self-gravitating and partially and/or fully ionized flo ws. Explicit methods may be viewed as a very special class of so-called 'de fect-correction iterations' for solving an implicit discretization. Within this context one can design various implicit methods, ranging from weakly t o fully implicit, which allow to follow evolutionary phases on much longer time scales than the dynamic one. We particularly present a new three-stage s implicit numerical method for searching strongly time-dependent, quasi-st ationary and steady-state solutions for the above-mentioned equations. Prec onditioned Krylov-space and multilevel techniques are employed for enhancin g the efficiency and robustness of the computation. The spatial discretizat ion is on highly nonuniform tenser-product meshes and uses cartesian, cylin der or spherical coordinates depending on the geometrical structure of the problem. The accuracy is of second-order in space and time and can easily b e increased without modifying the structure of the scheme. The algebraic so lver consists of a pre-conditioned transpose-free Krylov iteration for the conservation equations and optimized multigrid algorithms for solving the P oisson equation for the gravitational potential and the transport-diffusion equation for the radiation density in flow regions with dynamically varyin g optical depths. The existing implementation of the proposed method employ s axis-symmetry in order to reduce the problem to two dimensions. (C) 2001 Elsevier Science B.V. All rights reserved.