Near axis constraints

This tutorial shows how to find equilibrium solutions in DESC which are constrained to have the same axis and near-axis behavior as a NAE solution from the pyQSC or pyQIC codes

Creating a DESC Equilibrium from a pyQSC or pyQIC Near-Axis Equilibrium

Note that you must have pyQSC and/or pyQIC installed in order to make use of the Equilibrium.from_near_axis method. You can do so with pip install qsc.

[2]:
# must have installed pyQsc with `pip install qsc` in order to use this!
from qsc import Qsc
import numpy as np
import matplotlib.pyplot as plt
from desc.equilibrium import Equilibrium
from desc.objectives import get_fixed_boundary_constraints
from desc.plotting import (
    plot_comparison,
    plot_fsa,
    plot_section,
    plot_surfaces,
    plot_qs_error,
)
DESC version 0.10.2+351.g7a44e147.dirty,using JAX backend, jax version=0.4.13, jaxlib version=0.4.13, dtype=float64
Using device: CPU, with 14.42 GB available memory

DESC is able to create an equilibrium based off of a pyQsc NAE equilibrium object. First, we’ll make the NAE equilibrium using pyQsc

[3]:
qsc_eq = Qsc.from_paper("precise QA")

Then, to make the DESC equilibrium, the Equilibrium class has a method Equilibrium.from_near_axis. This method creates a DESC Equilibrium based off of the pyQsc equilibrium. It requires as input the desired DESC Fourier-Zernike resolutions, as well as the radius at which you want to evaluate the qsc equilibrium at to make the DESC equilibrium’s boundary. The equilibrium’s initial R_lmn, Z_lmn Fourier-Zernike coefficients are fit to the R,Z evaluated from the Qsc equilibrium, and the initial L_lmn are 0 (because the Qsc equilibrium uses Boozer angles, so there is no poloidal stream function)

[4]:
# ntheta is the number of points in theta to sample
# when fitting the DESC Fourier-Zernike basis to
# the NAE surfaces.
# The more shaped the NAE surfaces, the higher this should be,
# 75 is a good minimum value.
ntheta = 75
# r is the finite radius to evaluate the NAE surfaces at, this will
# roughly set the aspect ratio of the resulting solution.
# However if this is too large, the surfaces from the NAE may become
# self-intersecting! The r at which this occurs varies for different NAE
# solutions, r=0.35 was found to be a reasonable radius for precise QA
r = 0.35
desc_eq = Equilibrium.from_near_axis(
    qsc_eq,  # the Qsc equilibrium object
    r=r,  # the finite radius (m) at which to evaluate the Qsc surface to use as the DESC boundary
    L=8,  # DESC radial resolution
    M=8,  # DESC poloidal resolution
    N=8,  # DESC toroidal resolution
    ntheta=ntheta,
)
eq_fit = desc_eq.copy()  # copy so we can see the original Qsc surfaces later

Now we solve the equilibrium as normal in DESC

[5]:
# get the fixed-boundary constraints, which include also fixing the pressure and fixing the current profile (iota=False flag means fix current)
constraints = get_fixed_boundary_constraints(eq=desc_eq)
for c in constraints:
    print(c)

# solve the equilibrium
desc_eq.solve(
    verbose=3,
    ftol=1e-2,
    objective="force",
    maxiter=100,
    xtol=1e-6,
    constraints=constraints,
)

