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

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.     

A full-wave solver of the Maxwell's equations in 3D cold plasmas

A new solver for Maxwell's equations in three-dimensional (3D) plasma configurations is presented. The new code LEMan (Low-frequency ElectroMagnetic wave propagation) determines a global solution of the wave equation in a realistic stellarator geometry at low frequencies. The code is aimed at the applications with relatively small computational resources and is very efficient in the Alfvén frequency range. In the present work, the cold plasma model is implemented. Finite elements are applied for the radial discretization and the spectral representation is used for the poloidal and toroidal angles. Special care is taken to avoid the numerical pollution of the spectrum as well as to ensure the energy conservation. The numerical scheme and the convergence properties are discussed. Several benchmarks and results in different geometries are presented.     

Two-dimensional full-electromagnetic Vlasov code with conservative scheme and its application to magnetic reconnection

A detailed procedure of full-electromagnetic Vlasov simulation technique is presented. Our new unsplitting conservative scheme exactly satisfies the continuity equation for charge. The implicit Finite Difference Time Domain method is also adopted for computation of electromagnetic fields, which is not restricted by the CFL condition for light. The Geospace Environment Modeling magnetic reconnection challenge problem is adopted as a benchmark test. The characteristics of the present Vlasov code are studied by varying the resolution in configuration space.     

A transport simulation code for inertial confinement fusion relevant laser-plasma interaction

We present a code for the simulation of laser-plasma interaction processes relevant for applications in inertial confinement fusion. The code consists of a fully nonlinear hydrodynamics in two spatial dimensions using a Lagrangian, discontinuous Galerkin-type approach, a paraxial treatment of the laser field and a spectral treatment of the dominant non-local transport terms. The code is fully parallelized using MPI in order to be able to simulate macroscopic plasmas.One example of a fully nonlinear evolution of a laser beam in an underdense plasma is presented for the conditions previewed for the future MegaJoule laser project.     

A FV-TD electromagnetic solver using adaptive Cartesian grids

A second-order finite-volume (FV) method has been developed to solve the time-domain (TD) Maxwell equations, which govern the dynamics of electromagnetic waves. The computational electromagnetic (CEM) solver is capable of handling arbitrary grids, including structured, unstructured, and adaptive Cartesian grids, which are topologically arbitrary. It is argued in this paper that the adaptive Cartesian grid is better than a tetrahedral grid for complex geometries considering both efficiency and accuracy. A cell-wise linear reconstruction scheme is employed to achieve second-order spatial accuracy. Second-order time accuracy is obtained through a two-step Runge-Kutta scheme. Issues on automatic adaptive Cartesian grid generation such as cell-cutting and cell-merging are discussed. A multi-dimensional characteristic absorbing boundary condition (MDC-ABC) is developed at the truncated far-field boundary to reduce reflected waves from this artificial boundary. The CEM solver is demonstrated with several test cases with analytical solutions.     

Numerical solution of the Vlasov-Poisson system using generalized Hermite functions

Two different spectral approaches for solving the nonlinear Vlasov-Poisson equations are presented and discussed. The first approach is based on a standard spectral Galerkin method (SGM) using Hermite functions in the velocity space. The second method which belongs to the family of pseudospectral methods (SCM) uses Gauss-Hermite collocation points for the velocity discretization. The high-dimensional feature of these equations and the suspected presence of small scales in the solution suggested us to employ these methods that provide high order accuracy while considering a “small” number of ad hoc basis functions. The scaled Hermite functions allow us to treat the case of unbounded domains and to properly recover Gaussian-type distributions. Some numerical results on usual test cases are shown and prove the good agreement with the theory.     

Constructing plasma response models from full toroidal magnetohydrodynamic computations

We explore accurate and efficient algorithms for constructing plasma response models, based on the computed data using a full toroidal MHD stability code MARS-F. These response models are used to study feedback stabilization of resistive wall modes for fusion plasmas. Three approaches are discussed and compared. A direct full-model computation offers the most accurate response, unfortunately without producing analytical expressions for the response model. The pole-residue expansion methods yield analytical and asymptotically rigorous response models. A low-order Padé approximation serves as a model reduction technique that simplifies the controller design, while keeping a reasonable accuracy for the response models. From the computational viewpoint, the most efficient approaches are the pole-residue expansion based on eigenfunction projection, and the low-order Padé approximation.     

A parallel algorithm for global magnetic reconnection studies

A new algorithm is presented for the computation of two-dimensional magnetic reconnection in plasmas. Both resistive magnetohydrodynamic (MHD) and two-fluid models are considered. It has been implemented on several parallel platforms and shows good scalability up to 32 CPUs for reasonable problem sizes. A fixed, non-uniform rectangular mesh is used to resolve the different spatial scales in the reconnection problem. The resistive MHD version uses an implicit/explicit hybrid method, while the two-fluid version uses an alternating-direction implicit (ADI) method with high-order artificial dissipation. The technique has proven useful for comparing several different theories of collisional and collisionless reconnection.     

