Poincare Boundary Condition

[YigitElma]

trafficstars

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.

• Adds Rp_lmn, Zp_lmn and Lp_lmn to Equilibrium class These are not optimizable_parameter by default, instead when you call build method of one of the FixSection... constraints, they become optimizable and Rb_lmn and Zb_lmn become un-optimizable. This behavior can be changed by remove_optimizables argument.

• Adds make_optimizable and make_nonoptimizable for this conversion. These are called automatically in FixSection.. and FixBoundary.. constraints.

• Alters ProximalProjection to solve inner equilibrium with Poincare boundary condition. This is done automatically is a user gives FixSection.. constraint (they shouldn't mix any FixBoundary..)

• get_surface_at method of Equilibrium can now return any arbitrary zeta cross-section including the lambda values. If the asked cross-section is not 0 or pi/NFP, then the return ZernikeRZTorodialSection object doesn't have symmetry.

• Adds tutorial on how to use Poincare BC in different contexts. Also, adds an example to use in optimization without running it for long.

• Removes bdry_mode attribute

• Adds update_all_equilibria_backwards_compatibility.py script to devtools that load and then saves (almost) every equilibria in desc/examples and tests/inputs to run _set_up method.

[github-actions[bot]]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     +2.07 +/- 7.50     | +1.11e-02 +/- 4.00e-02 |  5.44e-01 +/- 3.9e-02  |  5.33e-01 +/- 1.1e-02  |
-test_equilibrium_init_medres            |    +19.10 +/- 1.40     | +7.29e-01 +/- 5.36e-02 |  4.54e+00 +/- 3.8e-02  |  3.82e+00 +/- 3.8e-02  |
-test_equilibrium_init_highres           |    +51.57 +/- 2.05     | +2.45e+00 +/- 9.72e-02 |  7.20e+00 +/- 8.6e-02  |  4.75e+00 +/- 4.6e-02  |
 test_objective_compile_dshape_current   |     +1.56 +/- 6.27     | +4.55e-02 +/- 1.83e-01 |  2.96e+00 +/- 1.4e-01  |  2.91e+00 +/- 1.2e-01  |
 test_objective_compute_dshape_current   |     +2.70 +/- 1.61     | +9.50e-05 +/- 5.67e-05 |  3.61e-03 +/- 5.3e-05  |  3.51e-03 +/- 1.9e-05  |
 test_objective_jac_dshape_current       |     -0.04 +/- 6.45     | -1.70e-05 +/- 2.48e-03 |  3.85e-02 +/- 1.7e-03  |  3.85e-02 +/- 1.8e-03  |
 test_perturb_2                          |     +0.88 +/- 2.48     | +1.33e-01 +/- 3.75e-01 |  1.52e+01 +/- 3.1e-01  |  1.51e+01 +/- 2.1e-01  |
 test_proximal_jac_atf_with_eq_update    |     +0.30 +/- 1.19     | +4.98e-02 +/- 2.00e-01 |  1.68e+01 +/- 1.6e-01  |  1.68e+01 +/- 1.2e-01  |
 test_proximal_freeb_jac                 |     +0.08 +/- 0.92     | +4.53e-03 +/- 5.30e-02 |  5.74e+00 +/- 3.9e-02  |  5.74e+00 +/- 3.6e-02  |
 test_solve_fixed_iter_compiled          |     -1.34 +/- 3.46     | -2.39e-01 +/- 6.15e-01 |  1.75e+01 +/- 3.4e-01  |  1.78e+01 +/- 5.1e-01  |
 test_LinearConstraintProjection_build   |     +0.17 +/- 4.32     | +1.34e-02 +/- 3.47e-01 |  8.05e+00 +/- 2.6e-01  |  8.04e+00 +/- 2.3e-01  |
 test_objective_compute_ripple_spline    |     -0.51 +/- 2.88     | -1.62e-03 +/- 9.22e-03 |  3.19e-01 +/- 3.6e-03  |  3.20e-01 +/- 8.5e-03  |
 test_objective_grad_ripple_spline       |     -0.15 +/- 2.40     | -2.18e-03 +/- 3.47e-02 |  1.44e+00 +/- 2.6e-02  |  1.45e+00 +/- 2.3e-02  |
 test_build_transform_fft_midres         |     -0.20 +/- 5.51     | -1.10e-03 +/- 3.07e-02 |  5.57e-01 +/- 2.9e-02  |  5.58e-01 +/- 9.8e-03  |
 test_build_transform_fft_highres        |     -0.18 +/- 3.51     | -1.52e-03 +/- 2.96e-02 |  8.42e-01 +/- 2.7e-02  |  8.44e-01 +/- 1.3e-02  |
-test_equilibrium_init_lowres            |     +6.03 +/- 1.31     | +2.20e-01 +/- 4.77e-02 |  3.86e+00 +/- 3.9e-02  |  3.65e+00 +/- 2.7e-02  |
 test_objective_compile_atf              |     +0.16 +/- 2.55     | +9.89e-03 +/- 1.58e-01 |  6.20e+00 +/- 9.6e-02  |  6.19e+00 +/- 1.3e-01  |
 test_objective_compute_atf              |     +1.14 +/- 2.80     | +1.04e-04 +/- 2.56e-04 |  9.23e-03 +/- 1.7e-04  |  9.13e-03 +/- 1.9e-04  |
 test_objective_jac_atf                  |     -0.07 +/- 2.55     | -1.36e-03 +/- 4.65e-02 |  1.83e+00 +/- 3.8e-02  |  1.83e+00 +/- 2.7e-02  |
 test_perturb_1                          |     +0.28 +/- 2.47     | +3.49e-02 +/- 3.04e-01 |  1.24e+01 +/- 2.5e-01  |  1.23e+01 +/- 1.8e-01  |
 test_proximal_jac_atf                   |     -0.14 +/- 0.85     | -1.10e-02 +/- 6.81e-02 |  7.98e+00 +/- 6.0e-02  |  7.99e+00 +/- 3.3e-02  |
 test_proximal_freeb_compute             |     -1.27 +/- 2.92     | -2.10e-03 +/- 4.82e-03 |  1.63e-01 +/- 3.4e-03  |  1.65e-01 +/- 3.5e-03  |
 test_solve_fixed_iter                   |     -0.36 +/- 3.06     | -1.03e-01 +/- 8.76e-01 |  2.85e+01 +/- 6.1e-01  |  2.86e+01 +/- 6.3e-01  |
 test_objective_compute_ripple           |     -0.63 +/- 1.39     | -1.78e-02 +/- 3.93e-02 |  2.82e+00 +/- 3.2e-02  |  2.83e+00 +/- 2.3e-02  |
 test_objective_grad_ripple              |     -0.10 +/- 1.02     | -7.12e-03 +/- 7.54e-02 |  7.35e+00 +/- 5.7e-02  |  7.35e+00 +/- 4.9e-02  |

[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

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 88.11189% with 51 lines in your changes missing coverage. Please review.

Project coverage is 95.57%. Comparing base (3b5e928) to head (6a626f4).

Files with missing lines	Patch %	Lines
desc/objectives/linear_objectives.py	84.21%	27 Missing :warning:
desc/optimize/_constraint_wrappers.py	41.93%	18 Missing :warning:
desc/geometry/surface.py	95.74%	2 Missing :warning:
desc/optimizable.py	86.66%	2 Missing :warning:
desc/continuation.py	50.00%	1 Missing :warning:
desc/equilibrium/equilibrium.py	98.68%	1 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #782      +/-   ##
==========================================
- Coverage   95.71%   95.57%   -0.15%     
==========================================
  Files         101      101              
  Lines       26370    26720     +350     
==========================================
+ Hits        25240    25537     +297     
- Misses       1130     1183      +53     

Files with missing lines	Coverage Δ
desc/compat.py	86.95% <100.00%> (+0.24%)	:arrow_up:
desc/equilibrium/initial_guess.py	96.59% <100.00%> (-0.61%)	:arrow_down:
desc/equilibrium/utils.py	100.00% <100.00%> (ø)
desc/input_reader.py	91.96% <ø> (-0.05%)	:arrow_down:
desc/objectives/__init__.py	100.00% <ø> (ø)
desc/objectives/getters.py	96.29% <100.00%> (+0.74%)	:arrow_up:
desc/objectives/utils.py	100.00% <100.00%> (ø)
desc/optimize/optimizer.py	96.90% <100.00%> (+0.07%)	:arrow_up:
desc/perturbations.py	89.83% <100.00%> (+0.78%)	:arrow_up:
desc/vmec.py	91.84% <100.00%> (+0.01%)	:arrow_up:
... and 6 more

... and 2 files with indirect coverage changes

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[YigitElma]

Only failing regression tests are Near axis ones. @dpanici will look into those. Probably related to #839

[dpanici]

Only failing regression tests are Near axis ones. @dpanici will look into those. Probably related to #839

Some plot tests seem to fail too?

[YigitElma]

Only failing regression tests are Near axis ones. @dpanici will look into those. Probably related to #839

Some plot tests seem to fail too?

Yeah, I mean among regression tests😁

[dpanici]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     -2.28 +/- 1.42     | -2.12e-04 +/- 1.32e-04 |  9.08e-03 +/- 1.2e-04  |  9.30e-03 +/- 5.1e-05  |
 test_build_transform_fft_midres         |     -0.25 +/- 1.26     | -2.16e-04 +/- 1.10e-03 |  8.74e-02 +/- 9.1e-04  |  8.76e-02 +/- 6.1e-04  |
 test_build_transform_fft_highres        |     +1.58 +/- 1.19     | +7.22e-03 +/- 5.41e-03 |  4.63e-01 +/- 4.6e-03  |  4.56e-01 +/- 2.9e-03  |
 test_equilibrium_init_lowres            |     +7.39 +/- 2.51     | +1.97e-02 +/- 6.70e-03 |  2.87e-01 +/- 6.6e-03  |  2.67e-01 +/- 1.2e-03  |
-test_equilibrium_init_medres            |    +53.88 +/- 1.66     | +3.26e-01 +/- 1.01e-02 |  9.31e-01 +/- 8.7e-03  |  6.05e-01 +/- 5.1e-03  |
-test_equilibrium_init_highres           |    +68.07 +/- 1.70     | +1.48e+00 +/- 3.71e-02 |  3.67e+00 +/- 3.4e-02  |  2.18e+00 +/- 1.5e-02  |
 test_objective_compile_dshape_current   |     +5.67 +/- 8.54     | +2.20e-01 +/- 3.31e-01 |  4.10e+00 +/- 2.3e-01  |  3.88e+00 +/- 2.3e-01  |
 test_objective_compile_atf              |    +10.54 +/- 5.19     | +8.04e-01 +/- 3.96e-01 |  8.44e+00 +/- 2.8e-01  |  7.63e+00 +/- 2.8e-01  |
 test_objective_compute_dshape_current   |     +8.62 +/- 6.90     | +4.57e-04 +/- 3.66e-04 |  5.76e-03 +/- 3.6e-04  |  5.30e-03 +/- 6.4e-05  |
 test_objective_compute_atf              |     +6.38 +/- 3.93     | +7.20e-04 +/- 4.44e-04 |  1.20e-02 +/- 4.1e-04  |  1.13e-02 +/- 1.7e-04  |
-test_objective_jac_dshape_current       |   +110.41 +/- 10.48    | +5.20e-02 +/- 4.93e-03 |  9.90e-02 +/- 4.4e-03  |  4.71e-02 +/- 2.2e-03  |
-test_objective_jac_atf                  |    +16.31 +/- 3.40     | +3.48e-01 +/- 7.24e-02 |  2.48e+00 +/- 4.3e-02  |  2.13e+00 +/- 5.8e-02  |
-test_perturb_1                          |    +22.90 +/- 6.79     | +3.62e+00 +/- 1.07e+00 |  1.94e+01 +/- 9.7e-01  |  1.58e+01 +/- 4.5e-01  |
-test_perturb_2                          |    +20.06 +/- 1.61     | +4.19e+00 +/- 3.37e-01 |  2.51e+01 +/- 9.8e-02  |  2.09e+01 +/- 3.2e-01  |
-test_proximal_jac_atf                   |     +7.35 +/- 1.07     | +5.28e-01 +/- 7.67e-02 |  7.72e+00 +/- 6.5e-02  |  7.19e+00 +/- 4.1e-02  |
 test_proximal_freeb_compute             |     +1.89 +/- 1.07     | +2.34e-03 +/- 1.33e-03 |  1.26e-01 +/- 6.2e-04  |  1.24e-01 +/- 1.2e-03  |
 test_proximal_freeb_jac                 |     -0.16 +/- 1.24     | -1.12e-02 +/- 8.82e-02 |  7.08e+00 +/- 6.0e-02  |  7.09e+00 +/- 6.5e-02  |

These are some pretty large slowdowns, mainly for the init due to poincare surface being innitialized (OK), but for the rest likely just due to the addition of Zp Rp Lp attributes which make the jacobian larger for the constraint even if it then gets factored out in the null space. Might be worth thinking of ways to improve this before we merge

[f0uriest]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     -2.28 +/- 1.42     | -2.12e-04 +/- 1.32e-04 |  9.08e-03 +/- 1.2e-04  |  9.30e-03 +/- 5.1e-05  |
 test_build_transform_fft_midres         |     -0.25 +/- 1.26     | -2.16e-04 +/- 1.10e-03 |  8.74e-02 +/- 9.1e-04  |  8.76e-02 +/- 6.1e-04  |
 test_build_transform_fft_highres        |     +1.58 +/- 1.19     | +7.22e-03 +/- 5.41e-03 |  4.63e-01 +/- 4.6e-03  |  4.56e-01 +/- 2.9e-03  |
 test_equilibrium_init_lowres            |     +7.39 +/- 2.51     | +1.97e-02 +/- 6.70e-03 |  2.87e-01 +/- 6.6e-03  |  2.67e-01 +/- 1.2e-03  |
-test_equilibrium_init_medres            |    +53.88 +/- 1.66     | +3.26e-01 +/- 1.01e-02 |  9.31e-01 +/- 8.7e-03  |  6.05e-01 +/- 5.1e-03  |
-test_equilibrium_init_highres           |    +68.07 +/- 1.70     | +1.48e+00 +/- 3.71e-02 |  3.67e+00 +/- 3.4e-02  |  2.18e+00 +/- 1.5e-02  |
 test_objective_compile_dshape_current   |     +5.67 +/- 8.54     | +2.20e-01 +/- 3.31e-01 |  4.10e+00 +/- 2.3e-01  |  3.88e+00 +/- 2.3e-01  |
 test_objective_compile_atf              |    +10.54 +/- 5.19     | +8.04e-01 +/- 3.96e-01 |  8.44e+00 +/- 2.8e-01  |  7.63e+00 +/- 2.8e-01  |
 test_objective_compute_dshape_current   |     +8.62 +/- 6.90     | +4.57e-04 +/- 3.66e-04 |  5.76e-03 +/- 3.6e-04  |  5.30e-03 +/- 6.4e-05  |
 test_objective_compute_atf              |     +6.38 +/- 3.93     | +7.20e-04 +/- 4.44e-04 |  1.20e-02 +/- 4.1e-04  |  1.13e-02 +/- 1.7e-04  |
-test_objective_jac_dshape_current       |   +110.41 +/- 10.48    | +5.20e-02 +/- 4.93e-03 |  9.90e-02 +/- 4.4e-03  |  4.71e-02 +/- 2.2e-03  |
-test_objective_jac_atf                  |    +16.31 +/- 3.40     | +3.48e-01 +/- 7.24e-02 |  2.48e+00 +/- 4.3e-02  |  2.13e+00 +/- 5.8e-02  |
-test_perturb_1                          |    +22.90 +/- 6.79     | +3.62e+00 +/- 1.07e+00 |  1.94e+01 +/- 9.7e-01  |  1.58e+01 +/- 4.5e-01  |
-test_perturb_2                          |    +20.06 +/- 1.61     | +4.19e+00 +/- 3.37e-01 |  2.51e+01 +/- 9.8e-02  |  2.09e+01 +/- 3.2e-01  |
-test_proximal_jac_atf                   |     +7.35 +/- 1.07     | +5.28e-01 +/- 7.67e-02 |  7.72e+00 +/- 6.5e-02  |  7.19e+00 +/- 4.1e-02  |
 test_proximal_freeb_compute             |     +1.89 +/- 1.07     | +2.34e-03 +/- 1.33e-03 |  1.26e-01 +/- 6.2e-04  |  1.24e-01 +/- 1.2e-03  |
 test_proximal_freeb_jac                 |     -0.16 +/- 1.24     | -1.12e-02 +/- 8.82e-02 |  7.08e+00 +/- 6.0e-02  |  7.09e+00 +/- 6.5e-02  |

These are some pretty large slowdowns, mainly for the init due to poincare surface being innitialized (OK), but for the rest likely just due to the addition of Zp Rp Lp attributes which make the jacobian larger for the constraint even if it then gets factored out in the null space. Might be worth thinking of ways to improve this before we merge

I think it's because the benchmarks aren't including the self-consistency constraints that remove the extra DoFs, will fix

[YigitElma]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     -4.52 +/- 3.72     | -8.46e-02 +/- 6.96e-02 |  1.79e+00 +/- 1.8e-02  |  1.87e+00 +/- 6.7e-02  |
 test_build_transform_fft_midres         |     -5.43 +/- 3.27     | -1.13e-01 +/- 6.80e-02 |  1.97e+00 +/- 3.9e-02  |  2.08e+00 +/- 5.6e-02  |
 test_build_transform_fft_highres        |     -0.87 +/- 2.38     | -2.07e-02 +/- 5.66e-02 |  2.35e+00 +/- 3.8e-02  |  2.37e+00 +/- 4.2e-02  |
 test_equilibrium_init_lowres            |     +1.18 +/- 3.98     | +1.15e-01 +/- 3.88e-01 |  9.87e+00 +/- 2.9e-02  |  9.75e+00 +/- 3.9e-01  |
 test_equilibrium_init_medres            |     +3.07 +/- 3.56     | +3.25e-01 +/- 3.77e-01 |  1.09e+01 +/- 3.5e-02  |  1.06e+01 +/- 3.8e-01  |
 test_equilibrium_init_highres           |     +3.79 +/- 8.28     | +5.04e-01 +/- 1.10e+00 |  1.38e+01 +/- 1.1e-01  |  1.33e+01 +/- 1.1e+00  |
 test_objective_compile_dshape_current   |    -13.66 +/- 7.65     | -5.87e-01 +/- 3.28e-01 |  3.71e+00 +/- 5.5e-02  |  4.29e+00 +/- 3.2e-01  |
 test_objective_compile_atf              |     -7.34 +/- 3.81     | -5.93e-01 +/- 3.08e-01 |  7.49e+00 +/- 1.1e-01  |  8.08e+00 +/- 2.9e-01  |
+test_objective_compute_dshape_current   |     -8.81 +/- 2.68     | -5.40e-04 +/- 1.64e-04 |  5.59e-03 +/- 4.8e-05  |  6.13e-03 +/- 1.6e-04  |
 test_objective_compute_atf              |     -2.81 +/- 1.92     | -6.16e-04 +/- 4.22e-04 |  2.13e-02 +/- 2.4e-04  |  2.19e-02 +/- 3.4e-04  |
 test_objective_jac_dshape_current       |     -3.28 +/- 6.43     | -1.38e-03 +/- 2.71e-03 |  4.07e-02 +/- 1.7e-03  |  4.21e-02 +/- 2.1e-03  |
 test_objective_jac_atf                  |     -0.56 +/- 4.60     | -1.05e-02 +/- 8.66e-02 |  1.87e+00 +/- 5.0e-02  |  1.88e+00 +/- 7.1e-02  |
 test_perturb_1                          |    +12.61 +/- 5.12     | +2.07e+00 +/- 8.38e-01 |  1.84e+01 +/- 3.0e-01  |  1.64e+01 +/- 7.8e-01  |
 test_perturb_2                          |     +8.79 +/- 4.30     | +1.91e+00 +/- 9.36e-01 |  2.37e+01 +/- 2.6e-01  |  2.18e+01 +/- 9.0e-01  |
-test_proximal_jac_atf                   |     +5.41 +/- 1.45     | +3.95e-01 +/- 1.05e-01 |  7.69e+00 +/- 9.2e-02  |  7.30e+00 +/- 5.2e-02  |
 test_proximal_freeb_compute             |     +1.51 +/- 0.58     | +1.89e-03 +/- 7.20e-04 |  1.27e-01 +/- 3.7e-04  |  1.25e-01 +/- 6.2e-04  |
 test_proximal_freeb_jac                 |     +0.44 +/- 3.16     | +3.23e-02 +/- 2.31e-01 |  7.36e+00 +/- 2.2e-01  |  7.33e+00 +/- 6.9e-02  |

It is much better now with the new tests!!!

[YigitElma]

#931 changed the error in factorize_linear_constraints to a warning but still pytest fails some test due to the warning, i.e.

=========================== short test summary info ============================
FAILED tests/test_equilibrium.py::test_poincare_as_bc - UserWarning: 
Not equal to tolerance rtol=3e-14, atol=3e-14
Incompatible constraints detected, cannot satisfy constraint <desc.objectives.linear_objectives.BoundaryRSelfConsistency object at 0x7f5982199220>.
Mismatched elements: 2 / 38 (5.26%)
Max absolute difference: 6.2394534e-14
Max relative difference: inf
 x: array([ 3.063705e-16, -1.948655e-16, -1.146194e-17, -4.111364e-18,
        4.678379e-17,  2.974076e-17, -1.001699e-15,  2.379133e-16,
        1.181742e-16,  3.608323e-17,  7.501427e-16, -1.843074e-16,...
 y: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
 This may indicate incompatible constraints, or be due to floating point error.

Either we should add manual flags like @pytest.mark.filterwarnings("ignore::UserWarning") to relevant tests or we can ignore all UserWarnings in setup.cfg.

I found out that,

with pytest.raises(UserWarning):
...

is not a proper solution due to numerical tolerances of different computers. For example, Github machines throw warning but my laptop doesn't. So, running test with with clause cause test to fail on my laptop, because the warning is not thrown.

[dpanici]

• check master and maybe raise tols for the factorize assertion
• if needed, throw an ignore simple filter on the test that throws this
• play with ordering of params

[YigitElma]

I tried to change the order of the arguments,

arg_order = (
            "R_lmn",
            "Z_lmn",
            "L_lmn",
            "p_l",
            "i_l",
            "c_l",
            "Psi",
            "Te_l",
            "ne_l",
            "Ti_l",
            "Zeff_l",
            "a_lmn",
            "Ra_n",
            "Za_n",
            "Rb_lmn",
            "Zb_lmn",
            "Rp_lmn",
            "Zp_lmn",
            "Lp_lmn",
            "I",
            "G",
            "Phi_mn",
        )

I put *p_lmn parameters before axis parameters, after axis and to the very end (after Phi_mn). Only putting it to the end changed the results a little bit (in a bad way). I will try some other orders too, but I think we should find some other solution to this, since it affects the accuracy noticably.

[YigitElma]

The numerical tolerance issue is probably caused by this. Here is the dimensions of matrix A for both master and Poincare branch,

constraints = get_fixed_boundary_constraints(eq)
constraints = maybe_add_self_consistency(eq, constraints)
constraints = ObjectiveFunction(constraints)

objective = ObjectiveFunction(ForceBalance(eq))
objective = LinearConstraintProjection(objective, constraints)
objective.build(verbose=3)
print(objective._A.shape)

master:     (650, 6356)
poincare:   (897, 6603)

So, it is roughly 1.43 times bigger than master in terms of number of matrix elements.

[dpanici]

change tol to be eps * size(A) or something

[YigitElma]

This branch is literally healing itself lol It still has some of the tests fail but the error tolerances are significantly reduced.

[YigitElma]

Failing Heliotron vacuum test results, I increased the maxiter but it seems to have no effect.

[YigitElma]

This is ready for review. Only the Heliotron vacuum test fails. I guess it is some numerical tolerance issue that is somehow more prominent for that equilibrium. For reference, that test again solves the fixed boundary problem, and the newly added variables Rp, Zp and Lp are supposed to be taken out by the self-consistency constraints. My hunch is that for that specific problem, factorize_linear_constraints cannot do a good job of finding a particular solution or null space. I checked the errors in self-consistency constraints and they were around 1e-15 and 1e-16 for some of them. I would assume that should be fine, but somehow the final result is not the same as master even if I increase the maxiter. I will be expecting your feedback for the PR and try to find some bugs or improve the accuracy of factorize_linear_constrains.

[dpanici]

write what the actual test is checking (area diff), check if iota goes thru rational, also can make coiils and check against field tracing, one must be correct and one must be incorrect

[YigitElma]

I have finally got to see the actual coverage of this PR 🥲🥲🥲 How I manage to do this a bit controversial but a good example of the problem with null space accuracy!

[YigitElma]

Do we ever want to fix the cross-section at zeta=Pi/NFP ? I don't test them and those were mainly for fixing 2 cross-sections together. I am pretty sure if I delete them the coverage will pass too. @f0uriest @dpanici

[dpanici]

dont do pi/nfp otion for zeta, make issue to make compat function to change zeta=0 definition

[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 https://github.com/PlasmaControl/DESC/blob/9ef0b386c2bbffa2fc52ee9dfa31da9f7faa7ee0/desc/optimize/_constraint_wrappers.py#L580

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]

OUTDATED

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?

[YigitElma]

I changed the hacky way with another approach where I declare some parameters non-optimizable. To do so, I added a new function inside desc.optimizable called make_nonoptimizable and a setter for cls.optimizable_params. If this works, it should solve the speed issue.

[YigitElma]

@f0uriest @dpanici, one final offer about this. I am pretty convinced that a user cannot solve or optimize an equilibrium (with cold-start) with just Poincare (minor changes in the section cause big changes in the overall equilibrium). Given this, I don't think we need to have Rp_lmn etc as an optimizable variable. What I suggest is to enforce the constraint by computing the Rp_lmn values during the build method of the objective. This way you can prevent the cross-section from changing during the optimization. In this PR, I had added some logic to different parts of the code like perturb etc, I can remove all of them if you are ok to this? I assume we do something similar in NAE.

If I do this, there will only be SectionRSelfConsistency no FixSectionR (actually I can switch the names) such that when you calculate the constraint matrix, there won't be 00000100000 rows but just 0101010111000000 rows.