comparing available energy

[unalmis]

NCSX

W7-X

[unalmis]

Source code

• https://github.com/unalmis/DESC/blob/master/desc/compute/_turbulence.py
• https://github.com/unalmis/DESC/blob/master/desc/objectives/_turbulence.py

[unalmis]

@RalfMackenbach

Benchmark setup

• Environment: conda run -n desc-env
• Equilibrium: desc.examples.get("DSHAPE")
• Flux surface: rho = 0.5
• Gradients: omn = 1.0, omt = 1.0
• Tokamak grid: grid.num_zeta == 1 via LinearGrid(..., N=0, ...)
• Bounce interpolation: nufft_eps = 0.0
• Field-line interval: one poloidal transit between adjacent global maxima of |B|
• Bounce wells: Bounce2D.points, not AEpy/BAD root finding
• Energy quadrature: same generalized Laguerre nodes/weights for DESC and AEpy

Normalization

On the same wells and bounce integrals,

DESC _ae = 4 * sqrt(pi) * AEpy.AE_per_lam_per_z

For the final output

desc available energy = (2/3) * AE_AEpy_ftvol

as expected.

Results

rho=0.5
omn=1.0, omt=1.0
iota=8.325000000000e-01
theta_bmax=3.141592653590e+00
alpha_start=5.491651220631e+00
num_field_periods=2
grid_num_zeta=1
Y_B=64
expected_global_maxima_zeta=(4.724476564731e+00, 1.227184630309e+01)
observed_global_maxima_zeta=(4.724854095907e+00, 1.227183431740e+01)
observed_global_maxima_alpha=(5.491651220631e+00, 4.439217681678e+00)
observed_global_maxima_B=(2.450379781484e-01, 2.450379795792e-01)
num_pitch=17, num_quad=16, num_energy=12
valid_wells=17
fieldline_weight_dl_over_B=2.170494463726e+02
fieldline_weight_from_V_psi_relation=2.132621216357e+02
max_rel_DESC_ae_vs_AEpy_integrand=4.642715e-16
max_abs_DESC_ae_vs_AEpy_integrand=1.065120214250e-15
DESC_eq_compute_available_energy=4.198763958449e-03
DESC_manual_scalar=4.186846415809e-03
DESC_compute_minus_manual=1.191754263960e-05
AEpy_flux_tube_scalar=6.280269623714e-03
AEpy_to_DESC_conversion_factor=6.666666666667e-01
AEpy_DESC_normalized=4.186846415809e-03
DESC_minus_AEpy_DESC_normalized=1.191754263960e-05
DESC_over_AEpy_DESC_normalized=1.002846424601e+00

Conclusion

The local per-pitch/per-well formula agrees to roundoff:

max_rel_DESC_ae_vs_AEpy_integrand = 4.642715e-16

The registered DESC compute result and the AEpy result converted to DESC
normalization agree within about 2.85e-3 relative for this low-resolution
benchmark configuration.

benchmark_available_energy_desc_vs_aepy.py

[RalfMackenbach]

Looks nice -- I quickly tried to reproduce the value of total AE on my local installation of AEpy, but failed to get a similar value [Total AE for DSHAPE: 2.109322e-02 (normalized units)]. Pipeline is

Get DSHAPE VMEC input namelist -> Generate VMEC equilibrium -> Generate GIST at ρ=0.5 -> Calculate AE using AEpy's GIST functionality.

Before going further, could I check the following:

1. Is rho=sqrt(s) in VMEC (i.e. is setting s=1/4 equivalent to your rho=1/2)?
2. Could you verify independence of npol (i.e., does including two wells keep the end result unchanged), making sure the normalisation is correct.
3. Could you generate a bounce-well plot, as below?

An additional exercise are the asymptotes discussed here: https://www.cambridge.org/core/journals/journal-of-plasma-physics/article/trappedparticle-precession-and-modes-in-quasisymmetric-stellarators-and-tokamaks-a-nearaxis-perspective/7E9CFB3C4C7DE910F43F3D1AB023C7EB . We can talk about it in a meeting perhaps?