Search

Currently, the dominating constraint in many high performance computing applications is data capacity and bandwidth, in both inter-node communications and even more-so in on-node data motion. A new approach to address this limitation is to make use of data compression in the form of a compressed data array. Storing data in a compressed data array and converting to standard IEEE-754 types as needed during a computation can reduce the pressure on bandwidth and storage. However, repeated conversions (lossy compression and decompression) introduce additional approximation errors, which need to be shown to not significantly affect the simulation results. We extend recent work [J. Diffenderfer, et al., Error Analysis of ZFP Compression for Floating-Point Data, SIAM Journal on Scientific Computing, 2019] that analyzed the error of a single use of compression and decompression of the ZFP compressed data array representation [P. Lindstrom, Fixed-rate compressed floating-point arrays, IEEE Transactions on Visualization and Computer Graphics, 2014] to the case of time-stepping and iterative schemes, where an advancement operator is repeatedly applied in addition to the conversions. We show that the accumulated error for iterative methods involving fixed-point and time evolving iterations is bounded under standard constraints. An upper bound is established on the number of additional iterations required for the convergence of stationary fixed-point iterations. An additional analysis of traditional forward and backward error of stationary iterative methods using ZFP compressed arrays is also presented. The results of several 1D, 2D, and 3D test problems are provided to demonstrate the correctness of the theoretical bounds.

As scientific applications strive towards increasingly realistic modeling of complex phenomena, they are integrating multiple models and simulations into complex, coupled scientific workflows. As a result, ensuring that existing codes can be combined and recombined correctly and flexibly as part of these workflows is essential. In this paper, we propose Benesh, a programming system for creating in-situ scientific workflows. Benesh provides a domain-specific abstraction that enables a programmer to instrument an existing simulation code to be used as a building block in defining complex workflows. Using Benesh, developers define a workflow-level shared specification of data objects over common or partitioned data domains. This permits dependency-based execution to be specified at the workflow level, distinct from the independent operation of the component simulations. We additionally describe features of a scalable runtime that builds on a distributed data services layer to implement the Benesh programming system.

It has been conjectured that higher-order discretizations for partial differential equations will have advantages over the lower-order counterparts commonly used today. The reasoning is that the increase in arithmetic operations will be more than offset by the reduction in data transfers and the increase in concurrent floating-point units. To evaluate this conjecture, the arithmetic intensity of a class of high-order finite-volume discretizations for hyperbolic systems of conservation laws is theoretically analyzed for spatial discretizations from orders three through eight in arbitrary dimensions. Three cache models are considered: the limiting cases of no cache and infinite cache as well as a finite-sized cache model. Models are validated experimentally by measuring floating-point operations and data transfers on an IBM Blue Gene/Q node. Theory and experiments demonstrate that high-order finite-volume methods will be able to provide increases in arithmetic intensity that will be necessary to make better utilization of on-node floating-point capability.

HPC applications use floating point arithmetic operations extensively to solve computational problems. Mixed-precision computing seeks to use the lowest precision data type that is sufficient to achieve a desired accuracy, improving performance and reducing power consumption. Manually optimizing a program to use mixed precision is challenging as it not only requires extensive knowledge about the numerical behavior of the algorithm but also estimates of the rounding errors. In this work, we present ADAPT, a scalable approach for mixed-precision analysis on HPC workloads using algorithmic differentiation to provide accurate estimates about the final output error. ADAPT provides a floating-point precision sensitivity profile while incurring an overhead of only a constant multiple of the original computation irrespective of the number of variables analyzed. The sensitivity profile can be used to make algorithmic choices and to develop mixed-precision configurations of a program. We evaluate ADAPT on six benchmarks and a proxy application (LULESH) and show that we are able to achieve a speedup of 1.2× on the proxy application.

We present a scheme that spatially couples two gyrokinetic codes using first-principles. Coupled equations are presented and a necessary and sufficient condition for ensuring accuracy is derived. This new scheme couples both the field and the particle distribution function. The coupling of the distribution function is only performed once every few time-steps, using a five-dimensional (5D) grid to communicate the distribution function between the two codes. This 5D grid interface enables the coupling of different types of codes and models, such as particle and continuum codes, or delta-f and total-f models. Transferring information from the 5D grid to the marker particle weights is achieved using a new resampling technique. Demonstration of the coupling scheme is shown using two XGC gyrokinetic simulations for both the core and edge. We also apply the coupling scheme to two continuum simulations for a one-dimensional advection–diffusion problem.

