Poloidal Transit Time discrepancy between full orbit and guiding center orbit.

[lstagner]

Notice the full orbit code should be integrating for the same amount of time but they don't match up.

using Equilibrium, GuidingCenterOrbits, PyPlot
δ = 0.33
ϵ = 0.32
κ = 1.7
B0=2 #5.3
R0=6.2
qstar = 1.57
alpha = -0.155
S = solovev(B0, R0, ϵ, δ, κ, alpha, qstar)

gcp = GCAlpha(350,0.5,8.0,0.0)
o = get_orbit(S,gcp,maxiter=0)
path, stat = get_full_orbit(S,gcp,tmax=1.0*o.tau_p/cyclotron_period(S,gcp));
plot(bdry.r,bdry.z,color="gray",ls="--")
plot(path.r,path.z)
plot(o.path.r,o.path.z)
grid(true)
gca().set_aspect("equal")

It gets worse at higher energies

gcp = GCAlpha(3500,0.5,8.0,0.0)
...

There is also a mass dependence

If I were to hazard a guess then maybe its a relativistic effect but I cant be sure. It could also be an error in setting the initial conditions for the full orbit integrator

[henrikjaerleblad]

Doesn't get_full_orbit() already handle relativistic effects? Line 143 in fullorbit.jl.

I tried with a 35 keV, still get mismatch. Should not see relativistic effects for that low energies?
gcp = GCAlpha(35,0.5,8.0,0.0)

