Discussion: Foot-point search behavior when invariant returns baddata

[lkncl]

Description:

I would like to open a discussion about the behavior of the foot-point search routine used in our field-line–based calculations (make_lstar1, drift_shell1, drift_bounce_orbit1).

The foot-point search proceeds by scanning different foot-point values of $$\theta{}$$, and computing the normalized second adiabatic invariant between the mirror points $$B_{mir}$$, using:

$$ \hat{J} = \int_{s_{m1}}^{s_{m2}} \sqrt{1-\frac{B(s)}{B_{mir}}} ds $$

If the computed $$\hat{J}$$ is smaller than the target invariant $$\hat{J}_0$$, the algorithm continues toward the north (step dtet=pi/Ntet), otherwise, it moves toward the south.
The search stops when the bracketing criterion is satisfied:

$$ (\hat{J}_1 - \hat{J}_0) * (\hat{J}_2 - \hat{J}_0) < 0 $$

i.e., when $$\hat{J}_0$$ lies between the two most recent values of $$\hat{J}_1$$ and $$\hat{J}_2$$, meaning that the corresponding $$\theta$$ are one dtet apart. The returned value for $$\theta{}$$ is taken as the midpoint of the bracket.

Observed behavior / potential issue

For small target values $$\hat{J}_0$$, this means that when sweeping from low to high $$\theta{}$$, the algorithm encounters a sequence of baddatavalues until it reaches the first $$\theta{}$$ for which the invariant computation succeeds. However, with the fixed step size dtet, this first valid $$\theta{}$$ may already correspond to a value of $$\hat{J}$$ that is significantly larger that the target. In practice, the routine therefore tends to return the first $$\theta{}$$ that does not fail rather than a $$\theta{}$$ that truly brackets the target invariant.

Example input that triggers this behavior:

For the inputs conditions detailed in the file below, many invariant values are baddata, but no error is triggered.

preview :

	iphi			phi				          theta				        leI (J1)		                 leI1 (J2)
           2 -0.34503019249923372       0.42694916024616159        1.5063455022609258E-002  -9.9999999999999996E+030
           3  -9.3702780212050252E-002  0.39422423677126789       -9.9999999999999996E+030  0.25531051056935800     
           4  0.15762463207513322       0.36804429799135296       -9.9999999999999996E+030  0.13395514395823421     
           5  0.40895204436231669       0.34840934390641676       -9.9999999999999996E+030   6.1391347097208136E-002
...

(see here for input/output details : output_std.txt)

Here, $$\hat{J}_0 $$ is equal to I = 0.00164629 and the resulting bounds are quite far apart, giving Lstar = 7.70959085.

Additional tests

I tried increasing the theta search depth with a geometrically decreasing step when one bound is baddata (division by 2, max depth 10). The final L* changes slightly (Lstar = 7.69832289), but for this particular point, is allowed to have true bracketing for the invariant.

preview :

	iphi			phi				          theta					  leI (J1)		                   leI1 (J2)		     depth
           2 -0.34503019249923372       0.42494646310642209        3.0106465265108888E-003   4.7356439528083226E-004         8
           3  -9.3702780212050252E-002  0.39308227329584211        3.3587277748371405E-004   1.7011950836060802E-002         4
           4  0.15762463207513322       0.36667223739774446       0.15294992340812927        3.4763669178272276E-004         1
           5  0.40895204436231669       0.34711398235220259        1.2619646341398830E-004   2.6830852320370675E-003         6
...

(see here for input/output details : output_refined.txt)

Additional remark
The computation of L* may also fail if the invariant evaluation returns baddata for a value of $$\theta{}$$ below the solution. This situation can arise at high latitude, where $$B_{mir}$$ can become very large. In such cases, the southern mirror point can fall inside the Earth, triggering this condition:

c Pourquoi?  "why?"
IF (rr2.LT.1.D0) THEN
         leI = baddata
ENDIF

(This can happen, for example, because of the South Atlantic Anomaly.)

Since baddata is negative, it is interpreted as being below the target invariant. As a result, the search algorithm moves to the next $$\theta{}$$ in the northward direction, even though the true solution may actually lie south of the failing point. This can cause the routine to skip over the correct bracketing region entirely.

Discussion points

Should this behavior remain as-is, given that the algorithm can terminate with one side of the bracket equal to baddata ?
Should the routine return some explicit failure indicator ?
Moreoever, when one of the bracketing invariant values is baddata because the southern mirror point lies into the Earth, should it return a negative L* value, as it is the expected behavior ?

Alternatively, would it be preferable to handle baddatamore robustly, for instance by ignoring invalid bounds in the bracketing test, or by retrying the search until two valid invariant evaluations are obtained ? I have already experimented with a fail-tolerant dichotomy based method for root finding that could be adapted here and with the adaptative dtet step if we want to keep the current step-by-step search.

I can also run a larger parameter sweeps to quantify how often this issue occurs and estimate the resulting L* error.

I am opening this issue to gather feedback on whether the current behavior is intended (bug or feature?) and whether a more robust approach would be desirable.