Add three new C-Mod physics methods with H-mode related variables

[zapatace]

This PR adds three new physics methods to disruption_py/machine/cmod/physics.py, along with their attribute definitions (description, units, validity limits) in disruption_py/machine/cmod/config.toml.

Column	Method	Units	Validity limits
h98	get_h98	dimensionless	[0.05, 1.75]
h_alpha	get_h_alpha	W/(m²·sr)	[0.0, 161.920898]
pow_thr_LH_Martin	get_itpa_pow_thr	W	[0, 0.1e8]

I used the following shot as reference example.

In this example, there are 6 alternating LH transitions followed by HL transition. Notice how H98 approach one, while H-alpha diminish just after Pinput -Prad crosses P_thr from Martin scaling (although not always).

Validation / error analysis

I run the workflow on scale using disruption warning table. I also run disruption-errors script and no tracebacks or errors originate from these 3 methods, 99.99% of fata was retrieved at ~0.015s per shot.

Next, we have the definition and value distributions for the 3 methods on the disruption warning table.

get_h_alpha → h_alpha

H-alpha line emission intensity, useful as a marker of ELMs, radiative events, and confinement-regime transitions. Reads \SPECTROSCOPY::HA_2_BRIGHT [mW/(cm²·sr)] and interpolates onto params.times, converting to SI brightness units W/(m²·sr). Falls back to NaNs (with a warning) if the signal is unavailable. Be adviced that for ELM detection the native signal time base is necessary to avoid losing fast transient events.

Multimodal distribution with median ≈ 14 W/(m²·sr), IQR ≈ [5.7, 32]. Some shots show small negative values (min ≈ −104), and there are discrete spikes/saturation lines at 80.523056 and 161.920898 — the latter is used as the upper validity bound.

Negative and saturated shots were inspected, here you have two examples

I checked with Bob Granetz, and all of this appears to be a unknown diagnostic errors (like wrong offset or maybe some gain adjust that caused the distortion).

get_h98 → h98

H98 energy-confinement enhancement factor, h98 = τ_E / τ_98, where the measured confinement time is τ_E = W_mhd / (P_input − dW_mhd/dt) and τ_98 is the ITER IPB98(y,2) scaling (eq. 20, ITER Physics Basis Ch. 2), evaluated with atomic mass A = 2. Non-physical non-positive values are clipped to 0.

Skewed Gaussian peaked around ~0.45 with a long tail toward 1.0 and a large spike at 0 (clipped non-positive values). Raw stats: median ≈ 0.45, IQR ≈ [0.32, 0.57]; mean/max are inf because P_input − dW_mhd/dt → 0 produces divide-by-zero spikes. The [0.05, 1.75] validity window in config.toml filters both the zero spike and the inf outliers.

get_itpa_pow_thr → pow_thr_LH_Martin

L–H transition power threshold from the Martin 2008 scaling (Y. R. Martin et al 2008 J. Phys.: Conf. Ser. 123 012033):

P_thr = 0.0488 · n_e^0.717 · B_t^0.803 · S^0.941   [MW]  →  ×1e6  [W]

with n_e in 10²⁰ m⁻³ (get_densities), B_t the baseline-subtracted absolute toroidal field, and S the plasma surface area from get_kappa_area. The np.sign(x)·|x|^p form guards against complex results from negative inputs (check the code).

Roughly log-normal distribution, peaked near ~0.28 MW with median ≈ 2.86e5 W and IQR ≈ [2.12e5, 3.64e5] W. A small number of unphysical outliers exist (min ≈ −2.4e7, max ≈ 1.7e8 W) driven by extreme density/area inputs; the [0, 0.1e8] validity range (to be removed by the user).

Warnings

• h98 produces inf when P_input ≈ dW_mhd/dt; relies on user postprocessing.
• h_alpha can be negative and saturates at discrete values on some shots.
• Tine baes for h_alpha matters.

[gtrevisan]

thanks, Enrique! it looks great.

I've left some minor comments after a quick look, I'll review again in the next few days.

[gtrevisan]

imas?

[nbarbour13]

There may be a way to use the spectrometer_visible category for this. I'm uncertain.

[zapatace]

We can use

