gray/input/equil_a.txt

2.96 0.0 1.25  ! R₀,z₀,a,  where ρ_p(R, z) = √[(R - R₀)² + (z - z₀)²]/a  (m)
6              ! B₀,       where B_φ(R) = B₀ R₀/R  (T)
3.5 10 2       ! q₀,q₁,α,  where q(ρ_p) = q₀ + (q₁-q₀)ρ_p^α
0              ! number of points in the limiter contourn, (R,z) pairs

! Notes:
! 1. use B₀>0 for clockwise B_φ
! 2. use q₀,q₁<0 for clockwise I_p,B_φ
-												rework the analytical model

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).

											
										
										
											2023-10-11 17:29:24 +02:00
+.96 0.0 1.25  ! R₀,z₀,a,  where ρ_p(R, z) = √[(R - R₀)² + (z - z₀)²]/a  (m)
 ! B₀,       where B_φ(R) = B₀ R₀/R  (T)
 .5 10 2       ! q₀,q₁,α,  where q(ρ_p) = q₀ + (q₁-q₀)ρ_p^α
 ! number of points in the limiter contourn, (R,z) pairs
 								! Notes:
 								! 1. use B₀>0 for clockwise B_φ
 								! 2. use q₀,q₁<0 for clockwise I_p,B_φ