Poorly conditioned initial field calculation

If the perturbed distribution function, f_1, is given by

f_1,s = (q_s phi/m_s) d f_0,s / dE + h_s exp(-ik.rho_s) [1]

The we can write our gyrokinetic Poisson’s equation as

-phi Sum_Species( n_s q_s^2 / T_s) + Sum_species(q_s * velocity_integral(J0 h_s)) = 0 [2]

However, in GS2 we don’t evolve the non-adiabatic contribution, h_s, rather we evolve (assuming an electrostatic simulation).

g_s = h_s - (q_s phi / T_s) J0 f_0,s

This means in practice that our Poisson’s equation looks like

-phi Sum_Species( n_s q_s^2 / T_s - q_s^2/T_s velocity_integral(J0^2 f_0,s)) + Sum_species(q_s * velocity_integral(J0 g_s)) = 0 [3]

‌

During initialisation we default to setting g_s and hence calculate the initial field by inverting [3]:

phi = Sum_species(q_s * velocity_integral(J0 g_s)) / Sum_Species( n_s q_s^2 / T_s - q_s^2/T_s velocity_integral(J0^2 f_0,s)) [4]

Suppose J0 → 1 (i.e. the argument goes to 0), clearly the denominator in [4] goes to 0 and this initial phi is not well defined. This can occur when kperp → 0 (e.g. very low ky or kx). In practice this means that as kperp → 0 the initial phi can potentially become very large (before our check for small denominator kicks in and sets phi to 0). As we set all species to have the same initial g_s I think we would in fact expect phi → 0 as kperp → 0 (expand J0 ~ 1 - kperp^2/4 + …, the leading order terms then cancel in the numerator and we’re left with some thing which decays to zero as kperp is reduced).

If we instead initialise h_s our initial field is found by inverting [2]:

phi = Sum_species(q_s * velocity_integral(J0 h_s)) / Sum_Species( n_s q_s^2 / T_s) [5]

which should be much better behaved.

‌

There is already a flag (https://gyrokinetics.gitlab.io/gs2/page/namelists/index.html#init_g_knobs-initial_condition_is_nonadiabatic_dfn) which can be set to true to tell GS2 to interpret the initial value as h_s rather than g_s, however we don’t yet calculate the initial fields directly from this. Instead we calculate the fields from h_s, use this to form g_s and then recalculate the fields using the “usual” approach (i.e. eqn 4).

Here’s a plot illustrating this issue. This shows the initial phi2 as a function of ky for a number of approaches. Blue and orange are high and low pitch angle resolution runs where initial_condition_is_nonadiabatic_dfn is true and we follow the process outlined above to first set g_s and to then calculate the fields. Red is where initial_condition_is_nonadiabatic_dfn is false (i.e. the default) and green is where initial_condition_is_nonadiabatic_dfn is true and we keep the initial fields calculated directly from this. These runs are electromagnetic.

‌

Actions arising:

Work out how to get GS2 to keep the initial fields calculated from h_s if we set h_s as the initial condition
Consider if https://gyrokinetics.gitlab.io/gs2/page/namelists/index.html#init_g_knobs-initial_condition_is_nonadiabatic_dfn should be the default.

Areas of possible impact:

Generally results should be insensitive to the initial condition, however there are a few reasons this might matter

By default, on restart (or timestep change) we recalculate the fields directly from g_s using the same routine as to find the initial fields. If this is poorly conditioned could this cause problems (e.g. if the denominator is sufficiently close to zero then we set phi = 0 here – this has the potential to zero our fields for low kx/ky on time step change/restart).
If the initial fields are not fully consistent might this introduce unanticipated structure in the fields (e.g. a smooth g_s giving rise to a noisy phi). For growing modes we may not care as the final phi will be dominated by that coming from the time advance, but for damped modes might we be more sensitive to this?

‌

Comments (3)