{'fraction_non_alphanumeric': 0.02873239436619718, 'fraction_numerical': 0.011936619718309859, 'mean_word_length': 4.493423597678917, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three Dimensional Pseudo-Spectral Compressible Magnetohydrodynamic GPU Code for Astrophysical Plasma Simulation

This paper presents the benchmarking and scaling studies of a GPU accelerated three dimensional compressible magnetohydrodynamic code. The code is developed keeping an eye to explain the large and intermediate scale magnetic field generation is cosmos as well as in nuclear fusion reactors in the light of the theory given by Eugene Newman Parker. The spatial derivatives of the code are pseudo-spectral method based and the time solvers are explicit. GPU acceleration is achieved with minimal code changes through OpenACC parallelization and use of NVIDIA CUDA Fast Fourier Transform library (cuFFT). NVIDIAs unified memory is leveraged to enable over-subscription of the GPU device memory for seamless out-of-core processing of large grids. Our experimental results indicate that the GPU accelerated code is able to achieve upto two orders of magnitude speedup over a corresponding OpenMP parallel, FFTW library based code, on a NVIDIA Tesla P100 GPU. For large grids that require out-of-core processing on the GPU, we see a 7x speedup over the OpenMP, FFTW based code, on the Tesla P100 GPU. We also present performance analysis of the GPU accelerated code on different GPU architectures - Kepler, Pascal and Volta.


I. INTRODUCTION
Plasmas can be viewed as a charged fluid. In stellar atmosphere, because of strong thermal radiation, the neutral atoms get ionised and form an ionised gas maintaining a quasi-neutrality in space. Thus the state of matter that dominates the structure of a young star is plasma. Because of very high kinetic energy the motion of the charged particles generate a magnetic field and thus the self-consistent motion of the plasma-fluid element gets primarily governed by the Maxwell's equations coupled with the equations of hydrodynamics. The subject that studies the self-consistent evolution of such a magnetised plasma is known as magneto-hydrodynamics (MHD). In laboratories around the globe, there are several Tokamak experiments being carried out to mimic the fusion process in turbulent plasmas taking place in stars. The fundamental motive of such experiments are to harness the binding energy of atoms and transform it to the form of daily usable electricity for the service of mankind. Thus the three dimensional magnetised plasma turbulence is crucial to understand the fundamentals of astro-plasmas present in sun or other young stars as well as the careful operation of complicated Tokamak reactors.
The dynamics of a fully developed turbulent plasma present in stellar objects are fundamentally different from that of a neutral fluid and hence needs a separate careful treatment. The inherent features of this plasma turbulence are needed to be understood and properly analysed to predict any unnatural behavior taking place in our nearest star "The Sun" or in the Tokamak reactors run over the several laboratories.
The large and intermediate scale magnetic field generation in sun as well as in accretion disks around any compact object (for example a neutron star or a black hole) are called "Dynamo Effect" first identified by Eugene Newman Parker. These large or intermediate scale magnetic field primarily governs the dynamics of the charged fluid (plasma) through a time dependent Lorentz force term added in the Navier-Stokes equation, thereby adding an extra dimension (and allowing the exchange of energy between kinetic and magnetic modes), to the already known and well-explored subject "fluid dynamics".
Addressing some of these issues of magnetic field generation and plasma dynamics require, efficient, well parallelized numerical solver capable of handling compressible MHD physics as well as numerical complexities. In order to simulate such a physical process, the finite difference DNS schemes falls short to cope up with the requirement of the high speed performance of the code. Pseudo-spectral technique is observed to be much more accurate and faster method than finite difference methods. We have developed a Three dimensional compressible MHD solver G-MHD3D, which considers continuity, momentum and energy equations for fluid and magnetic variables with a thermodynamic closure for pressure. Code uses conservative form using pseudo-spectral method in cartesian coordinates with GPU parallelization. The numerical solver we have developed captures similar Dynamo Effect for chaotic plasma flows. The code is expected to contribute even more for better understanding of such turbulent phenomena in cosmos.
In this paper, extensive physics benchmarking studies as well as numerical scaling studies is presented along with GPU parallelization, scalability and numerical as well as physics details. It is envisaged that this new code with its extensive physics diagnostics and GPU scalability should be able to address some of the fundamental physics issues in areas such as reconnection (a phenomena of transforming back the magnetic energy to kinetic energy of the plasma), dynamos as well as MHD turbulence.
The OpenMP parallel, FFTW based pseudo-spectral DNS code, MHD3D, has been developed in-house at IPR. The Graphics Processing Unit (GPU) version of the code (G-MHD3D) has been developed with NVIDIA and the performence results have been tested at NVIDIA clusters.

II. CODE G-MHD3D
G-MHD3D is an OpenACC [17] parallel three dimensional compressible, viscous, resistive magnetohydrodynamic GPU code using pseudo-spectral technique to address sponteneous magnetic field generation related problems. The pesudospectral technique uses the cuFFT libraries [20] which are one of the fastest Fourier Transform GPU libraries available till today. This technique is applied to calculate the spatial derivatives and to evaluate the non-linear terms involved in the basic underlying equations. The time derivative is solved using multiple time solvers viz. Adams-Bashforth, Runge-Kutta 4 and Predictor-Correcter algorithms and all the solvers has been checked to provide identical results. The velocity and the magnetic field strength has been normalised with the sound speed and the Alfven speed in the system. The salient features of MHD3D and G-MHD3D are listed in Table I. We develop a new code C-FD-2D, capable of handling compressible fluids, that uses pseudo-spectral technique for spatial discretisation with a standard de-aliasing by zero padding using the 3/2 rule to simulate the above set of equations. The equations are evaluated in two dimensions with periodic boundary conditions in cartesian co-ordinates. We use the multi-dimensional FFTW libraries to evaluate the fourier transforms, whenever needed. The time evolution of the code is performed with Adams-Bashforth, Predictor-Corrector and Runge-Kutta 4 techniques individually and all the three techniques are found to agree in the primary benchmark results. The results in this paper are calculated using Adams-Bashforth method for temporal evolution calculation.
Following we provide the specifications of the basic details and the salient features of our 3D MHD code both for CPU and GPU versions.

III. GOVERNING EQUATIONS
The basic equations that are evolved in the code with different specific initial conditions are as follows:  Mass density: Shear viscosity: Dilation: θ = ∇ · u (4) Equation of state: Ampere's law: Faraday's law: Ohm's law: No magnetic monopole: Putting 6 into 2 we get, Putting 8 into 7 and using 6 & 9 we get, The internal energy is evaluated by the time evolution of the following equation: where E = ρ int + 1 2 ρ| u| 2 + 1 8π | B| 2 and P tot = P + 1 8π | B| 2 .
Thus the complete set of MHD equtations are Hence the complete set of equations in component form for a three dimensional fluid becomes, The forcing determinates a spatial averaged velocity of the system (<| U (x, y, z, t)|> x,y,z ) in the presence of dissipation (ν = 0). We determine a reference scale of pressure (P ref ) and density (ρ ref ) to determine the reference sound speed (C s ref ). This provides a sonic Mach Number at steady state We also evaluate a dynamic Thus we evaluate the dynamic Alfven Mach number as .

A. Normalisation of the MHD equations
Below we provide the normalisation of MHD equations for a special case. We consider the equation of state as P = C 2 s ρ. Here we deal with the equations of unforced turbulence i.e. F = 0. For simplicity, we also assume, the off-diagonal components of the viscosity term to be zero. For the sake of convenience, we consider the Lorentz force term is absent in the momentum equation. The reduced set of equations are Define ρ = ρ 0ρ , t = t 0t and r = L r where, L is the characteristic length-scale and t 0 is a characteristic timescale. Eq. 14 can be written as Defining, u 0 = L t0 , we get, Eq. 15 can be written as, Cs and Re = ρ0u0L ν Eq. 16 can be written as where, R m = Lu0 η . Removing the "˜" sign from the variables, we write the equations as: This normalisation is followed throughout the code.

A. Incompressible hydrodynamic flow
In two dimensions, the Navier-Stokes equation in vorticity formalism become, The Navier-Stokes equation in velocity formalism with compressibility effects become, In order to simulate the two dimensional hydrodynamic fluid, we time evolve the equations in velocity formalism and reproduce the growth rate of Kelvin-Helmholtz instability for a broken jet in the incompressible limit analytically calculated earlier by Drazin [4] and later obtained by Ashwin et al [12] using molecular dynamics simulation. The analytical formula obtained by Drazin [4] in our notation becomes, is the growth rate, k x is the mode number of perturbation, R E is the Raynold's number, U 0 is the maximum fluid velocity, d is the shear width and ν is the coefficient of viscosity. We run our simulation with parameters N x = N y = 512; L x = L y = 2π; ω 0 = 25; d = 3π 128 ; ν = 1.5 × 10 −3 . The initial velocity profile is obtained from the vorticity profile given by ω = ω0 We obtain the identical growth rate from this code with M = 0.05 and another simple code that simulates the Navier-Stokes equation in vorticity formalism with M = 0 [Fig 1c]. We also obtain the same growth rate from another existing, wellbenchmarked code AG-Spect [13] at the limit of negligible visco-elasticity (τ m = 10 −4 ) [Fig 1d]. For this benchmarking the initial condition we have used is ω = We also reproduce Kelvin-Helmholtz instability with a velocity gradient and found the same growth rate given by Ray [5].
A. for hydrodynamic flow In order to simulate hydrodynamic flows we turn off Eq. 19 in the above set of equations. We start with an initial profile E(k) = Ak 4 exp(−2k 2 /k 2 0 ). The RMS (root-mean-square) of the velocity divergence is defined as θ = (∂u i /∂x i ) 2 1/2 and the Skewness is defined as S 3 ≡ <(∂u1/∂x1) 3 > <(∂u1/∂x1) 2 > 3/2 . We reproduce the time evolution of θ and S 3 as previously obtained by Samtaney et al .

B. for magnetohydrodynamic flow
First we turn off Eq. 17 and 18 and start with the initial condition ρ = ρ 0 as a uniform density fluid, the initial velocity profile as u and the initial magnetic field as B x = B y = B z = B 0 . We choose, ρ 0 = 1, U 0 = 1, k f = 1 and A = B = C = 1 as previously taken by Galloway et al. We reproduce the growth rate of the magnetic energy for the specified ABC flow for Reynold's numbers [ Fig. 3] Re = 120, 200 and 450. The exponential growth of magnetic energy of several orders of magnitude is known as "dynamo effect".  In the presence of external forcing the whole system acts as a forced relaxed system. Next, we time evolve all the terms viz. Eq. 17, 18 and 19, in addition with an external forcing defined as [ Fig. 3(h)]: In the absence of forcing, the initial conditions we start with are the following: Fig.4 Fig.4

(d)
We choose M A = 1 such that, the initial kinetic and magnetic energies are equal. The kinetic and magnetic energy oscillates by transferring energies between kinetic and magnetic modes keeping the total energy conserved. In Fig.4 the fall of total energy is solely due to viscous and resistive effects. For further check of the accuracy we increase the Alfven Mach number (M A ) such that the frequency of the Alfven waves increase leaving its signature in the oscillations of the kinetic and magnetic energy. We choose M A = 0.1 and still found consistant decay of total energy, though some oscillations of tiny amplitude appears.
We also reproduce the kinetic and magnetic energy spectra for an incompressible three dimensional magnetohydrodynamic ABC flow where the backreaction of magnetic field on the velocity fields are not considered. The results were obtained earlier by Sadek et al [10].
In order to develop some diagnostics of the 3D MHD code, we calculate the density and velocity structure functions for a initially random velocity and magnetic field. We get back similar structure functions as earlier found by Yang et al [11].
We have developed modules for wavenumber based forcing and gaussian white noise implementation algorithms. We have also incorporated a dynamical wave-number dependent turbulent viscosity model in our code. The coupling in the above set of equations can be altered easily to suit specific purposes.

VI. TRACER PARTICLES
Tracer particles act as a very important diagnostics for any fluid code. In MHD3D code also, tracer particles are added to discover the nature of turbulence that appears in a MHD plasma. The code evaluates the density, velocity and magnetic fields at the grid points using the MHD description of a plasma. The tracer particles are randomly sprinkled in the simulation domain at the beginning of the simulation. They can be placed anywhere in between the grids or at the grid location. These particles are considered to be passive elements of the code. Thus, they do not contribute to the time evolution of any of the fields that are calculated at the grid location from the MHD equations mentioned in the previous section. But the fields do affect the tracers and determine their motion within the simulation box. The initial velocity of the tracers are evaluated from the velocity at the grid. Next, the time evolution of the tracers are fully governed by the magnetic fields at the grids. The grid magnetic field is interpolated at the position of the tracers in between the grids and the Lorentz force acting on them is evaluated. This force pushes the particles to their new positions. The interpolation of the grid magnetic field is done again at the new positions of the tracers and the new force is evaluated. This procedure is repeated to follow the trajectory of the tracer particles in the plasma. Thus it gives us the information on the motion of the fluid elements in a turbulent plasma and one can follow the deterministic path of a tracer which may reveal more information about the underlying symmetries of the plasma flow. Following we describe the interpolation scheme and the particle pusher algorithm in detail.

A. Interpolation Scheme
In order to calculate the magnetic force acting on a tracer particle, we use the following Cloud-In-Cell (CIC) scheme. In this scheme, the magnetic field of eight grid points are interpolated at the particle position. In Fig. 5, we show the volume weighting scheme. The magnetic field is known at all the grid points at the simulaiton domain. The tracer particle divides the unit cube, into eight subshells. One of them (the bottom-left-front) are shown in dark colors (Fig. 5) with blue filled and hollow circles at the corners for better identification. The magnetic field at eight corner-grids (blue and black hollow circles at Fig. 5) are then interpolated at the particle locaton ( B p ) (yellow filled circle at Fig. 5) using the below scheme. Fig. 5: (Color online) CIC scheme for volume weighting for the evaluation of magnetic fields on the location of the tracer particles. The tracer particle (yellow filled circle) is located within the unit shell. It divides the unit shell into 8 subshells, one of them are shown in dark colors with blue circles at the edges. The magnetic field at the tracer particle location ( B p ) is interpolated from the grid magnetic field ( B g ) of the outer corner grids (black hollow circles). The weights of interpolation is given in Eq. 20 If the m th particle is placed within an unit cube whose bottom-left-front corner has a grid index (i, j, k), the magnetic field acting on the m th particle (= B m p ) will be given by, (20) where, the subscript p denotes the particle and subscript g denotes the grid. x i g , y j g , z k g denote the coordinate of the (i, j, k) th grid. Periodic boundary condition (PBC) has been used such that the tracer particle present at one boundary can feel the magnetic field from the opposite boundary.

B. Particle Pusher Scheme
We use Boris Algorithm for pushing the particles from its old to new position due to the effect of Lorentz force acting on the particles. This Lorentz force is evaluated by interpolating the grid magnetic fields to the tracer particle location using the CIC scheme mentioned above. If x m t and v m t are the position and velocity of m th tracer particle respectively at time t, we use the following steps to calculate x m t+dt and v m t+dt at time t + dt.
where, q and M are the charge and mass of the tracer particle respectively. In this code we work in the normalisations where we assume q = M = 1. The particles are guranteed to stay inside the simulation domain by enforcing periodic boundary condition in the position of the tracers.
Thus the following steps are executed to determine the dynamics of the tracers in the code.
1) Determine the grid index (i, j, k) of the bottom-left-front corner grid point of the unit cell in which the m th tracer particle is situated. 2) Use Eq. 20 to evaluate the magnetic field at the particle position. 3) Calculate the Lorentz force on the m th tracer particle using Eq. 21. 4) Determine the new position of the tracer particle using Eq. 22