# Save equilibrium as .h5 file
desc_eq.save("DESC_from_NAE_precise_QA_output.h5")
<desc.objectives.linear_objectives.FixBoundaryR object at 0x7f7f4d0d4520>
<desc.objectives.linear_objectives.FixBoundaryZ object at 0x7f7f4d0d4730>
<desc.objectives.linear_objectives.FixPsi object at 0x7f7f4d0d4580>
<desc.objectives.linear_objectives.FixPressure object at 0x7f7f4d0d4640>
<desc.objectives.linear_objectives.FixCurrent object at 0x7f7f3e6b2040>
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 2.87 sec
Timer: Objective build = 8.61 sec
Timer: Linear constraint projection build = 10.6 sec
Compiling objective function and derivatives: ['force']
Timer: Objective compilation time = 3.76 sec
Timer: Jacobian compilation time = 11.4 sec
Timer: Total compilation time = 15.2 sec
Number of parameters: 856
Number of objectives: 5346
Starting optimization
Using method: lsq-exact
   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality
       0              1          2.079e+01                                    1.241e+04
       1              2          2.355e+00      1.844e+01      3.138e-01      1.297e+03
       2              3          5.777e-02      2.297e+00      1.514e-01      3.875e+02
       3              4          9.046e-03      4.872e-02      1.865e-01      5.217e+01
       4              6          1.395e-03      7.652e-03      4.553e-02      2.322e+01
       5              8          1.059e-03      3.357e-04      2.327e-02      4.387e+00
       6             10          1.033e-03      2.653e-05      1.057e-02      5.663e-01
       7             11          1.011e-03      2.195e-05      1.761e-02      2.928e+00
       8             13          9.892e-04      2.155e-05      8.471e-03      5.916e-01
       9             14          9.722e-04      1.707e-05      1.537e-02      2.905e+00
      10             16          9.537e-04      1.846e-05      8.131e-03      6.297e-01
      11             17          9.391e-04      1.457e-05      1.428e-02      2.182e+00
      12             19          9.267e-04      1.239e-05      7.196e-03      3.560e-01
      13             20          9.151e-04      1.161e-05      1.330e-02      1.063e+00
      14             21          9.134e-04      1.742e-06      2.693e-02      2.125e+00
      15             22          8.891e-04      2.431e-05      6.961e-03      1.243e-01
      16             23          8.837e-04      5.384e-06      1.327e-02      4.947e-01
Optimization terminated successfully.
`ftol` condition satisfied.
         Current function value: 8.837e-04
         Total delta_x: 3.221e-01
         Iterations: 16
         Function evaluations: 23
         Jacobian evaluations: 17
Timer: Solution time = 1.97 min
Timer: Avg time per step = 6.95 sec
Start of solver
Total (sum of squares):  2.079e+01,
Maximum absolute Force error:  1.044e+07 (N)
Minimum absolute Force error:  7.686e-01 (N)
Average absolute Force error:  7.492e+04 (N)
Maximum absolute Force error:  6.465e+00 (normalized)
Minimum absolute Force error:  4.758e-07 (normalized)
Average absolute Force error:  4.638e-02 (normalized)
R boundary error:  0.000e+00 (m)
Z boundary error:  0.000e+00 (m)
Fixed-Psi error:  0.000e+00 (Wb)
Fixed-pressure profile error:  0.000e+00 (Pa)
Fixed-current profile error:  0.000e+00 (A)
End of solver
Total (sum of squares):  8.837e-04,
Maximum absolute Force error:  3.671e+04 (N)
Minimum absolute Force error:  2.581e-01 (N)
Average absolute Force error:  8.472e+02 (N)
Maximum absolute Force error:  2.273e-02 (normalized)
Minimum absolute Force error:  1.598e-07 (normalized)
Average absolute Force error:  5.245e-04 (normalized)
R boundary error:  2.003e-17 (m)
Z boundary error:  8.080e-17 (m)
Fixed-Psi error:  0.000e+00 (Wb)
Fixed-pressure profile error:  0.000e+00 (Pa)
Fixed-current profile error:  0.000e+00 (A)

Now we have a DESC equilibrium solved with the boundary from Qsc. It has zero toroidal current as its profile constraint along with zero pressure since the original Qsc equilibrium had no pressure or current.

However, if we plot the surfaces, we see that while the boundary matches, the interior deviates slightly from the NAE, especially near the core. This means we may have lost some of the optimized properties from the QSC equilibrium.

[6]:
plot_comparison(
    eqs=[desc_eq, eq_fit],
    labels=["DESC-solved equilibrium", "NAE surfaces"],
    figsize=(12, 12),
    theta=0,
    color=["k", "r"],
    ls=["-", ":"],
    lw=[1, 2],
);
../../_images/notebooks_tutorials_nae_constraint_14_0.png

Instead of evaluating the NAE at a finite radius, fixing that boundary, and hoping the axis behavior stays the same, we can instead directly fix the axis and the \(O(\rho)\) asymptotic behavior.

Solving Equilibria with Fixed Axis and Fixed NAE \(O(\rho)\) Behavior in DESC

[7]:
# utility functions for getting the NAE constraints
from desc.objectives import get_equilibrium_objective, get_NAE_constraints

eq_NAE = eq_fit.copy()
# this has all the constraints we need
constraints = get_NAE_constraints(eq_NAE, qsc_eq, order=1)

