Search

Maintaining conservation laws in the fully discrete setting is critical for accurate long-time behavior of numerical simulations and requires accounting for discrete conservation properties in both space and time. This paper derives arbitrary order finite element exterior calculus spatial discretizations for the two-dimensional (2D) Navier-Stokes and drift-reduced magnetohydrodynamic equations that conserve both energy and enstrophy to machine precision when coupled with generally symplectic time-integration methods. Both continuous and discontinuous-Galerkin (DG) weak formulations can ensure conservation, but only generally symplectic time integration methods, such as the implicit midpoint method, permit exact conservation in time. Moreover, the symplectic implicit midpoint method yields an order of magnitude speedup over explicit schemes. The methods are implemented using the MFEM library and the solutions are verified for an extensive suite of 2D neutral fluid turbulence test problems. Numerical solutions are verified via comparison to a semi-analytic linear eigensolver as well as to the finite difference Global Drift Ballooning (GDB) code. However, it is found that turbulent simulations that conserve both energy and enstrophy tend to have too much power at high wavenumber and that this part of the spectrum should be controlled by reintroducing artificial dissipation. The DG formulation allows upwinding of the advection operator which dissipates enstrophy while still maintaining conservation of energy. Coupling upwinded DG with implicit symplectic integration appears to offer the best compromise of allowing mid-range wavenumbers to reach the appropriate amplitude while still controlling the high-wavenumber part of the spectrum.

We present an extension of a linearized Coulomb collision operator, previously used in several Eulerian kinetic codes for like-species collisions and unlike-species collisions in the case where the backgrounds about which the linearization is made all are in collisional equilibrium, to the situation of interpenetrating plasma streams. In the latter case, the backgrounds cannot be taken to be in equilibrium and a significant generalization is required. Our development is targeted toward the Eulerian kinetic plasma code LOKI, which evolves the Vlasov-Poisson or Vlasov-Maxwell system in a Cartesian "2 + 2-dimensional" phase space. The extended operator has been implemented in a test code, and results of both quantitative verification and qualitative "realizability" tests are presented.

The interpenetration of counterstreaming plasmas is an important phenomenon in several application areas, such as astrophysical flows, design of controlled fusion devices, and laser–induced plasma experiments. Multispecies “single-fluid” codes are unable to model this phenomenon due to the single velocity representation for all the species/fluids. Kinetic codes, though capable of modeling interpenetration, are computationally prohibitive for at-scale simulations. In this paper, we propose a multifluid model that solves the fluid equations for each ion fluid or stream. This allows distinct flows that interact with each other through electrostatic and collisional forces. We introduce and describe our code, EUCLID, that uses a conservative finite-difference formulation to discretize the governing equations in space. The 5th-order Monotonicity-Preserving WENO scheme is used for the upwind approximation of the hyperbolic flux, and the explicit 4th-order Runge–Kutta scheme is used for time integration. The code is verified for several benchmark cases and manufactured solutions. We simulate one- and two-dimensional interactions of counterstreaming plasmas in vacuum as well as in the presence of gas fill, where the setups are representative of laser-induced plasma experiments.

Helicon waves have been recently proposed as an off-axis current drive actuator due to their expected high current drive efficiency in the mid-radius region in high beta tokamaks. This paper focuses on a numerical study to better understand effects of scrape-off-layer (SOL) turbulence on helicon wave propagation and absorption on the DIII-D tokamak using a recently developed helicon full-wave model with turbulent density inputs from synthetic single wavelength SOL turbulence and first-principles HERMES multi-wavelength turbulence models. With both input turbulence models, three key effects are observed: the helicon wave can scatter to undesirable locations in the plasma, large helicon wave electric fields can form in localized regions near the SOL turbulence, and the helicon wave can mode convert to slow waves in the SOL. This is shown to cause helicon wave refraction to undesirable locations and strong helicon wave absorption in the SOL resulting in significantly less helicon wave power in the core plasma. Using synthetic SOL turbulence, the simulations additionally show that high amplitudes and long wavelengths greater than a few cm on average have the largest effect on modifying the helicon wave propagation and absorption; the modeling predicts, for example, that approximately 60% of helicon power can be absorbed in the SOL for ñ/n ∼ 0.8 and lambda_perp ∼ 0.05 m. Several potential physical mechanisms that may explain the interaction of helicon waves with SOL turbulence in these simulations are discussed.

The nonlinear evolution of the ion-ion streaming instability (IISI) is studied using numerical techniques novel to this problem that afford direct insight into the evolution of the particle distributions of each species. During the linear phase of the instability, we demonstrate quantitative agreement with linear kinetic theory. Subsequently, the electrostatic field generated by the IISI causes ring-like velocity distributions of ions to form that are both heated and slowed to varying degrees relative to their initial flows. Due to variation in the trapping conditions for ion species of differing charge-to-mass ratio, when flows of multiple species interact, the nonlinear evolution of each species can be starkly different: we show a case where a lighter ion species is completely stopped by a heavier ion species via the IISI alone (i.e., without collisions) and, for the first time, demonstrate how the IISI can introduce a relative flow between ion species that initially have the same flow velocities, thereby separating them.

