From be8fec597298c99e6200a9bbec689eaf4a1ae3ad Mon Sep 17 00:00:00 2001 From: zapata_enrique Date: Wed, 6 May 2026 10:50:26 -0400 Subject: [PATCH 1/4] Added 3 methods C-Mod (H alpha line, LH power thr from Martin Scaling, H98) --- disruption_py/machine/cmod/config.toml | 15 +++ disruption_py/machine/cmod/physics.py | 164 +++++++++++++++++++++++++ 2 files changed, 179 insertions(+) diff --git a/disruption_py/machine/cmod/config.toml b/disruption_py/machine/cmod/config.toml index e5d92c811..bf89c0589 100644 --- a/disruption_py/machine/cmod/config.toml +++ b/disruption_py/machine/cmod/config.toml @@ -77,6 +77,16 @@ imas = "/summary/global_quantities/greenwald_fraction/value" units = "dimensionless" validity = [0, 1.5] +[cmod.physics.attributes.h98] +description = "H98 confinement enhancement factor." +units = "dimensionless" +validity = [0.05, 1.75] + +[cmod.physics.attributes.h_alpha] +description = "H-alpha line emission intensity." +units = "W/(m2*sr)" +validity = [0.0, 161.920898] + [cmod.physics.attributes.i_efc] description = "Error field correction (EFC) coil current." units = "A" @@ -184,6 +194,11 @@ units = "dimensionless" description = "Peaking factor of the pressure profile measured by the Thomson scattering diagnostic." units = "dimensionless" +[cmod.physics.attributes.pow_thr_LH_Martin] +description = "Martin 2008 L-H transition power threshold scaling." +units = "W" +validity = [0, 0.1e8] + [cmod.physics.attributes.q0] description = "Safety factor at magnetic axis (EFIT)." imas = "/equilibrium/time_slice(itime)/global_quantities/q_axis" diff --git a/disruption_py/machine/cmod/physics.py b/disruption_py/machine/cmod/physics.py index 23d5cbaaa..28721cf5a 100644 --- a/disruption_py/machine/cmod/physics.py +++ b/disruption_py/machine/cmod/physics.py @@ -22,6 +22,7 @@ from disruption_py.inout.mds import mdsExceptions from disruption_py.machine.cmod.thomson import CmodThomsonDensityMeasure from disruption_py.machine.tokamak import Tokamak +from disruption_py.machine.cmod.efit import CmodEfitMethods class CmodPhysicsMethods: @@ -2128,3 +2129,166 @@ def _is_on_blacklist(shot_id: int) -> bool: or 1150000000 < shot_id < 1150610000 or 1160000000 < shot_id < 1160303000 ) + + @staticmethod + @physics_method(columns=["h_alpha"], tokamak=Tokamak.CMOD) + def get_h_alpha(params: PhysicsMethodParams): + """ + Get the H_alpha line emission intensity. + + The intensity of H-alpha radiance indicates presence of ELMs, and/or + radiative events, and changes of confinement regimes. + + 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)]. + + Last major update by Enrique Zapata zapata_e@mit.edu + """ + output = { + "h_alpha": [np.nan], + } + # 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 + + @staticmethod + @physics_method(columns=["h98"], tokamak=Tokamak.CMOD) + def get_h98(params: PhysicsMethodParams): + """ + Compute H98 by getting tau_E + + Parameters + ---------- + params : PhysicsMethodParams + The parameters containing the MDSplus connection, shot id and more. + + Returns + ------- + dict + A dictionary with the H98 confinement enhancement factor (`h98`), which is dimensionless + + References + ---------- + Scaling from eq. 20, ITER Physics Basis Chapter 2 + https://iopscience.iop.org/article/10.1088/0029-5515/39/12/302/pdf + + Original Author: Andrew Maris (maris@mit.edu) + Last major update by Enrique Zapata (zapata_e@mit.edu) + """ + + # Get parameters for calculating confinement time + 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) + + # Get btor + 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)) + + # Get signals + 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] + wmhd = efit_df.get("wmhd") / 1.0e6 # [J] -> [MJ] + r0 = efit_df.get("rmagx") # [m] + a_minor = efit_df.get("a_minor") # [m] + epsilon = a_minor / r0 + kappa = efit_df.get("kappa") + tau = wmhd / (p_input - dwmhd_dt) + + # Compute 1998 tau_E scaling, taking A (atomic mass) = 2. + tau_98 = ( + 0.0562 + * (np.sign(n_e) * np.abs(n_e) ** 0.41) + * (2**0.19) + * (np.sign(ip) * np.abs(ip) ** 0.93) + * (np.sign(r0) * np.abs(r0) ** 1.97) + * (np.sign(epsilon) * np.abs(epsilon) ** 0.58) + * (np.sign(kappa) * np.abs(kappa) ** 0.78) + * (np.sign(btor) * np.abs(btor) ** 0.15) + * (np.sign((p_input - dwmhd_dt)) * np.abs((p_input - dwmhd_dt)) ** -0.69) + ) + h98 = tau / tau_98 + h98[h98 <= 0] = 0 + output = { + "h98": h98, + } + return output + + @staticmethod + @physics_method(columns=["pow_thr_LH_Martin"], tokamak=Tokamak.CMOD) + def get_itpa_pow_thr(params: PhysicsMethodParams): + """ + Power threshold for L-H transition. + + Parameters + ---------- + params : PhysicsMethodParams + The parameters containing the MDSplus connection, shot id and more. + + Returns + ------- + dict + A dictionary with the power threshold for L-H transition (`pow_thr_LH_Martin`), in Watts [W] + + References + ---------- + 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 + """ + + 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] + + # Estimate power threshold. + itpa_power_thr = ( + 0.0488 + * (np.sign(n_e) * np.abs(n_e) ** 0.717) + * (np.sign(btor) * np.abs(btor) ** 0.803) + * (np.sign(areao) * np.abs(areao) ** 0.941) + ) + output = { + "pow_thr_LH_Martin": 1e6 * itpa_power_thr, + } + return output From d936a9815778ba5fe03dce136a556d023f550489 Mon Sep 17 00:00:00 2001 From: zapata_enrique Date: Wed, 6 May 2026 14:44:50 -0400 Subject: [PATCH 2/4] Fixes imports physics.py (isort) Co-authored-by: Copilot --- disruption_py/machine/cmod/physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/disruption_py/machine/cmod/physics.py b/disruption_py/machine/cmod/physics.py index 28721cf5a..1bd79e92f 100644 --- a/disruption_py/machine/cmod/physics.py +++ b/disruption_py/machine/cmod/physics.py @@ -20,9 +20,9 @@ interp1, ) from disruption_py.inout.mds import mdsExceptions +from disruption_py.machine.cmod.efit import CmodEfitMethods from disruption_py.machine.cmod.thomson import CmodThomsonDensityMeasure from disruption_py.machine.tokamak import Tokamak -from disruption_py.machine.cmod.efit import CmodEfitMethods class CmodPhysicsMethods: From 636d1af536ef1e1b1b53fad42d0f18b7998cff1c Mon Sep 17 00:00:00 2001 From: zapata_enrique Date: Wed, 6 May 2026 14:54:35 -0400 Subject: [PATCH 3/4] Fixed toml linting --- disruption_py/machine/cmod/config.toml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/disruption_py/machine/cmod/config.toml b/disruption_py/machine/cmod/config.toml index bf89c0589..7274567c8 100644 --- a/disruption_py/machine/cmod/config.toml +++ b/disruption_py/machine/cmod/config.toml @@ -186,6 +186,11 @@ imas = "/bolometer/power_radiated_total" units = "W" validity = [0, inf] +[cmod.physics.attributes.pow_thr_LH_Martin] +description = "Martin 2008 L-H transition power threshold scaling." +units = "W" +validity = [0, 0.1e8] + [cmod.physics.attributes.prad_peaking] description = "Peaking factor of the radiated power." units = "dimensionless" @@ -194,11 +199,6 @@ units = "dimensionless" description = "Peaking factor of the pressure profile measured by the Thomson scattering diagnostic." units = "dimensionless" -[cmod.physics.attributes.pow_thr_LH_Martin] -description = "Martin 2008 L-H transition power threshold scaling." -units = "W" -validity = [0, 0.1e8] - [cmod.physics.attributes.q0] description = "Safety factor at magnetic axis (EFIT)." imas = "/equilibrium/time_slice(itime)/global_quantities/q_axis" From ce027e47bd9fad7cf223d1a4bf2279b824e0d510 Mon Sep 17 00:00:00 2001 From: zapata_enrique Date: Wed, 6 May 2026 15:05:29 -0400 Subject: [PATCH 4/4] Fixed pylint Co-authored-by: Copilot --- disruption_py/machine/cmod/physics.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/disruption_py/machine/cmod/physics.py b/disruption_py/machine/cmod/physics.py index 1bd79e92f..22f809949 100644 --- a/disruption_py/machine/cmod/physics.py +++ b/disruption_py/machine/cmod/physics.py @@ -2256,7 +2256,8 @@ def get_itpa_pow_thr(params: PhysicsMethodParams): Returns ------- dict - A dictionary with the power threshold for L-H transition (`pow_thr_LH_Martin`), in Watts [W] + A dictionary with the power threshold for L-H transition + (`pow_thr_LH_Martin`), in Watts [W] References ----------