[cmod.physics.attributes.h_alpha]
description = "H-alpha line emission intensity."
imas = "/spectrometer_visible/channel(i1)/grating_spectrometer/processed_line(i2)/radiance"
units = "W/(m^2*sr)"
validity = [0.0, 161.920898]

However, IMAS units are photons·m⁻²·s⁻¹·sr⁻¹ Do you want to use these units in the method (they are not SI)?

[gtrevisan]

imas?

[nbarbour13]

I did not find a related IMAS entry.

[zapatace]

Me neither. The only L–H-related IDSs are global_quantities/h_mode (confinement state label) and global_quantities/power_loss (actual P_sep), neither of which is the Martin-scaling threshold.

[gtrevisan]

don't you need a comma, here? please double check dimensions.
your slicing might work in the following line, but this might not be what you want -- a list of indices.

[nbarbour13]

I believe this line is ok.
np.where() returns indices if no optional array arguments are provided.
I also checked the dimensions, but please confirm that it behaves as intended.

[yumouwei]

This is interesting. It looks like arr[idx] automatically unpacks the tuple, so btor[baseline_indices] automatically becomes btor[baseline_indices[0]]. I actually didn't know this, which makes me wonder why we chose to do idx, = np.where(...) in the first place...

[zapatace]

You can see that behaviour is the same with this code snippet.

import numpy as np
t_mag = np.linspace(-2.0, 0.0, 11)
btor  = np.arange(11, dtype=float)

baseline_indices = np.where(t_mag <= -1.8)   # tuple -> (array([0, 1]),)
idx, = np.where(t_mag <= -1.8)                # comma idiom -> array([0, 1])

np.array_equal(btor[baseline_indices], btor[idx])   # True
np.mean(btor[baseline_indices]) == np.mean(btor[idx])  # True (both 0.5)

[gtrevisan]

same as above.

[gtrevisan]

same as above.

[Copilot]

Pull request overview

This PR adds three new CMOD physics retrieval methods and registers their metadata (descriptions/units/validity ranges) in the CMOD config, targeting H-mode related quantities for downstream analysis.

Changes:

• Add get_h_alpha to retrieve and interpolate H-alpha brightness onto params.times.
• Add get_h98 to compute the H98 confinement enhancement factor from EFIT, density, power, and magnetics signals.
• Add get_itpa_pow_thr to compute the Martin (2008) L–H power threshold scaling and expose it as pow_thr_LH_Martin, plus new attribute entries in config.toml.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.

File	Description
disruption_py/machine/cmod/physics.py	Adds the three new physics methods (get_h_alpha, get_h98, get_itpa_pow_thr).
disruption_py/machine/cmod/config.toml	Registers new attribute metadata for h98, h_alpha, and pow_thr_LH_Martin.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[nbarbour13]

A significant portion of the data used to compute the IPB98(y,2) scaling were published in the second version of the H mode confinement database ITERHDB.2 (DB2). An ITER H Mode Database Working Group paper detailing updates to this database was published in Nuclear Fusion in 1994. The units listed for line averaged density are 10^19 m^-3, and the paper lists two scaling laws for energy confinement time.
For reference, the IPB98(y,2) scaling can be found in Table 5 of ITER Physics Expert Group on Confinement and Transport et al 1999 Nucl. Fusion 39 2175. This paper also lists 10^19 m^-3 as units of density.

I'll leave this thread open in case a reviewer wants to take a second look.

[yumouwei]

@zapatace Can you make sure this is correct?

[zapatace]

I remember wrestling with this. Thank you, @nbarbour13, for tracking down the references again! The papers state it pretty clearly. This H98 scaling is formulated with cental line averaged densities in [10¹⁹ m⁻³]. Here is some more evidence: updated ITPA DB. Copilot is hallucinating here. While it is true that in power threshold scaling the density is [10²⁰ m⁻³], this is not the case in H98.

[yumouwei]

We might want to distill the btor computation in get_n_equal_1_amplitude, get_h98 and get_itpa_pow_thr into a single physics method similar to D3DPhysicsMethods.get_btor, or just call get_n_equal_1_amplitude to get the btor signal.

I do realize this suggestion is contradictory with my suggestion in R2180 so take this with a grain of salt.

[nbarbour13]

No action was taken on this.

[zapatace]

I can make a helper function for the 3 methods. However, since I need to modify get_n_equal_1_amplitude can we leave this for a follow-up small PR? Up to you!