1. Introduces enumerations (and some booleans) intended to replace all
the magic numbers used throughout the code to represent multiple
choices.
2. Replace the gray_params.sh script a new one that automatically
generates code for all the GRAY parameters by parsing
gray_params.f90.
3. Also generate extra code to accept the enum identifiers as valid
values in the configuration files and command line arguments.
4. Set sensible default values for parameters that are rarely changes.
This change modifies the analytical equilibrium in order to simplify the
computation of the poloidal flux normalization and the derivatives.
In the power law parametrisation of the safety factor, ρ_t is replaced
with ρ_p and, similarly, the normalised poloidal radius is now
identified with ρ_p, instead of ρ_t.
With the same parameters (q₀,q₁,α...), this choice slightly changes the
plasma current distribution, but enables us to obtain a closed form for
ψ_a = ψ(r=a) and the relation ρ_t(ρ_p). In fact, both expressions are
now obtained by integrating the q(ρ_p), instead of 1/q(ρ_t), which has
no elementary antiderivative.
As the normalisation is now computed exactly, the values of the
normalised flux ψ_n = ψ/ψ_a and the gradient ∇ψ (entering the raytracing
equations in X and ∇X, respectively) are computed to the same precision.
Previously, ψ_n was computed to a lower precision due to the use of a
simple trapezoid integration of 1/q(ρ_p) for ψ_a, while ∇ψ was computed
up to machine precision using an exact formula.
This error effectively caused a very slight decoupling between X=ω_p²/ω²
and ∇X that introduced a systematic error in the numerical solution of
the raytracing equations.
The error manifests itself as a bias with a weak dependency on X in the
values taken by the dispersion function Λ(r̅, n̅) on the phase-space
points generated by the integrator. More specifically,
lim h→0 Λ(r̅_i, n̅_i) = -kX(r̅_i)
where h is the integrator step size;
r̅_i is the position at the i-th step;
k ≈ -3.258⋅10⁻⁵ and depends only on the number of points used to
perform the trapedoid integral for ψ_a (as ~ 1/n²).
After this change Λ behaves consistently with being a conserved quantity
(zero) up to the cumulative integration error of the 4° order
Runge-Kutta method. In fact we now have that:
Λ(r̅_i, n̅_i) ∝ - h⁴ ‖∂⁴X(r̅_i)/∂r̅⁴‖
It must be said that within this model the relation ρ_p(ρ_t) can't be
computed analytically (inverting ρ_t(ρ_p) produces a trascendental
equation of the form b = x + c x^α). However, this relation is not
necessary for raytracing and is easily solved, up to machine
precision, using minpack.
In addition, this change also makes the model consistetly use the
cocos=3 and fully implements the ability to force the signs of I_p, B_φ
(via equilibrium.sgni,sgnb) and rescaling the field (via
equilibrium.factb).
1. Fix the mismatch between the psnbnd in coreprofiles and gray_core.
This happens whenever gray overrides the externally provided one
(i.e. the density tail would become negative before psnbnd and is so
rescaled to end exactly on the zero).
2. Make psnbnd no longer required by always computing it as in 1.
It hasn't been removed, because gray_params.data is sacrosant,
but it no longer has any effect.
3. Cleanup: mark public functions, restructure the global variables into
three categories; add comments explaining the analytical profiles
format, formulae and how the polynomial tail is computed.
This adds a new configuration file based on the INI format.
The new format will allow adding GRAY parameters without breaking
compatibility with existing configurations, unlike as of the old
gray_params.data.
