Add three new C-Mod physics methods with H-mode related variables by zapatace · Pull Request #562 · MIT-PSFC/disruption-py

zapatace · 2026-05-29T14:58:05Z

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.

…, H98)

Co-authored-by: Copilot <copilot@github.com>

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 · 2026-05-29T15:21:02Z

+[cmod.physics.attributes.h98]
+description = "H98 confinement enhancement factor."
+units = "dimensionless"
+validity = [0.05, 1.75]


gtrevisan · 2026-05-29T15:21:08Z

+[cmod.physics.attributes.h_alpha]
+description = "H-alpha line emission intensity."
+units = "W/(m2*sr)"
+validity = [0.0, 161.920898]


gtrevisan · 2026-05-29T15:21:16Z

+[cmod.physics.attributes.pow_thr_LH_Martin]
+description = "Martin 2008 L-H transition power threshold scaling."
+units = "W"
+validity = [0, 0.1e8]


gtrevisan · 2026-05-29T15:21:45Z

+        dict
+            A dictionary with the H-alpha signal (`h_alpha`). In IS brightness units [W/(m2*sr)].
+
+        Last major update by Enrique Zapata zapata_e@mit.edu


you might want to reference this very PR -- take a look at the other docstrings for inspiration/formatting.

I would remove the "Last major update..." line and add the "References" block like other physics methods.

gtrevisan · 2026-05-29T15:22:26Z

