`plot_surfaces` vartheta contours can be discontinuous at boundary, likely due to rootfind fails

[dpanici]

For the following input file (attached), the resulting asymmetric equilibrium, when plot_surfaces is called, has the vartheta contours make a sharp change at the boundary, when if you just do plot_section(eq,"theta_PEST"), you see none of this behavior, indicating something is happening with the rootfind for theta_PEST contours during the plot_surfaces call.

The equilibrium is perfectly nested so it is not an equilibrium issue. Maybe we can just do a grid calculation of theta_PEST and make contours based off of that instead of doing the rootfind for the plot_surfaces?

input.asym_tokamak.txt

[unalmis]

Check if removing period argument from fixes it

1. https://github.com/PlasmaControl/DESC/blob/89cedb4fe06cdecaed0f4ef66fe89e16246b897f/desc/plotting.py#L1710 2.https://github.com/PlasmaControl/DESC/blob/89cedb4fe06cdecaed0f4ef66fe89e16246b897f/desc/plotting.py#L1688

[dpanici]

check lambda derivative wrt theta, if it is >1 or not

[dpanici]

d lambda d theta > 1 at some points in the plasma, so yea this is messing up the rootfinding convergence.

[dpanici]

Try re-fitting boundary with better poloidal angle and see if this issue becomes resolved

[dpanici]

Spectral condensation of the boundary helps with this issue immensely and lowers the value of $d\lambda/d\theta$ we see. Closing this as resolved but this is motivation to make a spectral condensation routine for surface fitting available in the code

[dpanici]

I see this even in less strongly shaped cases, this can also happen due to poor initial guesses for the vartheta contour. In master currently the general rootfinding map_coordinates is used with inbasis=["rho","theta_PEST","phi"] which does NOT use the _map_PEST... method inside map_coordinates (as it is allowing for omega nonzero and so performing the rootfind in all the coordinates).

Changing on master for now to inbasis=["rho","theta_PEST","zeta"] helps it be more robust #1951

[unalmis]

You can use the closed-form expression that is guaranteed to converge the unique solution I give in the comments of #1154 . You can set that as the initial guess if you prefer too. The 3D coordinate mapping convergence would be improved by resolving #1207 .