Skip to content

Utilities - FEATURE - Write the response, tearing growth rates and NTV torque to IMAS - #679

Open
logan-nc wants to merge 4 commits into
developfrom
feature/imas-key-outputs
Open

logan-nc wants to merge 4 commits into
developfrom
feature/imas-key-outputs

Conversation

@logan-nc

Copy link
Copy Markdown
Collaborator

Release note

  • Audience: users
  • Numerical impact: none (harness @ de4e338)
  • Migration: none

write_imas now writes more than the least-stable energy per toroidal mode number: the perturbed-equilibrium response (the normal displacement and the normal field by poloidal harmonic on the code's own ψ grid), tearing growth rates, and the NTV torque as a neoclassical momentum source with the torque inside each surface. A stage that did not run writes nothing, and the existing entries are unchanged.

Regression report

================================================================
Case: diiid_n1 — DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium
================================================================

Regression Report: diiid_n1
==========================================================================================================
Ref 1: develop  @ f11d75a0a (2026-11-12)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 6f3e1b5f (pinned), 16 threads/16 BLAS, jit build, warm
Ref 2: feature/imas-key-outputs  @ de4e33826 (2026-11-14)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 6f3e1b5f (pinned), 16 threads/16 BLAS, jit build, warm
----------------------------------------------------------------------------------------------------------
Quantity                                      develop          feature/imas-key-outputs  Diff       Status
----------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        7.979258e-01     7.979258e-01              0.0e+00    OK    
total energy Im(et[1])                        1.122433e-06     1.122433e-06              0.0e+00    OK    
plasma energy Re(ep[1])                       -1.377303e+00    -1.377303e+00             0.0e+00    OK    
vacuum energy Re(ev[1])                       2.175228e+00     2.175228e+00              0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873975e-01     1.873975e-01              0.0e+00    OK    
plasma energy (all)                           [35 elem]        [35 elem]                 0.0e+00    OK    
vacuum energy (all)                           [35 elem]        [35 elem]                 0.0e+00    OK    
total energy (all)                            [35 elem]        [35 elem]                 0.0e+00    OK    
ODE steps (saved)                             1177             1177                      0.0e+00    OK    
ODE steps (total)                             1172             1172                      0.0e+00    OK    
q0                                            1.204212e+00     1.204212e+00              0.0e+00    OK    
q95                                           4.781723e+00     4.781723e+00              0.0e+00    OK    
beta_t                                        1.327024e-02     1.327024e-02              0.0e+00    OK    
beta_n                                        1.372511e+00     1.372511e+00              0.0e+00    OK    
internal inductance li1                       8.842269e-01     8.842269e-01              0.0e+00    OK    
internal inductance li2                       7.080757e-01     7.080757e-01              0.0e+00    OK    
internal inductance li3                       7.304340e-01     7.304340e-01              0.0e+00    OK    
poloidal beta betap1                          6.680739e-01     6.680739e-01              0.0e+00    OK    
poloidal beta betap2                          5.349836e-01     5.349836e-01              0.0e+00    OK    
poloidal beta betap3                          5.518764e-01     5.518764e-01              0.0e+00    OK    
# singular surfaces                           5                5                         0.0e+00    OK    
singular psi locations                        [5 elem]         [5 elem]                  0.0e+00    OK    
singular q values                             [5 elem]         [5 elem]                  0.0e+00    OK    
current beta betaj                            4.236480e-01     4.236480e-01              0.0e+00    OK    
plasma volume                                 1.829472e+01     1.829472e+01              0.0e+00    OK    
plasma current                                1.152130e+00     1.152130e+00              0.0e+00    OK    
mpert                                         35               35                        0.0e+00    OK    
npert                                         1                1                         0.0e+00    OK    
toroidal field bt0                            2.006573e+00     2.006573e+00              0.0e+00    OK    
wall field bwall                              3.880145e-01     3.880145e-01              0.0e+00    OK    
aspect ratio                                  2.845746e+00     2.845746e+00              0.0e+00    OK    
elongation kappa                              1.708322e+00     1.708322e+00              0.0e+00    OK    
q profile (checksum)                          b1a9e91f9145...  b1a9e91f9145...           identical  OK    
pressure profile (checksum)                   3bf10a935abe...  3bf10a935abe...           identical  OK    
Mercier D_I profile (checksum)                b72acb76916b...  b72acb76916b...           identical  OK    
resistive interchange D_R profile (checksum)  01cd20728b25...  01cd20728b25...           identical  OK    
ballooning Delta' profile (checksum)          34a4f4fac3bb...  34a4f4fac3bb...           identical  OK    
island half-widths                            [5 elem]         [5 elem]                  0.0e+00    OK    
Chirikov parameter                            [5 elem]         [5 elem]                  0.0e+00    OK    
||resonant area-weighted field||              5.286415e-04     5.286415e-04              0.0e+00    OK    
dominant-coupling singular values             [3 elem]         [3 elem]                  0.0e+00    OK    
|forcing overlap with dominant mode|          1.419632e-04     1.419632e-04              0.0e+00    OK    
|delta_nominal| of coil set 1                 7.074909e-05     7.074909e-05              0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.246839e-04     1.246839e-04              0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.175435e-07     4.175435e-07              0.0e+00    OK    
PE plasma energy                              3.321081e+00     3.321081e+00              0.0e+00    OK    
PE vacuum energy                              3.174319e+00     3.174319e+00              0.0e+00    OK    
PE surface energy                             5.521990e+00     5.521990e+00              0.0e+00    OK    
PE toroidal torque                            9.083813e-07     9.083813e-07              0.0e+00    OK    
NTV torque FGAR [N·m]                         5.539325e-01     5.539325e-01              0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                8.224109e-02     8.224109e-02              0.0e+00    OK    
Runtime (s)                                   239.1s           247.2s                               --    
||forcing b~|| (root-area-weighted)           4.683171e-04     4.683171e-04              0.0e+00    OK    
resonant area-weighted field b^r              [5 elem]         [5 elem]                  0.0e+00    OK    
==========================================================================================================
Summary: 53 unchanged

Notes for reviewers

This is a minimal first pass. It writes only to slots the current data dictionary already has, and leaves the rest for the restructuring of IMAS for this kind of data.

Result Where it goes
least-stable energy of each n (unchanged) mhd_linear.time_slice[].toroidal_mode[j].energy_perturbed, j = 1..npert
response ξ_n [m] and b_n [T] of each n, by poloidal harmonic same entry, plasma.displacement_perpendicular, plasma.b_field_perturbed.coordinate1; grid type 24 (ρ_pol,N × straight-field-line Fourier), grid.dim1 = √ψ_N, dim2 = m
tearing roots (one per surface uncoupled, one coupled) further toroidal_mode entries after the npert ideal ones, perturbation_type.name = "tearing", growthrate [1/s], frequency [Hz], m_pol_dominant for an uncoupled root
NTV torque, real part, per kinetic method core_sources.source[], identifier neoclassical (401), description names the method: momentum_tor = (dT/dψ_N)/(dV/dψ_N) [N m⁻²], torque_tor_inside [N m], on grid.psi_norm

Choices a reviewer may want to change:

  • Δ' and the resonant fields at the rational surfaces are not written: there is no slot for either.
  • Tearing roots are appended after the ideal entries, so toroidal_mode[1:npert] still holds exactly the ideal energies the FUSE actor reads.
  • The torque coordinate: IMAS wants grid.rho_tor_norm as the primary coordinate of a core profile. When the dd holds the equilibrium (the IMAS-fed case), it is interpolated from dd.equilibrium. For a run from a file there is no exact value, because the domain stops short of the separatrix, so the profiles are written on psi_norm with IMASdd's coordinate check turned off for those fields.
  • The response displacement written is the physical normal displacement ξ_n. It is not ξ^ψ, which has the wrong units for displacement_perpendicular.

Checks:

  • test/runtests_imas.jl: 54/54 pass. New testset: the n = 2 response block lands in toroidal_mode[2]; a uniform torque density comes back uniform, and its cumulative value matches the total; ρ_tor,N is interpolated when the dd has an equilibrium; tearing roots follow the ideal entry with their m.
  • On real runs (scratch script, not committed):
    • DIII-D-like ideal example: displacement_perpendicular equals pe.xi_n_modes; ∫(dV/dψ_N) dψ_N = 18.29476 against params.volume 18.29472, which confirms the volume derivative is per ψ_N.
    • fgar: torque_tor_inside[end] = 0.55448 N m and ∫ momentum_tor dV = 0.55448 N m, against the adaptive-quadrature total of 0.55393 N m (trapezoid on the profile grid, 0.1 %).
    • DIII-D-like SLAYER example: six uncoupled roots written with m = 2..7, equal to tearing.gamma.
  • No harness case calls write_imas, so the regression report cannot move. An IMAS round-trip case could be added once the output settles.

🤖 Generated with Claude Code

https://claude.ai/code/session_01BgNxCG7pGc8HRXGZ3r3psj

logan-nc and others added 4 commits November 14, 2026 22:34
…V torque to IMAS

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BgNxCG7pGc8HRXGZ3r3psj
…ordinate when the dd holds the equilibrium

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BgNxCG7pGc8HRXGZ3r3psj
…rho_tor_norm is unavailable

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BgNxCG7pGc8HRXGZ3r3psj
@logan-nc
logan-nc requested a review from matt-pharr November 15, 2026 04:11
@github-actions github-actions Bot added the feature New capability label Nov 15, 2026
@matt-pharr

Copy link
Copy Markdown
Collaborator

@logan-nc i will review this but we can also make a list of quantities that don't fit properly into IMAS right now and open a pull request upstream to add fields for them!

@logan-nc

Copy link
Copy Markdown
Collaborator Author

@matt-pharr agreed, here is a first list for the upstream pull request.
It is checked against the data dictionary shipped with IMASdd 8.6.3 (mhd_linear, ntms, core_sources, equilibrium, summary) and against the groups of a full gpec.h5.

No slot at all

GPEC quantity gpec.h5 Why it does not fit
Dimensionless (power-normalized) mode energy ForceFreeStates/FreeBoundaryStability/eigenmode_energies energy_perturbed is in J; the FUSE extension stability_metric (units -) is in data_structures_extra only, not upstream
Plasma and vacuum parts of each mode energy eigenmode_plasma_energies, eigenmode_vacuum_energies only one energy per toroidal_mode
Energy matrices W_plasma, W_vacuum, W_freeboundary and the eigenvectors in harmonic space FreeBoundaryStability/W_* no m × m matrix anywhere in mhd_linear
Least vacuum eigenvalue vacuum_eigenvalue none
Energy against the edge surface (edge scan) ForceFreeStates/EdgeScan/* none; the result depends on the truncation, so this matters
Truncation: ψ_lim, q_lim, harmonic range m_low..m_high Info/psilim, qlim, mlow, mhigh only as free text in code.parameters; needed to interpret every stability result
Mercier D_I, resistive interchange D_R, ballooning α and α_crit, ballooning Δ′ profiles LocalStability/* no profile slot; summary.pedestal_fits.*.stability covers only the pedestal α
Rational-surface coefficients E, F, G, H, K, M and D_I at each surface SingularSurfaces/* none
Δ′ matrix (coupled, all surfaces), raw outer-region matrix, coil Δ SingularSurfaces/Delta_prime_matrix, Delta_prime_raw, Delta_coil ntms.mode.deltaw[].value holds one m⁻¹ Rutherford term per mode, not a dimensionless matrix that couples surfaces
Resonant field, penetrated field and resonant current at each rational surface SingularCoupling/resonant_area_weighted_field, penetrated_area_weighted_field, resonant_current none
Chirikov parameter SingularCoupling/chirikov_parameter none
Coupling matrices from boundary field to resonant quantities, and their dominant-mode SVD SingularCoupling/C_*, DominantMode/* none
Response matrices: permeability, reluctance, plasma and surface inductance ResponseMatrices/* em_coupling is axisymmetric by its own definition
Applied (vacuum) and total boundary field spectra forcing_b*, response_b* vacuum.b_field_perturbed is a field on a 2D grid, not a boundary spectrum for each kind of applied field
Perturbed-equilibrium energies (plasma, vacuum, surface) and the ideal toroidal torque PerturbedEquilibrium/Energies/* none
Coil error-field sensitivities: field of each coil as designed, its derivatives with respect to shift and tilt, overlap with the dominant mode ErrorFields/CoilSensitivities/* none
Kinetic δW (energy) of the NTV calculation KineticForces/<method>/total_energy none
Imaginary part of the kinetic torque KineticForces/<method>/T, dTdpsi core_sources is real; this PR drops it
Torque by bounce harmonic ℓ KineticForces/<method>/dTdpsi_ell none
Validity of the kinetic model: banana, potato and Larmor widths, scale lengths, valid region KineticForces/Validity/* none

Fits in part

  • Island width: ntms.mode.width is a full width in m, and GPEC writes a half width in ψ_N. This needs a conversion and one entry per surface under ntms rather than mhd_linear.
  • NTV torque: mhd_linear.toroidal_mode.plasma.ntv exists, but it is a stress tensor (N m⁻²) on the mode grid, not a torque density on a flux grid. core_sources as used here is the better home. It needs the imaginary part and the method as a typed field instead of the description.
  • Other response components (ξ^ψ, tangential, R/Z/φ components): displacement_parallel and b_field_perturbed.coordinate2/3 exist, but the basis depends on the coordinate_system blocks, which are not filled.

Points on this PR found while checking

  1. Units: energy_perturbed is in J, but et is the power-normalized generalized eigenvalue (ForceFreeStates/Free.jl). This predates the PR (the line is unchanged).
  2. Fourier kernel: the dictionary describes each mode as exp(i(n_tor·φ − m_pol·θ − 2π f t − phase)). GPEC uses exp(i(mθ − nζ)) (docs/src/conventions.md). The coefficients written look like the complex conjugate of what the dictionary expects. The handedness of φ (COCOS 11) and ζ against φ also need checking before the response is read by other codes.
  3. Grid type 24 (straight-field-line Fourier) is right for Hamada, PEST and Boozer, but not for the equal-arclength Jacobian, which is type 25.
  4. Tearing type: the dictionary has Tearing as index 13, and eKINK as index 12 could label the ideal entries. Setting perturbation_type.index lets readers match on it.

🤖 Generated with Claude Code

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

feature New capability

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants