analysis

Check the (kinetic) Chodura condition

Chodura condition is: ∫d^3v F/vpa^2 ≤ mi ne/Te

Return a tuple (whose first entry is the result for the lower boundary and second for the upper) of the ratio which is 1 if the Chodura condition is satisfied (with equality): Te/(mi ne) * ∫d^3v F/vpa^2

Currently only evaluates condition for the first species: is=1

2D2V

In normalised form (normalised variables suffixed with 'N'): vpa = cref vpaN vperp = cref vperpN ne = nref neN Te = Tref TeN F = FN nref / cref^3 pi^3/2 cref = sqrt(2 Tref / mi)

cref^3 ∫d^3vN FN nref / cref^3 pi^3/2 cref^2 vpaN^2 ≤ mi nref neN / Tref TeN nref / (pi^3/2 cref^2) * ∫d^3vN FN / vpaN^2 ≤ mi nref neN / Tref TeN mi nref / (pi^3/2 2 Tref) * ∫d^3vN FN / vpaN^2 ≤ mi nref neN / Tref TeN 1 / (2 pi^3/2) * ∫d^3vN FN / vpaN^2 ≤ neN / TeN 1 / (2 pi^3/2) * ∫d^3vN FN / vpaN^2 ≤ neN / TeN TeN / (2 neN pi^3/2) * ∫d^3vN FN / vpaN^2 ≤ 1

Note that integrate_over_vspace() includes the 1/pi^3/2 factor already.

1D1V

The 1D1V code evolves the marginalised distribution function f = ∫d^2vperp F so the Chodura condition becomes ∫dvpa f/vpa^2 ≤ mi ne/Te

In normalised form (normalised variables suffixed with 'N'): vpa = cref vpaN ne = nref neN Te = Tref TeN f = fN nref / cref sqrt(pi) cref = sqrt(2 Tref / mi)

cref ∫dvpaN fN nref / cref sqrt(pi) cref^2 vpaN^2 ≤ mi nref neN / Tref TeN nref / (sqrt(pi) cref^2) * ∫dvpaN fN / vpaN^2 ≤ mi nref neN / Tref TeN mi nref / (sqrt(pi) 2 Tref) * ∫dvpaN fN / vpaN^2 ≤ mi nref neN / Tref TeN 1 / (2 sqrt(pi)) * ∫dvpaN fN / vpaN^2 ≤ neN / TeN 1 / (2 sqrt(pi)) * ∫dvpaN fN / vpaN^2 ≤ neN / TeN TeN / (2 neN sqrt(pi)) * ∫dvpaN fN / vpaN^2 ≤ 1

Note that integrate_over_vspace() includes the 1/sqrt(pi) factor already.

If ir0 is passed, only load the data for as single r-point (to save memory).

If find_extra_offset=true is passed, calculates how many entries of f_lower/f_upper adjacent to $v_∥=0$ would need to be zero-ed out in order for the condition to be satisfied.

Fit a cosine to a 1d array

Fit function is Acos(2πn(z + δ)/L)

The domain z is taken to be periodic, with the first and last points identified, so L=z[end]-z[begin]

Arguments

z : Array 1d array with positions of the grid points - should have the same length as data data : Array 1d array of the data to be fit amplitudeguess : Float Initial guess for the amplitude (the value from the previous time point might be a good choice) offsetguess : Float Initial guess for the offset (the value from the previous time point might be a good choice) n : Int, default 1 The periodicity used for the fit

Returns

amplitude : Float The amplitude A of the cosine fit offset : Float The offset δ of the cosine fit error : Float The RMS of the difference between data and the fit

Fit delta_phi to get the frequency and growth rate.

Note, expect the input to be a standing wave (as simulations are initialised with just a density perturbation), so need to extract both frequency and growth rate from the time-variation of the amplitude.

The function assumes that if the amplitude does not cross zero, then the mode is non-oscillatory and so fits just an exponential, not exp*cos. The simulation used as input should be long enough to contain at least ~1 period of oscillation if the mode is oscillatory or the fit will not work.

Arguments

z : Array{mkfloat, 1} 1d array of the grid point positions t : Array{mkfloat, 1} 1d array of the time points deltaphi : Array{mkfloat, 2} 2d array of the values of delta_phi(z, t)

Returns

phifitresult struct whose fields are: growthrate : mkfloat Fitted growth rate of the mode amplitude0 : mkfloat Fitted amplitude at t=0 frequency : mkfloat Fitted frequency of the mode offset0 : mkfloat Fitted offset at t=0 amplitudefiterror : mkfloat RMS error in fit to ln(amplitude) - i.e. ln(A) offsetfiterror : mkfloat RMS error in fit to offset - i.e. δ cosinefiterror : mkfloat Maximum of the RMS errors of the cosine fits at each time point amplitude : Array{mkfloat, 1} Values of amplitude from which growthrate fit was calculated offset : Array{mk_float, 1} Values of offset from which frequency fit was calculated

Get 1D Fourier transform (in r) of nonuniformdata

First interpolates to uniform grid, then uses FFT.

If zind is not given, find the zind where mode seems to be growing most strongly.

Get 2D Fourier transform (in r and z) of nonuniformdata

First interpolates to uniform grid, then uses FFT

Return (v - mean(v, dims=2))

Get the unnormalised distribution function and unnormalised ('lab space') coordinates.

Inputs should depend only on z and vpa.

Get the unnormalised distribution function and unnormalised ('lab space') dzdt coordinate at a point in space.

Inputs should depend only on vpa.

Calculate a moving average

result[i] = mean(v[i-n:i+n])

Except near the ends of the array where indices outside the range of v are skipped.

steady_state_residuals(variable, variable_at_previous_time, dt;
                       epsilon=0.0001, use_mpi=false)

Calculate how close a variable is to steady state.

Calculates several quantities. Define the 'squared absolute residual' $r_\mathrm{abs}(t)^2$ for a quantity $a(t,x)$ as

$r_\mathrm{abs}(t)^2 = \left( a(t,x) - a(t - \delta t,x) \right)$

and the 'squared relative residual' $r_\mathrm{rel}(t)^2$

$r_\mathrm{rel}(t)^2 = \left( \frac{a(t,x) - a(t - \delta t,x)}{\delta t \left| a(t,x) + \epsilon \max_x(a(t,x)) \right|} \right)$

where $x$ stands for any spatial and velocity coordinates, and the offset $\epsilon \max_x(a(t,x))$ is used to avoid points where $a(t,x)$ happens to be very close to zero from dominating the result in the 'squared relative residual', with $max_x$ being the maximum over the $x$ coordinate(s). Returns an OrderedDict containing: the maximum 'absolute residual' $\max_x\left( \sqrt{r_\mathrm{abs}(t)^2} \right)$ ("RMS absolute residual"); the root-mean-square (RMS) 'absolute residual' $\left< \sqrt{r_\mathrm{abs}(t)^2} \right>_x$ ("max absolute residual"); the maximum 'relative residual' $\max_x\left( \sqrt{r_\mathrm{rel}(t)^2} \right)$ ("RMS relative residual"); the root-mean-square (RMS) 'relative residual' $\left< \sqrt{r_\mathrm{rel}(t)^2} \right>_x$ ("max relative residual").

variable gives the value of $a(t,x)$ at the current time, variable_at_previous_time the value $a(t - \delta t, x)$ at a previous time and dt gives the difference in times $\delta t$. All three can be arrays with a time dimension of the same length, or have no time dimension.

By default runs in serial, but if use_mpi=true is passed, assume MPI has been initialised, and that variable has r and z dimensions but no species dimension, and use @loop_* macros. In this case the result is returned only on global rank 0. When using distributed-memory MPI, this routine will double-count the points on block boundaries.

If only_max_abs=true is passed, then only calculate the 'maxium absolute residual'. In this case the OrderedDict returned will have only one entry, for "max absolute residual".

steady_state_square_residuals(variable, variable_at_previous_time, dt;
                              variable_max=nothing, epsilon=1.0e-4,
                              use_mpi=false, only_max_abs=false)

Used to calculate the mean square residual for steady_state_residuals.

Useful to define this separately as it can be called on (equally-sized) chunks of the variable and then combined appropriately. If this is done, the global maximum of abs.(variable) should be passed to variable_max.

See steady_state_residuals for documenation of the other arguments. The return values of steady_state_residuals are the square-root of the return values of this function.