Collision models in Gkeyll¶
In Gkeyll we currently have two different collision operators for use in kinetic models: the Bhatnagar–Gross–Krook (BGK) and the Dougherty operators. We referred to the latter as the LBO for the legacy of Lenard-Bernstein. Its implementation in Gkeyll is detailed in [Hakim2020] [Francisquez2020]. See [Francisquez2022] for definitions of the Dougherty LBO coefficients and alternative choices.
Contents
BGK collisions¶
The BGK operator [Gross1956] for the effect of collisions on the distribution of species \(f_s\) is
where the sum is over all the species. The distribution functon \(f_{Msr}\) is the Maxwellian
with the primitive moments \(\mathbf{u}_{sr}\) and \(v_{tsr}^2\) properly defined to preserve some properties (such as conservation), and \(d_v\) is the number of velocity-space dimensions. For self-species collisions \(\mathbf{u}_{sr}=\mathbf{u}_s\) and \(v_{tsr}^2=v_{ts}^2\). For multi-species collisions we follow an approach similar to [Greene1973] and define the cross-species primitive moments as
but contrary to Greene’s definition of \(\alpha_E\), we currently use in Gkeyll the following expression
Little guidance is provided by Greene as to how to choose \(\beta\), although it seems clear that \(-1<\beta\). In Gkeyll the default value is \(\beta=0\), but the user can specify it in the input file (explained below). We have introduced the additional quantity \(\delta_s\) (which Greene indirectly assumed to equal 1) defined as
The BGK operator can be used with both the Vlasov-Maxwell solver and the gyrokinetic solver.
Dougherty collisions¶
The Doughery (LBO) model for collisions [Dougherty1964] in Gkeyll is given by
In this case we compute the cross-primitive moments by a process analogous to Greene’s with the BGK operator, yielding the following formulas for the cross flow velocity and thermal speed:
with \(\alpha_E\) defined in the BGK section above. The LBO used by the gyrokinetic solver is
Collisions in Gkeyll input files¶
Users can specify collisions in input files by adding an additional Lua table within each species one wishes to add collisions to. The collision frequency can be constant, varying in space and time, or it can also have a user-defined profile.
Constant collisionality¶
An example of adding
LBO collisions (for BGK collisions simply replace LBOcollisions
with
BGKCollisions
) to a species named ‘elc’ is
elc = Plasma.Species {
charge = q_e, mass = m_e,
-- Velocity space grid.
...
-- Initial conditions.
...
evolve = true,
-- Collisions.
coll = Plasma.LBOCollisions {
collideWith = { "elc" },
frequencies = { nu_ee },
},
},
If there were another species, say one named ‘ion’, this ‘elc’ species could
be made to collide with ‘ion’ by adding ‘ion’ to the collideWidth
table:
coll = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
frequencies = { nu_ee, nu_ei },
},
The constant collision frequencies nu_ee
and nu_ei
need to be previously
computed/specified in the input file. The user can specify the value of \(\beta\)
in the above formulas for the cross-species primitive moments (\(\mathbf{u}_{sr}\)
and \(v_{tsr}^2\)) by specifying the variable betaGreene
in the collisions
table (if the user does not specify it, betaGreene=0.0
is assumed) like
coll = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
frequencies = { nu_ee, nu_ei },
betaGreene = 0.9
},
In some cases the user may be interested in colliding species ‘elc’ with species ‘ion’, but not collide species ‘ion’ with species ‘elc’. Gkeyll supports this combination, but since the formulas for cross-species primitive moments involve both \(\nu_{ei}\) and \(\nu_{ie}\), the code will default to assuming \(\nu_{ie}=m_e\nu_{ei}/m_i\). Note however that this scenario is not energy conserving: for exact energy conservation, one must include the effect of binary collisions on both species.
It is also possible to specify both LBO and BGK collisions between different binary pairs in a single input file. For example, if there are three species ‘elc’, ‘ion’ and ‘neut’, the ‘elc’ species could be made collide with both ‘ion’ and ‘neut’ as follows:
cColl = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
frequencies = { nu_ee, nu_ei },
},
nColl = Plasma.BGKCollisions {
collideWith = { "neut" },
frequencies = { nu_en },
},
If no collisionality is specified in the input file, it is assumed that the user desires Gkeyll to build a spatially-varying collisionality from scratch using a Spitzer-like formula for \(\nu_{sr}\) (explained below).
Spatially varying collisionality¶
Currently there are three ways to run simulations with a spatially varying collisionality. All of these options lead to a spatially varying, cell-wise constant collisionality. We will be adding support for variation of the collisionality within a cell in the future.
Option A¶
The simplest way to run with spatially varying collisionality is to not specify
the table frequencies
. In this case the code computes \(\nu_{sr}\)
according to
where \(\nu_{\mathrm{frac}}\) is a scaling factor, the Coulomb logarithm is defined as
and the \(\alpha\)-sum is over all the species. For Vlasov-Maxwell simulations we do not add the correction due to gyromotion (\(\omega_{c\alpha}=0\) here). The relative velocity here is computed as \(u^2=3v_{tr}^2+3v_{ts}^2\), the reduced mass is \(m_{sr} = m_sm_r/\left(m_s+m_r\right)\), and \(\omega_{p\alpha}\) is the plasma frequency computed with the density and mass of species \(\alpha\). Simpler formulas for the Coulomb logarithm can be easily generated by developers if necessary.
The formulas above assume all the plasma quantities and universal constants are in SI units. The user can provide a different value for these variables by passing them to the collisions table in the input files, as shown here:
coll = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
epsilon0 = 1.0, -- Vacuum permitivity.
elemCharge = 1.0, -- Elementary charge value.
hBar = 1.0, -- Planck's constant h/2pi.
},
Additionally the user can pass the scaling factor \(\nu_{\mathrm{frac}}\) by
specifying nuFrac
in the collisions table.
Option B¶
Another way to use a spatially varying collisionality is to pass a reference
collisionality normalized to a combination of the density and thermal speed of the
colliding species. This normalized collisionality, is defined as
\(\nu_{srN}=\nu_{sr0}\left(v_{ts0}^2+v_{tr0}^2\right)^{3/2}/n_{r0}\) and one
provides through normNu
in the collisions table as shown below:
elc = Plasma.Species {
...
coll = Plasma.LBOCollisions {
collideWith = { "ion" },
normNu = { nu_ei*((vte^2+vti^2)^(3/2))/n_i0 }
},
},
where nu_ei
, vte
, vti
, n_e0
are computed in the Preamble of the
input file and it is up to the user to ensure that these all have consistent units.
Then, in each time step, the collisions will be applied with the following collisionality
Note that if one is using the normNu
feature for self-species collisions, one must
still use these formulas. In this case one would specify electron-electron collisions like
elc = Plasma.Species {
...
coll = Plasma.LBOCollisions {
collideWith = { "elc" },
normNu = { nu_ee*((2*(vte^2))^(3/2))/n_e0 }
},
},
Option C¶
The user may also wish to specify their own collisionality profile, so for this purpose
one can pass functions into the frequencies
table in the collisions table.
For example, suppose that one would like to run a simulation with a collisionality that decays exponentially in x. In this case we could create a exponentially decaying function in the preamble and pass it as the collision frequency as follows:
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell
local Constants = require "Lib.Constants"
eps0 = Constants.EPSILON0
eV = Constants.ELEMENTARY_CHARGE
me = Constants.ELECTRON_MASS
n0 = 7e19 -- Number density [1/m^3].
Te0 = 100*eV -- Electron temperature [J].
-- Reference electron collision frequency (at x=0).
logLambdaElc = 24.0 - 0.5*math.log(n0/1e6) + math.log(Te0/eV)
nu_ee = logLambdaElc*(eV^4)*n0
/(12*math.sqrt(2)*(math.pi^(3/2))*(eps0^2)*math.sqrt(me)*(Te0^(3/2)))
local function nu_eeProfile(t, xn)
local x = xn[1]
return nu_ee*math.exp(-x)
end
vlasovApp = Plasma.App {
...
elc = Plasma.Species {
...
-- Collisions.
coll = Plasma.LBOCollisions {
collideWith = { "elc" },
frequencies = { nu_eeProfile },
}
}
...
}
-- Run application.
vlasovApp:run()
At present all the frequencies
must either be constant numbers or functions. We do not
yet support having a combination of the two in the same collisions table.
Comments on stability¶
The are known issues with the implementation of the collision operators in Gkeyll. One of them, for example, is that we do not have a positivy preseving algorithm for the LBO. Positivity issues are often accompanied by large flows or negative temperatures and/or densities. For this reason we have taken three precautions:
Calculation of primitive moments \(\mathbf{u}_{sr}\) and \(v_{tsr}^2\) is carried out using cell-average values if the number density is non-positive at one of the corners of that cell.
The collision term is turned off locally if the flow velocity \(\mathbf{u}_{sr}\) is greater than the velocity limits of the domain, or if \(v_{tsr}^2\) is negative.
The collision frequency \(\nu_{sr}\) is locally set to zero if the cell-average values of \(n_r\) or \(v_{tsr}^2\) are negative.
We track the number of cells in which precaution 2 is used, and for stable simulations this is typically small (a few percent or less). Further discussion of why these precautions are necessary appears in [Hakim2020].
Examples¶
We offer two full examples of the use of collisions. One in Vlasov-Maxwell and one in Gyrokinetics.
Example 1: 1x1v collisional relaxation¶
Consider an initial distribution function in 1x1v phase space given by a Maxwellian and a large bump in its tail
Suppose we wish to collisionally relax this initial state, without the influence of collisionless terms. That is, we wish to evolve this distribution function according to equation (1). In this case our input file will use the VlasovMaxwell App (for 1x1v it would be equivalent to use the Gyrokinetic App), and we define the distribution in equation (2) in the Preamble via the function
-- Maxwellian with a Maxwellian bump in the tail.
local function bumpMaxwell(x,vx,n,u,vth,bN,bU,bVth,bL,bS)
local vSq = ((vx-u)/(math.sqrt(2.0)*vth))^2
local vbSq = ((vx-bU)/(math.sqrt(2.0)*bVth))^2
return (n/math.sqrt(2.0*math.pi*vth))*math.exp(-vSq)
+(bN/math.sqrt(2.0*math.pi*bVth))*math.exp(-vbSq)/((vx-bL)^2+bS^2)
end
In this case we chose constants for all densities, flow speed and temperatures. We
also set the charge to 0. Under these conditions the collisionless terms have no effect,
but we can explicitly turn them off with the evolveCollisionless
flag. We will also
request the total integrated bulk flow energy (intM2Flow
) and the total thermal
energy (intM2Thermal
) as diagnostics.
plasmaApp = Plasma.App {
tEnd = 80, -- End time.
nFrame = 80, -- Number of frames to write.
lower = {0.0}, -- Configuration space lower coordinate.
upper = {1.0}, -- Configuration space upper coordinate.
cells = {8}, -- Configuration space cells.
polyOrder = 2, -- Polynomial order.
periodicDirs = {1}, -- Periodic directions.
-- Neutral species with a bump in the tail.
bump = Plasma.Species {
charge = 0.0, mass = 1.0,
-- Velocity space grid.
lower = {-8.0*vt0}, upper = { 8.0*vt0},
cells = {32},
-- Initial conditions.
init = function (t, xn)
local x, v = xn[1], xn[2]
return bumpMaxwell(x,v,n0,u0,vt0,nb,ub,vtb,uL,sb)
end,
evolve = true, -- Evolve species?
evolveCollisionless = false, -- Evolve collisionless terms?
diagnosticIntegratedMoments = { "intM2Flow", "intM2Thermal" },
-- Collisions.
coll = Plasma.LBOCollisions {
collideWith = {'bump'},
frequencies = {nu},
},
},
}
We run this input file with the call
gkyl lboRelax.lua
On a 2015 MacBookPro this ran in 1.5 seconds and produced a screen output like this one.
We can start looking at the data by first, for example, making a movie
of the distribution function as function of time with pgkyl
:
pgkyl "lboRelax_bump_[0-9]*.bp" interp sel --z0 0. anim -x '$v$' -y '$f(x=0,v,t)$'
(note that postgkyl allows abbreviations,
so interp
= interpolate, sel
= select,
anim
= animate) This command produces the movie given below. We can see that from the
initial, bump-in-tail state the distribution relaxes to a Maxwellian.
The Maxwellian by the way is the analytic steady state of this operator.
Such relaxation should also take place without breaking momentum or
energy conservation. We can examine the evolution of the total energy
in the system by adding intM2Flow
and intM2Thermal
and plotting
it as a function of time. This is achieved in pgkyl
via:
pgkyl lboRelax_bump_intM2Flow.bp lboRelax_bump_intM2Thermal.bp ev 'f[0] f[1] +' pl -x 'time' -y 'energy'
As we can see in the figure below, and in particular in the \(10^{-14}\) scale of it, the total particle energy is conserved very well. The changes in energy over a collisional period are of the order of machine precision.
Example 2: 1x2v collisional Landau damping¶
We now explore the modification of Landau damping by inclusion of Dougherty collisions. Specifically, we will consider ion acoustic waves with adiabatic electrons. This means that the electron number density simply follows
and our gyrokinetic Poisson equation is simply replaced by the quasineutrality
So there is no need to evolve the electron distribution function. In the Gyrokinetic
App we can specify an adiabatic species using Plasma.AdiabaticSpecies
:
plasmaApp = Plasma.App {
...
adiabaticElectron = Plasma.AdiabaticSpecies {
charge = -1.0, mass = mElc,
temp = Te,
-- Initial conditions.. use ion background so that background is exactly neutral.
init = function (t, xn)
return nElc
end,
evolve = false, -- Evolve species?
},
...
}
This simulation then only needs to solve the electrostatic gyrokinetic equations for ions
and we do so with an initial condition that contains a sinusoidal perturbation (wavenumber \(k=0.5\)) in the ion density:
If the right side of this equation (4) were zero, this ion acoustic wave would damp at the collisionless rate calculated by Landau (well he did electron Langmuir waves). But collisions will change the picture and we wish to numerically find out how.
This simulation is setup in the ionSound.lua input file. This input file calls for discretizing the ion phase space \([-\pi/k,\pi/k]\times[-6v_t,6v_t]\times[0,m_i(5v_t^2)/(2B_0)]\) using \(64\times128\times16\) cells and a piecewise linear basis. With a collisionality of \(\nu_=0.005\), the simulation ran on a 2015 MacbookPro in 41 minutes, while a collisionality of \(\nu=0.05\) required 1.35 hours. They were run the command
gkyl ionSound.lua
and produced this screen output.
Note that this is really a linear problem, that is, one can sufficiently model it
with a linearized version of equation (4), using \(f_i=f_{i0}+f_{i1}\),
where the fluctuation \(f_{i1}\) is small compared to the equilibrium
(Maxwellian) \(f_{i0}\). Users may wish to output this fluctuation in time: in
order to to this specify the background with the initBackground
table:
ion = Plasma.Species {
...
-- Specify background so that we can plot perturbed distribution and moments.
initBackground = {"maxwellian",
density = function (t, xn)
return nIon
end,
temperature = function (t, xn)
return Ti
end,
},
...
},
This will output the fluctuation to a file with the name format
<simulation>_<species>_f1_#.bp
, where #
stands for the frame number. So for
example, in this ionSound.lua
case it creates files named ionSound_ion_f1_#.bp
.
We can plot this fluctuation along \(v_\parallel\) at \($t=5$\) with
pgkyl "ionSound_ion_f1_10.bp" interp sel --z0 0.0 --z2 0.0 pl -x '$v_\parallel$' -y '$f_{i1}(x=0,v_\parallel,\mu=0,t=5)$'
(note that postgkyl allows abbreviations,
so interp
= interpolate, sel
= select,
pl
= plot) which produces the following image
Perhaps most valuable to the physics of this simulation is to see a signature of the
decay of the ion acoustic wave. This simulation produced the integrated squared
electrostatic potential, \(\int\mathrm{d}x\,|\phi|^2\), which we take as a measure
of the wave energy. It is stored in a file with the name format <simulation>_phiSq.bp
.
If we had run two simulations, ionSound.lua with \(\nu=0.005\)
and ionSoundH.lua with \(\nu=0.05\), we could plot both electrostatic
energies in time with the following pgkyl
command:
pgkyl ionSound_phi2.bp -l '$\nu=0.005$' ionSoundH_phi2.bp -l '$\nu=0.05$' pl --logy -f0 -x 'time' -y 'Integrated $|\phi|^2$'
Notice that we are giving each file a label to use in the plot with the -l
flag. Postgkyl
then produces the following figure
We thus see that the wave energy is decaying as a function of time (the envelope of the curve is going down), and that the rate at which this happens decreases with collisionality. That is, for this case increasing collisionality decreased the damping rate. From this curve we can also read the period of the wave, using the spacing between the dips.
References¶
E. P. Gross & M. Krook. Model for collision precesses in gases: small-amplitude oscillations of charged two-component systems. Physical Review, 102(3), 593–604 (1956).
J. M. Greene. Improved Bhatnagar-Gross-Krook model of electron-ion collisions. Physics of Fluids, 16(11), 2022–2023 (1973).
J. P. Dougherty. Model Fokker-Planck Equation for a Plasma and Its Solution. Physics of Fluids, 7(11), 1788–1799 (1964).
A. Hakim, et al. (2020). Conservative Discontinuous Galerkin Schemes for Nonlinear Fokker-Planck Collision Operators. Journal of Plasma Physics Vol 86 No. 4, 905860403 (2020), arXiv:1903.08062.
M. Francisquez, et al. (2020). Conservative discontinuous Galerkin scheme of a gyro-averaged Dougherty collision operator. Nuclear Fusion 60 No. 9, 096021 (2020), arxiv:2009.06660.
M. Francisquez, et al. (2022). Improved multispecies Dougherty collisions. Journal of Plasma Physics Vol 88 No. 3, 905880303 (2022), arxiv:2109.10381.