Compression of floating-point data will play an important role in high-performance computing as data bandwidth and storage become dominant costs. Lossy compression of floating-point data is powerful, but theoretical results are needed to bound its errors when used to store look-up tables, simulation results, or even the solution state during the computation. \black{In this paper, we analyze the round-off error introduced by ZFP, a %state-of-the-art lossy compression algorithm.} The stopping criteria for ZFP depends on the compression mode specified by the user; either fixed rate, fixed accuracy, or fixed precision [P. Lindstrom, Fixed-rate compressed floating-point arrays, IEEE Transactions on Visualization and Computer Graphics, 2014]. While most of our discussion is focused on the fixed precision mode of ZFP, we establish a bound on the error introduced by all three compression modes. In order to tightly capture the error, we first introduce a vector space that allows us to work with binary representations of components. Under this vector space, we define operators that implement each step of the ZFP compression and decompression to establish a bound on the error caused by ZFP. To conclude, numerical tests are provided to demonstrate the accuracy of the established bounds.

We propose a modular framework for representing the real numbers that generalizes ieee, posits, and related floating-point number systems, and which has its roots in universal codes for the positive integers such as the Elias codes. This framework unifies several known but seemingly unrelated representations within a single schema while also introducing new representations. We particularly focus on variable-length encoding of the binary exponent and on the manner in which fraction bits are mapped to values. Our framework builds upon and shares many of the attractive properties of posits but allows for independent experimentation with exponent codes, fraction mappings, reciprocal closure, rounding modes, handling of under- and overflow, and underlying precision.

The vast separation dividing the characteristic times of energy confinement and turbulence in the core of toroidal plasmas makes first-principles prediction on long timescales extremely challenging. Here we report the demonstration of a multiple-timescale method that enables coupling global gyrokinetic simulations with a transport solver to calculate the evolution of the self-consistent temperature profile. This method, which exhibits resiliency to the intrinsic fluctuations arising in turbulence simulations, holds potential for integrating nonlocal gyrokinetic turbulence simulations into predictive, whole-device models., 7 pages, 3 figures

Here, COGENT is a continuum gyrokinetic code for edge plasma simulations being developed by the Edge Simulation Laboratory collaboration. The code is distinguished by application of a fourth-order finite-volume (conservative) discretization, and mapped multiblock grid technology to handle the geometric complexity of the tokamak edge. The distribution function F is discretized in v∥ – μ (parallel velocity – magnetic moment) velocity coordinates, and the code presently solves an axisymmetric full-f gyro-kinetic equation coupled to the long-wavelength limit of the gyro-Poisson equation. COGENT capabilities are extended by implementing the fully nonlinear Fokker-Plank operator to model Coulomb collisions in magnetized edge plasmas. The corresponding Rosenbluth potentials are computed by making use of a finite-difference scheme and multipole-expansion boundary conditions. Details of the numerical algorithms and results of the initial verification studies are discussed. (© 2014 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)

Multiphysics simulations frequently are composed from highly-optimized solvers for physical subprocesses through the use of temporal operator splitting. The subphysics are evolved sequentially, passing information between solver components as needed. It is often useful yet difficult to determine how the simulation error depends on the temporal splitting method, the time step size and the discretization parameters in the solver components. This paper proposes a framework to decompose the total error in a quantity of interest for an advection–diffusion simulation into two primary contributions: that due to the operator splitting, as implied by the structure of the code, and that due to the discretization errors from the component solvers. The method is applied to the advection–diffusion equation with boundaries. The advection and diffusion operators require separate boundary conditions upon splitting; the specification of these impacts the splitting error contribution. Computational examples demonstrate that the proposed method successfully identifies the splitting error contribution, including that induced by a suboptimal imposition of boundary conditions, and the decrease in the splitting error contribution in moving to a higher-order splitting method. The discretization error contribution is decomposed further into contributions attributed to the separate advection and diffusion solvers.

