Implementation of an implicit shear Alfvén operator in the M3D code

A new pair of linear operators is introduced into the partially implicit but largely explicit time advance algorithm for the MHD equations in the M3D code. These operators jointly render the algorithm implicit with respect to propagation of the shear Alfvén wave to leading order. As a result, the maximum stable time step is increased by up to a factor of three for many practical applications. The additional computational cost is trivial for the linear mode of operation and adds only approximately 30% additional execution time per step for nonlinear calculations.     

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.     

Vlasov simulation of kinetic shear Alfvén waves

The treatment of kinetic shear Alfvén waves in homogeneous magnetized plasmas by means of Vlasov simulation is examined. To this end, the driftkinetic version of the Vlasov-Maxwell equations is solved via various numerical schemes, all employing a grid in (1+1)D phase space. Since kinetic shear Alfvén waves are Landau damped, the use of an equidistant grid in velocity space leads to a recurrence problem. The latter can be circumvented, however, by damping the finest velocity space scales through higher-order collision operators. Of particular interest is the question if and under which circumstances the magnetohydrodynamic limit (small perpendicular wavenumber) can be recovered.     

Multi-fluid MHD model and calculations of Alfvén wave spectrum and dissipation in tokamaks

The numerical method is developed for calculations of wave excitation and dissipation in Alfvén and in Ion Cyclotron Range of Frequency (ICRF) in axisymmetric tokamaks. Multi-fluid magneto-hydrodynamic plasma model is used and two-dimensional inhomogeneity of plasma parameters with arbitrary cross section of magnetic surfaces is considered. The difference scheme for the wave equation is not connected to magnetic field geometry and is suitable for the method extensions to nonlinear and three-dimensional case. Special care is taken to avoid the spectrum distortion and pollution. Relevant benchmark cases are presented. Finally, the results of numerical calculations of Alfvén wave absorption are presented for the experimental conditions foreseen for the Tokamak Chauffage Alfvén wave experiment in Brazil (TCABR) [Nucl. Fusion 30 (1996) 503]. In particular, the effect of toroidal mode coupling on the power deposition of Global Alfvén Wave (GAW) eigenmodes is demonstrated.     

A global collisionless PIC code in magnetic coordinates

A global plasma turbulence simulation code, ORB5, is presented. It solves the gyrokinetic electrostatic equations including zonal flows in axisymmetric magnetic geometry. The present version of the code assumes a Boltzmann electron response on magnetic surfaces. It uses a Particle-In-Cell (PIC), δf scheme, 3D cubic B-splines finite elements for the field solver and several numerical noise reduction techniques. A particular feature is the use of straight-field-line magnetic coordinates and a field-aligned Fourier filtering technique that dramatically improves the performance of the code in terms of both the numerical noise reduction and the maximum time step allowed. Another feature is the capability to treat arbitrary axisymmetric ideal MHD equilibrium configurations. The code is heavily parallelized, with scalability demonstrated up to 4096 processors and 10⁹ marker particles. Various numerical convergence tests are performed. The code is validated against an analytical theory of zonal flow residual, geodesic acoustic oscillations and damping, and against other codes for a selection of linear and nonlinear tests.     

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.     

Parallelizaton of visual magneto-hydrodynamics code based on fluctuation distribution scheme on triangular grids

