On a semi-implicit scheme for spectral simulations of dispersive magnetohydrodynamics

https://doi.org/10.1016/j.cpc.2009.05.018Get rights and content

Abstract

The effectiveness of a semi-implicit (SI) temporal scheme is discussed in the context of the dispersive magnetohydrodynamics where, due to the whistler modes, stability of explicit algorithms requires a time step decreasing quadratically as the resolution is linearly increased. After analyzing the effects of this scheme on the Alfvén-wave dispersion relation, spectral simulations of nonlinear initial value problems where small-scale dispersion has a main effect on the global dynamics are presented. Permitting a moderate, albeit significant, increase of the time step for a minor additional cost relatively to explicit schemes, the SI algorithm provides an efficient tool in situations, such as turbulent regimes, where the time steps making fully implicit schemes efficient are too large to ensure a satisfactory accuracy.

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 L=Ridi where di 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 Δt<cΔx2 (or Δt<C/kmax2 in terms of the maximum wavenumber kmax), 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 getstv=xb,tb=xvixxb, 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 K<k [15]. At the level of the transverse velocity and magnetic components, this corresponds to the amplification of sideband modes of wavenumbers k±K. 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 Ri=4 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]tb+αξ(|b|2b)+iδξξb=0, for a complex magnetic field b=by+ibz, in the Alfvén wave moving frame defined by ξ=xt. In physical units, α=14(1β) 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 bSI(x,t) 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 bSI(x,t)) by performing n time steps of magnitude Δt 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)

  • D.A. Knoll et al.

    Jacobian-Free Newton–Krylov methods for the accurate time integration of stiff wave systems

    J. Comp. Phys.

    (2005)
  • L. Chacón

    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.

    View full text