LinearMaxwellVlasov.jl Documentation

The function parallel gives the same answers for a given set of inputs. Calculate the key given the inputs and the (identity) operator for storing and retrieving values

Identity operator for storing and retreiving values from CacheDict

ColdSpecies <: AbstractFluidSpecies

Cold plasma species. Fields:

• Π Plasma frequency [rad / s]
• Ω Cyclotron Frequency [rad / s]

(c::ConcertinaSinpi)(t)

Evaluate and sum the ConcertinaSinpi function at location (t, v⊥) where t ∈ (0, 1). t is mapped to (0, 1) + az[1]:az[2], so that it takes into account the sign of the denominator, sinpi(t + ax[1]:azp[2]), before finally dividing by |sinpi(t)|.

...

Arguments

• c::ConcertinaSinpi(tv⊥):

...

Configuration(args...)

This holds the frequency, wavenumber and options object of the calculation i.e. everything except the plasma.

...

Arguments

• args...:

...

Example

Kinetic plasma species defined by one coupled distribution function in momentum space such that the relativistic dielectric tensor can be calculated. Fields:

• Π Plasma frequency [rad / s]
• Ω Cyclotron Frequency [rad / s]
• mass Species particle mass [kg]
• F :: AbstractFRelativisticMomentum Distribution function in momentum space parallel and perpendicular to the background magnetic field (normalised)

CoupledRelativisticSpecies(Π,Ω,m,pthz::Number,pth⊥=pthz,pzdrift=0)

...

Arguments

• Π: classical plasma frequency [rad/s]
• Ω: classical cyclotron frequency [rad/s]
• m: mass [kg]
• pthz::Number: parallel thermal momentum [kg m/s]
• pth⊥=pthz: perpendicular thermal moment [kg m/s]
• pzdrift=0: parallel bulk momentum [kg m/s]

...

Kinetic plasma species defined by one coupled distribution function in velocity space, parallel and perpendicular to the background magnetic field Fields:

• Π Plasma frequency [rad / s]
• Ω Cyclotron Frequency [rad / s]
• F :: AbstractCoupledVelocity Distribution function in velocity space parallel and perpendicular to the background magnetic field (normalised)

CoupledVelocitySpecies(Π::Float64,Ω::Float64,vthz::Float64,vth⊥::Float64=vthz,vzdrift::Float64=0.0,v⊥drift::Float64=0.0)

...

Arguments

• Π::Float64: plasma frequency [rad/s]
• Ω::Float64: cyclotron frequency [rad/s]
• vthz::Float64: parallel thermal speed [m/s]
• vth⊥::Float64=vthz: perpendicular thermal speed [m/s]
• vzdrift::Float64=0.0: parallel bulk speed [m/s]
• v⊥drift::Float64=0.0: perpendicular bulk speed [m/s]

...

FBeam(vth::T,vd::U=0.0)where{T<:Number,U<:Number}

Create a suitably normalised parallel distribution function proportational to exp(-(vz-vd)^2/vthz^2)

...

Arguments

• vth::T: thermal velocity [m/s]
• vd::U=0.0: bulk parallel velocity [m/s]

...

FCoupledVelocityNumerical

A disribution function where vz and v⊥ are coupled, i.e. non-separable.

...

Arguments

• F::T: the distrubtion function
• normalisation::Tuple{U,U}: the speeds used for normalisation in parallel and perp [m/s]
• lower::Float64: minimum speed for integration bounds [m/s]
• upper::Float64: maximum speed for integration bounds [m/s]

...

Example

vth = 1e4
vshell = 1e5
fshell = FShell(vth, vshell)
lower = vshell - 12 * vth # where the shell is zero
upper = vshell + 12 * vth # where the shell is zero
FCoupledVelocityNumerical(fshell, (vshell, vshell), lower, upper)

The Dirac-delta function parallel distribution functions

Type for an arbitrary parallel distribution function that holds, a function for the distribution function, it's derivative and the limits of integral

The Dirac-delta function parallel distribution functions

Type for an arbitrary perpendicular distribution function that holds, a function for the distribution function, it's derivative and the limits of integral