dHybrid: A massively parallel code for hybrid simulations of space plasmas

A massively parallel simulation code, called dHybrid, has been developed to perform global scale studies of space plasma interactions. This code is based on an explicit hybrid model; the numerical stability and parallel scalability of the code are studied. A stabilization method for the explicit algorithm, for regions of near zero density, is proposed. Three-dimensional hybrid simulations of the interaction of the solar wind with unmagnetized artificial objects are presented, with a focus on the expansion of a plasma cloud into the solar wind, which creates a diamagnetic cavity and drives the Interplanetary Magnetic Field out of the expansion region. The dynamics of this system can provide insights into other similar scenarios, such as the interaction of the solar wind with unmagnetized planets.     

A finite-difference eigenvalue algorithm for calculating the band structure of a photonic crystal

A new method based on a finite difference of the governing differential equation for the eigenvalue problem is introduced to calculate the band structure of a two-dimensional photonic crystal. The effective medium technique is also used in the method. The problem is reduced to a standard matrix eigenvalue problem. Compared to the conventional plane wave expansion method, the present method improves the convergence of the solution and thus is a fast and accurate algorithm for calculating the band structure of a photonic crystal.     

A fourth-order symplectic exponentially fitted integrator

A numerical method for ordinary differential equations is called symplectic if, when applied to Hamiltonian problems, it preserves the symplectic structure in phase space, thus reproducing the main qualitative property of solutions of Hamiltonian systems. In a previous paper [G. Vanden Berghe, M. Van Daele, H. Van de Vyver, Exponential fitted Runge-Kutta methods of collocation type: fixed or variable knot points?, J. Comput. Appl. Math. 159 (2003) 217-239] some exponentially fitted RK methods of collocation type are proposed. In particular, three different versions of fourth-order exponentially fitted Gauss methods are described. It is well known that classical Gauss methods are symplectic. In contrast, the exponentially fitted versions given in [G. Vanden Berghe, M. Van Daele, H. Van de Vyver, Exponential fitted Runge-Kutta methods of collocation type: fixed or variable knot points?, J. Comput. Appl. Math. 159 (2003) 217-239] do not share this property. This paper deals with the construction of a fourth-order symplectic exponentially fitted modified Gauss method. The RK method is modified in the sense that two free parameters are added to the Buthcher tableau in order to retain symplecticity.     

A Runge-Kutta-Nyström pair for the numerical integration of perturbed oscillators

New Runge-Kutta-Nyström methods especially designed for the numerical integration of perturbed oscillators are presented in this paper. They are capable of exactly integrating the harmonic or unperturbed oscillator. We construct an embedded 4(3) RKN pair that is based on the FSAL property. The new method is much more efficient than previously derived RKN methods for some subclasses of problems.     

BOUT++: A framework for parallel plasma fluid simulations

A new modular code called BOUT++ is presented, which simulates 3D fluid equations in curvilinear coordinates. Although aimed at simulating Edge Localised Modes (ELMs) in tokamak x-point geometry, the code is able to simulate a wide range of fluid models (magnetised and unmagnetised) involving an arbitrary number of scalar and vector fields, in a wide range of geometries. Time evolution is fully implicit, and 3rd-order WENO schemes are implemented. Benchmarks are presented for linear and non-linear problems (the Orszag-Tang vortex) showing good agreement. Performance of the code is tested by scaling with problem size and processor number, showing efficient scaling to thousands of processors.Linear initial-value simulations of ELMs using reduced ideal MHD are presented, and the results compared to the ELITE linear MHD eigenvalue code. The resulting mode-structures and growth-rate are found to be in good agreement (γBOUT++=0.245ωA, γELITE=0.239ωA, with Alfvénic timescale 1/ωA=R/VA). To our knowledge, this is the first time dissipationless, initial-value simulations of ELMs have been successfully demonstrated.     

A Mathematica program for the two-step twelfth-order method with multi-derivative for the numerical solution of a one-dimensional Schrödinger equation

In this paper, we present the detailed Mathematica symbolic derivation and the program which is used to integrate a one-dimensional Schrödinger equation by a new two-step numerical method. We add the fourth- and sixth-order derivatives to raise the precision of the traditional Numerov's method from fourth order to twelfth order, and to expand the interval of periodicity from (0,6) to the one of (0,9.7954) and (9.94792,55.6062). In the program we use an efficient algorithm to calculate the first-order derivative and avoid unnecessarily repeated calculation resulting from the multi-derivatives. We use the well-known Woods-Saxon's potential to test our method. The numerical test shows that the new method is not only superior to the previous lower order ones in accuracy, but also in the efficiency. This program is specially applied to the problem where a high accuracy or a larger step size is required.

