Fix reactor_QA script to account for ScaledProfiles being returned now by rescale

[dpanici]

In our reactor QA script, we use FixCurrent assuming the [0,1] indices are the constant and rho^1 modes of the profile. However after the changes in #1871, rescaling an eq now returns ScaledProfile object which means now that [0,1] indices are the scale and the constant mode of the power series.

this resulted in letting the linear term be free, and a nonzero derivative of current at the axis which is unphysical (this also leads to unphysical iota profiles due to the way the axis limit of iota works compared to how iota is computed away from rho=0). This fixes that unintentional bug and update the reactor_QA output accordingly.

The reactor_QA example .h5 file had this issue since #1907 updated it (my bad...).

[github-actions]

Memory benchmark result

|               Test Name                |      %Δ      |    Master (MB)     |      PR (MB)       |    Δ (MB)    |    Time PR (s)     |  Time Master (s)   |
| -------------------------------------- | ------------ | ------------------ | ------------------ | ------------ | ------------------ | ------------------ |
  test_objective_jac_w7x                 |   -2.75 %    |     3.931e+03      |     3.823e+03      |   -107.95    |       22.19        |       20.18        |
  test_proximal_jac_w7x_with_eq_update   |    0.44 %    |     6.564e+03      |     6.593e+03      |    29.11     |       89.28        |       90.20        |
  test_proximal_freeb_jac                |   -0.43 %    |     1.333e+04      |     1.328e+04      |    -57.09    |       52.90        |       54.40        |
  test_proximal_freeb_jac_blocked        |   -0.59 %    |     7.673e+03      |     7.628e+03      |    -45.00    |       49.08        |       49.59        |
  test_proximal_freeb_jac_batched        |    0.42 %    |     7.567e+03      |     7.599e+03      |    31.97     |       47.88        |       48.77        |
  test_proximal_jac_ripple               |   -0.95 %    |     3.451e+03      |     3.418e+03      |    -32.71    |       35.56        |       35.55        |
  test_proximal_jac_ripple_bounce1d      |    1.35 %    |     3.702e+03      |     3.752e+03      |    49.98     |       44.79        |       46.19        |
  test_eq_solve                          |    0.60 %    |     1.892e+03      |     1.903e+03      |    11.30     |       35.26        |       36.76        |

For the memory plots, go to the summary of Memory Benchmarks workflow and download the artifact.

[dpanici]

Test was failing bc the iota profile is not well-represented by an even power series for reactorQA (see below, eq1 is the real eq, eq2 is after loading from the asciiIO, where the iota was obtained with get_profile as a power series with default L=eq.L=8 and even symmetry), it was surprising it passed before. I loosened the tols even more now since it is more an exercist of the geometry and the logic for kinettics

[codecov]

Codecov Report

❌ Patch coverage is 45.45455% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.31%. Comparing base (7841b52) to head (88b17d2).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
desc/compat.py	25.00%	3 Missing ⚠️
desc/objectives/_bootstrap.py	57.14%	3 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #2268   +/-   ##
=======================================
  Coverage   94.31%   94.31%           
=======================================
  Files         101      101           
  Lines       28870    28880   +10     
=======================================
+ Hits        27228    27239   +11     
+ Misses       1642     1641    -1     

Files with missing lines	Coverage Δ
desc/io/ascii_io.py	91.86% <ø> (ø)
desc/objectives/_power_balance.py	91.75% <ø> (ø)
desc/compat.py	87.80% <25.00%> (-0.81%)	⬇️
desc/objectives/_bootstrap.py	95.18% <57.14%> (-3.51%)	⬇️

... and 4 files with indirect coverage changes

🚀 New features to boost your workflow:

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

[YigitElma]

I don't know if it is really necessary, but I think the conversion is subtle enough to have more detailed explanation in docstring

[ddudt]

While reviewing this PR, I noticed that there were a few references to Equilibrium.ion_density missing after it was introduced in #2090. I added that fix to this PR and clarified some of the documentation. (Sorry for the scope creep.)

FYI I always do the following after calling rescale to avoid this issue (assuming the current profile was originally a power series):

eq = rescale(eq, L=("R0", R0), B=("B0", B0), scale_pressure=False, verbose=1)
eq.current = PowerSeriesProfile(params=eq.c_l[1:] * eq.c_l[0])

[YigitElma]

Looks good to me. I would prefer an even more elaborate explanation but non-blocking.

[ddudt]

Since I made some of the changes this approval doesn't really count.

[lkadz]

Before approving, I just had these minor changes / questions

[lkadz]

Suggested change

	errorif(
	eq.ion_density is None and eq.atomic_number is None,
	RuntimeError,
	"Bootstrap current calculation requires an ion density profile.",
	)
	errorif(
	eq.atomic_number is None,
	RuntimeError,
	"Bootstrap current calculation requires an atomic number (Zeff) profile.",
	)
	errorif(
	eq.ion_density is None and eq.electron_density is None,
	RuntimeError,
	"Bootstrap current calculation requires an ion density or electron density profile.",
	)

[lkadz]

Previously it only guards against both being missing. It does not require atomic_number when ion_density is supplied directly. But the actual Redl formula requires atomic_number in desc/compute/_bootstrap.py#L148:

Zeff = jnp.maximum(1 + 1.0e-14, profile_data["Zeff"])

And per desc/compute/_profiles.py#L472:

def _Zeff(...):
    if profiles["atomic_number"] is not None:
        data["Zeff"] = profiles["atomic_number"].compute(...)
    else:
        data["Zeff"] = jnp.nan * data["0"]

[dpanici]

DESC/desc/equilibrium/equilibrium.py

Line 374 in 7841b52

atomic_number = 1

this is set to 1 if not given, so it would be hard to have kinetic profiles without the atomic number existing. @ddudt also said we want ion density to exist going forward so it is a conscious decision on the user to set

[ddudt]

The Redl formula requires all kinetic profiles: Te, Ti, ne, ni, and Zeff. The only exception is that if ni is not given, it is assumed to be ne/Zeff. atomic_number = 1 as a default in Equilibrium.__init__ but you could initialize an eq with a pressure profile then switch it to kinetic profiles later and forget to set the atomic number.