Also, I wanted to try it with the 96100 equilibrium but it failed and produced a lost orbit that I know should be co-passing. Also, the pitch seems to be handled differently. It is lost downwards, but should be lost upwards (if it weren't co-passing but lost).

With old package versions
using Equilibrium, GuidingCenterOrbits, Plots
filename_eqdsk = "g96100_0-53.0012.eqdsk"
M, wall = read_geqdsk(filename_eqdsk)
gcp = GCDeuteron(35.0,0.7,3.2,M.axis[2])
o = get_orbit(M, gcp, wall=wall)
Plots.plot(o.path.r, o.path.z, aspect_ratio=:equal, color=:blue)
Plots.plot!(wall.r, wall.z, color=:black)

With new package versions
using Equilibrium, GuidingCenterOrbits, Plots, EFIT
filename_eqdsk = "g96100_0-53.0012.eqdsk"
M = efit(readg(filename_eqdsk);clockwise_phi=true)
wall = Wall(readg(filename_eqdsk))
gcp = GCDeuteron(35.0,0.7,3.2,M.axis[2])
o = get_orbit(M,gcp, wall=wall)
Plots.plot(o.path.r,o.path.z, aspect_ratio=:equal)
Plots.plot!(wall.r,wall.z,color=:black)

[lstagner]

M = efit(readg(filename_eqdsk);clockwise_phi=true) Are you sure the positive phi direction is in the clockwise direction? Most European and American machines define the positive toroidal direction to be counter-clockwise.

[lstagner]

Yeah changing clockwise_phi=true to clockwise_phi=false seems to fix the issue.

[henrikjaerleblad]

In the email that I forwarded from Jacob, he wrote that JET has it in the clockwise direction.

"
Hi Henrik,

As we discussed with Luke and Mirko last time, I have now extracted the larger version of the FBM file for 94701V01.
It is located on Heimdall in the folder / home / jeriks / to / Henrik-J / Full-FBM.
I have checked and it contains the following data:

nsjdotb = 1 # so j and Bt are parallel to each other

nsnccw = -1 # so Bt is not counterclockwise, ie clockwise, seen from above

So this should mean that plasma current and magnetic field are parallel, and directed in a negative toroidal direction (ie clockwise, seen from above). This agrees with my picture of how it "usually" is on JET.

Regards // Jacob
"

I will double-check the nsnccw variable myself tomorrow. Probably misunderstood something.

[lstagner]

The clockwise_phi has nothing to do with the plasma current or b field. It is just the coordinate system used.

[lstagner]

I added a function called orbit_projection(M, gcp) that calculates the rz projection of the orbit using the contouring method. It agrees with the guiding center code. So what I am thinking is that we have entered a regime where the magnetic moment conservation is no longer conserved. If you plot the standard mu = vperp^2/(2B) using the full orbit path you see that it is not conserved. In order to make it conserved you have to use higher order corrections for the gyromotion. See Physics of Plasmas 10, 3240 (2003); https://doi.org/10.1063/1.1592155

[henrikjaerleblad]

Alright. So we would need to implement higher-order corrections such as those in R. G. Littlejohn, Plasma Phys.29, 111(1983) and equation (14) of https://doi.org/10.1063/1.1592155 ?

[henrikjaerleblad]

And (10)

[lstagner]

Maybe....it might be easier just to use the full orbit code and do all the callbacks on the guiding center position determined at each time step. But I'm not sure how much the jacobian calculation would like that.

[henrikjaerleblad]

There could be a low-aspect ratio check/flag?

[lstagner]

Its not just a low-aspect ratio problem though. The problem seems to be that the gyro-radius is too large breaking the low-order approximations used in the formulations of the guiding center equations of motion.

I think the ASCOT code does something to switch between guiding center and full orbit based on some criterion it may be a good idea to look into what they did. Actually, I think David Pfefferle worked on that.

[henrikjaerleblad]

Alright, I think there is an issue with the full-orbit computation. Because, since the formula for the gyroradius is

r_g = p_perp / |q|*B

I examined the gyroradii for an alpha particle and a deuteron at the same E,p,R,Z point (in this case, a starting point for a banana orbit in JET equilibrium 96100).

As you can see, the full orbit code predicts a larger gyroradius for the alpha particle compared to the deuteron. At first glance, this seems to make sense. However, for the E,p,R,Z (orbit starting point) value in the figure above, we have

p_perp_alpha = 1.94e-20
p_perp_deuteron = 1.38e-20
q_alpha = 2
q_deuteron = 1

so the gyroradius should in fact be smaller for the alpha particle than for a deuteron (assuming same E,p,R,Z,B point), since

p_perp_alpha / q_alpha = 0.97e-20
p_perp_deuteron / q_deuteron = 1.94e-20

This is not what the full-orbit code predicts. Hence, I would suspect that the problem lies with the full-orbit code?

[lstagner]

I think I figured out the issue. We are not scaling by the charge number in the boris push

 q_m_half_dt = (0.5*dt*e0/m) # <-------- where is the q factor? 

[lstagner]

Including the charge number seems to fix most of the issue. There still is the problem that the poloidal transit times do not match.

gcp = GCAlpha(3500,0.5,8.0,0.0)
o = get_orbit(S,gcp,maxiter=0)
path, stat = get_full_orbit(S,gcp,tmax=1.0*o.tau_p/cyclotron_period(S,gcp));
#oproj = orbit_projection(S,gcp)
plot(bdry.r,bdry.z,color="gray",ls="--")
plot(path.r,path.z)
plot(o.path.r,o.path.z)
#plot(oproj.r, oproj.z)
grid(true)
gca().set_aspect("equal")

[lstagner]

There seems to be a dependence on the initial gyro-angle

gcp = GCAlpha(3500,0.5,8.0,0.0)
o = get_orbit(S,gcp,maxiter=0)
path1, stat = get_full_orbit(S,gcp,tmax=1.0*o.tau_p/cyclotron_period(S,gcp),gamma=0.0);
path2, stat = get_full_orbit(S,gcp,tmax=1.0*o.tau_p/cyclotron_period(S,gcp),gamma=pi/2);
path3, stat = get_full_orbit(S,gcp,tmax=1.0*o.tau_p/cyclotron_period(S,gcp),gamma=pi);
path4, stat = get_full_orbit(S,gcp,tmax=1.0*o.tau_p/cyclotron_period(S,gcp),gamma=3pi/2);

#oproj = orbit_projection(S,gcp)
fig,ax = plt.subplots(ncols=4)
fig.set_size_inches(12,5)
ax[1].plot(path1.r,path1.z)
ax[2].plot(path2.r,path2.z)
ax[3].plot(path3.r,path3.z)
ax[4].plot(path4.r,path4.z)
gamma = ["0","pi/2","pi","3pi/2"]
for i=1:4
    ax[i].plot(bdry.r,bdry.z,color="gray",ls="--")
    ax[i].plot(o.path.r,o.path.z)
    ax[i].set_aspect("equal")
    ax[i].set_title("gamma = $(gamma[i])")
end
fig.tight_layout()

[henrikjaerleblad]

I think to get the end points to (approximately) match for the full-orbit and gc-orbit integrations, one might need to do something else than to just normalize tmax by the cyclotron period at the starting point since:

δ = 0.33
ϵ = 0.32
κ = 1.7
B0=2
R0=6.2
qstar = 1.57
alpha = -0.155
S = solovev(B0, R0, ϵ, δ, κ, alpha, qstar)

E = 3500.0
p = 0.5
R = 8.0
Z = 0.0
gcp = GCAlpha(E,p,R,Z)
o = get_orbit(S,gcp,maxiter=0)

gamma = lorentz_factor(gcp)
cp_max = 0.0
cp_min = Inf64
for i=1:length(o.path.r)
    ri = (o.path.r)[i]
    zi = (o.path.z)[i]
    Babs = norm(Bfield(S, ri, zi))

    cf = abs(gcp.q*GuidingCenterOrbits.e0)*Babs/(gamma*gcp.m)
    cp = 2*pi/cf
    if cp > cp_max
        cp_max = cp
    end
    if cp < cp_min
        cp_min = cp
    end
end

cp_max # 8.099362532408623e-8
cp_min # 6.074522663284002e-8
cp_default = cyclotron_period(S,gcp) # 8.099362530886745e-8

[lstagner]

The poloidal transit times are only off by about 2%. I wonder if its caused by the lack of second order corrections in the guiding center equations. The GCDE check may not trigger but the transit times may be off

[lstagner]

Could also be how the starting position is determined. We aren't using the second order corrections to the larmor radius when we take a step away from the gc location

[lstagner]

Well it wasn't that. It helped a bit but the poloidal transit time discrepancy remains.

[henrikjaerleblad]

How does the discrepancy change as function of particle energy? I was thinking if we need the lorentz factor factor when computing q_m_half_dt in the borispush()? Or if it is a relativistic issue of some sort

[lstagner]

The lorentz factor is like 1.001 which isn't really relativistic 1.05 and above is.

[lstagner]

Using your gcde_check the orbit fails

gcde_check(S,o,verbose=true)

Evaluating gcde criterion for orbit path position 1 of 135... 
Evaluating gcde criterion for orbit path position 2 of 135... 
Evaluating gcde criterion for orbit path position 3 of 135... 
--- Criterion violated! ---
Criterion > Threshold: 0.08300099551478635 > 0.073
- Gyroradius (r_g) at bad position: 0.14437908238201114 m
- |B| at bad position: 1.6225121821580957 T
- √(λmax) at bad position: 0.9327537211911904
√(λmax) * r_g / |B| > Threshold

[henrikjaerleblad]

Alright, so we are walking on thin ice anyhow. Does the dependence on the initial gryo-angle persist with the second-order correction?

[lstagner]

The second order correction to they gyro step? Yeah it persists. The discrepancy also exists for orbits that GC valid but its much less pronounced.