This paper is concerned with the computation of numerical discretization error for uncertainty quantification. An a posteriori error formula is described for a functional measurement of the solution to a scalar advection equation that is estimated by finite volume approximations. An exact error formula and computable error estimate are derived based on an abstractly defined approximation of the adjoint solution. The adjoint problem is divorced from the finite volume method used to approximate the forward solution variables and may be approximated using a low-order finite volume method. The accuracy of the computable error estimate provably satisfies an a priori error bound for sufficiently smooth solutions of the forward and adjoint problems. Computational examples are provided that show support of the theory for smooth solutions. The application to problems with discontinuities is also investigated computationally.

This paper summarizes the outcomes from the Turbulent Flow Simulation at the Exascale: Opportunities and Challenges Workshop, which was held 4-5 August 2015, and was sponsored by the U.S. Department of Energy Office of Advanced Scientific Computing Research. The workshop objective was to define and describe the challenges and opportunities that computing at the exascale will bring to turbulent-flow simulations in applied science and technology. The need for accurate simulation of turbulent flows is evident across the U.S. Department of Energy applied-science and engineering portfolios, including combustion, plasma physics, nuclear-reactor physics, wind energy, and atmospheric science. The workshop brought together experts in turbulent-flow simulation, computational mathematics, and high-performance computing. Building upon previous ASCR workshops on exascale computing, participants defined a research agenda and path forward that will enable scientists and engineers to continually leverage, engage, and direct advances in computational systems on the path to exascale computing.

The estimation of discretization error in numerical simulations is a key component in the development of uncertainty quantification. In particular, there exists a need for reliable, robust estimators for finite volume and finite difference discretizations of hyperbolic partial differential equations. The approach espoused here, often called the error transport approach in the literature, is to solve an auxiliary error equation concurrently with the primal governing equation to obtain a point-wise (cell-wise) estimate of the discretization error. Nonlinear, time-dependent problems are considered. In contrast to previous work, fully nonlinear error equations are advanced, and potential benefits are identified. A systematic approach to approximate the local residual for both method-of-lines and space–time discretizations is developed. Behavior of the error estimates on problems that include weak solutions demonstrates the positive properties of nonlinear error transport.

We present an approach for constructing finite-volume methods for flux-divergence forms to any order of accuracy defined as the image of a smooth mapping from a rectangular discretization of an abstract coordinate space. Our approach is based on two ideas. The first is that of using higher-order quadrature rules to compute the flux averages over faces that generalize a method developed for Cartesian grids to the case of mapped grids. The second is a method for computing the averages of the metric terms on faces such that freestream preservation is automatically satisfied. We derive detailed formulas for the cases of fourth-order accurate discretizations of linear elliptic and hyperbolic partial differential equations. For the latter case, we combine the method so derived with Runge-Kutta time discretization and demonstrate how to incorporate a high-order accurate limiter with the goal of obtaining a method that is robust in the presence of discontinuities and underresolved gradients. For both elliptic and hyperbolic problems, we demonstrate that the resulting methods are fourth-order accurate for smooth solutions.

With the next generation of high-power laser facilities for inertial fusion coming online1,2, ensuring laser beam propagation through centimetre-scale plasmas is a key physics issue for reaching ignition. Existing experimental results3,4,5 including the most recent one6 are limited to small laser spots, low-interaction laser beam energies and small plasma volumes of 1–2 mm. Here, we demonstrate the propagation of an intense, high-energy, ignition-size laser beam through fusion-size plasmas on the National Ignition Facility (NIF) and find the experimental measurements to agree with full-scale modelling. Previous attempts to apply computer modelling as a predictive capability have been limited by the inherently multiscale description of the full laser–plasma interaction processes7,8,9,10,11. The findings of this study validate supercomputer modelling as an essential tool for the design of future ignition experiments.