A two-dimensional (2D) magneto-hydrodynamics (MHD) code which has visualization and parallel processing capability is presented in this paper. The code utilizes a fluctuation splitting (FS) scheme that runs on structured or unstructured triangular meshes. First FS scheme which included the wave model: Model-A had been developed by Roe [Roe PL. Discrete models for the numerical analysis of time-dependent multi-dimensional gas dynamics. J Comp Phys 1986;63:458-76.] for the solutions of Euler’s equations. The first 2D-MHD wave model: MHD-A, was then developed by Balci and Aslan [Balci ?. The numerical solutions of two dimensional MHD equations by fluctuation splitting scheme on triangular meshes, Ph.D. Thesis, University of Marmara, Science-Art Faculty, Physics Dept Istanbul, Turkey; 2000; Aslan N. MHD-A: A fluctuation splitting wave model for planar magnetohydrodynamics. J Comp Phys 1999;153:437-66.] to solve MHD problems including shocks and discontinuities. It was then shown in [Balci S, Aslan N. Two dimensional MHD solver by fluctuation splitting and dual time stepping. Int J Numer Meth Fluids, in press.] that this code was capable of producing reliable results in compressible and nearly incompressible limits and under the effect of gravitational fields and that it was able to identically reduce to model-A of Roe in Euler limit with no sonic problems at rarefaction fans (Balci and Aslan, in press). An important feature of this code is its ability to run time dependent or steady problems on structured or unstructured triangular meshes that can be generated automatically by the code for specified domains. In order to use the parallel processing capability of the code, the triangular meshes are decomposed into different blocks in order to share the workload among a number of processors (here personal computers) which are connected by Ethernet. Due to the compact nature of the FS scheme, only one set of data transfer is required between neighbor processors. As it will be shown, this phenomenon results in minimum amount of communication loss and makes the scheme rather robust for parallel processing. The other important feature of the new code is its visual capability. As the code is running, colorful images of scalar quantities (density, pressure, Mach number, etc.) or vector graphics of vectoral quantities (velocity, magnetic field, etc.) can be followed on the screen. The extended code, called PV-MHDA, also allows following the trajectories of the particles in time by means of a recently included particle in cell (PIC) algorithm. Because the numerical dissipation embedded in its wave model reflects real physical viscosity and resistivity, it is able to run accurately for compressible flows (including shocks) as well as nearly incompressible flows (e.g., Kelvin-Helmholtz instability). The user-friendly visual and large-scale computation capability of the code allow the user more thorough analysis of MHD problems in two-dimensional complex domains.     

Numerical model of a two-dimensional,non-plane transient magnetohydrodynamic flow

The equations describing two-dimensional three-component magnetohydrodynamic (MHD) transient flows are formulated for a system of spherical coordinates. With the numerical code based on Implicit Continuous Fluid Eulerian (ICE) scheme, MHD flows resulting from a sudden energy release in a stratified medium are examined. Because of the inclusion of out-of-plane components of velocity and magnetic fields, MHD transverse waves are observed in addition to fast, slow and entropy waves. Numerical results for compressible MHD shocks are found in satisfactory agreement with the theoretical predictions.     

Equilibrium and stability studies with the 3D MHD code tube

TUBE is a 3D code for computing ideal MHD equilibrium and stability by minimizing the potential energy of the system in successive iteration steps over all points of the spatial grid. The independent variables are chosen such that s = const describes a magnetic surface and, together with a linear relation between the periodic variables u and v, specifies a line of force on this surface. These variables are Lagrangian in the sense that the line of force through a certain triple u, v, s stays the same throughout the iteration procedure. The dependent variables describe the position corresponding to the grid points of u, v, s in real space. The rate of convergence of the iteration schemes implemented and the extrapolation of the results to the limit of an infinitely fine grid are extensively studied. TUBE results are compared with those of other codes and with analytical results. The TUBE code is described in some detail.     

Efficient magnetohydrodynamic simulations on distributed multi-GPU systems using a novel GPU Direct–MPI hybrid approach

