Efficient magnetohydrodynamic simulations on graphics processing units with CUDA

Magnetohydrodynamic (MHD) simulations based on the ideal MHD equations have become a powerful tool for modeling phenomena in a wide range of applications including laboratory, astrophysical, and space plasmas. In general, high-resolution methods for solving the ideal MHD equations are computationally expensive and Beowulf clusters or even supercomputers are often used to run the codes that implemented these methods. With the advent of the Compute Unified Device Architecture (CUDA), modern graphics processing units (GPUs) provide an alternative approach to parallel computing for scientific simulations. In this paper we present, to the best of the author?s knowledge, the first implementation of MHD simulations entirely on GPUs with CUDA, named GPU-MHD, to accelerate the simulation process. GPU-MHD supports both single and double precision computations. A series of numerical tests have been performed to validate the correctness of our code. Accuracy evaluation by comparing single and double precision computation results is also given. Performance measurements of both single and double precision are conducted on both the NVIDIA GeForce GTX 295 (GT200 architecture) and GTX 480 (Fermi architecture) graphics cards. These measurements show that our GPU-based implementation achieves between one and two orders of magnitude of improvement depending on the graphics card used, the problem size, and the precision when comparing to the original serial CPU MHD implementation. In addition, we extend GPU-MHD to support the visualization of the simulation results and thus the whole MHD simulation and visualization process can be performed entirely on GPUs.     

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.     

MINERVA: Ideal MHD stability code for toroidally rotating tokamak plasmas

A new linear MHD stability code MINERVA is developed for investigating a toroidal rotation effect on the stability of ideal MHD modes in tokamak plasmas. This code solves the Frieman-Rotenberg equation as not only the generalized eigenvalue problem but also the initial value problem. The parallel computing method used in this code realizes the stability analysis of both long and short wavelength MHD modes in short time. The results of some benchmarking tests show the validity of this MINERVA code. The numerical study with MINERVA about the toroidal rotation effect on the edge MHD stability shows that the rotation shear destabilizes the intermediate wavelength modes but stabilizes the short wavelength edge localized MHD modes, though the rotation frequency destabilizes both the long and the short wavelength MHD modes.     

High-β and toroidal effects on the internal kink mode in tokamaks

The inclusion of high-β and first-order toroidal terms in the reduced set of (resistive) MHD equations affords the possibility of improving the study of tokamak plasma behavior by three-dimensional numerical simulation. A new code, GALA, based on the reduced equations is developed. It is used to analyse the linear and nonlinear behavior of the internal kink mode in equilibria which are generated by a simple relaxation procedure. We find that the inclusion of toroidal effects in high-β equilibria provides considerable stabilization.     

A.I.K.E.F.: Adaptive hybrid model for space plasma simulations

The interaction between magnetized space plasmas and obstacles like comets, asteroids or planets is determined by a variety of physical processes that occur simultaneously on significantly different length and time scales. Frequently the dynamics of individual ions play a key role for the shape of the interaction region: strong velocity shear between light and heavy plasma constituents, non-Maxwellian particle distributions due to pick up and asymmetries in the magnetic field topology are crucial in determining this type of interaction. Covering these processes is beyond the scope of any Magnetohydrodynamic (MHD) model. In order to account for these effects we have developed a new adaptive hybrid code A.I.K.E.F. (Adaptive Ion-Kinetic Electron-Fluid). The code operates on Cartesian meshes that can adapt to the physical structures in both, space and time. To the authors' knowledge, there is no other adaptive hybrid simulation code in space plasma physics to the present day. Adaptivity is implemented by means of Hybrid-Block-AMR, that is individual octs are refined rather than entire blocks, where an oct is one eighth of a block. In order to account for a reasonable number of particles in each cell, particles are refined via splitting and merging. Both procedures conserve mass, momentum and kinetic energy. The code is implemented in C++ and efficiently parallelized for distributed systems by means of the Message Passing Interface (MPI). In order to demonstrate the validity of our newly developed code we have applied it to a series of fundamental test scenarios. On the one hand we demonstrate that the dispersion relation as well as the propagation characteristics of MHD and whistler mode waves are quantitatively reproduced by our simulation code. Wave propagation remains unaffected when traveling through regions that include different refinement levels. On the other hand we verify that the results obtained on high resolution uniform meshes are identical to the results from adaptive simulations that use coarse base meshes but include various levels of refinement. A remarkable speedup could be observed: the adaptive simulations required 71 times less CPU-hours than the uniform mesh simulations. Finally, we present a first series of global, three-dimensional simulations for the interaction of Mercury with the solar wind and a real time study of Titan's plasma interaction during a magnetosheath excursion.     

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.     

Extension of the Newcomb equation into the vacuum for the stability analysis of tokamak edge plasmas