The following results are presented from the development and application of TEMPEST, a fully nonlinear (full-f) five-dimensional (3d2v) gyrokinetic continuum edge-plasma code. (1) As a test of the interaction of collisions and parallel streaming, TEMPEST is compared with published analytic and numerical results for endloss of particles confined by combined electrostatic and magnetic wells. Good agreement is found over a wide range of collisionality, confining potential and mirror ratio, and the required velocity space resolution is modest. (2) In a large-aspect-ratio circular geometry, excellent agreement is found for a neoclassical equilibrium with parallel ion flow in the banana regime with zero temperature gradient and radial electric field. (3) The four-dimensional (2d2v) version of the code produces the first self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal flow (Rosenbluth–Hinton residual) with Boltzmann electrons using a full-f code. The electric field is also found to agree with the standard neoclassical expression for steep density and ion temperature gradients in the plateau regime. In divertor geometry, it is found that the endloss of particles and energy induces parallel flow stronger than the core neoclassical predictions in the SOL.

Structured-grid PDE solver frameworks parallelize over boxes, which are rectangular domains of cells or faces in a structured grid. In the Chombo framework, the box sizes are typically 163 or 323, but larger box sizes such as 1283 would result in less surface area and therefore less storage, copying, and/or ghost cells communication overhead. Unfortunately, current on node parallelization schemes perform poorly for these larger box sizes. In this paper, we investigate 30 different inter-loop optimization strategies and demonstrate the parallel scaling advantages of some of these variants on NUMA multicore nodes. Shifted, fused, and communication-avoiding variants for 1283 boxes result in close to ideal parallel scaling and come close to matching the performance of 163 boxes on three different multicore systems for a benchmark that is a proxy for program idioms found in Computational Fluid Dynamic (CFD) codes.

A coupled mode system is derived to investigate a three-wave parametric instability leading to energy transfer between co-propagating laser beams crossing in a plasma flow. The model includes beams of finite width refracting in a prescribed transverse plasma flow with spatial and temporal gradients in velocity and density. The resulting paraxial light equations are discretized spatially with a Crank-Nicholson-type scheme, and these algebraic constraints are nonlinearly coupled with ordinary differential equations in time that describe the ion acoustic response. The entire nonlinear differential-algebraic system is solved using an adaptive, backward-differencing method coupled with Newton's method. A numerical study is conducted in two dimensions that compares the intensity gain of the fully time-dependent coupled mode system with the gain computed under the further assumption of a strongly damped ion acoustic response. The results demonstrate a time-dependent gain suppression when the beam diameter is commensurate with the velocity gradient scale length. The gain suppression is shown to depend on time-dependent beam refraction and is interpreted as a time-dependent frequency shift.

The discontinuous Riemann initial value problem is considered for dissipative, constant coefficient hyperbolic systems with relaxation source terms. A short-time asymptotic expansion, a specialization of the results of Le Floch and Raviart [11], is constructed for the general linear system. Exact rates for the decay of the initial discontinuities are found to be in agreement with recent results from other approaches. A multiple scales analysis is used to identify the long-time asymptotic behavior. Many of the results can be exemplified in a simple model problem for which the exact solution can be found.

An analysis of the effects of ion trapping on ion acoustic waves excited by the stimulated Brillouin scattering of crossing intense laser beams is presented. Ion trapping alters the dispersion of ion acoustic waves by nonlinearly shifting the normal mode frequency and by reducing the ion Landau damping. This in turn can influence the energy transfer between two crossing laser beams in the presence of plasma flows such that stimulated Brillouin scattering (SBS) occurs. The same ion trapping physics can influence the saturation of SBS in other circumstances. A one-dimensional analytical model is presented along with reasonably successful comparisons of the theory to results from particle simulations and laboratory experiments. An analysis of the vulnerability of the National Ignition Facility Inertial Confinement Fusion point design [S. W. Haan et al., Fusion Sci. Technol. 41, 164 (2002)] is also presented.

We investigate the use of adaptive mesh refinement in the simulation of laser plasma filamentation. A numerical algorithm is constructed to solve model equations consisting of a fluid approximation of a quasineutral plasma combined with a paraxial light propagation model. The algorithm involves high-resolution plasma and light model discretizations on a block-structured, locally refined grid hierarchy, which is dynamically modified during the time integration to follow evolving fine-scale solution features. Comparisons of the efficiency of this approach to that of uniform grid calculations are presented.