eq_NAE.solve(
    verbose=3,
    ftol=1e-2,
    objective="force",
    maxiter=50,
    xtol=1e-6,
    constraints=constraints,
);
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 386 ms
Timer: Objective build = 1.09 sec
Timer: Linear constraint projection build = 14.7 sec
Compiling objective function and derivatives: ['force']
Timer: Objective compilation time = 3.69 sec
Timer: Jacobian compilation time = 11.0 sec
Timer: Total compilation time = 14.7 sec
Number of parameters: 1094
Number of objectives: 5346
Starting optimization
Using method: lsq-exact
   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality
       0              1          2.084e+01                                    8.055e+03
       1              3          1.269e+01      8.153e+00      2.157e-01      1.983e+04
       2              4          6.743e-01      1.202e+01      1.194e-01      1.470e+03
       3              5          1.259e-01      5.484e-01      2.337e-01      1.815e+02
       4              7          9.821e-03      1.161e-01      7.469e-02      4.620e+01
       5              9          1.694e-03      8.128e-03      4.481e-02      1.334e+01
       6             11          1.346e-04      1.559e-03      3.373e-02      3.294e+00
       7             13          2.017e-05      1.144e-04      2.348e-02      1.125e+00
       8             15          5.537e-06      1.463e-05      1.260e-02      2.340e-01
       9             17          3.431e-06      2.106e-06      6.067e-03      6.759e-02
      10             18          3.107e-06      3.244e-07      1.022e-02      1.998e-01
      11             19          2.496e-06      6.109e-07      8.939e-03      1.795e-01
      12             20          2.026e-06      4.700e-07      7.844e-03      1.338e-01
      13             21          1.717e-06      3.093e-07      7.130e-03      1.135e-01
      14             22          1.496e-06      2.208e-07      6.533e-03      1.353e-01
      15             23          1.327e-06      1.685e-07      6.036e-03      1.732e-01
      16             24          1.193e-06      1.348e-07      5.662e-03      1.980e-01
      17             25          1.082e-06      1.103e-07      5.424e-03      2.061e-01
      18             26          9.917e-07      9.056e-08      5.313e-03      2.057e-01
      19             27          9.157e-07      7.598e-08      5.247e-03      1.975e-01
      20             28          8.325e-07      8.321e-08      5.078e-03      1.678e-01
      21             29          7.252e-07      1.073e-07      4.754e-03      1.457e-01
      22             30          6.246e-07      1.006e-07      4.475e-03      1.301e-01
      23             31          5.468e-07      7.783e-08      4.324e-03      1.276e-01
      24             32          4.895e-07      5.724e-08      4.227e-03      1.396e-01
      25             33          4.481e-07      4.138e-08      4.135e-03      1.629e-01
      26             34          4.187e-07      2.946e-08      4.031e-03      1.912e-01
      27             35          2.990e-07      1.197e-07      1.016e-03      1.077e-02
      28             36          2.914e-07      7.646e-09      2.001e-03      5.302e-02
      29             37          2.775e-07      1.389e-08      1.968e-03      5.630e-02
      30             38          2.645e-07      1.299e-08      1.940e-03      5.930e-02
      31             39          2.522e-07      1.228e-08      1.913e-03      6.198e-02
      32             40          2.406e-07      1.166e-08      1.887e-03      6.399e-02
      33             41          2.295e-07      1.106e-08      1.862e-03      6.550e-02
      34             42          2.190e-07      1.049e-08      1.839e-03      6.648e-02
      35             43          2.091e-07      9.945e-09      1.817e-03      6.694e-02
      36             44          1.996e-07      9.423e-09      1.797e-03      6.691e-02
      37             45          1.907e-07      8.923e-09      1.779e-03      6.648e-02
      38             46          1.823e-07      8.444e-09      1.763e-03      6.577e-02
      39             47          1.743e-07      7.988e-09      1.749e-03      6.495e-02
      40             48          1.667e-07      7.556e-09      1.736e-03      6.413e-02
      41             49          1.596e-07      7.150e-09      1.725e-03      6.343e-02
      42             50          1.528e-07      6.771e-09      1.715e-03      6.302e-02
      43             51          1.464e-07      6.420e-09      1.706e-03      6.382e-02
      44             52          1.403e-07      6.094e-09      1.698e-03      6.472e-02
      45             53          1.345e-07      5.794e-09      1.691e-03      6.570e-02
      46             54          1.290e-07      5.518e-09      1.685e-03      6.673e-02
      47             55          1.237e-07      5.267e-09      1.680e-03      6.778e-02
      48             56          1.187e-07      5.039e-09      1.676e-03      6.878e-02
      49             57          1.138e-07      4.833e-09      1.672e-03      6.967e-02
      50             58          1.092e-07      4.649e-09      1.670e-03      7.041e-02
