Add <1/R> to geometry and calulate spr from <1/R>

[theo-brown]

• Add to geometry
• Add to CHEASE loader
• Add to IMAS loader
• Add to LY loader
• Add to EQDSK loader
• Handle division by zero
• Rerun tests and sense check the changes to outputs

[jcitrin]

@felicif , @theo-brown is carrying out the an extension of the geometry inputs to include the 1/ term. Where can this be found in the LY structure?

EDIT: nvm, it's Q0Q

[theo-brown]

Effect of changes on test cases

TLDR: Minor differences in j profiles (and hence Ibs), De, Ve, Qalpha (and hence P alpha). Differences can be pronounced in chi, but only for one or two timesteps before returning to better agreement. These differences tend to not make a difference to overall behaviour of the simulation. Sawteeth behaviour has been altered, but I need a second pair of eyes to ascertain whether it's ok.

iterhybrid_predictor_corrector cases

Observe differences in chi, D_e, j, Qalpha, but without any real impact on the wider behaviour.

iterhybrid_rampup cases

Observe how De is very different for one timestep but then returns to agreement for the next timestep.

test_iterhybrid_rampup_sawtooth

More pronounced flattening effect in the core. @jcitrin is this within the margin of acceptable changes to tests?

[theo-brown]

• A few tests are failing because they compare against torax._src.test_utils.torax_refs. What's the correct way of updating this file?
• I also don't quite understand the following check that is now failing:
  

  torax/torax/_src/core_profiles/tests/initialization_test.py

  Lines 321 to 325 in ec1988b

  	np.testing.assert_allclose(
  	np.mean(j_total1),
  	np.mean(j_ohmic2 + j_bootstrap2),
  	rtol=_TOL,
  	)

  
  Should this hold because the total plasma current is an externally imposed constraint?
• Finally, as there's no LY file bundled with TORAX I can't fully check that the loading is ok. Is Q0Q always present, or should I add a check like in the IMAS case where gm9 is optional?

[jcitrin]

More pronounced flattening effect in the core. @jcitrin is this within the margin of acceptable changes to tests?

Regarding the sawteeth, it's probably just because the sawtooth crash times are now slightly different, and in that specific time slice you see one simulation post-crash and the other is more recovered. Here what matters more is the average over time rather than the difference at one time slice. e.g. if you plot Te[0] of simulations over time, as well as the average difference over multiple sawtooth crashes. If the compare_sim_test scripts are showing small differences for the averaging, then it should be fine.

[jcitrin]

Observe how De is very different for one timestep but then returns to agreement for the next timestep.

Likely related to stiffness and very slight profile differences triggering for one time slice a seemingly large difference. Over time the averages should be quite similar

[jcitrin]

A few tests are failing because they compare against torax._src.test_utils.torax_refs. What's the correct way of updating this file?

We have an internal notebook to regenerate those, which I think should be exposed if you need it. However, the quickest would be to just leave it for now and we will regenerate those refs during internal review. In due course we'll make the regeneration script public and provide instructions on its use

[jcitrin]

I also don't quite understand the following check that is now failing:
Should this hold because the total plasma current is an externally imposed constraint?

Strange that it's failing. By how much? However look again I think that comparing the mean is a bit weird. We should be comparing all elements the cumulative integral (using spr as the metric coefficient), i.e. the plasma current up to a given flux surface. That should hold since the only components of plasma current in the test case is Ohmic and bootstrap. Feel free to change it, otherwise can leave it and we'll look at it during internal review

[jcitrin]

Finally, as there's no LY file bundled with TORAX I can't fully check that the loading is ok. Is Q0Q always present, or should I add a check like in the IMAS case where gm9 is optional?

Q0Q should always be present

[theo-brown]

Observe how De is very different for one timestep but then returns to agreement for the next timestep.

Likely related to stiffness and very slight profile differences triggering for one time slice a seemingly large difference. Over time the averages should be quite similar

Verified with taking a moving average in time over the D_e profiles, the averaged profiles end up matching very well except for two timesteps (a difference which decreases as you increase the width of the averaging window).

from scipy.signal import savgol_filter
import matplotlib.cm as cm

colors = cm.viridis(np.linspace(0, 1, len(original.time)))

# Average over time
smoothed_modified_D_turb_e = savgol_filter(
    modified.profiles.D_turb_e, window_length=5, polyorder=0, axis=0
)
smoothed_original_D_turb_e = savgol_filter(
    original.profiles.D_turb_e, window_length=5, polyorder=0, axis=0
)
for t in range(41):
  plt.plot(smoothed_modified_D_turb_e[t], color=colors[t])
  plt.plot(smoothed_original_D_turb_e[t], "--", color=colors[t])

plt.plot([], [], "k", label="Modified")
plt.plot([], [], "k--", label="Original")
plt.legend()

[theo-brown]

Regarding the sawteeth, it's probably just because the sawtooth crash times are now slightly different

Absolutely correct, verified using same method as above and more careful manual inspection of the nonsmoothed timetraces

[theo-brown]

I also don't quite understand the following check that is now failing:
Should this hold because the total plasma current is an externally imposed constraint?

Strange that it's failing. By how much? However look again I think that comparing the mean is a bit weird. We should be comparing all elements the cumulative integral (using spr as the metric coefficient), i.e. the plasma current up to a given flux surface. That should hold since the only components of plasma current in the test case is Ohmic and bootstrap. Feel free to change it, otherwise can leave it and we'll look at it during internal review

1. Ah, it's failing with a 0.0305 relative difference, and the tolerance is set to 0.03. That's not really worth losing sleep over imo, I think we can just bump the tolerance to 0.04?

2. There's a comment justifying the choice of comparing the mean in test_compare_initial_currents_with_different_initial_j_ohmic that says

torax/torax/_src/core_profiles/tests/initialization_test.py

Lines 170 to 172 in 4b3c09c

	# Due to approximation errors in psi-->j_total conversions, as well as
	# modifications to j_total on axis, we only compare the current profile
	# mean values up to relatively loose tolerance.

To keep this PR well-scoped I might just boost the tolerance, and leave a note.

[theo-brown]

A few tests are failing because they compare against torax._src.test_utils.torax_refs. What's the correct way of updating this file?

We have an internal notebook to regenerate those, which I think should be exposed if you need it. However, the quickest would be to just leave it for now and we will regenerate those refs during internal review. In due course we'll make the regeneration script public and provide instructions on its use

Remaining failing CI is only the tests from this file.

[jcitrin]

Looks great! Happy to bring it in after the one requested change.

[theo-brown]

Thanks for the review @jcitrin. Made the change and updated the necessary test files (changes were order 1e-4). Reminder that the remaining tests failing are due to torax/_src/physics/tests/psi_calculations_test.py and the reference values in torax._src.test_utils.torax_refs.