DESC Effective ripple ε

Resolves #519 . Automatically differentiable

Apr 21 '24 06:04 unalmis

| -------------------------------------- test_build_transform_fft_midres test_build_transform_fft_highres test_equilibrium_init_lowres test_objective_compile_atf test_objective_compute_atf test_objective_jac_atf test_perturb_1 test_proximal_jac_atf test_proximal_freeb_compute test_solve_fixed_iter_compiled test_build_transform_fft_lowres test_equilibrium_init_medres test_equilibrium_init_highres test_objective_compile_dshape_current test_objective_compute_dshape_current test_objective_jac_dshape_current test_perturb_2 test_proximal_freeb_jac test_solve_fixed_iter test_LinearConstraintProjection_build

Apr 21 '24 08:04 github-actions[bot]

Just sharing some outputs from the tests I added. We know that the bounce integral stuff works correctly, so it should be easier to debug this if necessary. ripple

Apr 30 '24 04:04 unalmis

Transcribing the discussion in meeting:

check the |B| calculated from the eq, as you increase the number of points you evaluate along zeta, does this data show the wiggle? or is the wiggle just after you spline and is in the splined values?
We can bound the fourier resolution of |B| on a surface by thinking of the resolution of R,Z, and |B| is something like a cubic poly in R,Z

May 01 '24 19:05 dpanici

Check STELLOPT and NCSX Paper for flux surface integral form of effective ripple

May 08 '24 19:05 dpanici

Figure_2_adaptive Ok after the changes made in the last few commits, we can go to higher resolution, and I can now get curves with features we expect. This one looks similar to what neo gives, besides a normalization factor. I'd be surprised if we should expect tight agreement between what neo and desc due to implementation differences in code. If anyone thinks otherwise let me know.

Many of these neoclassical quantities are surface integrals over phase space, which involves summing a set of bounce integrals throughout the plasma volume and then a simple integration over a velocity space coordinate, i.e. the pitch angle lambda. The hard part of doing all the bounce integrals is tested to be correct. The velocity integration is reduced to either calling jax or quadax, whose algorithms are also well-tested algorithmically. In theory if these two operations 1. bounce integral and 2. velocity integral are implemented correctly then we're done.

In addition to the already passing algorithmic unit tests, there's benefit in adding physics-based tests for these two operations so that we know the full end-to-end computation produced by DESC gives the correct result. For the bounce integration, we match the analytic approximation for the bounce averaged binormal drift $\langle \omega_D \rangle$. I would like to write a similar test for the velocity integration. Ideally we'd have an analytic expression for the effective ripple or one of the gammas.

@rahulgaur104 Do you know of any existing analytical work to obtain an expressions for the effective ripple for some simplified geometry? Something that's not axisymmetric / tokamak like? (Since effective ripple ~ cvdrift0 and the radial component of drift vanishes for omnigenous geometry, things like low beta tokamak won't be the best test).

If not, I think it would be just as effective to write a test where we replace cvdrift0 in the effective ripple computation with cvdrift. Let's call this modified effective ripple. I can compute the velocity space integration to compute $\int d k \langle \omega_D \rangle$ analytically for the low beta tokamak example, and probably this modified effective ripple as well. Then our physics-based end-to-end test for would be comparing numerical modified ripple to this analytical approximation.

May 27 '24 23:05 unalmis

Check STELLOPT and NCSX Paper for flux surface integral form of effective ripple

I couldn't find a paper like this. I just wrote out the math to do it myself.

I should mention the current implementation of the code is agnostic to whether you use a long field line integral or attempt an average. It's just a question of reducing memory since long field line integral is less efficient speed and memory wise.

Edited to make more clear: Also the equivalence for converting a long field line integral to a flux surface integral as is done in Velasco's gamma alpha paper is clear so long as we assume irrational surface. You can always write it as a surface integral of some function, in this case the function is a weighted sum over all the wells on the surface.

So it is fine to truncate the length of the field line after one rotation, assuming the field line is long enough to capture the endpoints of the well. That is pick a set of field lines (rho, alpha), detect and integrate bounce integrals over one toroidal rotation, and then integrate over alpha. Since everything is continuous, as the number of nodes tend to infinity these should converge to the same thing, again assuming irrational surface. However, the order at which all the bounce integrals in the plasma volume are added up is changed, so the convergence rate can differ.

May 28 '24 00:05 unalmis

This is great progress Kaya! Which equilibrium did you use to obtain the blue plots? I agree that the match may not be perfect but at least the trends should match.