Warning: Maximum number of iterations has been exceeded.
         Current function value: 1.092e-07
         Total delta_x: 3.322e-01
         Iterations: 50
         Function evaluations: 58
         Jacobian evaluations: 51
Timer: Solution time = 5.32 min
Timer: Avg time per step = 6.26 sec
Start of solver
Total (sum of squares):  2.079e+01,
Maximum absolute Force error:  1.044e+07 (N)
Minimum absolute Force error:  7.686e-01 (N)
Average absolute Force error:  7.492e+04 (N)
Maximum absolute Force error:  6.465e+00 (normalized)
Minimum absolute Force error:  4.758e-07 (normalized)
Average absolute Force error:  4.638e-02 (normalized)
R axis error:  0.000e+00 (m)
Z axis error:  0.000e+00 (m)
Fixed-Psi error:  0.000e+00 (Wb)
Fixed-pressure profile error:  0.000e+00 (Pa)
Fixed-current profile error:  0.000e+00 (A)
Fixed-R sum modes error:  6.029e-06 (m)
Fixed-R sum modes error:  3.393e-06 (m)
Fixed-R sum modes error:  7.190e-06 (m)
Fixed-R sum modes error:  1.582e-05 (m)
Fixed-R sum modes error:  1.702e-05 (m)
Fixed-R sum modes error:  1.628e-05 (m)
Fixed-R sum modes error:  1.248e-05 (m)
Fixed-R sum modes error:  3.286e-06 (m)
Fixed-R sum modes error:  9.086e-06 (m)
Fixed-R sum modes error:  3.605e-06 (m)
Fixed-R sum modes error:  5.182e-06 (m)
Fixed-R sum modes error:  1.026e-05 (m)
Fixed-R sum modes error:  1.490e-05 (m)
Fixed-R sum modes error:  1.185e-05 (m)
Fixed-R sum modes error:  7.719e-07 (m)
Fixed-R sum modes error:  1.415e-05 (m)
Fixed-R sum modes error:  1.513e-05 (m)
Fixed-Z sum modes error:  1.058e-05 (m)
Fixed-Z sum modes error:  1.938e-05 (m)
Fixed-Z sum modes error:  6.087e-06 (m)
Fixed-Z sum modes error:  9.129e-06 (m)
Fixed-Z sum modes error:  1.831e-05 (m)
Fixed-Z sum modes error:  1.872e-05 (m)
Fixed-Z sum modes error:  1.293e-05 (m)
Fixed-Z sum modes error:  5.624e-06 (m)
Fixed-Z sum modes error:  8.495e-07 (m)
Fixed-Z sum modes error:  3.280e-07 (m)
Fixed-Z sum modes error:  3.646e-06 (m)
Fixed-Z sum modes error:  7.968e-06 (m)
Fixed-Z sum modes error:  1.521e-05 (m)
Fixed-Z sum modes error:  2.188e-05 (m)
Fixed-Z sum modes error:  2.505e-05 (m)
Fixed-Z sum modes error:  2.206e-05 (m)
Fixed-Z sum modes error:  9.957e-06 (m)
End of solver
Total (sum of squares):  1.092e-07,
Maximum absolute Force error:  3.574e+02 (N)
Minimum absolute Force error:  3.939e-03 (N)
Average absolute Force error:  1.052e+01 (N)
Maximum absolute Force error:  2.212e-04 (normalized)
Minimum absolute Force error:  2.439e-09 (normalized)
Average absolute Force error:  6.513e-06 (normalized)
R axis error:  3.469e-18 (m)
Z axis error:  5.294e-23 (m)
Fixed-Psi error:  0.000e+00 (Wb)
Fixed-pressure profile error:  0.000e+00 (Pa)
Fixed-current profile error:  0.000e+00 (A)
Fixed-R sum modes error:  1.665e-16 (m)
Fixed-R sum modes error:  8.327e-17 (m)
Fixed-R sum modes error:  6.106e-16 (m)
Fixed-R sum modes error:  3.348e-16 (m)
Fixed-R sum modes error:  5.421e-17 (m)
Fixed-R sum modes error:  6.235e-16 (m)
Fixed-R sum modes error:  1.077e-16 (m)
Fixed-R sum modes error:  2.573e-16 (m)
Fixed-R sum modes error:  4.539e-16 (m)
Fixed-R sum modes error:  3.447e-17 (m)
Fixed-R sum modes error:  2.148e-16 (m)
Fixed-R sum modes error:  6.201e-16 (m)
Fixed-R sum modes error:  1.113e-16 (m)
Fixed-R sum modes error:  3.686e-17 (m)
Fixed-R sum modes error:  1.110e-16 (m)
Fixed-R sum modes error:  2.776e-17 (m)
Fixed-R sum modes error:  3.886e-16 (m)
Fixed-Z sum modes error:  2.220e-16 (m)
Fixed-Z sum modes error:  2.220e-16 (m)
Fixed-Z sum modes error:  1.110e-16 (m)
Fixed-Z sum modes error:  1.197e-16 (m)
Fixed-Z sum modes error:  1.214e-17 (m)
Fixed-Z sum modes error:  3.746e-17 (m)
Fixed-Z sum modes error:  5.954e-17 (m)
Fixed-Z sum modes error:  6.917e-17 (m)
Fixed-Z sum modes error:  2.271e-16 (m)
Fixed-Z sum modes error:  4.564e-17 (m)
Fixed-Z sum modes error:  3.716e-18 (m)
Fixed-Z sum modes error:  2.215e-16 (m)
Fixed-Z sum modes error:  2.982e-19 (m)
Fixed-Z sum modes error:  1.429e-16 (m)
Fixed-Z sum modes error:  1.752e-16 (m)
Fixed-Z sum modes error:  1.874e-16 (m)
Fixed-Z sum modes error:  8.327e-17 (m)

