- it is important to make sure the profiles being used (pressure and either iota or net toroidal current) are correct and physical
- specifically, at the axis if using net toroidal current, need zero net current there, and also need it to be ~ rho^2 (pressure also must be ~rho^2 but it is less important than current as iota is affected by current)
-the way I prefer to do this (with succcess) in my own comparisons with DIII/EFIT is to take what the input profile data are from EFIT and fit them as a PowerSeriesProfile.from_values(x=rho, y=data, order=L, ..., sym=True) which will enforce the correct ~rho^2 dependence,
- then take the current profile object and set the first parameter to zero (like current_profile.params[0] = 0.0 should work, can confirm it by checking current_profile(0.0))
- Another thing is when doing eq.solve (and I recommend eq.solve over solve_continuation_automatic for tokamak, that continuation method is meant for stellarators), you need to not use the default grids and tolerances.
- Something like this is what I usually use, maybe with higher maxiter if needed. Basically bc it is a tokamak, you know the solution should exist and get to low error (bc 2D, so has nested surfaces for certain), so you want low ftol, xtol, gtol, and also want to solve on a QuadratureGrid which has more points than a ConcentricGrid (important for making sure you compute iota accurately everywhere including near axis, and for getting the force error as low as it can be)
eq.solve(
ftol=1e-12,
maxiter=150,
gtol=0,
verbose=3,
xtol=0,
objective=ObjectiveFunction(
ForceBalance(eq, grid=QuadratureGrid(L=eq.L_grid, M=eq.M, N=0))
),
)another note: when you use VMECIO.load(vmecfile), it by default loads spline fits of the pressure and iota, so if you ever want the current profile, you need to pass profile="current" , though I recommend getting the current profile directly from the g file (or the VMEC input file if you have that instead of the g file), instead of from VMEC wout, as VMEC does not save the current profile directly, but rather indirectly as buco so it makes it slightly different than the current profile input sometimes
-the way I prefer to do this (with succcess) in my own comparisons with DIII/EFIT is to take what the input profile data are from EFIT and fit them as a PowerSeriesProfile.from_values(x=rho, y=data, order=L, ..., sym=True) which will enforce the correct ~rho^2 dependence,
another note: when you use VMECIO.load(vmecfile), it by default loads spline fits of the pressure and iota, so if you ever want the current profile, you need to pass profile="current" , though I recommend getting the current profile directly from the g file (or the VMEC input file if you have that instead of the g file), instead of from VMEC wout, as VMEC does not save the current profile directly, but rather indirectly as buco so it makes it slightly different than the current profile input sometimes