I can try and work out and (semi)-analytical setup for a rotating ellipse, just like Hegna's paper but it will be extra tedious (due to lack of toroidal symmetry). A better way to create a test is to try near-axis expansions and calculate these quantities near the magnetic axis.

What about comparing with NEO? Let's say we can scale our ripple by a constant, how does it look compared to the result from NEO.

I don't think cvdrift(binormal drift) zero is any measure of the ripple or a good indicator of the actual ripple (coming from cvdrift0(radial drift)).

Let me think more about it.

May 28 '24 00:05 rahulgaur104

To clarify the top plot I have shown is an intermediate quantity ∫ db ∑ⱼ Hⱼ² / Iⱼ that is essentially the unnormalized effective ripple. The bottom plot is effective ripple. I'll play around with normalization with what neo gives, and see how they differ then. When I said the bottom plot looks similar to neo in the previous post I mean the overall trend.

May 28 '24 00:05 unalmis

I don't think cvdrift(binormal drift) zero is any measure of the ripple or a good indicator of the actual ripple (coming from cvdrift0(radial drift)).

Let me think more about it.

I'm thinking we could replace the dependency data["cvdrift0"] = data["cvdrift"] as we pass in this dependency to compute the effective ripple. Then we'll get an output quantity that's computed in the code in an identical manner to the true effective ripple. If this output matches it's analytic expression well, and we also know separately that cvdrift0 can be calculated correctly in DESC when needed, then we can likely conclude that the code will compute the effective ripple correctly when we actually give it cvdrift0

May 28 '24 00:05 unalmis

I'm thinking we could replace the dependency data["cvdrift0"] = data["cvdrift"] as we pass in this dependency to compute the effective ripple. Then we'll get an output quantity that's computed in the code in an identical manner to the true effective ripple. If this output matches it's analytic expression well, and we also know separately that cvdrift0 can be calculated correctly in DESC when needed, then we can likely conclude that the code will compute the effective ripple correctly when we actually give it cvdrift0

Ah, ok. I understand now. Yes, that's a good test, just to check the calculations done by the code

May 28 '24 01:05 rahulgaur104

cvdrift0 is not the same as kappa_g, and epsilon effective should definitely use kappa_g. neemov's kappa_g is the same as desc kappa_g.

May 29 '24 07:05 unalmis

You don't have to use cvdrift0 is you don't want to but cvdrift0/kappa_g can at most be a constant.

May 29 '24 07:05 rahulgaur104

https://github.com/PlasmaControl/DESC/tree/master/publications/dudt2024 benchmark with these cases and data from Daniel

Jul 03 '24 19:07 dpanici

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

Jul 11 '24 23:07 review-notebook-app[bot]

I think the final ripple that people care about is epsilon_eff not epsilon_eff^(3/2). The rule of thumb is to have a ripple(epsilon_eff) < 1-2%. So, I would suggest doing that for the EffectiveRipple objective.

Sep 15 '24 09:09 rahulgaur104

Codecov Report

Attention: Patch coverage is 96.77419% with 4 lines in your changes missing coverage. Please review.

Project coverage is 95.55%. Comparing base (0e1fdb9) to head (3d34a23). Report is 1450 commits behind head on master.

Files with missing lines	Patch %	Lines
desc/objectives/_neoclassical.py	93.33%	3 Missing :warning:
desc/plotting.py	83.33%	1 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1003      +/-   ##
==========================================
+ Coverage   95.53%   95.55%   +0.02%     
==========================================
  Files          96       98       +2     
  Lines       25044    25115      +71     
==========================================
+ Hits        23925    23999      +74     
+ Misses       1119     1116       -3

Files with missing lines	Coverage Δ
desc/compute/__init__.py	`100.00% <ø> (ø)`
desc/compute/_equil.py	`100.00% <ø> (ø)`
desc/compute/_neoclassical.py	`100.00% <100.00%> (ø)`
desc/compute/_omnigenity.py	`99.37% <ø> (ø)`
desc/compute/_profiles.py	`99.19% <100.00%> (ø)`
desc/compute/_stability.py	`100.00% <ø> (ø)`
desc/equilibrium/equilibrium.py	`96.07% <100.00%> (ø)`
desc/integrals/_bounce_utils.py	`84.37% <100.00%> (ø)`
desc/integrals/_interp_utils.py	`100.00% <100.00%> (ø)`
desc/integrals/basis.py	`93.04% <100.00%> (ø)`
... and 8 more

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

Sep 16 '24 06:09 codecov[bot]