VII. GPU ACCELERATION
Over the last decade, there has been an increasing use of graphics processing units (GPUs) as general-purpose highly parallel computing units. The GPU is a massively parallel computing device consisting of a large number of streaming multiprocessors (SM), each of which is a set of computing cores. Each SM has a fixed number of registers and a fast onchip memory that is shared among its cores. The different SMs share a slower off-chip memory called the device memory, which is much larger in capacity. For example, the Tesla P100 GPU has 56 SMs each with 64 cores, resulting in a total of 3584 cores. Each SM of the Tesla P100 has a 256 KB register file size and 64 KB of fast on-chip memory. All SMs share a 4 MB L2 cache and 16 GB of device memory.
GPGPU computing refers to general purpose computing on graphics processing units (GPUs). In this context, a GPU is typically called a device, while the CPU is called a host. A kernel refers to a program that is executed on the GPU device using a large number of threads.
OpenACC [17] is an approach for quick and easy porting of existing applications from CPU to GPU. The next section introduces the OpenACC programming paradigm and describes how the MHD3D code was accelerated using the OpenACC pragmas.

A. OpenACC Acceleration of MHD3D
OpenACC is a compiler directive-based programming model, designed to facilitate quick porting of existing applications to accelerators like GPUs, without significant programming effort. OpenACC directives, which are pragmas in C/C++ and specialized comments in Fortran, provide the compiler with additional information, enabling them to optimize the code within the directives for a specific accelerator. With openACC, the same code can be compiled for different accelerators, thereby providing performance portability.
OpenACC provides two types of directives broadly -data directives and compute directives. The data directives provide hints to the compiler on movement of data between the host memory and accelerator memory. The compute directives define regions of the host program that need to be compiled as accelerator code. More details on the OpenACC directives can be found in the OpenACC Programming and Best Practices Guide [18] and the OpenACC specifications document [19].
The MHD3D code has triply nested loops interspersed with fourier transforms. To port the triply nested loop regions to the GPU, we use the OpenACC parallel and loop constructs coupled with the collapse clause [19].
The compiler that supports OpenACC, compiles the loop region into a GPU kernel and distributes the iteration space of the triply nested loop between threads on the GPU. There are clauses that the programmer can use to tune how the loops are distributed amongst the threads. Figure 6 is a snapshot of the compile command and the compiler output. The -acc flag tells the compiler to process the openACC directives, -ta=tesla:XX specifies the architecture for which the accelerator code should be compiled and -Minfo=accel makes the compiler output information corresponding to all the accelerator regions in the code. The compiler output tells the user which loops are transformed into GPU kernels and how the threads are organized. Information about what data is copied into or out of the accelerator is also dumped.