The parallel electron heat transport in a weakly collisional regime can be represented in the framework of the Landau-fluid model (Hammett et al., 1990). Practical implementation of Landau-fluid transport has become possible due to the recent invention of an efficient non-spectral method for the non-local closure operators (Dimits et al., 2014). Here the implementation of a Landau-fluid based model for the parallel plasma transport is described, and the model is tested for different collisionality regimes against Fokker–Planck simulations. The new method appears to represent the weakly collisional electron transport more accurately than the conventional flux-limiter based models, on the other hand it is computationally efficient enough to be incorporated in comprehensive edge plasma simulations.

The coupled fluid plasma and kinetic neutral physics equations are analyzed through theory and simulation of benchmark cases. It is shown that coupling methods that do not treat the coupling rates implicitly are restricted to short time steps for stability. Fast charge exchange, ionization and recombination coupling rates exist, even after constraining the solution by requiring that the neutrals are at equilibrium. For explicit coupling, the present implementation of Monte Carlo correlated sampling techniques does not allow for complete convergence in slab geometry. For the benchmark case, residuals decay with particle number and increase with grid size, indicating that they scale in a manner that is similar to the theoretical prediction for nonlinear bias error. Progress is reported on implementation of a fully implicit Jacobian-free Newton–Krylov coupling scheme. The present block Jacobi preconditioning method is still sensitive to time step and methods that better precondition the coupled system are under investigation. Keywords: Divertor modeling, Charge exchange, Ionization, Recombination, Implicit, Newton–Krylov, UEDGE, DEGAS2, EIRENE, SOLPS, MSC: 00-01, 99-00

A recently introduced algorithm for performing integrated kinetic simulations of the fast ignition approach to laser fusion is reviewed. The integrated algorithm uses a conventional fully electromagnetic relativistic particle-in-cell algorithm in vacuum and plasma up to a density well above cutoff for the incident laser, above which density the plasma is sufficiently collisional so that light waves and electron plasma waves are unimportant and a simplified physics model is justified that leads to improved computational efficiencies. Integrated comprehensive kinetic simulation of fast ignition is rendered more practical with this two-region algorithm. Professor C. K. Birdsall was a pioneer in computational plasma physics and enthusiastically championed the expansion of its use for discovery science and to simulate plasma phenomena in a wall-to-wall manner.

In diverse areas of science and technology, including inertial confinement fusion (ICF), astrophysics, geophysics, and engineering processes, turbulent mixing induced by hydrodynamic instabilities is of scientific interest as well as practical significance. Because of the fundamental roles they often play in ICF and other applications, three classes of hydrodynamic instability-induced turbulent flows—those arising from the Rayleigh-Taylor, Richtmyer-Meshkov, and Kelvin-Helmholtz instabilities—have attracted much attention. ICF implosions, supernova explosions, and other applications illustrate that these phases of instability growth do not occur in isolation, but instead are connected so that growth in one phase feeds through to initiate growth in a later phase. Essentially, a description of these flows must encompass both the temporal and spatial evolution of the flows from their inception. Hydrodynamic instability will usually start from potentially infinitesimal spatial perturbations, will eventually transition to a turbulent flow, and then will reach a final state of a true multiscale problem. Indeed, this change in the spatial scales can be vast, with hydrodynamic instability evolving from just a few microns to thousands of kilometers in geophysical or astrophysical problems. These instabilities will evolve through different stages before transitioning to turbulence, experiencing linear, weakly, and highly nonlinear states. The challenges confronted by researchers are enormous. The inherent difficulties include characterizing the initial conditions of such flows and accurately predicting the transitional flows. Of course, fully developed turbulence, a focus of many studies because of its major impact on the mixing process, is a notoriously difficult problem in its own right. In this pedagogical review, we will survey challenges and progress, and also discuss outstanding issues and future directions.

The extension of Langevin-equation Monte-Carlo algorithms for Coulomb collisions from the conventional Euler-Maruyama time integration to the next higher order of accuracy, the Milstein scheme, has been developed, implemented, and tested. This extension proceeds via a formulation of the angular scattering directly as stochastic differential equations in the fixed-frame spherical-coordinate velocity variables. Results from the numerical implementation show the expected improvement [O(@Dt) vs. O(@Dt^1^/^2)] in the strong convergence rate both for the speed |v| and angular components of the scattering. An important result is that this improved convergence is achieved for the angular component of the scattering if and only if the ''area-integral'' terms in the Milstein scheme are included. The resulting Milstein scheme is of value as a step towards algorithms with both improved accuracy and efficiency. These include both algorithms with improved convergence in the averages (weak convergence) and multi-time-level schemes. The latter have been shown to give a greatly reduced cost for a given overall error level when compared with conventional Monte-Carlo schemes, and their performance is improved considerably when the Milstein algorithm is used for the underlying time advance versus the Euler-Maruyama algorithm. A new method for sampling the area integrals is given which is a simplification of an earlier direct method and which retains high accuracy. This method, while being useful in its own right because of its relative simplicity, is also expected to considerably reduce the computational requirements for the direct conditional sampling of the area integrals that is needed for adaptive strong integration.