+
+        Last major update by Enrique Zapata zapata_e@mit.edu
+        """
+        output = {


I don't think you need this: any unhandled exception will result in returning NaNs.

gtrevisan · 2026-05-29T15:25:38Z

+        # Get btor
+        btor, t_mag = params.data_conn.get_data_with_dims(
+            r"\btor", group="magnetics"
+        )  # tmag: [s]


at this point mention units for both quantities, here, I guess T and s.

gtrevisan · 2026-05-29T15:26:23Z

+        )  # tmag: [s]
+        # Toroidal power supply takes time to turn on, from ~ -1.8 and should be
+        # on by t=-1. So pick the time before that to calculate baseline
+        baseline_indices = np.where(t_mag <= -1.8)


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.

gtrevisan · 2026-05-29T15:27:13Z

+        Scaling from Y R Martin et al 2008 J. Phys.: Conf. Ser. 123 012033
+        DOI 10.1088/1742-6596/123/1/012033
+
+        Last major update by Enrique Zapata


same as above.

gtrevisan · 2026-05-29T15:27:33Z

+        # Get BT
+        btor, t_mag = params.data_conn.get_data_with_dims(
+            r"\btor", group="magnetics"
+        )  # tmag: [s]


same as above.

gtrevisan · 2026-05-29T15:27:36Z

+        )  # tmag: [s]
+        # Toroidal power supply takes time to turn on, from ~ -1.8 and should be
+        # on by t=-1. So pick the time before that to calculate baseline
+        baseline_indices = np.where(t_mag <= -1.8)


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.

+        # Get signals from SPECTROSCOPY tree
+        try:
+            h_alpha, time_halpha = params.data_conn.get_data_with_dims(
+                r"\SPECTROSCOPY::HA_2_BRIGHT", group="SPECTROSCOPY"
+            )  # [mW/(cm2*sr)], [s]
+            # Interpolate Halpha to params.times
+            h_alpha_interp = interp1(time_halpha, h_alpha, params.times)
+            output["h_alpha"] = 10 * h_alpha_interp  # [W/(m2*sr)]
+        except ValueError as e:
+            params.logger.warning("Failed to get H_alpha signal. Returning NaNs.")
+            params.logger.opt(exception=True).debug(e)
+        return output


+        In case of using this signal for ELM detection, it is recommended to
+        use the native time base of the signal to avoid loosing ELMs.
+
+        Parameters
+        ----------
+        params : PhysicsMethodParams
+            The parameters containing the MDSplus connection, shot id and more.
+
+        Returns
+        -------
+        dict
+            A dictionary with the H-alpha signal (`h_alpha`). In IS brightness units [W/(m2*sr)].
+


+        ip = np.abs(ip_df.get("ip")) / 1.0e6  # [A] -> [MA]
+        n_e = density_df.get("n_e") / 1.0e19  # [m^-3] -> [10^19 m^-3]
+        p_input = powers_df.get("p_input") / 1.0e6  # [W] -> [MW]
+        dwmhd_dt = efit_df.get("dwmhd_dt") / 1.0e6  # [W] -> [MW]


+        density_df = CmodPhysicsMethods.get_densities(params=params)
+        areao_df = CmodPhysicsMethods.get_kappa_area(params=params)
+
+        # Get BT
+        btor, t_mag = params.data_conn.get_data_with_dims(
+            r"\btor", group="magnetics"
+        )  # tmag: [s]
+        # Toroidal power supply takes time to turn on, from ~ -1.8 and should be
+        # on by t=-1. So pick the time before that to calculate baseline
+        baseline_indices = np.where(t_mag <= -1.8)
+        btor = btor - np.mean(btor[baseline_indices])
+        btor = np.abs(interp1(t_mag, btor, params.times))  # [T]
+        n_e = density_df.get("n_e") / 1.0e20  # [m^-3] -> [10^20 m^-3]
+        areao = areao_df.get("kappa_area")  # [m^2]
+


yumouwei · 2026-06-06T19:52:42Z

+        except ValueError as e:
+            params.logger.warning("Failed to get H_alpha signal. Returning NaNs.")
+            params.logger.opt(exception=True).debug(e)


I don't think this try-except block would catch any ValueError. The only thing that could fail is R2140-2142 params.data_conn.get_data_with_dims (which can now be shortened to params.get_data_with_dims) which could be automatically caught by the framework.

Since the method will just return output=[nan] if R2140-2142 fails, we can safely move R2140-2145 outside the try-except block, then remove R2135-2137 & R2146-2148.

yumouwei · 2026-06-06T20:06:18Z

+        powers_df = CmodPhysicsMethods.get_power(params=params)
+        efit_df = CmodEfitMethods.get_efit_parameters(params=params)
+        density_df = CmodPhysicsMethods.get_densities(params=params)
+        ip_df = CmodPhysicsMethods.get_ip_parameters(params=params)


All these methods return a dictionary, not a dataframe, so it'd be better to change the variable names to reflect their data type.

I'm always a little wary of calling other physics methods instead of fetching the raw signals directly from MDSplus, and then do some additional computation. The physics methods will return the signals already interpolated onto the requested timebase, so there could be a chance that this would introduce some errors when the original signals are sampled at a much lower frequency, or the timebase of the original signal is severely misaligned with the requested timebase. That being said I've done similar things before, and I don't think this is going to be a major problem here; it's just something to consider for future physics methods that involve using data from e.g. the Thomson system.

yumouwei · 2026-06-06T20:16:24Z

+        References
+        ----------
+        Scaling from eq. 20, ITER Physics Basis Chapter 2
+        https://iopscience.iop.org/article/10.1088/0029-5515/39/12/302/pdf


I would make sure the format here is consistent with the References section in other physics methods' docstring. We've never referenced materials outside the disruption-py repo before; I think we can put this as a standalone bullet point and put the URL to a hyperlink, so something like this:

References -------- - Scaling from eq. 20, [ITER Physics Basis Chapter 2](https://iopscience.iop.org/article/10.1088/0029-5515/39/12/302/pdf) - pull requests: #[562](<link_to_this_pr>)

yumouwei · 2026-06-06T20:18:30Z

+        return output
+
+    @staticmethod
+    @physics_method(columns=["pow_thr_LH_Martin"], tokamak=Tokamak.CMOD)


I'd use all lowercase for the signal name to align with the convention.

I also just found out that RetrievalSettings(run_columns=['pow_thr_LH_Martin'] doesn't call the method correctly unless I switch it to all lowercase

yumouwei · 2026-06-06T20:23:49Z

+        # Get BT
+        btor, t_mag = params.data_conn.get_data_with_dims(
+            r"\btor", group="magnetics"
+        )  # tmag: [s]
+        # Toroidal power supply takes time to turn on, from ~ -1.8 and should be
+        # on by t=-1. So pick the time before that to calculate baseline
+        baseline_indices = np.where(t_mag <= -1.8)
+        btor = btor - np.mean(btor[baseline_indices])
+        btor = np.abs(interp1(t_mag, btor, params.times))  # [T]


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.

yumouwei · 2026-06-06T20:24:22Z

+        output = {
+            "pow_thr_LH_Martin": 1e6 * itpa_power_thr,
+        }
+        return output


We can just do return {"pow_thr_lh_martin": 1e6 * itpa_power_thr}

zapatace and others added 5 commits May 6, 2026 10:50

Added 3 methods C-Mod (H alpha line, LH power thr from Martin Scaling…

be8fec5

…, H98)

Fixes imports physics.py (isort)

d936a98

Co-authored-by: Copilot <copilot@github.com>

Fixed toml linting

636d1af

Fixed pylint

ce027e4

Co-authored-by: Copilot <copilot@github.com>

Merge dev into ezc/h_mode_vars (ready to PR)

a203588

gtrevisan reviewed May 29, 2026

View reviewed changes

gtrevisan requested review from Copilot and yumouwei May 29, 2026 15:28

gtrevisan added the machine: C-MOD Related to the Alcator C-MOD tokamak label May 29, 2026

Copilot started reviewing on behalf of gtrevisan May 29, 2026 15:28 View session

Copilot AI reviewed May 29, 2026

View reviewed changes

yumouwei requested changes Jun 6, 2026

View reviewed changes

Conversation

zapatace commented May 29, 2026

Validation / error analysis

get_h_alpha → h_alpha

get_h98 → h98

get_itpa_pow_thr → pow_thr_LH_Martin

Warnings

Uh oh!

gtrevisan left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yumouwei Jun 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yumouwei Jun 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

`get_h_alpha` → `h_alpha`

`get_h98` → `h98`

`get_itpa_pow_thr` → `pow_thr_LH_Martin`

yumouwei Jun 6, 2026 •

edited

Loading

yumouwei Jun 6, 2026 •

edited

Loading