The formulation for solving numerically the two-dimensional Newcomb equation has been extended to calculate the vacuum energy integral by using a vector potential method. According to this extension, a stability code MARG2D has been adapted, and coded for parallel computing in order to reduce substantially the CPU time. The MARG2D code enables a fast stability analysis of ideal external MHD modes from low to high toroidal mode numbers on the basis of the single physical model, and then the code works as a powerful tool in an integrated simulation where it is combined with transport codes, and also in the analysis of tokamak edge plasma experiments.     

Large-eddy simulation of turbulent flow in a stirred tank driven by a Rushton turbine

Large-eddy simulation (LES) of mixing process in a baffled tank was presented. The impeller rotation was modeled using the sliding mesh technique. In this study the CFD code was used for simulation of a standard vessel agitated by a 6-blade Rushton turbine and results were evaluated in terms of the predicted flow field, power number, mean velocity components, mixing time, turbulent kinetic energy and turbulent dissipation rate using published experimental data. Subsequently, the effects of varying injection position of the passive scalar have been investigated. The results show that LES is a reliable tool to investigate the unsteady behavior of the turbulent flow in stirred tank.     

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.     

Application of Compact Schemes to Large Eddy Simulation of Turbulent Jets

We present 3-D large eddy simulation (LES) results for a turbulent Mach 0.9 isothermal round jet at a Reynolds number of 100,000 (based on jet nozzle exit conditions and nozzle diameter). Our LES code is part of a Computational Aeroacoustics (CAA) methodology that couples surface integral acoustics techniques such as Kirchhoff's method and the Ffowcs Williams– Hawkings method with LES for the far field noise estimation of turbulent jets. The LES code employs high-order accurate compact differencing together with implicit spatial filtering and state-of-the-art non-reflecting boundary conditions. A localized dynamic Smagorinsky subgrid-scale (SGS) model is used for representing the effects of the unresolved scales on the resolved scales. A computational grid consisting of 12 million points was used in the present simulation. Mean flow results obtained in our simulation are found to be in very good agreement with the available experimental data of jets at similar flow conditions. Furthermore, the near field data provided by the LES is coupled with the Ffowcs Williams–Hawkings method to compute the far field noise. Far field aeroacoustics results are also presented and comparisons are made with experimental measurements of jets at similar flow conditions. The aeroacoustics results are encouraging and suggest further investigation of the effects of inflow conditions on the jet acoustic field.     

A parallel,load-balanced MHD code on locally-adapted,unstructured grids in 3d

An efficient parallel code for the approximate solution of initial boundary value problems for hyperbolic balance laws is introduced. The method combines three modern numerical techniques: locally-adaptive upwind finite-volume methods on unstructured grids, parallelization based on non-overlapping domain decomposition, and dynamic load balancing. Key ingredient is a hierarchical mesh in three space dimensions.The proposed method is applied to the equations of compressible magnetohydrodynamics (MHD). Results for several testproblems with computable exact solution and for a realistic astrophysical simulation are shown.     

Modelling of a Utility Boiler Using Parallel Computing

A mathematical model for the simulation of the turbulent reactive flow and heat transfer in a power station boiler has been parallelized. The mathematical model is based on the numerical solution of the governing equations for mass, momentum, energy and transport equations for the scalar quantities. The k- model and the conserved scalar/prescribed probability density function formalism are employed. Radiative heat transfer is calculated using the discrete ordinates method. The code has been fully parallelized using the spatial domain decomposition approach and MPI. Calculations were performed using an IBM-SP2. It is shown that the computational requirements are reduced and the parallel efficiency increases if the mean temperature and density are calculated a priori, and stored. The role of the different parts of the code on the parallel performance is discussed. A speedup of 5.9 is achieved using 8 processors.     

Direct numerical simulation of turbulent flows in fuel rod bundles

Numerical methods for solving conservation equations using the DINUS code are presented. The DINUS code has been developed for direct numerical simulation of thermal-hydraulic phenomena in fuel rod bundles. To examine the methods, two test problems have been studied: turbulent flows between parallel plates and in a Triangular-Arrayed rod bundle.     

Using the dinus code for direct numerical simulation of hydrodynamic processes in VVER-440 fuel rod bundles

The paper presents results of the DINUS code verification against data from experiments conducted at a full-size test facility in Finland. The DINUS code is intended for direct numerical simulation of thermal-hydraulic processes in reactor fuel rod bundles. The experiments studied hydrodynamic processes in the VVER-440 fuel rod bundle at the Loviisa NPP and were carried out for the Reynolds number of 50000. A comparison is given for calculated and experiment data on the axial velocity and flow turbulence intensity distributions at the rod bundle cross sections at different distances from the spacer grid. The DINUS code was used to analyze the effect of the spacer grid on the intensity of inter-subchannel turbulent mixing and convective transfer.     

