On a semi-implicit scheme for spectral simulations of dispersive magnetohydrodynamics
Introduction
High-resolution numerical integration of the magnetohydrodynamic equations with a Hall term (Hall-MHD, hereafter referred to as HMHD) [1], for an electron–proton plasma permeated by a strong ambient field experiences specific difficulties due to the right-hand polarized Alfvén waves (whistler modes) whose frequency scales like the squared wavenumber, making explicit schemes very penalizing at high resolution [2]. At the opposite, implicit algorithms [3], [4], [5], [6] are unconditionally stable but, due to the computational cost of each time step (that often requires a preconditioner specific to the considered problem), they are only worthwhile when the latter is relatively large. This may result in an insufficient accuracy, for example in dispersive turbulent regimes where a broad range of scales (together with the associated frequencies) is to be resolved. This contrasts with dissipative problems where the modes whose decay time is much smaller than the time step, are essentially slaved on the time scale of the dynamically relevant phenomena, as expected on the basis of a central manifold theorem. No such argument exists in the dispersive case.
In the presence of dispersion, the stability of implicit schemes results in fact from a distortion of the dispersion relation (at a scale that increases with the time step), aimed to produce a saturation of the wave frequency at large wavenumbers. An alternative approach whose implementation is especially simple for spectral codes and is practically free of computational overhead is provided by the so-called semi-implicit (SI) schemes that, by also distorting the dispersion relation at small scales, allow to relax the temporal stability constraint. This is done by including an additional linear operator in such a way that consistency is not affected. Such algorithms were used for the HMHD equations [7], and also in ideal MHD in order to artificially slow down the highest compressional modes or both the compressional and shear Alfvén waves ([8], [9], [10] and references therein). A different approach (not discussed here) to improve stability, consists in “sub-cycling the Hall physics”, by using a multi-time step integration algorithm where the magnetic field is advanced with a time step much shorter than the ideal MHD time step used for the velocity and the density [11].
The aim of the present paper is to analyze the advantages and limitations of a SI scheme for spectral simulations of compressible HMHD equations in a periodic geometry, when small-scale dispersion plays a central role while viscous and Joule dissipations are absent or very weak. All the simulations are performed in periodic domains, using a Fourier spectral method. It turns out that for the problems addressed in this paper, aliasing must be removed. Because of the quadratic character of most of the nonlinearities involved in the HMHD equations, this is done by spectral truncation of the computed nonlinear terms outside an ellipsoid whose principal axes retain 2/3 of the maximal wavenumber in each direction. In cases where nonlinearities are cubic, as in the derivative nonlinear Schrödinger (DNLS) equation considered in Sections 7 Strong Alfvénic turbulence, 8 Monitoring the error, only 1/2 of the spectral domain is kept. The spatial resolutions given in the following sections are the effective ones, after aliasing has been suppressed. In Section 2, we introduce the HMHD equations and their one-dimensional reduction, Sections 3 defines a SI algorithm chosen for its efficiency. Its properties are compared in Section 4 with those of implicit schemes in the context of linear Alfvén waves for which the optimal choice for a free parameter involved by the SI scheme is also discussed. The other sections concentrate on the validation of this scheme for various nonlinear problems where dispersion plays a central role. This includes the modulational instabilities of Alfvén waves and their nonlinear developments (Section 5), the filamentation instability that leads to the formation of coherent structures (Section 6), the strong Alfvénic turbulence governed by the randomly driven DNLS equation (Section 7). A procedure to monitor the errors is discussed in Section 8. Section 9 is the conclusion.
Section snippets
The Hall-MHD description
HMHD can be viewed as a bi-fluid description of a plasma when the electron inertia is neglected. The presence of the Hall term in the generalized Ohm's law, and thus in the induction equation, allows a decoupling of the ion fluid from the electron one in which the magnetic field lines are frozen. Choosing as units the Alfvén speed, the amplitude of the ambient magnetic field (taken in the x direction), the equilibrium density and a typical length scale where is the ion inertial length
A semi-implicit scheme
In the absence of dissipation, second-order explicit temporal schemes of Runge–Kutta or Adams–Bashford type, are unconditionally unstable for advection-like equations ([12], pp. 137 and 150) and a fortiori for HMHD. Differently, third-order Adams–Bashford (AB3) or Runge–Kutta (RK3) are conditionally stable (and dissipative), the stability condition on the time step Δt and the mesh size Δx being of the type (or in terms of the maximum wavenumber ), because of the whistler
Efficiency of the SI scheme for linear Alfvén waves
When Eqs. (8) are linearized about a uniform ambient field, the equations for the transverse components of the velocity and magnetic field (that govern the Alfvén waves) decouple from the longitudinal ones which, together with the density fluctuations, describe the sonic waves. One then gets that contains all the difficulties of the original problem for what concerns the linear stability of temporal schemes. This system thus provides a convenient soluble model to address
Modulational instabilities in one dimension
As a simple one-dimensional evolution problem, we consider the modulation instabilities of Alfvén waves and their nonlinear developments. Denoting by k the pump wavenumber, the instability is said to be modulational when the wavenumber K of the amplified density and parallel velocity modes obeys [15]. At the level of the transverse velocity and magnetic components, this corresponds to the amplification of sideband modes of wavenumbers . When the range of amplified wavenumbers K is
Filamentation instability
As a first example of three-dimensional simulations of the HMHD equations, we use the SI scheme to reproduce a simulation of Alfvén wave filamentation in the presence of an initial low-density channel aligned with the ambient magnetic field, performed in [17] with an explicit RK3 time stepping.
In the simulations presented in this section, we take for easier comparison with previous publications, thus choosing a length unit equal to 4 ion inertial lengths. Using this unit, the size of the
Strong Alfvénic turbulence
In order to concentrate on the dynamics of nonlinear Alfvén waves propagating along the ambient field, it is convenient to perform a reductive perturbative expansion of inviscid Eqs. (5), (6), (7), (8) in the long-wavelength limit, assuming the wave amplitude small compared with the magnitude of the ambient field. This leads to the DNLS equation [21] for a complex magnetic field , in the Alfvén wave moving frame defined by . In physical units, and
Monitoring the error
In order to quantitatively evaluate the accuracy of the SI scheme at a given stage of the evolution, a relatively simple procedure was implemented. We consider the SI solution at a given time t and perform one more time step of size Δt with the SI scheme, taking off the random forcing if present. We also advance the solution of the same time interval (starting from the same initial condition ) by performing n time steps of magnitude with the explicit RK3 scheme, still in
Conclusion
We have discussed the implementation of a SI time stepping to large-resolution simulations of dispersive magnetohydrodynamics where small scale dispersion plays a main dynamical role. As well known, the stability conditions of fully explicit schemes that require the accurate resolution of the frequencies associated with all the retained spatial scales, prescribe in this regime tiny time steps and thus huge computational costs. On the other hand, implicit schemes that are computationally
Acknowledgements
The work was supported by “Programme National Terre Soleil” of INSU-CNRS. Three-dimensional simulations were performed on the “Mesocetre SIGAMM” machine, hosted by Observatoire de la Côte d'Azur. The work of DB was partly supported by the Euratom Communities under the contract of Association between EURATOM/ENEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
References (24)
- et al.
A 2D high-β Hall MHD implicit nonlinear solver
J. Comp. Phys.
(2003) - et al.
An iterative semi-implicit scheme with robust damping
J. Comp. Phys.
(2008) - et al.
Accurate semi-implicit treatment of the Hall effect in magnetohydrodynamic computations
J. Comp. Phys.
(1989) - et al.
Semi-implicit magnetohydrodynamic calculations
J. Comp. Phys.
(1987) Derivation of implicit difference schemes by the method of differential approximations
J. Comp. Phys.
(1991)- et al.
A new semi-implicit method for MHD computations
J. Comp. Phys.
(1991) Low storage Runge Kutta schemes
J. Comp. Phys.
(1980)- et al.
A semi-implicit Hall-MHD solver using whistler wave preconditioning
Comput. Phys. Comm.
(2008) Hall magnetohydrodynamics in space and laboratory plasmas
Phys. Plasmas
(1995)- et al.
Energy transfer in Hall-MHD turbulence; cascades, backscatter, and dynamo action
J. Plasma Phys.
(2007)
Jacobian-Free Newton–Krylov methods for the accurate time integration of stiff wave systems
J. Comp. Phys.
A fully implicit 3D extended magnetohydrodynamics algorithm
Cited by (0)
- 1
Present address: Burning Plasma Research Group, Dipartimento di Energetica, Politecnico di Torino, C.so Duca degli Abruzzi 24, 10129 Torino, Italy.