Add Shaing corrections to Angioni-Sauter model

[theo-brown]

Part of #1406.

Currently in draft form for development purposes. Before merging, will need to move many of the calculations to physics.collisions.

@fcasson @jcitrin for discussions.

Sources

Electrons: https://pubs.aip.org/aip/pop/article/4/5/1375/980556/Electron-transport-processes-close-to-magnetic
Ions: https://pubs.aip.org/aip/pop/article/4/3/771/263553/Ion-transport-process-around-magnetic-axis-in

[theo-brown]

In tandem with this PR I am producing a closer comparison of STEP-like plasmas in JETTO and TORAX; this is taking a while, so I thought I'd get the feature development under way as well.

Below is a plot of the initial state from a modified examples/step_flattop_bgb.py simulation comparing Shaing and Angioni-Sauter contributions for the full radius.

• Once I've got a more comparable JETTO run (requires turning off a bunch of models) I'll add JETTO NCLASS profiles as a reference.
• The chi produced by this model still don't seem to be looking like what we expect from NCLASS - namely, significant peaking in ion transport on the axis.
• However, there is an obvious improvement over Angioni-Sauter in that transport doesn't disappear on-axis.

I've added a smooth exponential transition between the two models, shown here:

[theo-brown]

Some notes from a discussion with Sarah Newton (neoclassical expert at UKAEA) and @fcasson

• NEO solves the force balance, solving for viscosities and transport simultaneously
• Anecdotally, NEO behaves similarly to NCLASS on axis in JET cases simulated in JETTO (chi neo i explodes at the axis)
• NCLASS uses fitted values for the viscosities as functions of collisionality (?), and then computes transport using those viscosities. The fit is not good everywhere.
• The Shaing 1997 model implemented here uses an asymptotic expansion for viscosity- effectively like evaluating NCLASS once and using that value always
• Angioni-Sauter is a fit to NEO data, but the fit is not good everywhere.

[theo-brown]

In integrated scenarios, the difference in value on-axis caused by this model actually makes a minimal difference, because the ion transport is still so tiny. Compare with the effect of adding in a patch to mimic NCLASS behaviour on the examples/step_flattop_bgb.py simulation:

(Solid=AS+Shaing, dashed=AS only)

(Solid=with on-axis patch, dashed=AS only)

(Solid=with on-axis patch, dashed=AS+Shaing)

As a recap of why we're interested in this, current holes at the core could be a huge problem for STEP. As you can see by comparing the q profiles, the ion transport in the core has a noticeable effect on the current hole.

[jcitrin]

Main comment is on the dimensionality of chi, and the need I believe to divide by (dpsi/dr)^2 which will significantly increase the output values and hopefully alleviate the "smallness" you are seeing

[theo-brown]

I've given it a crack. Here's the step_flattop_bgb example with and without Shaing corrections, using the self-consistent potato width. In this case, the potato orbit width for electrons is 6e-4 and for ions is ~0.15 (both in rho norm units).

On-axis dpsi/drhon has machine epsilon added to avoid division by zero, and then is clipped at chi max.

[jcitrin]

The width for ions looks less than 0.15 in the plot?

Also, since the on-axis values of chi is never actually used, maybe for the sake of plot clarity, we can just set chi[0]=chi[1] ?

[theo-brown]

The width looks less because even in the Shaing transport the ion neoclassical contribution is very small compared to the turbulent one, so you only really see the effect when it blows up to infinity at the axis.

With Shaing (L), without Shaing (R)

Comparison of the TORAX models with JETTO-NCLASS:

[jcitrin]

Looked into this a bit closer.

I think there's a missing 2*pi term in the [Wb] to [m] conversion of chi.

This is because the canonical momentum used in the Shaing paper is e\psi, which implies that it's actually divided by from [Wb] psi. The canonical momentum term is actually eRA , and the relation between A and psi comes from

psi = int (B * dS) = int ( curlA * dS) = int_contour(A * dl) = 2piAR . So we see that we must divide by 2pi to get the canonical momentum term as just psi.

So that means that there would be a (2*pi)^2 multiplication factor still needed in the formula, which I think would get us to O(1) of NCLASS, and the extra tuning coefficient default can be set according to that

[jcitrin]

A couple of other points for more simplification:

1. Only chi_i really matters here. I think for simplicity we can only modify chi_i with the Shaing term

2. We can maybe forget about the DeltaPsi? The Shaing terms drop off quickly anyway, likely on the same order at DeltaPsi. So maybe we can just do a blend with a given window width and that's it?

3. Regarding the on-axis point, we can make some decisions like spline extrapolation of Shaing for rho_norm < 0.05 or somesuch (should play with it a bit and see), using similar or the same functions used in the StandardGeometryIntermediates post_init

[theo-brown]

Great catch on the psi conversion, sorry for not being on top of that! STEP's chi profiles are now looking much more like what we expect from NCLASS. I haven't adjusted the scaling factors at all, need to see how this performs in ITER cases as well.

Regarding simplifications:

• Is electron neoclassical transport much smaller in all cases? If so, then yes we can switch to just doing Shaing correction for the ions
• I've dropped the DeltaPsi bit, reverting back to default window width of rho_norm = 0.2
• Currently I've just set chi[0] = chi[1] for simplicity as per your original suggestion. Since it's never used I reckon this or 0 makes most sense, rather than doing splines?

[theo-brown]

If this is good enough, I'll add some tests and move a bunch of the calculations to collisions.py. Lmk if you spot any other TODOs in the meantime.

[theo-brown]

Any particular requirements on tests for this? The Angioni-Sauter test is currently only a regression test against NEOS values. I can test against reference values from running standalone NCLASS if required, but would rather not if there's a simpler option as I would have to first learn how to run NCLASS.

[fcasson]

Any particular requirements on tests for this? The Angioni-Sauter test is currently only a regression test against NEOS values. I can test against reference values from running standalone NCLASS if required, but would rather not if there's a simpler option as I would have to first learn how to run NCLASS.

I don't see any reason the test should be different to the Angioni-Sauter case. NCLASS is probably not better than NEOS?

[jcitrin]

Is electron neoclassical transport much smaller in all cases? If so, then yes we can switch to just doing Shaing correction for the ions

Yes, electron neoclassical transport is always small since it is related to the Larmor radius

I've dropped the DeltaPsi bit, reverting back to default window width of rho_norm = 0.2

This is fine.

Currently I've just set chi[0] = chi[1] for simplicity as per your original suggestion. Since it's never used I reckon this or 0 makes most sense, rather than doing splines?

This is fine.

[jcitrin]

Any particular requirements on tests for this? The Angioni-Sauter test is currently only a regression test against NEOS values. I can test against reference values from running standalone NCLASS if required, but would rather not if there's a simpler option as I would have to first learn how to run NCLASS.

The Angioni-Sauter test is just run against Angioni-Sauter values, to check for regressions. It's not against NEOS (which was used to calibrate Angioni-Sauter).

Can do something similar here: can just update the test to run against the combined Angioni-Sauter+Shaing model.

[jcitrin]

Profiles look very good for STEP. How does it compare for ITER? Can set the prefactor default to be something that minimizes the difference to NCLASS for the 2 cases.