Large-Eddy Simulation of double-row compound-angle film cooling: Setup and validation

Film cooling is an important technique allowing to increase the thermal efficiency of gas turbines. By blowing cool air through an array of small holes in the turbine blades a thin fluid film is set up shielding the blades from the hot gas arriving from the combustion chamber.This work presents a Large-Eddy Simulation of a particular film-cooling configuration known to provide a high level of effectiveness. It incorporates spanwise rows of holes which, by pairwise combination, generate a so-called anti-kidney-vortex. The simulation setup employs the Navier–Stokes code NSMB for compressible flow including the Approximate Deconvolution Model for subgrid turbulence modeling, the Synthetic-Eddy Method for turbulent inflow generation and reflection-reducing outflow boundary conditions. The setup is first validated by simulations of standard flat-plate turbulent boundary-layer flow without blowing. Results of the film-cooling simulation are then compared with related experimental data. They show reasonable agreement of the cooling effectiveness and temperature distribution, thus confirming the validity of the simulation approach which will be used in future studies of film cooling.     

A high order finite-difference solver for investigation of disturbance development in turbulent boundary layers

This paper outlines a velocity–vorticity based numerical simulation method for modelling perturbation development in laminar and turbulent boundary layers at large Reynolds numbers. Particular attention is paid to the application of integral conditions for the vorticity. These provide constraints on the evolution of the vorticity that are fully equivalent to the usual no-slip conditions. The vorticity and velocity perturbation variables are divided into two distinct primary and secondary groups, allowing the number of governing equations and variables to be effectively halved. Compact finite differences are used to obtain a high-order spatial discretization of the equations. Some novel features of the discretization are highlighted: (i) the incorporation of the vorticity integral conditions and (ii) the related use of a co-ordinate transformation along the semi-infinite wall-normal direction. The viability of the numerical solution procedure is illustrated by a selection of test simulation results. We also indicate the intended application of the simulation code to parametric investigations of the effectiveness of spanwise-directed wall oscillations in inhibiting the growth of streaks within turbulent boundary layers.     

Solution of the time-dependent,three-dimensional resistive magnetohydrodynamic equations

Resistive magnetohydrodynamics (MHD) is described by a set of eight coupled, nonlinear, three-dimensional, time-dependent, partial differential equations. A computer code, IMP (Implicit MHD Program), has been developed to solve these equations numerically by the method of finite differences on an Eulerian mesh. In this model, the equations are expressed in orthogonal curvilinear coordinates, making the code applicable to a variety of coordinate systems. The Douglas-Gunn algorithm for Alternating-Direction Implicit (ADI) temporal advancement is used to avoid the limitations in timestep size imposed by explicit methods. The equations are solved simultaneously to avoid synchronization errors. While the continuity and magnetic field equations are expressed as conservation laws, the momentum and energy equations are nonconservative. This is to: (1) provide enhanced numerical stability by eliminating errors introduced by the nonvanishing of τ · B on the finite difference mesh; and, (2) allow the simulation of low β plasmas. The resulting finite difference equations are a coupled system of nonlinear algebraic equations which are solved by the Newton-Raphson iteration technique. We apply our model to a number of problems of importance in magnetic fusion research. Ideal and resistive internal kink instabilities are simulated in a Cartesian geometry. Growth rates and nonlinear saturation amplitudes are found to be in agreement with previous analytic and numerical predictions. We also simulate these instabilities in a torus, which demonstrates the versatility of the orthogonal curvilinear coordinate representation.     

Normal mode analysis for resistive cylindrical plasmas

The compressible, resistive MHD equations are linearized around an equilibrium with cylindrical symmetry and solved numerically as a complex eigenvalue problem. This normal mode code allows one to solve for very small resistivity γ≈10^-10. The scaling of growth rates and layer width agrees very well with analytical theory. Especially, the influence of both current and pressure on the instabilities is studied in detail; the effect of resistivity on the ideally unstable internal kink is analyzed.     

基于Hadoop分布式计算平台的磁流体动力学模型仿真研究

在Hadoop分布式云计算平台上进行科学计算仿真,具有节省软硬件投资、缩短模拟时间等研究意义。针对需要高计算能力的磁流体动力学(MHD)仿真问题,设计了一种基于Hadoop分布式计算平台的MHD仿真器。首先,将Spark和HAMA两种分布式并行计算模型整合到Hadoop生态系统中,分别用于支持内存计算和整体同步并行计算。然后,将Hadoop集群构建成Master-Slave对等结构,解决全局同步和局部同步问题。最后,在Hadoop集群上,利用有限体积法和黎曼问题来求解MHD方程。实验结果表明,该仿真器能够精确模拟MHD,同时大大缩短了仿真计算时间。