Comparing near axis behavior

Again we can plot the surfaces, and we see that in this case, the near axis behavior is preserved, while the outer surfaces differ. This is expected, as the NAE from QSC is only valid asymptotically near the axis. By constraining this behavior we are able to keep all the desireable properties of the NAE where they are valid without overly constraining the problem.

[8]:
fig, ax = plot_comparison(
    eqs=[eq_fit, eq_NAE],
    labels=[f"NAE Surfaces r={r}", f"DESC NAE constrained"],
    color=["k", "g"],
    ls=["-", "--"],
    figsize=(12, 12),
    theta=0,
    lw=[1, 2],
)
fig, ax = plot_comparison(
    eqs=[desc_eq, eq_NAE],
    labels=[f"DESC Fixed Surface", f"DESC NAE constrained"],
    color=["r", "g"],
    ls=[":", "--"],
    figsize=(12, 12),
    theta=0,
    lw=[1, 2],
);
../../_images/notebooks_tutorials_nae_constraint_20_0.png
../../_images/notebooks_tutorials_nae_constraint_20_1.png

As a final comparison, we can plot the QS error for the fixed boundary and fixed NAE solves. We see that by fixing the NAE behavior, we are able to preserve the QS from the original QSC equilibrium.

[12]:
fix, ax = plt.subplots()
plot_qs_error(
    desc_eq,
    fC=False,
    fT=False,
    log=True,
    ax=ax,
    color=["r"],
    legend=False,
    labels=["Fixed boundary"],
    figsize=(8, 8),
)
plot_qs_error(
    eq_NAE,
    fC=False,
    fT=False,
    log=True,
    ax=ax,
    color=["g"],
    legend=False,
    labels=["Fixed NAE"],
)
ax.legend()
ax.set_ylabel(
    "$f_B$", fontsize=16
)  # f_B = sqrt(sum(symmetry-breaking modes**2)) / sqrt(sum(all modes**2))
[12]:
Text(0, 0.5, '$f_B$')
../../_images/notebooks_tutorials_nae_constraint_22_1.png

The correct behaviour near the axis may be also checked by assessing directly other features of the equilibrium; most notably, the behaviour of |B| on axis, the rotational transform on axis, as well as the form of \(\lambda\). Starting from the rotational transform, we may compare the values obtained in the equilibrium to those of the near-axis expansion.

[23]:
import matplotlib.pyplot as plt
from desc.plotting import plot_fsa