Both binary and grid-based Langevin equations models for Coulomb collisions are used in particle simulation of plasmas. We introduce a variant of the conventional binary collision algorithm for performing Coulomb collisions. In this algorithm particles in a configuration space cell are not paired for collisions. Instead, for every test particle in the cell, a unique field particle is defined by randomly sampling a velocity distribution defined on the grid by accumulating moments of the particle distribution function(s). The test and field particle pair then undergoes a collision using the standard methodology for binary collisions. The performance of the new algorithm is illustrated in example computations and compared to a drag-diffusion Langevin equations algorithm. The grid-based algorithms do not conserve momentum and energy, although with good particle statistics the non-conservation is relatively small. Conservation can be restored after collisions using a shift and scaling of the momenta. The comparative merits of the new algorithm are discussed.

The accuracy of first-order Euler and higher-order time-integration algorithms for grid-based Langevin equations collision models in a specific relaxation test problem is assessed. We show that statistical noise errors can overshadow time-step errors and argue that statistical noise errors can be conflated with time-step effects. Using a higher-order integration scheme may not achieve any benefit in accuracy for examples of practical interest. We also investigate the collisional relaxation of an initial electron-ion relative drift and the collisional relaxation to a resistive steady-state in which a quasi-steady current is driven by a constant applied electric field, as functions of the time step used to resolve the collision processes using binary and grid-based, test-particle Langevin equations models. We compare results from two grid-based Langevin equations collision algorithms to results from a binary collision algorithm for modeling electron-ion collisions. Some guidance is provided on how large a time step can be used compared to the inverse of the characteristic collision frequency for specific relaxation processes.

We investigate the accuracy of and assumptions underlying the numerical binary Monte Carlo collision operator due to Nanbu [K. Nanbu, Phys. Rev. E 55 (1997) 4642]. The numerical experiments that resulted in the parameterization of the collision kernel used in Nanbu's operator are argued to be an approximate realization of the Coulomb-Lorentz pitch-angle scattering process, for which an analytical solution for the collision kernel is available. It is demonstrated empirically that Nanbu's collision operator quite accurately recovers the effects of Coulomb-Lorentz pitch-angle collisions, or processes that approximate these (such interspecies Coulomb collisions with very small mass ratio) even for very large values of the collisional time step. An investigation of the analytical solution shows that Nanbu's parameterized kernel is highly accurate for small values of the normalized collision time step, but loses some of its accuracy for larger values of the time step. Careful numerical and analytical investigations are presented, which show that the time dependence of the relaxation of a temperature anisotropy by Coulomb-Lorentz collisions has a richer structure than previously thought, and is not accurately represented by an exponential decay with a single decay rate. Finally, a practical collision algorithm is proposed that for small-mass-ratio interspecies Coulomb collisions improves on the accuracy of Nanbu's algorithm.

The interactions of charged particles in a plasma are governed by long-range Coulomb collision. We compare two widely used Monte Carlo models for Coulomb collisions. One was developed by Takizuka and Abe in 1977, the other was developed by Nanbu in 1997. We perform deterministic and statistical error analysis with respect to particle number and time step. The two models produce similar stochastic errors, but Nanbu's model gives smaller time step errors. Error comparisons between these two methods are presented.

The turbulent Rayleigh–Taylor instability is investigated in the limit of strong mode-coupling using a variety of high-resolution, multimode, three dimensional numerical simulations (NS). The perturbations are initialized with only short wavelength modes so that the self-similar evolution (i.e., bubble diameter Db∝amplitude hb) occurs solely by the nonlinear coupling (merger) of saturated modes. After an initial transient, it is found that hb∼αbAgt2, where A=Atwood number, g=acceleration, and t=time. The NS yield Db∼hb/3 in agreement with experiment but the simulation value αb∼0.025±0.003 is smaller than the experimental value αb∼0.057±0.008. By analyzing the dominant bubbles, it is found that the small value of αb can be attributed to a density dilution due to fine-scale mixing in our NS without interface reconstruction (IR) or an equivalent entrainment in our NS with IR. This may be characteristic of the mode coupling limit studied here and the associated αb may represent a lower bound that is insensiti...

A new electromagnetic kinetic electron simulation model that uses a generalized split-weight scheme, where the adiabatic part is adjustable, along with a parallel canonical momentum formulation has been developed in three-dimensional toroidal flux-tube geometry. This model includes electron–ion collisional effects and has been linearly benchmarked. It is found that for H-mode parameters, the nonadiabatic effects of kinetic electrons increase linear growth rates of the ion-temperature-gradient-driven (ITG) modes, mainly due to trapped-electron drive. The ion heat transport is also increased from that obtained with adiabatic electrons. The linear behaviour of the zonal flow is not significantly affected by kinetic electrons. The ion heat transport decreases to below the adiabatic electron level when finite plasma β is included due to finite-β stabilization of the ITG modes. This work is being carried out using the 'Summit Framework'.

Turbulent hydrodynamic mixing induced by the Rayleigh–Taylor (RT) and Richtmyer–Meshkov (RM) instabilities occurs in settings as varied as exploding stars (supernovae), inertial confinement fusion (ICF) capsule implosions, and macroscopic flows in fluid dynamics facilities such as shock tubes. Turbulence theory and modeling have been applied to RT and RM induced flows and developed into a quantitative description of turbulence from the onset to the asymptotic end-state. The treatment, based on a combined approach of theory, direct numerical simulation (DNS), and experimental data analysis, has broad generality. Three areas of progress will be reported. First, a robust, easy to apply criteria will be reported for the mixing transition in a time-dependent flow. This allows an assessment of whether flows, be they from supernova explosions or ICF experiments, should be mixed down to the molecular scale or not. Second, through DNS, the structure, scaling, and spectral evolution of the RT instability induced fl...