FRing(vth::T,vd::U=0.0)where{T<:Number,U<:Number}

Create a suitably normalised ring distribution in perpendicular velocity proportional to exp(-(v⊥-vd)^2/vthz^2)

...

Arguments

• vth::T: thermal velocity [m/s]
• vd::U=0.0: bulk perpendicular velocity [m/s]

...

The parallel integral of the Beam only requires an integral over a drifting Maxwellian subject to the relevant kernels. All this is calculated here.

Perpendicular integrals for a Maxwellian

NeutralPlasma checks that the plasma is neutral

Options(::Type{T}=Float64;kwargstuple...)where{T}

Create an Options object a set of kwargs. Best to read the code to see what counts as a duplicate. Copying that logic here will only create broken docs in the future.

...

Arguments

• ::Type{T}=Float64: the number type to express the tolerances in

Optional args:

• kwargstuple...: the keyword arguments that are preferred over the defaults

...

The kernel of the perpendicular integral

Plasma doesn't check if the plasma is neutral

Relativestic dielectric tensor as a function of cyclotron harmonic This is only used for the principal part.

RelativisticMaxwellian Return a drifting Maxwellian function.

Members

pthz::Real - thermal momentum parallel to magnetic field [kg m/s] pth⊥::Real - thermal momentum perpendicular to magnetic field [kg m/s] pzdrift::Real=0.0 - drift parallel to the magnetic field [kg m/s] lognormalisation::Real - the log of the normalisation constant

Kinetic plasma species with separable distribution functions parallel and perpendicular to the magnetic field. Fields:

• Π Plasma frequency [rad / s]
• Ω Cyclotron Frequency [rad / s]
• Fz :: AbstractFParallel Distribution function parallel to magnetic field (normalised)
• F⊥ :: AbstractFPerpendicular Distribution function perpendicular to magnetic field (normalised)

The probability density at parallel velocity, v, of a shifted maxwellian from thermal velocity vth, and drift velocity vd.

The probability density at perpendicular velocity, v, of a shifted maxwellian from thermal velocity vth, and drift velocity vd.

Transform a function from domain [-∞, ∞]ⁿ down to [-1, 1]ⁿ

Transform a function to transform an integrand from domain [-∞, 0]×[∞, ∞] down to [-1, -π/2]×[1, π/2].

Example:

julia> f(x) = exp(-(x[1]^2 + x[2]^2)/2) * cos(x[2])^2 * sin(x[1])^2
julia> hcubature(f, [-12.0, 0.0], [12.0, 12.0])
(0.7710130943379178, 1.1482318484139944e-8)
julia> hcubature(UnitSemicircleIntegrandTransform(f, 2.0), [0, -π/2], [1, π/2])
(0.7710130940919337, 1.1464196819507393e-8)

Warm plasma species, with speeds of sound parallel and perpendicular to the magnetic field. Fields:

• Π Plasma frequency [rad / s]
• Ω Cyclotron Frequency [rad / s]
• soundspeed Sound speed [m/s]

WarmSpecies(Π::Float64,Ω::Float64,thermalspeed::Float64,adiabiaticindex::Number)

Warm plasma species - accept thermalspeed and ratio of specific heats to get sound speed

...

Arguments

• Π: Plasma frequency [rad / s]
• Ω: Cyclotron Frequency [rad / s]
• thermalspeed: Thermal speed of Maxwellian distribution
• adiabiaticindex: Equation of state Gruneisen gamma (ratio of specific heats)

...