rho = np.linspace(0, 1e-1)
fig, ax, iota_nae = plot_fsa(
    eq_NAE, "iota", rho=rho, return_data=True, lw=6, linecolor="g"
)
fig, ax, iota_surf = plot_fsa(
    desc_eq, "iota", rho=rho, ax=ax, return_data=True, linecolor="r", lw=3
)
plt.plot(rho, np.ones(np.size(rho)) * qsc_eq.iota, linestyle="--", lw=3, c="b")
plt.legend(["Fixed NAE", "Fixed boundary", "NAE"])
print(
    "Relative error in the rotational transform (fixed NAE): ",
    (iota_nae["iota"][0] - qsc_eq.iota) / qsc_eq.iota,
)
print(
    "Relative error in the rotational transform (fixed surface): ",
    (iota_surf["iota"][0] - qsc_eq.iota) / qsc_eq.iota,
)
# print(iota[1].yaxis)
Relative error in the rotational transform (fixed NAE):  -2.6032466642813977e-05
Relative error in the rotational transform (fixed surface):  0.11412206021843516
../../_images/notebooks_tutorials_nae_constraint_24_1.png

This clearly shows that the constraint is working as it is intended. A similar expected behaviour can be found for other features of the equilibrium such as the magnetic field on axis.

[21]:
from desc.compute import data_index
from desc.grid import LinearGrid

grid = LinearGrid(theta=np.array(0.0), N=100, NFP=desc_eq.NFP, rho=np.array(0.0))
# Evaluate B modes near the axis
data_surf = desc_eq.compute(["|B|"], grid=grid)
data_nae = eq_NAE.compute(["|B|"], grid=grid)
phi = grid.nodes[:, 2]

plt.plot(phi, data_surf["|B|"], lw=3, c="r")
plt.plot(phi, data_nae["|B|"], lw=6, c="g")
plt.plot(phi, np.ones(np.size(phi)) * qsc_eq.B0, linestyle="--", lw=3, c="b")
plt.xlabel("$\phi$")
plt.ylabel("$|B|$")
plt.legend(["Fixed surface", "Fixed NAE", "NAE"])
plt.show()
../../_images/notebooks_tutorials_nae_constraint_26_0.png

Finally, we may check \(\lambda\) on axis and compare it to \(\nu\) from within the near-axis expansion.

\(\lambda\) is the poloidal stream function used in DESC which transforms the computational poloidal angle \(\theta\) into the straight-field-line (SFL) poloidal angle \(\theta_{PEST} = \theta + \lambda\) (So B is “straight” when plotted in this poloidal angle and the cylindrical toroidal angle \(\phi\)).

Now, the NAE describes flux surfaces using the Boozer poloidal angle \(\theta_B\), and we are constraining our solution to be asymptotically equivalent to the NAE, so our computational poloidal angle on-axis should approach \(\theta = \theta_B\).

We also know that \(\theta_{PEST} = \theta_B - \iota \nu\) where \(\nu = \phi_B - \phi\), with \(\phi_B\) being the Boozer toroidal angle.

Taken together with \(\theta_{PEST} = \theta + \lambda\), we find that on-axis we expect our DESC NAE-constrained solution to have \(\lambda = -\iota \nu\), where \(\nu\) is known already from the NAE.

[22]:
# Evaluate lambda near the axis
data_surf = desc_eq.compute("lambda", grid=grid)
data_nae = eq_NAE.compute("lambda", grid=grid)
lam_surf = data_surf["lambda"]
lam_nae = data_nae["lambda"]

print(
    "Deviation of theta from Boozer angle (fixed surface)",
    np.mean(np.abs(lam_surf + qsc_eq.iota * qsc_eq.nu_spline(phi))),
)
print(
    "Deviation of theta from Boozer angle (fixed nae)",
    np.mean(np.abs(lam_nae + qsc_eq.iota * qsc_eq.nu_spline(phi))),
)
plt.plot(phi, lam_surf, lw=3, c="r")
plt.plot(phi, lam_nae, lw=6, c="g")
plt.plot(phi, -qsc_eq.iota * qsc_eq.nu_spline(phi), linestyle="--", lw=3, c="b")
plt.xlabel("$\phi$")
plt.ylabel("$\lambda$")
plt.legend(["Fixed surface", "Fixed NAE", "NAE"])
plt.show()
Deviation of theta from Boozer angle (fixed surface) 0.0993022659320299
Deviation of theta from Boozer angle (fixed nae) 5.016654348199815e-05
../../_images/notebooks_tutorials_nae_constraint_28_1.png

The above comparison shows that the DESC equilibrium constructed with the near-axis constraint is consistent with the near-axis behaviour.