B. cuFFT Library
The MHD3D code, as mentioned earlier, has loop regions interspersed with fourier transforms. To accelerate the fourier transforms on the GPU, we call routines of the cuFFT library. The NVIDIA cuFFT library [20] APIs operate on data in the GPU device memory. Hence, before calling the cuFFT routines, we need to move the relevant data to the  GPU, if not already present, and pass the device pointers to the API calls. The cuFFT plan creation (this is needed for initializing some state information that is necessary before calling the fft routines) is done once in the beginning and destroyed at the end. Figure 7 shows the steps needed to call the cuFFT routines. To call the cuFFT library routines from our OpenACC Fortran application, we create C wrapper modules as shown in Figure 8. In-place transforms are used, i.e the output transforms are written to the input arrays. This reduces the memory footprint (amount of memory used by the application).

C. Out-of-core Processing
While computing the MHD3D equations for large 3D grids of size 512x512x512, we observe that the data operated on within the triply nested loops (aka working set) exceeds the device memory capacity of the GPU. Each variable array in this case is 1 GB in size. This results in the need for efficient out-of-core processing. Out-of-core processing refers to processing of data that does not completely fit within the memory of a computing system. One approach to out-of-core processing is to split the data into blocks (also called tiles), move one tile to the GPU memory and completely process it before moving it back to the host memory. However, since the solver code is interspersed with fourier transforms that need the complete data, it is not possible to do a tile-by-tile processing.
For seamless out-of-core processing without programmer intervention, in G-MHD3D, we leverage NVIDIA's unified memory [21]. Unified memory is a component of NVIDIA's GPU programming model that defines a managed memory space in which all processors (CPU and GPU) see a single coherent memory image with a common address space. This means that there is no need to explicitly move and manage data between the host memory and GPU. For example, when a variable is accessed on the GPU within a kernel and the page that holds the variable is not present on the GPU device, (a page is a contiguous block of virtual memory), a page fault occurs. The system transparently handles this and moves the page to the GPU device if necessary. This is called on-demand page migration.
Unified memory also facilitates over-subscription of the GPU memory. This means that you can use variables whose sizes exceed the physical memory capacity of the GPU and the system transparently manages moving parts of data between host and GPU memory as needed.
To enable unified memory, you need to compile and link your application with the -ta:tesla=managed flag.