Program summary

Title of program: ShdEq.nbCatalogue number: ADTTProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADTTProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputer for which the program is designed and others on which it has been tested: The program has been designed for the microcomputer and been tested on the microcomputer.Computers: IBM PCOperating systems under which the program has been tested: Windows XPProgramming language used: Mathematica 4.2Memory required to execute with typical data: 51 712 bytesNo. of bytes in distributed program, including test data, etc.: 45 381No. of lines in distributed program, including test data, etc.: 7311Distribution format: tar gzip fileCPC Program Library subprograms used: noNature of physical problem: Numerical integration of one-dimensional or radial Schrödinger equation to find the eigenvalues for a bound states and phase shift for a continuum state.Method of solution: Using a two-step method twelfth-order method to integrate a Schrödinger equation numerically from both two ends and the connecting conditions at the matching point, an eigenvalue for a bound state or a resonant state with a given phase shift can be found.Restrictions on the complexity of the problem: The analytic form of the potential function and its high-order derivatives must be known.Typical running time: Less than one second.Unusual features of the program: Take advantage of the high-order derivatives of the potential function and efficient algorithm, the program can provide all the numerical solution of a given Schrödinger equation, either a bound or a resonant state, with a very high precision and within a very short CPU time. The program can apply to a very broad range of problems because the method has a very large interval of periodicity.References: [1] T.E. Simos, Proc. Roy. Soc. London A 441 (1993) 283.[2] Z. Wang, Y. Dai, An eighth-order two-step formula for the numerical integration of the one-dimensional Schrödinger equation, Numer. Math. J. Chinese Univ. 12 (2003) 146.[3] Z. Wang, Y. Dai, An twelfth-order four-step formula for the numerical integration of the one-dimensional Schrödinger equation, Internat. J. Modern Phys. C 14 (2003) 1087.     

A semi-Lagrangian code for nonlinear global simulations of electrostatic drift-kinetic ITG modes

A semi-Lagrangian code for the solution of the electrostatic drift-kinetic equations in straight cylinder configuration is presented. The code, CYGNE, is part of a project with the long term aim of studying microturbulence in fusion devices. The code has been constructed in such a way as to preserve a good control of the constants of motion, possessed by the drift-kinetic equations, until the nonlinear saturation of the ion-temperature-gradient modes occurs. Studies of convergence with phase space resolution and time-step are presented and discussed. The code is benchmarked against electrostatic Particle-in-Cell codes.     

Qprop: A Schrödinger-solver for intense laser-atom interaction

The Qprop package is presented. Qprop has been developed to study laser-atom interaction in the nonperturbative regime where nonlinear phenomena such as above-threshold ionization, high order harmonic generation, and dynamic stabilization are known to occur. In the nonrelativistic regime and within the single active electron approximation, these phenomena can be studied with Qprop in the most rigorous way by solving the time-dependent Schrödinger equation in three spatial dimensions. Because Qprop is optimized for the study of quantum systems that are spherically symmetric in their initial, unperturbed configuration, all wavefunctions are expanded in spherical harmonics. Time-propagation of the wavefunctions is performed using a split-operator approach. Photoelectron spectra are calculated employing a window-operator technique. Besides the solution of the time-dependent Schrödinger equation in single active electron approximation, Qprop allows to study many-electron systems via the solution of the time-dependent Kohn-Sham equations.

Program summary

