Poincare Boundary Condition

[YigitElma]

Resolves #204
Resolves #773

This PR focuses on solving equilibrium by giving a Poincare section at zeta 0 surfaces instead of giving the last closed flux surface, LCFS.

This method gives exactly same results for high NFP >12, but different results for low NFP cases. Note: as far as my observation, the difference is caused by the NFP, but it might caused by something else that I am missing.

Tasks,

• Add tests (make them run in reasonable time)
• make sure API works with poincare (try to resolve inaccuracy issue) (kind of done)
• make sure we don't use eq.surface anywhere expecting explicitly a toroidal surface
• move utility functions for making poincare XS from an equilibrium to proper part (utils.py getters.py) Moved it to Equilibrium class i.e. eq.set_poincare_equilibrium() gives you separate equilibrium
• Get poincare to work with BoundaryRSelfConsistency (and Z), and remove poincareR and poincareZ objectives
• Add get_fixed_xsection_constraints() utility
• Add new fixed and self consistency objectives for Poincare BC

Previous tasks that I am not sure about,

• remove bdry_mode from Equilibrium and instead add new DOFs of the equilibrium Rp_lmn Zp_lmn Lp_lmn and xsection.
• When creating an Equilibrium with a Poincare section, check that initial axis is consistent with the equilibrium (it should just be a circular axis at the poincare section axis position) (I think it is like that but I don't know what Dario meant exactly)
• decide if LambdaGauge is needed or not with Poincare equilibria

[github-actions]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     -0.10 +/- 3.78     | -5.44e-04 +/- 1.97e-02 |  5.21e-01 +/- 1.8e-02  |  5.21e-01 +/- 9.0e-03  |
-test_equilibrium_init_medres            |    +12.48 +/- 1.37     | +5.09e-01 +/- 5.60e-02 |  4.59e+00 +/- 4.5e-02  |  4.08e+00 +/- 3.3e-02  |
-test_equilibrium_init_highres           |    +32.53 +/- 0.71     | +1.74e+00 +/- 3.79e-02 |  7.08e+00 +/- 2.2e-02  |  5.34e+00 +/- 3.1e-02  |
 test_objective_compile_dshape_current   |     +2.86 +/- 1.13     | +1.11e-01 +/- 4.40e-02 |  4.00e+00 +/- 3.9e-02  |  3.89e+00 +/- 2.0e-02  |
-test_objective_compute_dshape_current   |    +34.22 +/- 1.79     | +1.75e-03 +/- 9.12e-05 |  6.85e-03 +/- 5.4e-05  |  5.10e-03 +/- 7.4e-05  |
 test_objective_jac_dshape_current       |     +4.22 +/- 4.54     | +1.81e-03 +/- 1.94e-03 |  4.46e-02 +/- 1.3e-03  |  4.28e-02 +/- 1.4e-03  |
-test_perturb_2                          |    +24.39 +/- 1.85     | +4.75e+00 +/- 3.61e-01 |  2.42e+01 +/- 2.9e-01  |  1.95e+01 +/- 2.1e-01  |
 test_proximal_freeb_jac                 |     +1.20 +/- 1.15     | +8.85e-02 +/- 8.48e-02 |  7.46e+00 +/- 5.8e-02  |  7.37e+00 +/- 6.2e-02  |
-test_solve_fixed_iter                   |    +10.18 +/- 1.68     | +3.21e+00 +/- 5.31e-01 |  3.48e+01 +/- 4.0e-01  |  3.16e+01 +/- 3.5e-01  |
-test_LinearConstraintProjection_build   |    +27.03 +/- 2.41     | +2.74e+00 +/- 2.44e-01 |  1.29e+01 +/- 1.1e-01  |  1.02e+01 +/- 2.2e-01  |
 test_build_transform_fft_midres         |     -0.46 +/- 6.57     | -2.83e-03 +/- 4.02e-02 |  6.09e-01 +/- 3.2e-02  |  6.12e-01 +/- 2.5e-02  |
 test_build_transform_fft_highres        |     +0.40 +/- 3.23     | +3.83e-03 +/- 3.11e-02 |  9.66e-01 +/- 2.2e-02  |  9.62e-01 +/- 2.2e-02  |
 test_equilibrium_init_lowres            |     +5.61 +/- 3.13     | +2.17e-01 +/- 1.21e-01 |  4.08e+00 +/- 7.8e-02  |  3.86e+00 +/- 9.2e-02  |
 test_objective_compile_atf              |     -2.14 +/- 3.80     | -1.73e-01 +/- 3.07e-01 |  7.92e+00 +/- 2.8e-01  |  8.09e+00 +/- 1.2e-01  |
-test_objective_compute_atf              |    +11.41 +/- 3.19     | +1.80e-03 +/- 5.04e-04 |  1.76e-02 +/- 3.3e-04  |  1.58e-02 +/- 3.8e-04  |
 test_objective_jac_atf                  |     +0.71 +/- 2.81     | +1.37e-02 +/- 5.47e-02 |  1.96e+00 +/- 4.3e-02  |  1.94e+00 +/- 3.4e-02  |
-test_perturb_1                          |    +32.21 +/- 2.11     | +4.73e+00 +/- 3.10e-01 |  1.94e+01 +/- 2.9e-01  |  1.47e+01 +/- 1.2e-01  |
-test_proximal_jac_atf                   |     +7.34 +/- 1.17     | +5.99e-01 +/- 9.52e-02 |  8.77e+00 +/- 5.9e-02  |  8.17e+00 +/- 7.5e-02  |
-test_proximal_freeb_compute             |    +30.94 +/- 1.78     | +6.21e-02 +/- 3.58e-03 |  2.63e-01 +/- 2.7e-03  |  2.01e-01 +/- 2.4e-03  |
-test_solve_fixed_iter_compiled          |    +30.78 +/- 3.08     | +6.25e+00 +/- 6.26e-01 |  2.66e+01 +/- 6.1e-01  |  2.03e+01 +/- 1.3e-01  |

[f0uriest]

I think the better API would be to introduce new degrees of freedom Rp, Zp etc like we have for the LCFS and axis, and then constrain those to be self consistent.

This would eliminate the boundary_mode thing, and allow more flexibility for setting up constraints.

Each equilibrium, in addition to having a .surface and .axis attribute would have a .poincare or .section property with the ZernikeRZToroidalSurface object. This would then have an assigned zeta value so we wouldn't need to pass around zeta to different functions.

[YigitElma]

I think the better API would be to introduce new degrees of freedom Rp, Zp etc like we have for the LCFS and axis, and then constrain those to be self consistent.

This would eliminate the boundary_mode thing, and allow more flexibility for setting up constraints.

Each equilibrium, in addition to having a .surface and .axis attribute would have a .poincare or .section property with the ZernikeRZToroidalSurface object. This would then have an assigned zeta value so we wouldn't need to pass around zeta to different functions.

@f0uriest I made some changes to the API. First of all, I created a new class called PoincareSurface (which is an inherited class of ZernikeRzToroidalSection) to be able to store lambda basis and coefficients in the section. That class needs more cleaning but my main aim is to decide which type of objectives to use by checking the eq.surface type. This can be either FourierRZToroidalSurface or PoincareSurface. I changed the previous eq.bdry_mode stuff to that. Also, I am using the zeta of the surface now, so no need to pass zeta value to optimizer and objectives anymore.

[ddudt]

Let's get rid of support for the Poincare BC from a text input file in this PR

[YigitElma]

@dpanici @f0uriest I made the maybe_add_self_consistency function cleaner. For your point on why do we need to check extra cases, the problem starts at this line

DESC/desc/optimize/_constraint_wrappers.py

Line 580 in 9ef0b38

self._linear_constraints = maybe_add_self_consistency(

This is the set of constraints used for ProximalProjections equilibrium update. We add self consistency constraints at the beginning. Then, we pass these constraints to perturb and solve which again call maybe_add_self_consistency function. This is the main problem why we need the checks. My previous implementation was really bad. Now it is a bit cleaner.

The idea here is to add self consistency only if there is a corresponding fixed constraint, AND there is no opposing self consistency constraint. For example,
{FixBoundaryR} -> add BoundaryRSelfConsistency and FixSectionR
{FixBoundaryR, BoundaryRSelfConsistency, FixSectionR} -> don't add anything

Applying the similar thing to LCFS and axis stuff requires a bit more complex checks. For example, during optimization if none of the constraints is FixBoundary, FixAxis or FixSection, then it messes up the optimization by adding a FixBoundary constraint.

[YigitElma]

By the way, I am aware of the full red benchmarks, and don't want to affect 99% of the users which won't use this. I will take a closer look at it once the main code works as expected.

[YigitElma]

I think these are the all related cases.

[dpanici]

We have added lots of stuff since the first time we noticed things were different here, can we check now again to see if the equilibrium solve without any hacky removal of constraints results in the same answer?

[YigitElma]

We have added lots of stuff since the first time we noticed things were different here, can we check now again to see if the equilibrium solve without any hacky removal of constraints results in the same answer?

I think it is better than it used to be but the test_HELITRON_vac_results still fails. If I remember correctly, that was the most problematic one for some reason. I tried to change tolerances but putting lower tolerances cause the flux surfaces to be not nested after the first step of continuation and perturb. And it quits. I can achieve the intended result like this tho,

args = ["../tests/inputs/HELIOTRON_vacuum", "-vvv", "-o", "heli_vac.h5"]
main(args)

vmec_path = "..//tests//inputs//wout_HELIOTRON_vacuum.nc"
eq = Equilibrium.load("heli_vac.h5")[-1]
VMECIO.plot_vmec_comparison(eq, vmec_path);
plt.savefig("vmec_comp.png", dpi=1000)

eq.solve(maxiter=100, ftol=2e-5, gtol=0, verbose=3);

rho_err, theta_err = area_difference_vmec(eq, vmec_path)
np.testing.assert_allclose(rho_err.mean(), 0, atol=1e-2)
np.testing.assert_allclose(theta_err.mean(), 0, atol=2e-2)
curr = eq.get_profile("current")
np.testing.assert_allclose(curr(np.linspace(0, 1, 20)), 0, atol=1e-8)

Is there a way to give different stopping tolerance for different steps of continuation?