VIII. PERFORMANCE RESULTS
In this section, we present the performance analysis of the GPU accelerated 3D compressible magnetohydrodynamic code on different GPU architectures. The accelerated code, as specified in Section VII, uses OpenACC [17] for GPU parallelization. In-place fourier transforms are computed on the GPU using the cuFFT library [20]. Computing the transforms in-place reduces the memory footprint, resulting in reduced data transfers between the host and GPU and improved performance. For seamless out-of-core processing of large grids that do not fit in the GPU memory, we leverage NVIDIA's unified memory [21].
Experiments are run on a Intel Xeon E5-2698 v3 16 core, dual socket CPU @2.3 GHz, 256 GB RAM and the NVIDIA Tesla K80 (Kepler) GPU, NVIDIA Tesla P100 (Pascal) and NVIDIA Tesla V100 (Volta) GPUs. We use cuda toolkit 9.0.176 and pgi 17.9. Table II shows the execution time of the solver per timestep for grid sizes of 64X64X64 and 128X128X128 on the CPU and the Kepler, Pascal and Volta GPUs. The CPU version of the code is developed using OpenMP parallelization and the FFTW library. We see that the GPU code is about two orders of magnitude faster than the CPU code.
Across GPU architectures, we see that Tesla P100 is 2x to 3x faster than the Tesla K80. This is due to the difference in the memory bandwidth between the two architectures. Tesla P100 has a theoretical peak memory bandwidth of 720 GB/s which is 3x that of Tesla K80 whose peak memory bandwidth is 240 GB/s. As the MHD3D code is memory intensive, we can see 2x-3x performance difference between the two architectures. Performance results over the Pascal and Volta architecture are quite similar (within 30%) as there is not a great difference in their memory bandwidths (720 GB/s vs 900 GB/s). To verify the correctness of the GPU accelerated solver, we compared the CPU and GPU outputs and computed the mean square error. This was found to be of the order of 5.89E-15.
For larger grid sizes which do not fit in GPU memory, we use unified memory. On demand page migration and oversubscription features of unified memory are not supported on     Table III shows the execution time per timestep for large grids. For 256X256X256 grid size, the total memory footprint of the solver exceeds the GPU device memory of 16 GB. However, the working set of each kernel is within the GPU memory capacity. When the grid size becomes 512X512X512, the working set of the GPU kernels is in some cases 40GB and does not fit completely within the GPU memory. As mentioned in Section VII-C, tiling the data typically helps in such a scenario, but since the solver code is interspersed with fourier transforms that need the complete data, it is not possible to do a tile-by-tile processing.
Therefore, in this case, there is frequent movement of pages back and forth between the CPU and GPU and a large number of page faults on the GPU. To reduce the page faults on the GPU (GPU page fault handling is expensive), we ensure that the data arrays on the CPU are pre-mapped on the GPU. This avoids a page fault on the first access. Eventually the page may migrate to the GPU (based on the unified memory system heuristics), but the mappings on the GPU page tables are kept updated.
Even with the above optimizations, the out-of-core solver is an order magnitude slower than the in-core solver, largely due to movement of pages between host and GPU and handling of the GPU page faults. Even with this additional overhead, the GPU accelerated code is able to show a speedup of 7x over the OpenMP version.

IX. FUTURE SCOPE
The implementation of some diagnostics viz. field line tracer, poincare section, modo-mode energy transfer in a triad [14] is in progress. The performance of the code can be improved by implementing implicit time solvers. In order to further enhance the speed of the code, the GPGPU parallel cuda version of the code is being developed. At present the code can handle only periodic boundaries. The simulation of plasma flows in bounded domains with upgraded version of the code is under development.