Exascale computing systems are essential for the scientific fields that will transform the 21st century global economy, including energy, biotechnology, nanotechnology, and materials science. Progress in these fields is predicated on the ability to perform advanced scientific and engineering simulations, and analyze the deluge of data. On July 29, 2013, ASCAC was charged by Patricia Dehmer, the Acting Director of the Office of Science, to assemble a subcommittee to provide advice on exascale computing. This subcommittee was directed to return a list of no more than ten technical approaches (hardware and software) that will enable the development of a system that achieves the Department's goals for exascale computing. Numerous reports over the past few years have documented the technical challenges and the non¬-viability of simply scaling existing computer designs to reach exascale. The technical challenges revolve around energy consumption, memory performance, resilience, extreme concurrency, and big data. Drawing from these reports and more recent experience, this ASCAC subcommittee has identified the top ten computing technology advancements that are critical to making a capable, economically viable, exascale system.

The coupling of radiation transport and hydrodynamics is discussed for the Eulerian frame. The discussion is aimed at developing a suitable set of equations for nonrelativistic radiation hydrodynamics (RHD) that can be numerically integrated using high-resolution methods for conservation laws. We outline how numerical methods based on a wave decomposition may be developed, along with the importance of conservation, particularly in the equilibrium regime. The properties of the RHD equations are examined through asymptotic and dispersion analyses. The conditions required to obtain the classical equilibrium limit are rigorously studied. The results show that a simple coupling term developed recently by Morel, which retains a minimum of relativistic corrections, may be sufficient for nonrelativistic flows. We also give two constraints on the relativistic corrections that result in retaining terms on the order of the truncation. In addition, the dispersion results for the P1 approximation are studied in detail and are compared with both the exact-transport results and a full relativistic treatment. We also examine some nonintuitive behavior in the dispersion results.

A three-dimensional grain model, in which the grains are represented by regular truncated octahedra, has been developed to study probabilistic time-dependent intergranular failure in polycrystalline arrays. In this model, grain boundary facets are assumed to fail randomly in time, as a function of the facet normal stress. A simple approximate method for calculating the load shed by failed facets and a reasonable choice of failure criterion complete the model. This leads to a conceptually simple, but computationally complex, model capable of handling assemblages consisting of relatively large numbers (> 5000) of grains. The predicted scatter in the times-to-failure and the variation in number of failed facets with time are in quite reasonable agreement with available experimental data.

The first 4D (axisymmetric) high-order continuum gyrokinetic transport simulations that span the magnetic separatrix of a tokamak are presented. The modeling is performed with the COGENT code, which is distinguished by fourth-order finite-volume discretization combined with mapped multiblock grid technology to handle the strong anisotropy of plasma transport and the complex X-point divertor geometry with high accuracy. The calculations take into account the effects of fully nonlinear Fokker-Plank collisions, electrostatic potential variations, and anomalous radial transport. Topics discussed include: (a) ion orbit loss and the associated toroidal rotation and (b) edge plasma relaxation in the presence of anomalous radial transport.

Direct discretization of continuum kinetic equations, like the Vlasov equation, are under-utilized because the distribution function generally exists in a high-dimensional (>3D) space and computational cost increases geometrically with dimension. We propose to use high-order finite-volume techniques with block-structured adaptive mesh refinement (AMR) to reduce the computational cost. The primary complication comes from a solution state comprised of variables of different dimensions. We develop the algorithms required to extend standard single-dimension block structured AMR to the multi-dimension case. Specifically, algorithms for reduction and injection operations that transfer data between mesh hierarchies of different dimensions are explained in detail. In addition, modifications to the basic AMR algorithm that enable the use of high-order spatial and temporal discretizations are discussed. Preliminary results for a standard 1D+1V Vlasov-Poisson test problem are presented. Results indicate that there is potential for significant savings for some classes of Vlasov problems.