Modern graphics processing units (GPUs) have been widely utilized in magnetohydrodynamic (MHD) simulations in recent years. Due to the limited memory of a single GPU, distributed multi-GPU systems are needed to be explored for large-scale MHD simulations. However, the data transfer between GPUs bottlenecks the efficiency of the simulations on such systems. In this paper we propose a novel GPU Direct–MPI hybrid approach to address this problem for overall performance enhancement. Our approach consists of two strategies: (1) We exploit GPU Direct 2.0 to speedup the data transfers between multiple GPUs in a single node and reduce the total number of message passing interface (MPI) communications; (2) We design Compute Unified Device Architecture (CUDA) kernels instead of using memory copy to speedup the fragmented data exchange in the three-dimensional (3D) decomposition. 3D decomposition is usually not preferable for distributed multi-GPU systems due to its low efficiency of the fragmented data exchange. Our approach has made a breakthrough to make 3D decomposition available on distributed multi-GPU systems. As a result, it can reduce the memory usage and computation time of each partition of the computational domain. Experiment results show twice the FLOPS comparing to common 2D decomposition MPI-only implementation method. The proposed approach has been developed in an efficient implementation for MHD simulations on distributed multi-GPU systems, called MGPU–MHD code. The code realizes the GPU parallelization of a total variation diminishing (TVD) algorithm for solving the multidimensional ideal MHD equations, extending our work from single GPU computation (Wong et al., 2011) to multiple GPUs. Numerical tests and performance measurements are conducted on the TSUBAME 2.0 supercomputer at the Tokyo Institute of Technology. Our code achieves 2 TFLOPS in double precision for the problem with 1200³ grid points using 216 GPUs.     

The modified differential transform method for solving MHD boundary-layer equations

A new analytical method (DTM-Padé) was developed for solving magnetohydrodynamic boundary-layer equations. It was shown that differential transform method (DTM) solutions are only valid for small values of independent variable. Therefore the DTM is not applicable for solving MHD boundary-layer equations, because in the boundary-layer problem y→∞. Numerical comparisons between the DTM-Padé and numerical methods (by using a fourth-order Runge-Kutta and shooting method) revealed that the new technique is a powerful method for solving MHD boundary-layer equations.     

Multi-dimensional semi-Lagrangian scheme that guarantees exact conservation

A new numerical method that guarantees exact mass conservation is proposed to solve multi-dimensional hyperbolic equations in semi-Lagrangian form without directional splitting. The method is based on a concept of CIP scheme and keep the many good characteristics of the original CIP scheme. The CIP strategy is applied to the integral form of variable. Although the advection and non-advection terms are separately treated, the mass conservation is kept in a form of spatial profile inside a grid cell. Therefore, it retains various advantages of the semi-Lagrangian schemes with exact conservation that has been beyond the capability of conventional semi-Lagrangian schemes.     

Chebyshev wavelet collocation method for magnetohydrodynamic flow equations

This study proposes Chebyshev wavelet collocation method for partial differential equation and applies to solve magnetohydrodynamic (MHD) flow equations in a rectangular duct in the presence of transverse external oblique magnetic field. Approximate solutions of velocity and induced magnetic field are obtained for steady‐state, fully developed, incompressible flow for a conducting fluid inside the duct. Numerical results of the MHD flow problem show that the accuracy of proposed method is quite good even in the case of a small number of grid points. The results for velocity and induced magnetic field are visualized in terms of graphics for values of Hartmann number H_a ≤ 1000.

     

Convolution and iterative methods applied to low-frequency waves in 3D warm configurations

Kinetic effects have been implemented in the three-dimensional full-wave code LEMan. These allow modeling Landau damping as well as the Kinetic Alfvén Wave (KAW). Special attention has been taken on the exact computation of the parallel wave vector as it can have a strong impact on the results. Due to numerical specificities, two techniques have been used depending on the frequency domain under investigation. The damping of global modes like Toroidally-induced Alfvén Eigenmodes (TAE), which is a key issue for future reactors as they can be driven by fast ions, can then be determined. Multiple couplings are shown to lead to an increase of the damping value and thus to a stabilization of the mode. In the Ion Cyclotron Range of Frequencies (ICRF), a Large Helical Device (LHD) configuration has been investigated. For this specific case, the difference between the results obtained with an accurate determination of the parallel wave vector and the approximation k_∥=n/R has been shown to be insignificant.     

New numerical method for the solutions of the MCDF equations based on the CIP method

