Sign problem in vmec_fieldlines

[rahulgaur104]

The quantities cvdrift and gbdrift output by the vmec_fieldlines routine must always satisfy: cvdrift >= gbdrift These quantities are defined in Dr. Edmund Highcock's thesis in section 3.7 here.

Right now, this is not the case. In fact, cvdrift <= gbdrift. Below is a simple script demonstrating the problem:

#!/usr/bin/env python3

import pdb
import numpy as np
from simsopt.mhd.vmec import Vmec
from simsopt.mhd.vmec_diagnostics import vmec_fieldlines
from simsopt.mhd.vmec_diagnostics import vmec_splines

v         = Vmec("wout_d23p4_tm.nc")
theta_fac = int(11)
ntheta    = int(701)
theta     = np.linspace(-1*theta_fac*np.pi, theta_fac*np.pi, ntheta)
s         = 0.5

vs           = vmec_splines(v)
f1           = vmec_fieldlines(vs, 0.5, 0.0, theta1d=theta)

a_N         = f1.L_reference
B_N         = f1.B_reference
bmag        = f1.bmag[0][0]
gbdrift     = f1.gbdrift[0][0]
cvdrift     = f1.cvdrift[0][0] 

# Note that (cvdrift-gbdrift)*bmag**2 is a constant
dPdrho      = -1*np.mean((cvdrift-gbdrift)*bmag**2)

# The printed number should be negative
print(dPdrho)

pdb.set_trace()

The equilibrium I have used can be downloaded from here.

My hypothesis is that there is a sign error in the vmec_fieldlines routine at line 1174

   cvdrift = gbdrift + 2 * B_reference * L_reference * L_reference * sqrt_s[:, None , None] * mu_0 * d_pressure_d_s[:, None, None] \
              * toroidal_flux_sign / (edge_toroidal_flux_over_2pi * modB * modB)                                                               

specifically, the second term. A more important concern is the effect of the source of this problem on the rest of the geometric quantities. I hope we haven't made the same error somewhere else.

[rahulgaur104]

Upon comparing the geometry coefficients from Simsopt with the ones from the gyrokinetics code GX, I get the following results
where every geometric quantity except cvdrift and gbdrift matches well. Upon flipping the sign of cvdrift and gbdrift, I get this You can download and look at these figures closesly to make sure they match. I am also confident that this sign problem is only limited to these two geometric quantities. You can also find these figures here and here.

[garethrc]

I think Rahul is right here about cvdrift being >= gbdrift, on account of the beta term being positive definite (using bhat cross grad psi dot grad alpha = B_N = bmag). I think he's also right that (cvdrift - gbdrift)* bmag^2 is constant.

Here's an example where gbdrift (orange) is always greater than cvdrift (green) plotted versus arc length. It's taken from the Ku & Boozer (2011) NFP4 Aspect Ratio 8 QH configuration, at s=0.5 and alpha=0.

So I suspect he's right about a sign issue.

[rahulgaur104]

Thanks for doing that check, Gareth!

A few weeks ago, Matt suspected that there is another coefficient with a sign issue in vmec_fieldlines.

To understand this better for a simpler case, I modified vmec_fieldlines for axisymmetric equilibria and compared the geometric coefficients with my VMEC2GK interface. The VMEC2GK interface uses a different mechanism to calculate the coefficients but the final answer should be the same for a given equilibrium. I take a simple DIIID-like equilibrium and compare the coefficients at rho = 0.6. Here are the results:

As you can see, the problematic quantity is cvdrift0. For this case, we need to flip the sign of cvdrift0 as well. I don't have any recommendations right now. I'll post again as I learn more about this problem.

[rahulgaur104]

Ideal-ballooning marginal stability test

The marginal stability contour against ideal-ballooning modes has been plotted by many physicists before us. Therefore, I'll use that to clear some of the confusion surrounding this issue. First, I compare the two coefficients needed to solve the ballooning equation from a VMEC equilibrium with their respective analytical expression from Edmund Highcock's thesis:

The important point is to get the signs and the general trend right. Using these coefficients, we solve the ideal ballooning equation (code given here) and compare the marginal stability curve from our solver with the ground truth below:

On the other hand, if we use an older version of vmec_geometry (previously vmec_fieldlines), we get a sign discrepancy: and a fairly different marginal stability curve (absence of a second-stability boundary + fairly different first-stability boundary).

Please ignore the S in the lower right.

This analysis is done at theta0 = 0 for circular (up-down symmetric) equilibria. In other words, this test only checks the sign of cvdrift and gbdrift. We still need tests/further analysis for the two important sign-depends quantities: gds21 and cvdrift0.

Unrelated bug in vmec_geometry for 2D equilibria

This is an unrelated issue that one must be cognizant of. In a VMEC file, we parametrize the boundary shape using the following Fourier modes to get the geometric coefficients plotted in the figure at the beginning of this message

RBC(0,1) = 40.0 ZBS(0,0) = 0.0 RBC(0,1) = 2.0 ZBS(0,1) = 2.0

However, if you flip the sign of ZBS(0, 1) (which corresponds to the exact same boundary shape) you get the following coefficients: This is more than just a sign change. We have changed the sign as well as redefined theta, theta = -pi now being on the outboard side. The difference is even more pronounced for strongly-shaped equilibria. I have a new axisymmetric function that I have tested extensively, solves this problem, and is slightly faster than vmec_geometry for 2D equilibria. I'll create a pull request for that soon.

This message is already too long. I'll create another thread discussing the invariance of the geometric coefficients w.r.t the transformation psi -> -psi.