Three-dimensional high-resolution simulations (up to 8 billion zones) have been performed for a Richtmyer–Meshkov instability produced by passing a shock through a contact discontinuity with a two-scale initial perturbation. The setup approximates shock-tube experiments with a membrane pushed through a wire mesh. The simulation produces mixing-layer widths similar to those observed experimentally. Comparison of runs at various resolutions suggests a mixing transition from unstable to turbulent flow as the numerical Reynolds number is increased. At the highest resolutions, the spectrum exhibits a region of power-law decay, in which the spectral flux is approximately constant, suggestive of an inertial range, but with weaker wave number dependence than Kolmogorov scaling, about k−6/5. Analysis of structure functions at the end of the simulation indicates the persistence of structures with velocities largest in the stream-wise direction. Comparison of three-dimensional and two-dimensional runs illustrates th...

An electromagnetic hybrid scheme (fluid electrons and gyrokinetic ions) is elaborated in example calculations and extended to toroidal geometry. The scheme includes a kinetic electron closure valid for βe>me/mi (βe is the ratio of the plasma electron pressure to the magnetic field energy density). The new scheme incorporates partially linearized (δf ) drift-kinetic electrons whose pressure and number density moments are used to close the fluid momentum equation for the electron fluid (Ohm’s law). The test cases used are small-amplitude kinetic shear-Alfven waves with electron Landau damping, the ion-temperature-gradient instability, and the collisionless drift instability (universal mode) in an unsheared slab as a function of the plasma βe. Attention is given to resolution and convergence issues in simulations of turbulent steady states.