WarmSpecies(s::T,γ=5/3)where{T<:AbstractSeparableVelocitySpecies{

Create a warm plasma species from the AbstractSeparableVelocitySpecies

...

Arguments

• s::T:
• γ=5/3where{T<:AbstractSeparableVelocitySpecies{:

...

Wavenumber decomposed into parallel and perpendicular components Construct with parallel and wavenumber components, or by keyword arguement pairs

• parallel & perpendicular
• wavenumber & propagationangle

Regarding integrals with besselj(n, x)* besselj(n±1, x), e.g.: where m >= 0

besselj(-m,x) * besselj(-m-1,x) == -besselj(m,x) * besselj(m+1,x)

besselj(-m+1,x) * besselj(-m-1,x) == +besselj(m-1,x) * besselj(m+1,x)

Base.get!(f::F,data::Dict{UInt64,CacheDict{C}},species,config)where{F<:Function,C<:AbstractCacheType}

Return the correct cache dictionary for the given species and configuration.

This creates a CacheKeyOp with the arguments given so that it can create a CacheDict with the correct type, which then populates the CacheDict if it's empty. This process acts as a function barrier to get better type inference.

...

Arguments

• f::F: The function to cache
• data::Dict{UInt64: the dictionary of caches
• species: The species argument to f
• config: The configuration argument to f

...

Base.get!(f::F,cache::CacheDict{C,V},i...)::Vwhere{F<:Function,C<:AbstractCacheType,V}

description

...

Arguments

• f::F:
• cache::CacheDict:
• i...:

...

FShell(vth::Real,vshell::Real)

The shell distribution function, as though f is only non-zero on or close to the surface of a sphere.

...

Arguments

• vth::Real: the thermal velocity of the shell [m/s]
• vshell::Real: the speed of the shell [m/s]

...

Example

FSlowingDown(vbeam::Real,vcrit::Real,vcutoffwidth::Real)

The slowing down distribution

...

Arguments

• vbeam::Real: the beam speed [m/s]
• vcrit::Real: the critival velocity [m/s]
• vcutoffwidth::Real: the width of the error function used to smooth

the distribution function at vbeam [m/s] ...

Example

The non-integer argument to the BesselJ

MaxwellianSpecies(Π,Ω,vthb,vth⊥=vthb,vdb=0.0)

Kinetic Maxwellian Plasma species that can optionally have a drift along the magnetic field

...

Arguments

• Π: Plasma frequency [rad / s]
• Ω: Cyclotron Frequency [rad / s]
• vthb: parallel thermal speed [m/s]
• vth⊥=vthb: perpendicular thermal speed [m/s]
• vdb=0.0: parallel beam speed [m/s]

Returns

• SeparableVelocitySpecies(Π, Ω, FBeam(vthb, vdb), FPerpendicularMaxwellian(vth⊥))

...

Example

Create a kinetic plasma species with separable distribution functions parallel $f(v_\parallel)$ and perpendicular $f(v_\perpendicular)$ to the magnetic field, which are defined as a drifting beam and a ring respectively. ...

Arguments

• Π: Plasma frequency [rad / s]
• Ω: Cyclotron Frequency [rad / s]
• vthb: parallel thermal speed [m/s]
• vth⊥=vthb: perpendicular thermal speed [m/s]
• vdb=0.0: parallel beam speed [m/s]
• vd⊥=0.0: perpendicular ring speed [m/s]

Returns

• SeparableVelocitySpecies(Π, Ω, FBeam(vthb, vdb), FRing(vth⊥, vd⊥))

...

Example

Calculate the Alfven speed given Π_Ωs, which is an iterable container of the ratio of the plasma to the cyclotron frequency of all the species of the plasma

The conducivity tensor for a given species

Calculate the unitless susceptibility tensor contribution for a coupled velocity distribution function, without the need to sum over harmonics.

Calculate the unitless susceptibility tensor for a warm plasma species Swanson 3.63

Calculate the converged unitless susceptibility tensor contribution for a maxwellian distribution function, having summed over bessel indices n.

Calculate the converged unitless susceptibility tensor contribution for a maxwellian distribution function, having summed over bessel indices n.

Calculate the unitless susceptibility tensor for a cold plasma species

Calculate the unitless susceptibility tensor for a maxwellian distribution function given the bessel indices n.

converge(f::T,tol::Tolerance=Tolerance()) where {T}

Sum the output of f starting with argument 0 and (1,-1), (2,-2)... etc until the convergence criterion is met based on the tolerance

...

Arguments

• f::T: input function or functor
• tol::Tolerance=Tolerance: Tolerance object containing relative and absolute tolerances

...

Example

julia> using Symbolics julia> @syms kx ky kz; julia> k = [kx, ky, kz]; julia> curl = im * [0 -k[3] k[2]; k[3] 0 -k[1]; -k[2] k[1] 0]; julia> curl * curl 3×3 Matrix{Any}: ky^2 + kz^2 -kxky -kxkz -kxky kx^2 + kz^2 -kykz -kxkz -kykz kx^2 + ky^2

cyclotronfrequency(B::Number,mass::Number,Z::Number=1)

...

Arguments

• B::Number: magnetic field in T
• mass::Number: mass in kg
• Z::Number=1: charge number

...

defaults(::Type{T}=Float64)where{T}

The default settings for the Options object.

...

Arguments

• ::Type{T}=Float64where{T}: the number type to express the tolerances in

...

The dielectric tensor for a given plasma

The contribution to the dielectric tensor for a given species

discretefouriertransform(f::T,n::Int,N=512)where{T<:Function}

calculate the DFT of function f at Fourier mode n with N points in the interval [0,2π]

...

Arguments

• f::T: DFT of this function
• n::Int: DFT of this Fourier mode
• N=512: evaluate f at this many points

...

The electrostatic dielectric tensor for a given plasma, a zero valued determinant of which represents a solution to the linear poisson-vlasov system of equations

fastisapprox(n,x,y,atol,rtol,nans)

...

Arguments

• n: normalisation
• x: compare this
• y: with that
• atol: absolute tolerance
• rtol: relative tolerance
• nans: whether to treat NaNs as equal

...

Example

Takes a function that when integrated between -Inf and +Inf returns value x, and returns a new function that returns x when integrated between real(pole) and +Inf.

Integration of Abstract Arbitrary FParallel

The integrals over the distribution functions and associated integral kernel We avoid doing integrals involving the derivative of the Dirac delta function by doing integral by parts and knowing that f-> 0 and the integral limits

Calculate the parallel integral with principal part and residue

The integrals over the distribution functions and associated integral kernel We avoid doing integrals involving the derivative of the Dirac delta function by doing integral by parts and knowing that f-> 0 and the integral limits

Integrate over the parallel distribution function multiplied by various kernels. If the imaginary part of the pole is zero, then do a trick that folds over the integral from the left of the pole to right. We need to take into account whether the real part of the pole is negative, otherwise positive slopes at negative velocities give rise to instability and not damping.

integrate(f::AbstractFPerpendicular,kernel::T,∂F∂v::Bool,tol::Tolerance=Tolerance()

Integration of Abstract Arbitrary FPerpendicular with respect to kernel where the distribution function may be differentiated. If quadrature is involved it will adaptively calculate the integral until the tolerance is met.

...

Arguments

• f::AbstractFPerpendicular:
• kernel::T: the kernel e.g. v⊥^2 * BesselJ(n, v⊥^2 * β^2)^2 etc.
• ∂F∂v::Bool: differentiate with respect to velocity, or not.
• tol::Tolerance=Tolerance:

...

Example

isneutral(vs::AbstractVector{T},atol=100eps()

Calculate if a plasma is neutral

...

Arguments

• vs::AbstractVector{T}: A vector of species
• atol=100eps(): tolerance for deciding upon neutrality

...

Transform the limits of an integrand quadrature(foldnumeratoraboutpole(integrand, pole), limitsfolder(limits, pole)...)

Tool to normalize a function f between two integral limits a and b

Interface The parallel integral of the distribution function and the relevant kernel

Parallel integrals for the Beam - used for testing purposes

Compile the kernels of the parallel integral and fetch it from the DistributionFunctions module. If k parallel is zero then there is no Landau damping, so this case is made separate, as is easier to deal with.

parallel_integral(species::S,config::Configuration,parallelcaches=Dict{UInt64,CacheDict{ParallelCache}}()) where {S}

Memoise the inputs to the integrals; export the memoised function, and the dictionary of inputs -> outputs. This uses a parallel cache which holds the logic for the integrals and ensure that the values are only calculated once. It only matters what n-m is, not the values of n and m individually So change (n, m)->(n-m, 0) to do fewer calculations

...

Arguments

• species::S: the species to do parallel integrals over
• config::Configuration: the configuration holding frequency and wavenumber, and options.
• parallelcaches=Dict{UInt64,CacheDict{ParallelCache}}: The dict to store parallel integrals

...

Example

Interface The perpendicular integral of the distribution function and the relevant kernel

Calculate the definite integral kernels of the perpenedicular distribution function, or derivative thereof, multiplied by the kernel functions

Memoise the inputs to the integrals; export the memoised function, and the dictionary of inputs -> outputs besselj harmonic numbers n and l can be any order Use this to do only half the calculatons (n, l)->(min(n, l), max(n, l)) Also know that besselj(-n, x) = (-1)^n * besselj(n, x)

Return the value of the plasma dispersion function This implementation includes the residue, which is easy to verify because Z(0) = im sqrt(π).

• x is the argument to the plasma disperion function
• power is the moment of the integral

[1] S.D. Baalrud, Phys. Plasmas 20, 012118 (2013) and put ν = -Inf [2] M. Sampoorna et al., Generalized Voigt functions and their derivatives, Journal of Quantitative Spectroscopy & Radiative Transfer (2006), doi:10.1016/j.jqsrt.2006.08.011

plasmafrequency(n::Number,mass::Number,Z::Number=1)

...

Arguments

• n::Number: species number density in m^-3
• mass::Number: particle mass in kg
• Z::Number=1: Charge number

...

residue(principalpart,pole::Number)

...

Arguments

• principalpart:
• pole::Number:

...

residue(numerator::T,pole::Number)where{T<:Function}

Calculate the residue of a number function at a pole

...

Arguments

• numerator::T:
• pole::Number:

...

Expressing f(x) = ∑ᵢ aᵢ (x - p)ⁱ find a₋₁

residuepartadaptive(f::T,pole::Number,radius::Real=(isreal(pole)?abs(pole):imag(pole))*sqrt(eps()),

Expressing f(x) = ∑ᵢ aᵢ (x - p)ⁱ find a₋₁ by doing a Laurent transform via discrete Fourier transform by evaluating f on walk around the pole. The number of evaluations grows with each loop until the tolerance has been met.

...

Arguments

• f::T: function to calculate residue of...
• pole::Number: ... at this pole
• radius::Real=(isreal(pole) ? abs(pole) : imag(pole)) * sqrt(eps()): radius of the residue
• N::Int=64: the number of points to evaluate the residue initially
• tol::Tolerance=Tolerance(): Tolerance object

...

residuesigma(pole::Number) = imag(pole) < 0 ? 2 : imag(pole) == 0 ? 1 : 0

Calculate the "sigma" factor of the residue

Calculate the tensor representing the linear Maxwell-Vlasov set of equations. The determinant is zero when the wavenumber and frequency represent a solution to the linear Maxwell-Vlasov system of equations for these species. Note that the curlcurl operator keeps track of the factor of -1 from im^2.

thermalspeed(ϵ_V::Number, mass::Number)

Thermal speed defined as sqrt(2 * q₀ * ϵ_V / mass) ''' # Arguments

• ϵ_V::Number: energy in electron Volts
• mass::Number: particle mass in kg

'''

Concertina and rescale a function in the sections between its roots such that all roots lie at -1 or +1. Integration over -1..1 of the resulting function gives the same answer as integration over original -Inf..Inf domain

wavedirectionalityhandler(x::Number,kz::Number)

Take into account the size of kz when performing parallel integrals

Most codes, assume that kz is real, but this is not always the case. Usually, they take the absolute value of kz, but this does not respect complex kz in convetive instability calculations, for example

...

Arguments

• x::Number:
• kz::Number:

...

Eq. 25 of R. O. Dendy, C. N. Lashmore-Davies, and K. G. McClements, G. A. Cottrell, The excitation of obliquely propagating fast Alfven waves at fusion ion cyclotron harmonics, Phys. Plasmas 1 (6), June 1994