This changes all the print_* subroutines to not depend on q_spline,
either for the ψ data points or just their number.
In fact, in the new analytical model q_spline doesn't exist anymore.
Pro: This allows to mark q_spline as private.
This change combines `equian` and `equinum_psi` into a new `pol_flux`
subroutine that computes ψ(R,z) and derivatives for either numerical or
analytical equilibria.
Similarly, `equian` (the fpol and dfpol outputs) and `equinum_fpol` are
combined intro `pol_curr` that computes F(ψ) and its derivative.
Callers of these subroutines do not select a specific version based on
the value of `iequil` anymore.
- Hide the implementation of the re-normalisation of the ψ(R,z) spline
by adjusting the spline coefficients, instead of shifting and
rescaling after each evaluation.
- Correct the value of `psia` = ψ(X point) - ψ(O point) after the ψ(R,z)
spline has been re-normalised.
This fixes another instance of decoupling between the values of X and
∇X that introduce a systematic error in the numerical integration of
the raytracing equations.
Here the issue is caused by the different normalisations used,
specifically X=X(ψ_n') and ∇X=∇X(ψ_n), where ψ_n' is the re-normalised
spline of the normalised flux and ψ_n the spline of the normalised
flux.
See d0a5a9f for more details on problems resulting from this error.
- Change `equian` and `equinum_psi` to return the normalised values
for both the flux and its derivatives, to avoid confusions.
Callers that needed the unnormalised derivatives now multiply
explicitly by `psia`.
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).
The ρ_p/ρ_t mapping is 1:1, so the interpolation must always preserve
monotonicity, of which cubic splines generally make no guarantee.
Note: Linear interpolation does not provide even C¹ continuity, but
these data is not directly used in the numerical integration, so it
should be fine. Ideally this should be replaced with cubic splines
computed with the Fritsch–Carlson algorithm.
This adds a new `splines` module which implements a high-level interface
for creating and evaluating splines and rewrite almost all modules to
use it. Also, notably:
1. both `simplespline` and DIERCKX splines can now used with a uniform
interface
2. most complexity due to handling working space arrays is gone
3. memory management has been significantly simplified too
This change structures the arguments of most functions, in particular
gray_main, into well-defined categories using derived types.
All types are defined in the gray_params.f90 (location subject to
change) and are organised as follows:
gray_parameters (statically allocated data)
├── antenna_parameters
├── ecrh_cd_parameters
├── equilibrium_parameters
├── misc_parameters
├── output_parameters
├── profiles_parameters
└── raytracing_parameters
gray_data - inputs of gray_main (dynamically-allocated arrays)
├── equilibrium_data
└── profiles_data
gray_results - outputs of gray_main (dynamically-allocated arrays)