I think the final ripple that people care about is epsilon_eff not epsilon_eff^(3/2). The rule of thumb is to have a ripple(epsilon_eff) < 1-2%. So, I would suggest doing that for the EffectiveRipple objective.

It looks like NEO returns $\varepsilon_{eff}^{3/2}$, but I agree we should return the raw value

Sep 18 '24 19:09 ddudt

I think the final ripple that people care about is epsilon_eff not epsilon_eff^(3/2). The rule of thumb is to have a ripple(epsilon_eff) < 1-2%. So, I would suggest doing that for the EffectiveRipple objective.

It looks like NEO returns $\varepsilon_{eff}^{3/2}$, but I agree we should return the raw value

Are you sure? All the plots I see have the 3/2 power on it, and to me that makes sense since that's the form of it that occurs in the transport equations.

Sep 22 '24 11:09 dpanici

i haven't changed the 3/2 power on effective ripple because I don't know what to change the name to if I do that

Sep 24 '24 23:09 unalmis

i haven't changed the 3/2 power on effective ripple because I don't know what to change the name to if I do that

I think that having the 3/2 is fine, I have never seen anyone use just $\epsilon_{eff}$, it always has the 3/2 power on it. As long as the documentation in list of variables makes it clear what it is and that it is raised to the 3/2 power it is fine with me

Sep 25 '24 15:09 dpanici

i haven't changed the 3/2 power on effective ripple because I don't know what to change the name to if I do that

I think that having the 3/2 is fine, I have never seen anyone use just ϵ e f f , it always has the 3/2 power on it. As long as the documentation in list of variables makes it clear what it is and that it is raised to the 3/2 power it is fine with me

The thing without the 3/2 power is literally the effective ripple. Plenty of papers don't use the 3/2 factors. Especially experiment-based papers. But it's fine whatever you choose.

Sep 25 '24 16:09 rahulgaur104

i haven't changed the 3/2 power on effective ripple because I don't know what to change the name to if I do that

I think that having the 3/2 is fine, I have never seen anyone use just ϵ e f f , it always has the 3/2 power on it. As long as the documentation in list of variables makes it clear what it is and that it is raised to the 3/2 power it is fine with me

The thing without the 3/2 power is literally the effective ripple. Plenty of papers don't use the 3/2 factors. Especially experiment-based papers. But it's fine whatever you choose.

Sorry, I did not realize. I am fine with either

Sep 25 '24 16:09 dpanici

make second compute quantity for 3/2

have objective be on the one without the 3/2

Sep 25 '24 18:09 dpanici

Check against WISTELL

Oct 02 '24 17:10 dpanici

https://github.com/landreman/vmec_equilibria/tree/master/W7-X/HighMirror one more data pt, matt has a DKES run of epsilon effective apparently along with this vacuum W7X, I re-ran in DESC and ran eps eff, our min values match at least, this is not so useful without a DKES res scan as well though, just putting as I saw the data: (also, the x coordinate from the dkes data is "r/a" acc. to the text file? not sure if that means rho but I assume it does?)

Oct 02 '24 18:10 dpanici

Thank for this. How do these DESC plots change as the DESC equilibrium resolution varies? I'd expect accuracy in effective ripple etc. via a field line integration requires higher DESC equilibrium resolution than a typical objective. After all the integrals are over ripples, which will not appear if higher fourier harmonics are suppressed.

It i interesting that DESC converges with few toroidal transits

Oct 02 '24 20:10 unalmis

Note that DKES is generally not great for eps_eff. Since eps_eff is computed as a limit as collisionality goes to zero, and DKES has problems converging for low collisionality. It also usually comes with error bars?

Oct 02 '24 21:10 f0uriest

also need to add the new objective to docs (or do we want to wait until #1290 for that?)

Oct 31 '24 01:10 f0uriest

also need to add the new objective to docs (or do we want to wait until #1290 for that?)

That works. I suggest merging the PRs that are upstream to 1290 into master in the order I shared in the slack. If you want to make explicit to users that these methods will be replaced with 1290, then one can add a soft block to the compute function in this PR:

errorif(not kwargs.get("ignore", False), NotImplementedError, msg="PR #1290 will replace this.")

Oct 31 '24 03:10 unalmis

Did we ever run a comparison against WISTELL?

Ok, I spent an afternoon on this. Their code errors with singular integral exceptions. Maybe it's something simple. Waiting on their advise.

Oct 31 '24 05:10 unalmis

DESC DESC copied to clipboard

Effective ripple ε

Codecov Report

DESC
DESC copied to clipboard