The non-linear 3-D toroidal gyrokinetic simulation code PG3EQ is used to study toroidal ion temperature gradient (ITG) driven turbulence - a key cause of the anomalous transport that limits tokamak plasma performance. Surveys of ?i versus E ? B and toroidal shear show that (a)?the maximum growth rate is not a good predictor of the E ? B shear required to suppress turbulent transport, (b)?there is often a `plateau region' in which E ? B shear significantly reduces the maximum linear growth rate but not the transport and (c)?the parallel velocity shear component of toroidal velocity shear can negate much of the transport reduction by the E ? B shear. Simulations in which the ion temperature gradient Ti(r) initially varies with radial position evolve towards a state without strong radial variations in Ti(r) while approximately preserving the total (E ? B + diamagnetic) ion flow. If the electrostatic potential is initialized to zero, then sheared E ? B flows result which can significantly reduce or quench the transport. If, however, the total ion flow profile (or equivalently the radial force) is initialized to be radially uniform, then the initial radial variations in Ti(r) do not result in a significant reduction in the transport. Surveys of ?i versus magnetic shear S show a (often very sharp) peak at S 0.5 - the value at which the orientation of the ITG modes is in the major radius direction and is independent of poloidal angle in the region near the outer midplane. Surveys of ?i versus Ti show that in most cases the thermal flux Q (?i Ti) has a linear dependence on Ti, with an intercept value of Ti which is often significantly larger than than linear critical value.

Results are presented from non-linear gyrokinetic simulations of toroidal ion temperature gradient turbulence and transport. The ion thermal fluxes are found to have an offset linear dependence on the temperature gradient and are significantly lower than gyrofluid or IFS-PPPL model predictions. A new phenomenon of non-linear effective critical gradients larger than the linear instability threshold gradients is observed and is associated with undamped flux surface averaged shear flows. The non-linear gyrokinetic codes have passed extensive tests, including comparison against independent linear calculations, a series of non-linear convergence tests and a comparison between two independent non-linear gyrokinetic codes. The most realistic simulations to date used actual reconstructed equilibria from experiments and a model for dilution by impurity and beam ions. These simulations highlight the importance of both self-generated and external E × B flow shear as well as the need for still more physics to be included.

Implicit, perturbative algorithms [B. I. Cohen, A. M. Dimits, J. J. Stimson, and D. C. Barnes, Phys. Rev. E {bold 53}, 2708 (1996)] are extended to address electromagnetic particle simulation of drift-wave instabilities in electron and ion plasmas with a strong applied magnetic field. The algorithms use gyrokinetic ions and drift-kinetic electrons with a realistic mass ratio. Implicit moment and direct implicit solutions of the coupled Maxwell and plasma dynamical equations are demonstrated. Aspects of the numerical dispersion properties are confirmed in test simulations. Applications to ion-temperature-gradient and collisionless drift instabilities in a two-dimensional slab with magnetic shear are presented that illustrate the efficacy of the algorithms and the electromagnetic effects. {copyright} {ital 1997} {ital The American Physical Society}

We present a new, for plasma physics, highly efficient multilevel Monte Carlo numerical method for simulating Coulomb collisions. The method separates and optimally minimizes the finite-timestep and finite-sampling errors inherent in the Langevin representation of the Landau-Fokker-Planck equation. It does so by combining multiple solutions to the underlying equations with varying numbers of timesteps. For a desired level of accuracy epsilon, the computational cost of the method is order(epsilon^{-2}) or order(epsilon^{-2} (\ln epsilon)^2), depending on the underlying discretization, Milstein or Euler-Maruyama respectively. This is to be contrasted with a cost of order(epsilon^{-3}) for direct simulation Monte Carlo or binary collision methods. We successfully demonstrate the method with a classic beam diffusion test case in 2D, making use of the Levy area approximation for the correlated Milstein cross terms, and generating a computational saving of a factor of 100 for epsilon = 10^{-5}. We discuss the importance of the method for problems in which collisions constitute the computational rate limiting step, and its limitations., 32 pages

Summary form only given. Fully electromagnetic, relativistic particle-in-cell (PIC) simulation of simulation of laser-plasma interactions has a forty-year history. The simulation of the fast ignition approach to laser fusion presents extreme challenges to traditional explicit algorithms because of the multi-scale nature of the plasma - the plasma density ranges from vacuum to near solid density. Here we describe a PIC-fluid hybrid algorithm that has proven successful in modeling laserplasma interactions in fast-ignition laser fusion.1 In this hybrid algorithm a traditional fully electromagnetic, relativistic, explicit PIC method is used in the vacuum region and in the plasma up to a density much greater than the critical density for the laser, e.g., ne/nc~100. For higher plasma densities (as long as the plasma remains sufficiently cold and collisional) an Ohm's law is used to determine the electric field in the plasma, and inclusion of the displacement current in Ampere's law becomes optional, which eliminates light waves and electron plasma waves as justified by the high collisionality. This approach achieves efficiencies by eliminating requirements on resolving the electron Debye length, the skin depth and the plasma frequency in the high-density region. A particle representation of the plasma electrons and ions is retained throughout the simulation domain. An important aspect of this hybrid model is the inclusion of an accurate and efficient Coulomb collision operator that is compatible with unequally weighted particles, which can be achieved with a variant of the Takizuka-Abe (or Nanbu) collision algorithm extended to a grid-based model.

This paper reports on the theoretical and simulation results of a gyro-Landau-fluid extension of the BOUT++ code, which contributes to increasing the physics understanding of edge-localized-modes (ELMs). Large ELMs with low-to-intermediate-n peeling-ballooning (P-B) modes are significantly suppressed due to finite Larmor radius (FLR) effects when the ion temperature increases. For type-I ELMs, it is found from linear simulations that retaining complete first order FLR corrections as resulting from the incomplete "gyroviscous cancellation" in Braginskii's two-fluid model is necessary to obtain good agreement with gyro-fluid results for high ion temperature cases (Ti3 keV) when the ion density has a strong radial variation, which goes beyond the simple local model of ion diamagnetic stabilization of ideal ballooning modes. The maximum growth rate is inversely proportional to Tibecause the FLR effect is proportional to Ti. The FLR effect is also proportional to toroidal mode number n, so for high n cases, the P-B mode is stabilized by FLR effects. Nonlinear gyro-fluid simulations show results that are similar to those from the two-fluid model, namely that the P-B modes trigger magnetic reconnection, which drives the collapse of the pedestal pressure. Due to the additional FLR-corrected nonlinear E × B convection of the ion gyro-center density, for a ballooning-dominated equilibrium the gyro-fluid model further limits the radial spreading of ELMs. In six-field two fluid simulations, the parallel thermal diffusivity is found to prevent the ELM encroachment further into core plasmas and therefore leads to steady state L-mode profiles. The simulation results show that most energy is lost via ion channel during an ELM event, followed by particle loss and electron energy loss. Because edge plasmas have significant spatial inhomogeneities and complicated boundary conditions, we have developed a fast non-Fourier method for the computation of Landau-fluid closure terms based on an accurate and tunable approximation. The accuracy and the fast computational scaling of the method have been demonstrated. © 2013 AIP Publishing LLC.

We formulate and test a hybrid fluid-Monte Carlo scheme for the treatment of elastic collisions in gases and plasmas. While our primary focus and demonstrations of applicability are for moderately collisional plasmas, as described by the Landau-Fokker-Planck equation, the method is expected to be applicable also to collision processes described by the Boltzmann equation. This scheme is similar to the previously discussed velocity-based scheme [R. Caflisch, et. al, Multiscale Modeling & Simulation 7, 865, (2008)] and the scattering-angle-based scheme [A.M. Dimits, et. al, Bull. APS 55, no. 15 (2010, Abstract: XP9.00006)], but with a firmer theoretical basis and without the inherent limitation to the Landau-Fokker-Planck case. It gives a significant performance improvement (e.g., error for a given computational effort) over the velocity-based scheme. These features are achieved by assigning passive scalars to each simulated particle and tracking their evolution through collisions. The method permits a detailed error analysis that is confirmed by numerical results. The tests performed are for the evolution from anisotropic Maxwellian and a bump-on-tail distribution.

Results are presented from the first systematic nonlinear gyrokinetic simulation study of the scalings and parameter dependences of toroidal ion-temperature-gradient-driven (ITG) turbulence and transport in large magnetic-fusion-relevant tokamaks. Such simulations are made possible by significant advances in computational physics algorithms which are summarized. The simulated ion-thermal transport rates have clear ``gyroreduced Bohm'' scaling and are too low to account for the anomalous transport seen in the outer half of some TFTR $L$-mode discharges, indicating that additional mechanisms, such as nonadiabatic-electron effects, are operative.

An algorithm for kinetic simulation of electron and ion plasmas is introduced that combines both implicit and perturbative methods to deal efficiently with disparate time scales and signal-to-noise issues associated with particle discreteness. The algorithm is introduced for an electrostatic model in a strong magnetic field with gyrokinetic ions and drift-kinetic electrons. Some of the numerical dispersion properties are analyzed and confirmed in test simulations. Two-dimensional applications of the algorithm to drift-wave instabilities excited either by an ion temperature gradient (ion-temperature-gradient instability) or by electron inverse Landau damping in the presence of a density gradient (collisionless drift instability) are presented that illustrate the usefulness of the method. Clear gains in resolution and computational efficiency over previous methods are achieved. {copyright} {ital 1996 The American Physical Society.}

A series of DIII-D [J. L. Luxon, Nucl. Fusion 42, 614 (2002)] low torque quiescent H-mode experiments show that density gradient driven trapped electron mode (DGTEM) turbulence dominates the inner core of H-mode plasmas during strong electron cyclotron heating (ECH). Adding 3.4 MW ECH doubles Te/Ti from 0.5 to 1.0, which halves the linear DGTEM critical density gradient, locally reducing density peaking, while transport in all channels displays extreme stiffness in the density gradient. This suggests that fusion α-heating may degrade inner core confinement in H-mode plasmas with moderate density peaking and low collisionality, with equal electron and ion temperatures, key conditions expected in burning plasmas. Gyrokinetic simulations using GYRO [J. Candy and R. E. Waltz, J. Comput. Phys. 186, 545 (2003)] (and GENE [Jenko et al., Phys. Plasmas 7, 1904 (2000)]) closely match not only particle, energy, and momentum fluxes but also density fluctuation spectra from Doppler backscattering (DBS), with and witho...

In this paper, a numerical instability first observed in a 2D electrostatic gyrokinetic code is described. The instability should also be present in some form in many versons of particle-in-cell simulation codes that employ guiding center drifts. A perturbation analysis of the instability is given and its results agree quantitatively with the observations from the gyrokinetic code in all respects. The basic mechanism is a false divergence of the E × B flow caused by the interpolation between the grid and the particles as coupled with the specific numerical method for calculating E - ??. Stability or instability depends in detail on the specific choice of particle interpolation method and field method. One common interpolation method, subtracted dipole, is stable. Other commonly used interpolation methods, linear and quadratic, are unstable when combined with a finite difference for the electric field. Linear and quadratic interpolation can be rendered stable if combined with another method for the electric field, the analytic differential of the interpolated potential.

Algorithms are presented for energy- and momentum-conserving like-particle Coulomb collisions in partially linearized ([delta][ital f]) particle simulations. They are developed and implemented in particular for gyrokinetic simulation models of a strongly magnetized plasma. The collision operators include both drag and diffusion terms, are not restricted to a single or few Fourier modes, and approximately conserve both momentum and energy locally in space in a statistical sense. The first algorithm is a many-mode generalization of a test-particle--plus--source algorithm previously proposed. The second is easier to implement and improves upon the first significantly by not requiring many time steps for good conservation. Implementations for the case for ion-ion collisions are given and conservation properties are demonstrated, both directly with non-self-consistent test simulation runs and indirectly with self-consistent runs. The computational cost of particle pushing and solving for fields depends on the relative collisionality and can result in a tripling of the total computational costs if collisions are done at each time step, but typically will be a small fraction of the total simulation cost. It is also shown that binary-collision-based algorithms are unsuitable for partially linearized simulations.

A gyrokinetic theory is developed under a set of orderings applicable to the edge region of tokamaks and other magnetic confinement devices, as well as to internal transport barriers. The result is a practical set equations that is valid for large perturbation amplitudes [q{\delta}{\psi}/T = O(1), where {\delta}{\psi} = {\delta}{\phi} - v_par {\delta}A_par/c], which is straightforward to implement numerically, and which has straightforward expressions for its conservation properties. Here, q is the particle charge, {\delta}{\phi} and {\delta}A_par are the perturbed electrostatic and parallel magnetic potentials, v_par is the parallel velocity, c is the speed of light, and T is the temperature. The derivation is based on the quantity {\epsilon}:=({\rho}/{\lambda})q{\delta}{\psi}/T << 1 as the small expansion parameter, where {\rho} is the gyroradius and {\lambda} is the perpendicular wavelength. Physically, this ordering requires that the E\times B velocity and the component of the parallel velocity perpendicular to the equilibrium magnetic field are small compared to the thermal velocity. For nonlinear fluctuations saturated at "mixing-length" levels (i.e., at a level such that driving gradients in profile quantities are locally flattened), {\epsilon} is of order {\rho}/L, where L is the equilibrium profile scale length, for all scales {\lambda} ranging from {\rho} to L. This is true even though q{\delta}{\psi}/T = O(1) for {\lambda} ~ L. Significant additional simplifications result from ordering L/R =O({\epsilon}), where R is the spatial scale of variation of the magnetic field. We argue that these orderings are well satisfied in strong-gradient regions, such as edge and screapeoff layer regions and internal transport barriers in tokamaks, and anticipate that our equations will be useful as a basis for simulation models for these regions., Comment: Accepted for publication in the Physics of Plasmas, 12/30/2011

If the collisional time scale for Coulomb collisions is comparable to the characteristic time scales for a plasma, then simulation of Coulomb collisions may be important for computation of kinetic plasma dynamics. This can be a computational bottleneck because of the large number of simulated particles and collisions (or phase-space resolution requirements in continuum algorithms), as well as the wide range of collision rates over the velocity distribution function. This paper considers Monte Carlo simulation of Coulomb collisions using the binary collision models of Takizuka and Abe and of Nanbu. It presents a hybrid method for accelerating the computation of Coulomb collisions. The hybrid method represents the velocity distribution function as a combination of a thermal component (a Maxwellian distribution) and a kinetic component (a set of discrete particles). Collisions between particles from the thermal component preserve the Maxwellian; collisions between particles from the kinetic component are per...

Particle simulations are used to solve the fully nonlinear, collisional kinetic equation describing the interaction of a high-temperature, high-density, deuterium-tritium plasma with absorbing boundaries, a plasma source, and the influence of kinetic effects on fusion reaction rates. Both hydrodynamic and kinetic effects influence the end losses, and the simulations show departures of the ion velocity distributions from Maxwellian due to the reduction of the population of the highest energy ions (Knudsen-layer effects). The particle simulations show that the interplay between sources, plasma dynamics, and end losses results in temperature anisotropy, plasma cooling, and concomitant reductions in the fusion reaction rates. However, for the model problems and parameters considered, particle simulations show that Knudsen-layer modifications do not significantly affect the velocity distribution function for velocities most important in determining the fusion reaction rates, i.e., the thermal fusion reaction rates using the local densities and bulk temperatures give good estimates of the kinetic fusion reaction rates.

An efficient and versatile non-Fourier method for the computation of Landau-fluid (LF) closure operators [Hammett and Perkins, Phys. Rev. Lett. 64, 3019 (1990)] is presented, based on an approximation by a sum of modified-Helmholtz-equation solves (SMHS) in configuration space. This method can yield fast-Fourier-like scaling of the computational time requirements and also provides a very compact data representation of these operators, even for plasmas with large spatial nonuniformity. As a result, the method can give significant savings compared with direct application of “delocalization kernels” [e.g., Schurtz et al., Phys. Plasmas 7, 4238 (2000)], both in terms of computational cost and memory requirements. The method is of interest for the implementation of Landau-fluid models in situations where the spatial nonuniformity, particular geometry, or boundary conditions render a Fourier implementation difficult or impossible. Systematic procedures have been developed to optimize the resulting operators for...

The steadily increasing power of supercomputing systems is enabling very high resolution simulations of compressible, turbulent flows in the high Reynolds number limit, which is of interest in astrophysics as well as in several other fluid dynamical applications. This paper discusses two such simulations, using grids of up to 8 billion cells. In each type of flow, convergence in a statistical sense is observed as the mesh is refined. The behavior of the convergent sequences indicates how a subgrid-scale model of turbulence could improve the treatment of these flows by high-resolution Euler schemes like PPM. The best resolved case, a simulation of a Richtmyer-Meshkov mixing layer in a shock tube experiment, also points the way toward such a subgrid-scale model. Analysis of the results of that simulation indicates a proportionality relationship between the energy transfer rate from large to small motions and the determinant of the deviatoric symmetric strain as well as the divergence of the velocity for the large-scale field.

The primary research objective of the Numerical Tokamak Turbulence Project (NTTP) is to develop a predictive ability in modeling turbulent transport due to drift-type instabilities in the core of tokamak fusion experiments, through the use of three-dimensional kinetic and fluid simulations and the derivation of reduced models.

The primary research objective of the Numerical Tokamak Turbulence Project (NTTP) is to develop a predictive ability in modeling turbulent transport due to drift-type instabilities in the core of tokamak fusion experiments, through the use of three-dimensional kinetic and fluid simulations and the derivation of reduced models.

Understanding turbulence and mix in compressible flows is of fundamental importance to real-world applications such as chemical combustion and supernova evolution. The ability to run in three dimensions and at very high resolution is required for the simulation to accurately represent the interaction of the various length scales, and consequently, the reactivity of the intermixin species. Toward this end, we have carried out a very high resolution (over 8 billion zones) 3-D simulation of the Richtmyer-Meshkov instability and turbulent mixing on the IBM Sustained Stewardship TeraOp (SST) system, developed under the auspices of the Department of Energy (DOE) Accelerated Strategic Computing Initiative (ASCI) and located at Lawrence Livermore National Laboratory. We have also undertaken an even higher resolution proof-of-principle calculation (over 24 billion zones) on 5832 processors of the IBM system, which executed for over an hour at a sustained rate of 1.05 Tflop/s, as well as a short calculation with a modified algorithm that achieved a sustained rate of 1.18Tflop/s. The full production scientific simulation, using a further modified algorithm, ran for 27,000 timesteps in slightly over a week of wall time using 3840 processors of the IBM system, clockin a sustained throughput of roughly 0.6 teraflop per second (32-bit arithmetic). Nearly 300,000 graphics files comprising over three terabytes of data were produced and post-processed. The capability of running in 3-D at high resolution enabled us to get a more accurate and detailed picture of the fluid-flow structure - in particular, to simulate the development of fine scale structures from the interactions of long-and short-wavelength phenomena, to elucidate differences between two-dimensional and three-dimensional turbulence, to explore a conjecture regarding the transition from unstable flow to fully developed turbulence with increasing Reynolds number, and to ascertain convergence of the computed solution with respect to mesh resolution.

One of the programs of the Magnetic fusion Energy (MFE) Theory and computations Program is studying the anomalous transport of thermal energy across the field lines in the core of a tokamak. We use the method of gyrokinetic particle-in-cell simulation in this study. For this LDRD project we employed massively parallel processing, new algorithms, and new algorithms, and new formal techniques to improve this research. Specifically, we sought to take steps toward: researching experimentally-relevant parameters in our simulations, learning parallel computing to have as a resource for our group, and achieving a 100 {times} speedup over our starting-point Cray2 simulation code`s performance.