Program title:QPROPCatalogue number:ADXBProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXBProgram obtainable from:CPC Program Library, Queen's University of Belfast, N. IrelandComputer on which program has been tested:PC Pentium IV, AthlonOperating system:LinuxProgram language used:C++Memory required to execute with typical data:Memory requirements depend on the number of propagated orbitals and on the size of the orbitals. For instance, time-propagation of a hydrogenic wavefunction in the perturbative regime requires about 64 KB RAM (4 radial orbitals with 1000 grid points). Propagation in the strongly nonperturbative regime providing energy spectra up to high energies may need 60 radial orbitals, each with 30000 grid points, i.e. about 30 MB. Examples are given in the article.No. of bits in a word:Real and complex valued numbers of double precision are usedNo. of lines in distributed program, including test data, etc.:69 995No. of bytes in distributed program, including test data, etc.: 2 927 567Peripheral used:Disk for input-output, terminal for interaction with the userCPU time required to execute test data:Execution time depends on the size of the propagated orbitals and the number of time-stepsDistribution format:tar.gzNature of the physical problem:Atoms put into the strong field of modern lasers display a wealth of novel phenomena that are not accessible to conventional perturbation theory where the external field is considered small as compared to inneratomic forces. Hence, the full ab initio solution of the time-dependent Schrödinger equation is desirable but in full dimensionality only feasible for no more than two (active) electrons. If many-electron effects come into play or effective ground state potentials are needed, (time-dependent) density functional theory may be employed. Qprop aims at providing tools for (i) the time-propagation of the wavefunction according to the time-dependent Schrödinger equation, (ii) the time-propagation of Kohn-Sham orbitals according to the time-dependent Kohn-Sham equations, and (iii) the energy-analysis of the final one-electron wavefunction (or the Kohn-Sham orbitals).Method of solution:An expansion of the wavefunction in spherical harmonics leads to a coupled set of equations for the radial wavefunctions. These radial wavefunctions are propagated using a split-operator technique and the Crank-Nicolson approximation for the short-time propagator. The initial ground state is obtained via imaginary time-propagation for spherically symmetric (but otherwise arbitrary) effective potentials. Excited states can be obtained through the combination of imaginary time-propagation and orthogonalization. For the Kohn-Sham scheme a multipole expansion of the effective potential is employed. Wavefunctions can be analyzed using the window-operator technique, facilitating the calculation of electron spectra, either angular-resolved or integratedRestrictions onto the complexity of the problem:The coupling of the atom to the external field is treated in dipole approximation. The time-dependent Schrödinger solver is restricted to the treatment of a single active electron. As concerns the time-dependent density functional mode of Qprop, the Hartree-potential (accounting for the classical electron-electron repulsion) is expanded up to the quadrupole. Only the monopole term of the Krieger-Li-Iafrate exchange potential is currently implemented. As in any nontrivial optimization problem, convergence to the optimal many-electron state (i.e. the ground state) is not automatically guaranteedExternal routines/libraries used:The program uses the well established libraries blas, lapack, and f2c     

A new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation

In this paper we present a new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation. The scheme uses a symmetrical multi-point difference formula to represent the partial differentials of the two-dimensional variables, which can improve the accuracy of the numerical solutions to the order of Δx2Nq+2 when a (2Nq+1)-point formula is used for any positive integer Nq with Δx=Δy, while Nq=1 equivalent to the traditional scheme. On the other hand, the new scheme keeps the same form of the traditional matrix equation so that the standard algebraic eigenvalue algorithm with a real, symmetric, large sparse matrix is still applicable. Therefore, for the same dimension, only a little more CPU time than the traditional one should be used for diagonalizing the matrix. The numerical examples of the two-dimensional harmonic oscillator and the two-dimensional Henon-Heiles potential demonstrate that by using the new method, the error in the numerical solutions can be reduced steadily and extensively through the increase of Nq, which is more efficient than the traditional methods through the decrease of the step size.     

A numerical differentiation library exploiting parallel architectures

We present a software library for numerically estimating first and second order partial derivatives of a function by finite differencing. Various truncation schemes are offered resulting in corresponding formulas that are accurate to order O(h), O(h²), and O(h⁴), h being the differencing step. The derivatives are calculated via forward, backward and central differences. Care has been taken that only feasible points are used in the case where bound constraints are imposed on the variables. The Hessian may be approximated either from function or from gradient values. There are three versions of the software: a sequential version, an OpenMP version for shared memory architectures and an MPI version for distributed systems (clusters). The parallel versions exploit the multiprocessing capability offered by computer clusters, as well as modern multi-core systems and due to the independent character of the derivative computation, the speedup scales almost linearly with the number of available processors/cores.

Program summary

Program title: NDL (Numerical Differentiation Library)Catalogue identifier: AEDG_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEDG_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 73 030No. of bytes in distributed program, including test data, etc.: 630 876Distribution format: tar.gzProgramming language: ANSI FORTRAN-77, ANSI C, MPI, OPENMPComputer: Distributed systems (clusters), shared memory systemsOperating system: Linux, SolarisHas the code been vectorised or parallelized?: YesRAM: The library uses O(N) internal storage, N being the dimension of the problemClassification: 4.9, 4.14, 6.5Nature of problem: The numerical estimation of derivatives at several accuracy levels is a common requirement in many computational tasks, such as optimization, solution of nonlinear systems, etc. The parallel implementation that exploits systems with multiple CPUs is very important for large scale and computationally expensive problems.Solution method: Finite differencing is used with carefully chosen step that minimizes the sum of the truncation and round-off errors. The parallel versions employ both OpenMP and MPI libraries.Restrictions: The library uses only double precision arithmetic.Unusual features: The software takes into account bound constraints, in the sense that only feasible points are used to evaluate the derivatives, and given the level of the desired accuracy, the proper formula is automatically employed.Running time: Running time depends on the function's complexity. The test run took 15 ms for the serial distribution, 0.6 s for the OpenMP and 4.2 s for the MPI parallel distribution on 2 processors.