An a posteriori error formula is described when a statistical measurement of the solution to a hyperbolic conservation law in 1D is estimated by finite volume approximations. This is accomplished using adjoint error estimation. In contrast to previously studied methods, the adjoint problem is divorced from the finite volume method used to approximate the forward solution variables. An exact error formula and computable error estimate are derived based on an abstractly defined approximation of the adjoint solution. This framework allows the error to be computed to an arbitrary accuracy given a sufficiently well resolved approximation of the adjoint solution. The accuracy of the computable error estimate provably satisfies an a priori error bound for sufficiently smooth solutions of the forward and adjoint problems. The theory does not currently account for discontinuities. Computational examples are provided that show support of the theory for smooth solutions. The application to problems with discontinuities is also investigated computationally.

Experiments demonstrate energy and power transfer between copropagating, same frequency, beams crossing at a small angle in a plasma with a Mach 1 flow. The process is interpreted as amplification of the low intensity probe beam by the stimulated scatter of the high intensity pump beam. The observed probe amplification increases slowly with pump intensity and decreases with probe intensity, indicative of saturation limiting the energy and power transfer due to ion-wave nonlinearities and localized pump depletion. The results are consistent with numerical modeling including ion-wave nonlinearities.

We explore the issues involved in creating a Godonov-type method for flows with stiff source terms. We briefly survey the theory of such flows with particular emphasis on different asymptotic regimes and the solution of the Riemann problem. We then discuss practical issues, such as the choice of working variables, with particular reference to a set of eleven equations describing a non-equilibrium diatomic gas. We propose a uniformly-valid discretisation technique incorporating an approximate Riemann solver, but the details of the solver remain to be determined.

Geodesic acoustic modes (GAMs) are an important phenomenon in a tokamak edge plasma. They regulate turbulence in a low confinement (L-mode) regime and can play an important role in the low to high (L–H) mode transition. It is therefore of considerable importance to develop a detailed theoretical understanding of their dynamics and relaxation processes. The present work reports on the numerical modelling of collisionless GAM relaxation, including the effects of a strong radial electric field characteristic of a tokamak pedestal in a high confinement (H-mode) regime. The simulations demonstrate that the presence of a strong radial electric field enhances the GAM decay rate, and heuristic arguments elucidating this finding are provided. The numerical modelling is performed by making use of the continuum gyrokinetic code COGENT.

The development of the continuum gyrokinetic code COGENT for edge plasma simulations is reported. The present version of the code models a nonlinear axisymmetric 4D (R, v∥, μ) gyrokinetic equation coupled to the long-wavelength limit of the gyro-Poisson equation. Here, R is the particle gyrocenter coordinate in the poloidal plane, and v∥ and μ are the guiding center velocity parallel to the magnetic field and the magnetic moment, respectively. The COGENT code utilizes a fourth-order finite-volume (conservative) discretization combined with arbitrary mapped multiblock grid technology (nearly field-aligned on blocks) to handle the complexity of tokamak divertor geometry with high accuracy. Topics presented are the implementation of increasingly detailed model collision operators, and the results of neoclassical transport simulations including the effects of a strong radial electric field characteristic of a tokamak pedestal under H-mode conditions.

Two-dimensional Vlasov simulations of nonlinear electron plasma waves are presented, in which the interplay of linear and nonlinear kinetic effects is evident. The plasma wave is created with an external traveling wave potential with a transverse envelope of width Δy such that thermal electrons transit the wave in a “sideloss” time, tsl~Δy/ve. Here, ve is the electron thermal velocity. The quasisteady distribution of trapped electrons and its self-consistent plasma wave are studied after the external field is turned off. In cases of particular interest, the bounce frequency, ωbe=keϕ/me, satisfies the trapping condition ωbetsl>2π such that the wave frequency is nonlinearly downshifted by an amount proportional to the number of trapped electrons. Here, k is the wavenumber of the plasma wave and ϕ is its electric potential. For sufficiently short times, the magnitude of the negative frequency shift is a local function of ϕ. Because the trapping frequency shift is negative, the phase of the wave on axis lags ...

We are developing a new class of finite-volume methods on locally-refined and mapped grids, which are at least fourth-order accurate in regions where the solution is smooth. This paper discusses the implementation of such methods for time-dependent problems on both Cartesian and mapped grids with adaptive mesh refinement. We show 2D results with the Berger-Colella shock-ramp problem in Cartesian coordinates, and fourth-order accuracy of the solution of a Gaussian pulse problem in a polytropic gas in mapped coordinates.

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.