An error in (17) of B. I. Cohen, A. M. Dimits, A. Friedman, and R. E. Caflisch, IEEE Trans. Plasma Sci. 38, 2394 (2010), is corrected. The consequences of this error, which are minor, are discussed.

A set of coordinates for simulations in toroidal magnetic geometry, called quasiballooning coordinates, is proposed and implemented. Quasiballooning coordinates are straight-field-line coordinates in which one of the coordinate directions is as close as possible to that of the magnetic field consistent with the near-parallel grid lines meshing exactly (without interpolation) as they cross the boundaries of the simulation region. This allows the true periodicity conditions in the toroidal-poloidal plane to be satisfied in a straightforward and seamless way, even for sheared magnetic fields. Quasiballooning coordinates are useful in the simulation of instabilities and turbulence of interest in fusion plasmas since the number of grid cells needed to represent structures that are elongated along the magnetic field with a given resoluton is greatly reduced compared with toroidal-poloidal or other nontwisting coordinates. For explicit codes, they allow shorter time steps, and it is anticipated that for particle codes, their use will naturally minimize the numerical noise. The key details necessary for the implementation of quasiballooning coordinates, both in finite-difference and pseudospectral fluid codes are presented, and a fluid code has been written. The advantages of quasiballooning coordinates are demonstrated by applying this code to turbulence driven by the ${\mathrm{\ensuremath{\upsilon}}}_{\mathrm{\ensuremath{\Vert}}}^{\ensuremath{'}}$ instability.

A gyrokinetic theory has been developed in an extended ordering in which the small parameter is the ratio of the E×B shearing rate to the gyrofrequency. This allows for long wavelength E×B flows of order of the thermal velocity, instead of the more restrictive standard orderings which either require that the electrostatic potential or the E×B flow velocity be small compared with the thermal levels. This theory generalizes prior work to allow for time dependence in the large long-wavelength component of the electric field and a continuum of scales in the field components rather than just two distinct components. In the new theory, a significant part of the polarization drift now resides in the equations of motion. However, there is still an identifiable polarization density that can be used to solve for the electrostatic potential from a quasineutrality or vorticity equation. The present derivation is carried out for the case of electrostatic perturbations and a slab equilibrium magnetic field, as this is ...

We present edge gyrokinetic simulations of tokamak plasmas using the fully non-linear (full-f) continuum code TEMPEST. A non-linear Boltzmann model is used for the electrons. The electric field is obtained by solving the 2D gyrokinetic Poisson equation. We demonstrate the following. (1) High harmonic resonances (n > 2) significantly enhance geodesic-acoustic mode (GAM) damping at high q (tokamak safety factor), and are necessary to explain the damping observed in our TEMPEST q-scans and consistent with the experimental measurements of the scaling of the GAM amplitude with edge q 95 in the absence of obvious evidence that there is a strong q-dependence of the turbulent drive and damping of the GAM. (2) The kinetic GAM exists in the edge for steep density and temperature gradients in the form of outgoing waves, its radial scale is set by the ion temperature profile, and ion temperature inhomogeneity is necessary for GAM radial propagation. (3) The development of the neoclassical electric field evolves through different phases of relaxation, including GAMs, their radial propagation and their long-time collisional decay. (4) Natural consequences of orbits in the pedestal and scrape-off layer region in divertor geometry are substantial non-Maxwellian ion distributions and parallel flow characteristics qualitatively like those observed in experiments.

Recent gyrokinetic simulations of electron temperature gradient (ETG) turbulence with flux-tube continuum codes vs. the global particle-in-cell (PIC) code GTC yielded different results despite similar plasma parameters. Differences between the simulations results were attributed to insufficient phase-space resolution and novel physics associated with toroidicity and/or global simulations. We have reproduced the results of the global PIC code using the flux-tube PIC code PG3EQ, thereby eliminating global effects as the cause of the discrepancy. We show that the late-time decay of ETG turbulence and the steady-state heat transport observed in these PIC simulations results from discrete particle noise. Discrete particle noise is a numerical artifact, so both these PG3EQ simulations and the previous GTC simulations have nothing to say about steady-state ETG turbulence and the associated anomalous heat transport. In the course of this work we develop three diagnostics which can help to determine if a particular PIC simulation has become dominated by discrete particle noise.