A new numerical method based on the constrained interpolation profile (CIP) method to solve the Multiconfiguration Dirac-Fock (MCDF) equations is presented. The radial wave functions are represented by the values and the spatial derivatives on an arbitrary grid system, and approximated by cubic polynomials. Owing to this representation, the values and the spatial derivatives of the effective charge distribution and inhomogeneous term are calculated using the previous cycle's wave functions. Then the homogeneous MCDF equations are integrated to obtain two linearly independent solutions, which are used to construct the Green function, by the adaptive stepsize controlled Runge-Kutta method controlling the truncation errors within a prescribed accuracy. The radial wave function is improved by taking the convolution of the Green function and the inhomogeneous term. The effectiveness of this numerical procedure is investigated after implementing it into the relativistic atomic structure code GRASP92.     

MOMCON: A spectral code for obtaining three-dimensional magnetohydrodynamic equilibria

A new code, MOMCON (spectral moments code with constraints), is described that computes three-dimensional ideal magnetohydrodynamic (MHD) equilibria in a fixed toroidal domain using a Fourier expansion for the inverse coordinates (R, Z) representing nested magnetic surfaces. A set of nonlinear coupled ordinary differential equations for the spectral coefficients of (R, Z) is solved using an accelerated steepest descent method. A stream function, λ, is introduced to improve the mode convergence properties of the Fourier series for R and Z. The convergence rate of the R - Z spectra is optimized on each flux surface by solving nonlinear constraint equations relating the m ≥ 2 spectral coefficients of R and Z.     

Numerical experimentation on spherically symmetric one-dimensional magnetohydrodynamic (MHD) wave propagation

Radial propagation of one-dimensional magnetohydrodynamic (MHD) waves are analyzed numerically on the basis of the Implicit-Continuous-Fluid-Eulerian (ICE) scheme. Accuracy of the numerical method and other properties are tested through the study of MHD wave propagation. The three different modes of MHD waves (i.e. fast- slow- and Alfven (transverse) mode) are generated by applying physically consistent boundary perturbations derived from MHD compatibility relations. It is shown that the resulting flow following these waves depend upon the relative configurations of the initial magnetic field and boundary perturbations.     

Accurate numerical method for the solutions of the Schrödinger equation and the radial integrals based on the CIP method

A new accurate numerical method based on the constrained interpolation profile (CIP) method to solve the Schrödinger wave equation for bound and free states in central fields and to calculate radial integrals is presented. The radial wave equation is integrated on an arbitrary grid system by the adaptive stepsize controlled Runge-Kutta method controlling the truncation errors within a prescribed accuracy. For the continuum orbitals in the highly oscillating region, the non-linear radial wave equation in the phase-amplitude representation is used. In the evaluation of the derivatives of the radial wave function, the potential energy is approximated by the CIP method. In addition, the radial integrals encountered in the computation of various atomic process are accomplished with the CIP method using the values and their analytical derivatives at the grids. This numerical procedure can be extended in a straightforward way to solve the Dirac wave equation.     

A new MHD code with adaptive mesh refinement and parallelization for astrophysics

A new code, named MAP, is written in FORTRAN language for magnetohydrodynamics (MHD) simulations with the adaptive mesh refinement (AMR) and Message Passing Interface (MPI) parallelization. There are several optional numerical schemes for computing the MHD part, namely, modified Mac Cormack Scheme (MMC), Lax–Friedrichs scheme (LF), and weighted essentially non-oscillatory (WENO) scheme. All of them are second-order, two-step, component-wise schemes for hyperbolic conservative equations. The total variation diminishing (TVD) limiters and approximate Riemann solvers are also equipped. A high resolution can be achieved by the hierarchical block-structured AMR mesh. We use the extended generalized Lagrange multiplier (EGLM) MHD equations to reduce the non-divergence free error produced by the scheme in the magnetic induction equation. The numerical algorithms for the non-ideal terms, e.g., the resistivity and the thermal conduction, are also equipped in the code. The details of the AMR and MPI algorithms are described in the paper.