Gyrokinetic App: Electromagnetic gyrokinetic model for magnetized plasmas¶

The Gyrokinetic App solves the (electromagnetic) gyrokinetic system on a Cartesian grid.

Overall structure of app¶

To set up a gyrokinetic simulation, we first need to load the Gyrokinetic App package. This should be done at the top of the input file, via

local Plasma = (require "App.PlasmaOnCartGrid").Gyrokinetic()

This creates a table Plasma that loads the gyrokinetic species, fields, etc. packages.

The general structure of the input file is then

-------------------------------------------------------------------------------
-- App dependencies.
local Plasma = (require "App.PlasmaOnCartGrid").Gyrokinetic()
...

-------------------------------------------------------------------------------
-- Preamble.
...

-------------------------------------------------------------------------------
-- App initialization.
plasmaApp = Plasma.App {
  -----------------------------------------------------------------------------
  -- Common
  ...

  -----------------------------------------------------------------------------
  -- Species
  electron = Plasma.Species {
    -- GkSpecies parameters
    ...
  },

  -- other species, e.g. ions

  -----------------------------------------------------------------------------
  -- Fields
  field = Plasma.Field {
    -- GkField parameters
    ...
  },

  -----------------------------------------------------------------------------
  -- ExternalFields
  extField = Plasma.Geometry {
    -- GkGeometry parameters
    ...
  },
}
-------------------------------------------------------------------------------
-- App run.
plasmaApp:run()

Kinetic simulations may take additional parameters in the input file Common, such as

Species parameters¶

The Gyrokinetic App works with an arbitrary number of species.

Each species should be declared as

-----------------------------------------------------------------------------
-- Species
species_name = Plasma.Species {
  -- GkSpecies parameters
  ...
},

The species name (species_name here) is arbitrary, but will be used for naming in diagnostic files, so names like ion or electron are common.

Here we describe all possible parameters used to specify a gyrokinetic species. Parameters that have default values can be omitted. Units are arbitrary, but often SI units are used. In the following, VDIM refers to the velocity space dimension, and CDIM refers to the configuration space dimension. The gyrokinetic app works for 1X1V (CDIM=1, VDIM=1), 1X2V, 2X2V, and 3X2V. The velocity coordinates are \((v_\parallel, \mu)\). See [Shi2017] for details.

Diagnostics¶

There are species-specific diagnostics available, which mainly consist of moments of the distribution function and integrals (over configuration-space) of these moments. There are also additional species diagnostics which serve as metrics of positivity and collisions-related errors.

Currently there are four types of diagnostic moments, defined below. Note that in these definitions \(\mathrm{d}\mathbf{w}=\mathrm{d}v_\parallel\) or \(\mathrm{d}\mathbf{w}=(2\pi B_0/m)\mathrm{d}v_\parallel\mathrm{d}\mu\) depending on whether it is a 1V or a 2V simulation. We also use the notation \(d_v\) to signify the number of physical velocity-space dimensions included, i.e. \(d_v=1\) for 1V and \(d_v=3\) for 2V. Also, \(v^2=v_\parallel^2\) for 1V and \(v^2=v_\parallel^2+2\mu B_0/m\) for 2V.

• Velocity moments of the distribution function, written as functions of configuration-space position on each diagnostic frame. The options are

  • M0: number density, \(n = M_0 = \int\mathrm{d}\mathbf{w}~f\).

  • M1: parallel momentum density, \(M_1=\int\mathrm{d}\mathbf{w}~v_\parallel f\).

  • M2: energy density, \(M_2 = \int\mathrm{d}\mathbf{w}~v^2 f\).

  • Upar: parallel flow velocity, \(u_\parallel= M_1/n\).

  • Temp: temperature, \(T = (m/d_v)(M_2 - M_1 u_\parallel)/n\)

  • Beta: plasma beta, \(\beta = 2\mu_0 nT/B^2\)

  • Energy: particle energy density (kinetic + potential), \(\mathcal{E}_H = \int\mathrm{d}\mathbf{w}~H f\), where \(H = mv^2/2 + q\phi\) is the Hamiltonian.

• Velocity moments integrated over configuration-space, written as time-series. The options are

  • intM0: particle number, \(N = \int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~f\)

  • intM1: parallel momentum, \(U = \int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~v_\parallel f\)

  • intM2: \(\int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~v^2 f\)

  • intKE: kinetic energy, \(\mathcal{E}_K = ({m}/{2})\int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~v^2 f\)

  • intEnergy: total (kinetic + potential) energy, \(E_H = \int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~H f\), where \(H = mv^2/2 + q\phi\) is the Hamiltonian.

Boundary flux diagnostics¶

One can request diagnostics of the fluxes through non-periodic boundaries by providing a diagnostics = {} table to the boundary condition Apps. For example, sheath boundary conditions along z would do this via

bcz = {Plasma.SheathBC{diagnostics={"M0"}}, Plasma.SheathBC{diagnostics={"M0"}}}

in order to request the particle flux through the sheaths. Another example is provided in the gyrokinetic Quickstart page.

The boundary fluxes are computed via integrals of the time rates of change computed in the ghost cells. If we consider a simple phase-space advection equation in 2X2V without any forces

the weak form used by the algorithm is obtained by multiplying this equation by a test function \(\psi\) and integrating over phase space in a single cell. After an integration by parts one obtains

where the hat means that a numerical flux is constructed, and \(\mathrm{d}\mathbf{z}=\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{v}\). In ghost cells only the surface terms corresponding to fluxes through the physical domain boundaries are computed. This means tha in the ghost cell at the upper boundary along \(x\), for example

This is phase-space flux through the upper \(x\) boundary during a stage of the PDE solver. For Runge-Kutta steppers one must form a linear combination of these fluxes from every stage in the same manner as the time rates of change are combined for forward time stepping. For the sake of simplicity here we just assume a single forward Euler step, and define phase-space flux during a single time step through the upper \(x\) boundary as

where the volume factor \(V\) arises from the phase-space integral on the left side of equation (1). Note that these integrals are over a single cell, and that the quantity \(\Gamma_{\mathbf{z},x_+}\) is phase-space field, \(\Gamma_{\mathbf{z},x_+}=\Gamma_{\mathbf{z},x_+}(\mathbf{x},\mathbf{v})\).

With this boundary flux in mind, if one requests the particle density of the boundary flux through diagnosticBoundaryFluxMoments={GkM0} the diagnostic would be computed as

This yields the rate of number density crossing the upper \(x\) boundary (per cell-length in the \(x\) direction of the ghost cell). In order to compute the number of particles per unit time crossing the upper \(x\) boundary (diagnosticIntegratedBoundaryFluxMoments={intM0}) we simply integrate the above quantity over \(y\) (and multiply it by the \(x\)-cell length of the ghost cell)

The final detail is that the files created by these diagnostics contain the fluxes through the boundary accumulated since the last snapshot (frame), not since the beginning of the simulation.

References¶