From e4c43d103c26d739ded8ddb6b5c051487468e9ff Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 5 Mar 2025 03:27:08 +0800 Subject: [PATCH 01/28] Fix EAST get_mirnov_std method --- disruption_py/machine/east/physics.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 922a1afe9..7e827f3d7 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -1465,22 +1465,22 @@ def get_mirnov_std(params: PhysicsMethodParams): Last major update: 2014/11/22 by William Wei """ - mirnov_sensor = "cmp1t" # default one used in disruption_warning_database.m - mirnov_std = np.full(len(params.times), np.nan) mirnov_std_normalized = [np.nan] # Get the toroidal magnetic field. btor = EastPhysicsMethods.get_btor(params)["btor"] - # Get the Mirnov signal + # Get the Mirnov signal from \cmp1t (5 MHz) time_window = 0.001 bp_dot, bp_dot_time = params.mds_conn.get_data_with_dims( - r"\Mirnov_sensor" + mirnov_sensor, tree_name="east" - ) + r"\cmp1t", tree_name="east" + ) # [T/s], [s] for i, time in enumerate(params.times): - (indices,) = np.where((time - time_window) < bp_dot_time < time) - mirnov_std[i] = np.nanstd(bp_dot[indices]) + (indices,) = np.where( + ((time - time_window) < bp_dot_time) & (bp_dot_time < time) + ) + mirnov_std[i] = np.nanstd(bp_dot[indices], ddof=1) mirnov_std_normalized = mirnov_std / abs(btor) From 23b30364195cfeec22d3f48210fe61f9355cf7ff Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sat, 15 Mar 2025 03:34:51 +0800 Subject: [PATCH 02/28] Remove h98 from test column (output matches matlab script, wrong data in sql table) --- disruption_py/machine/east/config.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index e9bf40fe8..7c9ed1907 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -22,7 +22,6 @@ expected_failure_columns = [ "dli_dt", "dwmhd_dt", "greenwald_fraction", - "h98", "ip_error_rt", "kappa", "kappa_area", @@ -69,7 +68,6 @@ test_columns = [ "dli_dt", "dwmhd_dt", "greenwald_fraction", - "h98", "ip", "ip_error", "ip_error_normalized", From fc5143c5821aa9e96075fae29810c15857661fec Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sat, 15 Mar 2025 04:42:58 +0800 Subject: [PATCH 03/28] Remove parea, pchisq, pconvergence, pkappa_area, pq0, pqstar from testing (signals not in table/all nulls, match MATLAB) --- disruption_py/machine/east/config.toml | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 7c9ed1907..f6596da51 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -39,12 +39,6 @@ expected_failure_columns = [ "p_input", "p_oh", "p_rad", - "parea", - "pchisq", - "pconvergence", - "pkappa_area", - "pq0", - "pqstar", "q0", "q95", "qstar", @@ -97,17 +91,11 @@ test_columns = [ "p_rad", "p_rad_rt", "paminor", - "parea", "pbeta_n", "pbeta_p", - "pchisq", - "pconvergence", "pkappa", - "pkappa_area", "pli", - "pq0", "pq95", - "pqstar", "prad_peaking", "pwmhd", "q0", From e390e3d2562c52e60105aed4d09e09da34ba5c01 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 18 Mar 2025 02:47:21 +0800 Subject: [PATCH 04/28] Move AXUV calibration factors to EastUtilMethods.get_axuv_calib_factors() --- disruption_py/machine/east/physics.py | 119 ++---------------------- disruption_py/machine/east/util.py | 129 ++++++++++++++++++++++++++ 2 files changed, 136 insertions(+), 112 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 7e827f3d7..46e605d44 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -1272,112 +1272,7 @@ def get_prad_peaking(params: PhysicsMethodParams): # TODO: Move calibration factors to a separate settings file # Calibration factors - fac_1 = ( - np.array( - [ - 1.3681, - 1.3429, - 1.3215, - 1.3039, - 1.2898, - 1.2793, - 1.2723, - 1.2689, - 1.2689, - 1.2723, - 1.2793, - 1.2898, - 1.3039, - 1.3215, - 1.3429, - 1.3681, - ] - ) - * 1e4 - ) - fac_2 = np.array( - [ - [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], - [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], - ] - ) # factors of Amp.Gain - fac_3 = np.array([1, 1, 1, 1]) # cross calibration factors between arrays - fac_4 = 1 * 1e-3 # unit convert - # Fac5=2.5 # corrected factor by cross calibration with foil bolometer - maj_r = 1.85 - del_r = ( - np.array( - [ - 3.6, - 3.6, - 3.5, - 3.4, - 3.4, - 3.3, - 3.3, - 3.2, - 3.2, - 3.1, - 3.1, - 3.0, - 3.0, - 2.9, - 2.9, - 2.8, - 2.9, - 2.8, - 2.8, - 2.8, - 2.8, - 2.8, - 2.8, - 2.7, - 2.7, - 2.7, - 2.7, - 2.7, - 2.6, - 2.6, - 2.6, - 2.6, - 2.6, - 2.6, - 2.6, - 2.6, - 2.7, - 2.7, - 2.7, - 2.7, - 2.7, - 2.8, - 2.8, - 2.8, - 2.8, - 2.8, - 2.8, - 2.9, - 2.8, - 2.9, - 2.9, - 3.0, - 3.0, - 3.1, - 3.1, - 3.2, - 3.2, - 3.3, - 3.3, - 3.4, - 3.4, - 3.5, - 3.6, - 3.6, - ] - ) - * 0.01 - ) + calib_factors = EastUtilMethods.get_axuv_calib_factors() # Get XUV data (xuvtime,) = params.mds_conn.get_dims(r"\pxuv1", tree_name="east_1") @@ -1398,14 +1293,14 @@ def get_prad_peaking(params: PhysicsMethodParams): signal_smoothed = smooth(signal, smoothing_window) xuv[:, ichord] = ( signal_smoothed - * fac_1[ichan] - * fac_2[iarray][ichan] - * fac_3[iarray] - * fac_4 + * calib_factors["Fac1"][ichan] + * calib_factors["Fac2"][iarray][ichan] + * calib_factors["Fac3"][iarray] + * calib_factors["Fac4"] * 2 * np.pi - * maj_r - * del_r[ichan] + * calib_factors["Maj_R"] + * calib_factors["Del_r"][ichan] ) # from Duan Yanming's program # Correction for bad channels (from Duan Yanmin's program) diff --git a/disruption_py/machine/east/util.py b/disruption_py/machine/east/util.py index ce92c2102..157f54fc8 100644 --- a/disruption_py/machine/east/util.py +++ b/disruption_py/machine/east/util.py @@ -80,3 +80,132 @@ def retrieve_ip(mds_conn: MDSConnection, shot_id: int): # Subtract baseline offset ip = EastUtilMethods.subtract_ip_baseline_offset(ip, ip_time) return ip, ip_time + + @staticmethod + def get_axuv_calib_factors() -> dict: + """ + Return the calibration factors for the AXUV arrays + + Used by get_radiated_power and get_prad_peaking. + + - fac_1: (unknown) + - fac_2: factors of Amp.Gain + - fac_3: cross calibration factors between arrays + - fac_4: unit convert + - fac_5: corrected factor by cross calibration with foil bolometer + - maj_r: major radius (of the machine cross-section?) + - del_r: (unknown) + """ + fac_1 = ( + np.array( + [ + 1.3681, + 1.3429, + 1.3215, + 1.3039, + 1.2898, + 1.2793, + 1.2723, + 1.2689, + 1.2689, + 1.2723, + 1.2793, + 1.2898, + 1.3039, + 1.3215, + 1.3429, + 1.3681, + ] + ) + * 1e4 + ) + fac_2 = np.array( + [ + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + ] + ) + fac_3 = np.array([1, 1, 1, 1]) + del_r = ( + np.array( + [ + 3.6, + 3.6, + 3.5, + 3.4, + 3.4, + 3.3, + 3.3, + 3.2, + 3.2, + 3.1, + 3.1, + 3.0, + 3.0, + 2.9, + 2.9, + 2.8, + 2.9, + 2.8, + 2.8, + 2.8, + 2.8, + 2.8, + 2.8, + 2.7, + 2.7, + 2.7, + 2.7, + 2.7, + 2.6, + 2.6, + 2.6, + 2.6, + 2.6, + 2.6, + 2.6, + 2.6, + 2.7, + 2.7, + 2.7, + 2.7, + 2.7, + 2.8, + 2.8, + 2.8, + 2.8, + 2.8, + 2.8, + 2.9, + 2.8, + 2.9, + 2.9, + 3.0, + 3.0, + 3.1, + 3.1, + 3.2, + 3.2, + 3.3, + 3.3, + 3.4, + 3.4, + 3.5, + 3.6, + 3.6, + ] + ) + * 0.01 + ) + # Use original names in the MATLAB script + return { + "Fac1": fac_1, + "Fac2": fac_2, + "Fac3": fac_3, + "Fac4": 1e-3, + "Fac5": 2.5, + "Maj_R": 1.85, + "Del_r": del_r, + } From 1ad9d5aad2e499f37b616103f10fbaf25042f6b7 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:17:22 +0800 Subject: [PATCH 05/28] Add get_radiated_power method and get_raw_axuv_data function --- disruption_py/machine/east/physics.py | 161 +++++++++++++++++++------- disruption_py/machine/east/util.py | 1 + 2 files changed, 119 insertions(+), 43 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 46e605d44..91ba48ba9 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -441,10 +441,110 @@ def get_density_parameters(params: PhysicsMethodParams): return {"n_e": ne, "greenwald_fraction": greenwald_fraction, "dn_dt": dn_dt} + @staticmethod + def get_raw_axuv_data(params: PhysicsMethodParams): + """ + Get the raw (uncalibrated) data from the AXUV arrays for calculating + the total radiated power and prad_peaking. The AXUV diagnostic contains + 4 arrays of 16 channels each. + + Note that the original implementations in Prad_bulk_xuv2014_2016.m and + get_prad_peaking.m are slightly different in where the calibration factors + are applied. This difference isn't explained in the scripts so I keep + these functions different. + + TODO: should this function be put in EastUtilMethods? + """ + # Get XUV data + (xuvtime,) = params.mds_conn.get_dims(r"\pxuv1", tree_name="east_1") # [s] + # There are 64 AXUV chords, arranged in 4 arrays of 16 channels each + xuv = np.full((len(xuvtime), 64), np.nan) + + dt = xuvtime[1] - xuvtime[0] + smoothing_time = 1e-3 + smoothing_window = int(max([round(smoothing_time / dt, 1)])) + + for iarray in range(4): + for ichan in range(16): + ichord = 16 * iarray + ichan + # TODO: confirm the actual node of each chord + signal = params.mds_conn.get_data( + r"\pxuv" + str(ichord + 1), tree_name="east_1" + ) + # Subtract baseline + signal = signal - np.mean(signal[:100]) + # TODO: change this to causal smoothing + xuv[:, ichord] = smooth(signal, smoothing_window) * 1e3 # [kW] -> [W] + return xuv, xuvtime + + @staticmethod + @physics_method(columns=["p_rad"], tokamak=Tokamak.EAST) + def get_radiated_power(params: PhysicsMethodParams): + """ + Calculate total radiated power in bulk plasma from the AXUV arrays. + + Parameters + ---------- + params : PhysicsMethodParams + Parameters containing MDS connection and shot information + + Returns + ------- + dict + A dictionary containing the following keys: + - 'p_rad' : array + Radiated power [W]. + + References + ------- + - Prad_bulk_xuv2014_2016.m (Not in repo) + + TODO: Update docstring + """ + # Get the raw AXUV data + try: + xuv, xuvtime = EastPhysicsMethods.get_raw_axuv_data(params) # [W], [s] + except mdsExceptions.MdsException: + return {"p_rad": [np.nan]} + + # Get calibration factors + calib_factors = EastUtilMethods.get_axuv_calib_factors() + # Apply calibration factors Fac1 to Fac4 + for i in range(64): + xuv[:, i] *= ( + calib_factors["Fac1"][i % 16] + * calib_factors["Fac2"][i // 16][i % 16] + * calib_factors["Fac3"][i // 16] + * calib_factors["Fac4"] + ) + # Correction for bad channels (from Duan Yanmin's program) + xuv[:, 11] = 0.5 * (xuv[:, 10] + xuv[:, 12]) + xuv[:, 35] = 0.5 * (xuv[:, 34] + xuv[:, 35]) + # Apply geometric calibrations and Fac5 + for i in range(64): + xuv[:, i] *= ( + 2 + * np.pi + * calib_factors["Maj_R"] + * calib_factors["Del_r"][i] + * calib_factors["Fac5"] + ) + + # Select core measurements + core_indices = ( + list(range(7, 15)) + + list(range(20, 32)) + + list(range(32, 48)) + + list(range(53, 59)) + ) + p_rad = np.sum(xuv[:, core_indices], axis=1) # TODO: check nans + # Interpolate to requested time base + p_rad = interp1(xuvtime, p_rad, params.times, kind="linear", fill_value=0) + return {"p_rad": p_rad} + @staticmethod @physics_method( columns=[ - "p_rad", "p_ecrh", "p_lh", "p_icrf", @@ -584,12 +684,7 @@ def get_heating_power(nodes, tree): p_ohm = EastPhysicsMethods.get_p_ohm(params)["p_oh"] # [W] # Get radiated power - # tree = 'prad_east' - # TODO: find and implement `prad_bulk_xuv2014_2016` - # Original from Duan Yanming, and modified by RSG - # p_rad, rad_time = prad_bulk_xuv2014_2016(params.shot_id) - # p_rad = interp1(rad_time, p_rad, params.times, kind="linear", fill_value=0) - p_rad = [np.nan] * len(params.times) + p_rad = EastPhysicsMethods.get_radiated_power(params)["p_rad"] # [W] # Get Wmhd and calculate dWmhd_dt wmhd, efittime = params.mds_conn.get_data_with_dims( @@ -606,7 +701,6 @@ def get_heating_power(nodes, tree): rad_loss_frac = p_rad / (p_input - dwmhd_dt) output = { - "p_rad": p_rad, "p_ecrh": p_ecrh, "p_lh": p_lh, "p_icrf": p_icrf, @@ -1266,45 +1360,26 @@ def get_prad_peaking(params: PhysicsMethodParams): Last major update: 2014/11/22 by William Wei """ + xuv, xuvtime = EastPhysicsMethods.get_raw_axuv_data(params) - # The following section about calibration factors was copied-and-pasted - # directly from Duan Yanmin 's code. - - # TODO: Move calibration factors to a separate settings file - # Calibration factors + # Get calibration factors calib_factors = EastUtilMethods.get_axuv_calib_factors() - - # Get XUV data - (xuvtime,) = params.mds_conn.get_dims(r"\pxuv1", tree_name="east_1") - # There are 64 AXUV chords, arranged in 4 arrays of 16 channels each - xuv = np.full((len(xuvtime), 64), np.nan) - - dt = xuvtime[1] - xuvtime[0] - smoothing_time = 1e-3 - smoothing_window = int(max([round(smoothing_time / dt, 1)])) - for iarray in range(4): - for ichan in range(16): - ichord = 16 * iarray + ichan - # TODO: confirm the actual node of each chord - signal = params.mds_conn.get_data( - r"\pxuv" + str(ichord + 1), tree_name="east_1" - ) - signal = signal - np.mean(signal[:100]) # Subtract baseline - signal_smoothed = smooth(signal, smoothing_window) - xuv[:, ichord] = ( - signal_smoothed - * calib_factors["Fac1"][ichan] - * calib_factors["Fac2"][iarray][ichan] - * calib_factors["Fac3"][iarray] - * calib_factors["Fac4"] - * 2 - * np.pi - * calib_factors["Maj_R"] - * calib_factors["Del_r"][ichan] - ) # from Duan Yanming's program - + # Apply calibration factors + for i in range(64): + xuv[:, i] *= ( + calib_factors["Fac1"][i % 16] + * calib_factors["Fac2"][i // 16][i % 16] + * calib_factors["Fac3"][i // 16] + * calib_factors["Fac4"] + * 2 + * np.pi + * calib_factors["Maj_R"] + * calib_factors["Del_r"][i] + # * calib_factors["Fac5"] + ) # Correction for bad channels (from Duan Yanmin's program) xuv[:, 11] = 0.5 * (xuv[:, 10] + xuv[:, 12]) + # xuv[:, 35] = 0.5 * (xuv[:, 34] + xuv[:, 35]) # Define the core chords to be #28 to #37 (centermost 10 chords), and # define the non-divertor chords to be #09 to #56. (Chords #1-8 view the diff --git a/disruption_py/machine/east/util.py b/disruption_py/machine/east/util.py index 157f54fc8..1d2ea4729 100644 --- a/disruption_py/machine/east/util.py +++ b/disruption_py/machine/east/util.py @@ -6,6 +6,7 @@ import scipy from disruption_py.inout.mds import MDSConnection +from disruption_py.core.utils.math import smooth class EastUtilMethods: From e367f888576c738ca1cadf06bdbe52b32c1627ac Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:31:28 +0800 Subject: [PATCH 06/28] Add notes to highlight difference between get_prad_peaking and get_radiated_power --- disruption_py/machine/east/physics.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 91ba48ba9..eba10f2dd 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -518,6 +518,7 @@ def get_radiated_power(params: PhysicsMethodParams): * calib_factors["Fac4"] ) # Correction for bad channels (from Duan Yanmin's program) + # NOTE: why not [35] = ([34]+[36])/2 when [35] is the bad channel? xuv[:, 11] = 0.5 * (xuv[:, 10] + xuv[:, 12]) xuv[:, 35] = 0.5 * (xuv[:, 34] + xuv[:, 35]) # Apply geometric calibrations and Fac5 @@ -1366,6 +1367,7 @@ def get_prad_peaking(params: PhysicsMethodParams): calib_factors = EastUtilMethods.get_axuv_calib_factors() # Apply calibration factors for i in range(64): + # NOTE: Original script has Del_r(ichan) which should be a mistake xuv[:, i] *= ( calib_factors["Fac1"][i % 16] * calib_factors["Fac2"][i // 16][i % 16] @@ -1374,10 +1376,11 @@ def get_prad_peaking(params: PhysicsMethodParams): * 2 * np.pi * calib_factors["Maj_R"] - * calib_factors["Del_r"][i] + * calib_factors["Del_r"][i] # * calib_factors["Fac5"] ) # Correction for bad channels (from Duan Yanmin's program) + # NOTE: get_radiated_power does this first and then apply geometric correction. xuv[:, 11] = 0.5 * (xuv[:, 10] + xuv[:, 12]) # xuv[:, 35] = 0.5 * (xuv[:, 34] + xuv[:, 35]) From 25b8dc032d74be7a9368b216490731082a4d0977 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 19 Mar 2025 09:46:51 -0400 Subject: [PATCH 07/28] Update get_radiated_power comment --- disruption_py/machine/east/physics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index eba10f2dd..6b0e6ed2a 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -531,7 +531,7 @@ def get_radiated_power(params: PhysicsMethodParams): * calib_factors["Fac5"] ) - # Select core measurements + # Subtract divertor measurements and overlapping measurments core_indices = ( list(range(7, 15)) + list(range(20, 32)) @@ -1376,7 +1376,7 @@ def get_prad_peaking(params: PhysicsMethodParams): * 2 * np.pi * calib_factors["Maj_R"] - * calib_factors["Del_r"][i] + * calib_factors["Del_r"][i] # * calib_factors["Fac5"] ) # Correction for bad channels (from Duan Yanmin's program) From 884a001e03e8c094b1c6043522cf230222068a9e Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sat, 29 Mar 2025 03:42:08 +0800 Subject: [PATCH 08/28] Complete get_n_equal_1_data --- disruption_py/machine/east/config.toml | 5 - disruption_py/machine/east/physics.py | 34 +++--- disruption_py/machine/east/util.py | 156 ++++++++++++++++++++++++- 3 files changed, 171 insertions(+), 24 deletions(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index f6596da51..14d7d265b 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -33,9 +33,6 @@ expected_failure_columns = [ "n1rms_normalized", "n2rms", "n2rms_normalized", - "n_equal_1_mode", - "n_equal_1_normalized", - "n_equal_1_phase", "p_input", "p_oh", "p_rad", @@ -44,8 +41,6 @@ expected_failure_columns = [ "qstar", "rad_input_frac", "rad_loss_frac", - "rmp_n_equal_1", - "rmp_n_equal_1_phase", "wmhd", "wmhd_rt", ] diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 6b0e6ed2a..0e2806bd5 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -853,6 +853,7 @@ def get_n_equal_1_data(params: PhysicsMethodParams): rmp_n_equal_1_phase = [np.nan] # Get the rmp coil currents + # Translated from get_rmp_and_saddle_signals.m (rmptime,) = params.mds_conn.get_dims(r"\irmpu1", tree_name="east") rmp = np.full((len(rmptime), 16), np.nan) for i in range(8): @@ -874,11 +875,9 @@ def get_n_equal_1_data(params: PhysicsMethodParams): r"\sad_fg", r"\sad_hi", r"\sad_jk", - r"\sad_lk", - r"\sad_no", + r"\sad_lm", + # r"\sad_no", # not operational in 2015 ] - - # TODO: verify this! for i, node in enumerate(saddle_nodes[:7]): try: saddle[:, i] = params.mds_conn.get_data(node, tree_name="east") @@ -888,33 +887,34 @@ def get_n_equal_1_data(params: PhysicsMethodParams): sad_lo = params.mds_conn.get_data(r"\sad_lo", tree_name="east") sad_lm = params.mds_conn.get_data(r"\sad_lm", tree_name="east") saddle[:, 7] = sad_lo - sad_lm + # End of get_rmp_and_saddle_signals # Calculate RMP n=1 Fourier component amplitude and phase (on the timebase # of the saddle signals) - # TODO: Find 'rmp_sadle_coeff_matrix.mat' - # coeff_matrix = None - # rmp_pickup = np.transpose(coeff_matrix * np.transpose(rmp)) + coeff_matrix = EastUtilMethods.load_rmp_saddle_coeff_matrix() + rmp_pickup = np.transpose(np.matmul(coeff_matrix, np.transpose(rmp))) # Interpolate rmp_pickup onto the saddle time timebase - # rmp_pickup = interp1(rmptime, rmp_pickup, saddletime) + rmp_pickup = interp1(rmptime, rmp_pickup, saddletime, axis=0) # Calculate fast Fourier transforms of the RMP-induced signals, and get # mode amplitudes and phases # Take FFT along 2nd dimension (phi) - # rmp_fft_output = scipy.fft.fft(rmp_pickup, axis=1) - # amplitude = abs(rmp_fft_output) / rmp_fft_output.shape[1] - # amplitude[:, 1:] *= 2 # TODO: Why? - # phase = np.arctan2(np.imag(rmp_fft_output), np.real(rmp_fft_output)) + rmp_fft_output = scipy.fft.fft(rmp_pickup, axis=1) + amplitude = abs(rmp_fft_output) / len(rmp_fft_output[0]) + amplitude[:, 1:] *= 2 # TODO: Why? + phase = np.arctan2(np.imag(rmp_fft_output), np.real(rmp_fft_output)) # Only want n=1 Fourier component - # rmp_n_equal_1 = amplitude[:, 1] - # rmp_n_equal_1_phase = phase[:, 1] + rmp_n_equal_1 = amplitude[:, 1] + rmp_n_equal_1_phase = phase[:, 1] # Interpolate onto the requested timebase - # rmp_n_equal_1 = interp1(saddletime, rmp_n_equal_1, params.times) - # rmp_n_equal_1_phase = interp1(saddletime, rmp_n_equal_1_phase, params.times) + # TODO: interpolate phase may cause unexpected error + rmp_n_equal_1 = interp1(saddletime, rmp_n_equal_1, params.times) + rmp_n_equal_1_phase = interp1(saddletime, rmp_n_equal_1_phase, params.times) # The saddle signals can include direct pickup from the RMP coils, when the # RMP coils are active. This direct pickup must be subtracted from the # saddle signals to leave just the plasma response and baseline drifts. - # saddle -= rmp_pickup + saddle -= rmp_pickup # Calculate fast Fourier transforms and get mode amplitudes and phases # Take FFT along 2nd dimension (phi) saddle_fft_output = scipy.fft.fft(saddle, axis=1) diff --git a/disruption_py/machine/east/util.py b/disruption_py/machine/east/util.py index 1d2ea4729..6d61128ed 100644 --- a/disruption_py/machine/east/util.py +++ b/disruption_py/machine/east/util.py @@ -6,8 +6,6 @@ import scipy from disruption_py.inout.mds import MDSConnection -from disruption_py.core.utils.math import smooth - class EastUtilMethods: """ @@ -210,3 +208,157 @@ def get_axuv_calib_factors() -> dict: "Maj_R": 1.85, "Del_r": del_r, } + + @staticmethod + def load_rmp_saddle_coeff_matrix(): + """ + Load the rmp saddle coefficient matrix. + Source: rmp_saddle_coeff_matrix.mat + """ + coeff_matrix = [ + [ + -0.000325933950706964, + -5.34427695316486e-05, + -9.60845947944944e-06, + -5.04741218778060e-06, + -5.59732219661124e-06, + -5.91438750088332e-06, + -7.67820045826030e-06, + -5.20559984855693e-05, + -0.000302471162280883, + -5.35568562149143e-05, + -9.56803697152602e-06, + -5.88885912077860e-06, + -7.99784778858085e-06, + -6.80525414488422e-06, + -8.71081837989025e-06, + -5.76306168766297e-05, + ], + [ + -5.22709907162920e-05, + -0.000320133097404118, + -5.50071761972670e-05, + -7.37274061429450e-06, + -6.90578359290429e-06, + -5.02997540459485e-06, + -5.10697946272394e-06, + -8.99898260398151e-06, + -5.66360083538718e-05, + -0.000311386125033681, + -5.23092310075144e-05, + -8.27585100002495e-06, + -8.00023807246884e-06, + -6.44383893574099e-06, + -6.54598718347685e-06, + -1.00724151886658e-05, + ], + [ + -8.08294716760164e-06, + -5.55046044711059e-05, + -0.000333809266236893, + -5.60440355539882e-05, + -9.50511197149052e-06, + -5.91247287479164e-06, + -5.20599513677853e-06, + -6.16780683251725e-06, + -1.15430117601721e-05, + -5.51203237240259e-05, + -0.000316394172222160, + -5.48875912260367e-05, + -1.06537242928572e-05, + -6.47369008200614e-06, + -6.27085788458336e-06, + -6.55919579400570e-06, + ], + [ + -5.12806076473295e-06, + -8.49939605799587e-06, + -5.92476041414216e-05, + -0.000328925331027763, + -5.69684251594807e-05, + -8.70337525570129e-06, + -4.92573011225709e-06, + -7.82605756687085e-06, + -8.55723442676295e-06, + -8.33796262503122e-06, + -5.77746097768847e-05, + -0.000325606660148107, + -5.93152260417858e-05, + -8.66043101910145e-06, + -6.37549893634100e-06, + -6.85515607000913e-06, + ], + [ + -4.65502958467029e-06, + -4.85211383504998e-06, + -1.02244946378116e-05, + -4.93116559918267e-05, + -0.000332458254995670, + -5.51888967708278e-05, + -7.85611622708372e-06, + -6.92063279447296e-06, + -8.36497581142532e-06, + -5.98112870895933e-06, + -1.04474088707701e-05, + -5.44799256210205e-05, + -0.000326342886258557, + -5.27273498884329e-05, + -9.18704874796283e-06, + -6.78335745990761e-06, + ], + [ + -4.24256096566107e-06, + -4.36018675342360e-06, + -6.89150573826226e-06, + -6.63787474197891e-06, + -5.61987082604549e-05, + -0.000301238014687208, + -5.52453316551176e-05, + -9.13919157975075e-06, + -8.37854203786939e-06, + -5.93322139470238e-06, + -6.58104941723929e-06, + -8.33865462199638e-06, + -5.68131341626709e-05, + -0.000317650462941753, + -5.63542945554439e-05, + -9.17221707740544e-06, + ], + [ + -7.79176701405456e-06, + -4.90936998737812e-06, + -6.98520489721429e-06, + -4.40637734472976e-06, + -9.65641938226777e-06, + -5.33787503140721e-05, + -0.000318500478130682, + -5.72159746691176e-05, + -1.15055865602172e-05, + -5.31875434926984e-06, + -6.76325866931821e-06, + -5.68472690837507e-06, + -9.76815563125627e-06, + -5.75939865428040e-05, + -0.000325379112579523, + -5.60371003488067e-05, + ], + [ + -5.72759424001505e-05, + -8.13153861612048e-06, + -7.30722987804662e-06, + -5.64219477471097e-06, + -6.40689982720806e-06, + -1.06255695592301e-05, + -6.29862623902986e-05, + -0.000330999771429404, + -5.87929295332690e-05, + -7.60812324388172e-06, + -7.52054522920347e-06, + -5.45222900075003e-06, + -8.27237279295883e-06, + -9.78221302986545e-06, + -6.81254152056869e-05, + -0.000336854492397579, + ], + ] + return np.array(coeff_matrix) From 66cdad0d651ac36fc8675bca2779b241df890dd1 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sat, 29 Mar 2025 04:00:33 +0800 Subject: [PATCH 09/28] Remove p_rad, add prad_peaking to XFAIL columns --- disruption_py/machine/east/config.toml | 2 +- disruption_py/machine/east/util.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 14d7d265b..3576f1f74 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -35,7 +35,7 @@ expected_failure_columns = [ "n2rms_normalized", "p_input", "p_oh", - "p_rad", + "prad_peaking", "q0", "q95", "qstar", diff --git a/disruption_py/machine/east/util.py b/disruption_py/machine/east/util.py index 6d61128ed..33f08b79f 100644 --- a/disruption_py/machine/east/util.py +++ b/disruption_py/machine/east/util.py @@ -7,6 +7,7 @@ from disruption_py.inout.mds import MDSConnection + class EastUtilMethods: """ A class of helper methods that might fetch and compute data from MDSplus From 53432b8295c9a778eea5a938443504b6b24c394f Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sun, 30 Mar 2025 04:07:12 +0800 Subject: [PATCH 10/28] Fix get_n1rms_n2rms. Remove signals from test --- disruption_py/machine/east/config.toml | 8 -------- disruption_py/machine/east/physics.py | 18 +++++++++--------- 2 files changed, 9 insertions(+), 17 deletions(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 3576f1f74..09e93bf94 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -29,10 +29,6 @@ expected_failure_columns = [ "li_rt", "mirnov_std", "mirnov_std_normalized", - "n1rms", - "n1rms_normalized", - "n2rms", - "n2rms_normalized", "p_input", "p_oh", "prad_peaking", @@ -67,10 +63,6 @@ test_columns = [ "li_rt", "mirnov_std", "mirnov_std_normalized", - "n1rms", - "n1rms_normalized", - "n2rms", - "n2rms_normalized", "n_e", "n_equal_1_mode", "n_equal_1_normalized", diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 0e2806bd5..0fcb90e58 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -1515,16 +1515,16 @@ def get_n1rms_n2rms(params: PhysicsMethodParams): r"\mitde2", r"\mitef2", r"\mitfg2", - r"\mit2gh", - r"\mit2hi", + r"\mitgh2", + r"\mithi2", r"\mitij2", r"\mitjk2", r"\mitkl2", - r"\mit2lm", - r"\mit2mn", + r"\mitlm2", + r"\mitmn2", r"\mitno2", r"\mitop2", - r"\mit2pa", + r"\mitpa2", ] for i, node in enumerate(mir_nodes): try: @@ -1537,7 +1537,7 @@ def get_n1rms_n2rms(params: PhysicsMethodParams): # Compute the n=1 and n=2 mode signals from the Mirnov array signals # TODO: Verify the output structure of scipy.fft.fft; MATLAB one has index 3 (so 0, 1, 2?) - mir_fft_output = scipy.fft.fft(mir, axis=0) + mir_fft_output = scipy.fft.fft(mir, axis=1) amplitude = abs(mir_fft_output) / len(mir_fft_output[0]) amplitude[:, 1:] *= 2 # TODO: Why? phase = np.arctan2(np.imag(mir_fft_output), np.real(mir_fft_output)) @@ -1548,10 +1548,10 @@ def get_n1rms_n2rms(params: PhysicsMethodParams): time_window = 0.001 for i, time in enumerate(params.times): (indices,) = np.where( - (time - time_window <= mirtime) & (mirtime < time + time_window) + (time - time_window < mirtime) & (mirtime < time + time_window) ) - n1rms[i] = np.nanstd(n1[indices]) - n2rms[i] = np.nanstd(n2[indices]) + n1rms[i] = np.nanstd(n1[indices], ddof=1) + n2rms[i] = np.nanstd(n2[indices], ddof=1) # Calculate the normalized signals n1rms_normalized = n1rms / abs(btor) From 86a2133aeeac26e1b2ee4b07d0241b99655a8702 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sun, 30 Mar 2025 05:32:21 +0800 Subject: [PATCH 11/28] Add get_efit_gaps --- disruption_py/machine/east/config.toml | 2 + disruption_py/machine/east/physics.py | 115 +++++++++++++++++++++++++ 2 files changed, 117 insertions(+) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 09e93bf94..5da585f33 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -94,6 +94,8 @@ test_columns = [ "rmp_n_equal_1", "rmp_n_equal_1_phase", "time_until_disrupt", + "upper_gap", + "lower_gap", "v_loop", "wmhd", "wmhd_rt", diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 0fcb90e58..dcd01894a 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -1601,3 +1601,118 @@ def get_h98(params: PhysicsMethodParams): h98_y2 = interp1(h98_y2_time, h98_y2, params.times) return {"h98": h98_y2} + + @staticmethod + @physics_method( + columns=["upper_gap", "lower_gap"], + tokamak=Tokamak.EAST, + ) + def get_efit_gaps(params: PhysicsMethodParams): + """ + This script calculates upper and lower gaps from EFIT information + on plasma boundary and first wall geometry. This is needed because + the EAST version of EFIT does not output the upper and lower gaps + directly. + + Parameters + ---------- + params : PhysicsMethodParams + Parameters containing MDS connection and shot information + + Returns + ------- + dict + A dictionary containing the following keys: + - 'upper_gap' : array + Upper gap [m]. + - 'lower_gap' : array + Lower gap [m]. + + References + ------- + https://github.com/MIT-PSFC/disruption-py/blob/matlab/EAST/get_EFIT_gaps.m + + Original Authors + ---------------- + Robert Granetz + Jiaxiang Zhu + + Last major update: 2014/11/22 by William Wei + """ + upper_gap = [np.nan] + lower_gap = [np.nan] + + # Get plasma boundary data + data, efittime = params.mds_conn.get_data_with_dims( + r"\top.results.geqdsk:bdry", tree_name="_efit_tree" + ) + # Convert the order of indices to MATLAB order + # MATLAB (0, 1, 2) -> Python (2, 1, 0) + data = np.transpose(data, [2, 1, 0]) + xcoords = np.reshape(data[0, :, :], (-1, len(efittime))) + ycoords = np.reshape(data[1, :, :], (-1, len(efittime))) + + # Get first wall geometry data + xfirstwall = params.mds_conn.get_data( + r"\top.results.geqdsk:xlim", tree_name="_efit_tree" + ) + yfirstwall = params.mds_conn.get_data( + r"\top.results.geqdsk:ylim", tree_name="_efit_tree" + ) + seed = np.ones((len(xcoords), 1)) + xfirstwall = np.reshape(xfirstwall, (-1, 1)) + yfirstwall = np.reshape(yfirstwall, (-1, 1)) + xfirstwall_mat = np.tile(xfirstwall, (1, len(efittime))) + yfirstwall_mat = np.tile(yfirstwall, (1, len(efittime))) + + # Calculate upper & lower gaps + index_upperwall, _ = np.where(yfirstwall > 0.6) + index_lowerwall, _ = np.where(yfirstwall < -0.6) + + xupperwall = xfirstwall_mat[index_upperwall, :] + xupperwall_mat = np.reshape( + np.kron(xupperwall, seed), (len(seed), -1, len(efittime)) + ) + xlowerwall = xfirstwall_mat[index_lowerwall, :] + xlowerwall_mat = np.reshape( + np.kron(xlowerwall, seed), (len(seed), -1, len(efittime)) + ) + + yupperwall = yfirstwall_mat[index_upperwall, :] + yupperwall_mat = np.reshape( + np.kron(yupperwall, seed), (len(seed), -1, len(efittime)) + ) + ylowerwall = yfirstwall_mat[index_lowerwall, :] + ylowerwall_mat = np.reshape( + np.kron(ylowerwall, seed), (len(seed), -1, len(efittime)) + ) + + xupperplasma_mat = np.reshape( + np.tile(xcoords, (len(xupperwall), 1)), (-1, len(xupperwall), len(efittime)) + ) + yupperplasma_mat = np.reshape( + np.tile(ycoords, (len(yupperwall), 1)), (-1, len(yupperwall), len(efittime)) + ) + xlowerplasma_mat = np.reshape( + np.tile(xcoords, (len(xlowerwall), 1)), (-1, len(xlowerwall), len(efittime)) + ) + ylowerplasma_mat = np.reshape( + np.tile(ycoords, (len(ylowerwall), 1)), (-1, len(ylowerwall), len(efittime)) + ) + + uppergap_mat = np.sqrt( + np.square(xupperplasma_mat - xupperwall_mat) + + np.square(yupperplasma_mat - yupperwall_mat) + ) + lowergap_mat = np.sqrt( + np.square(xlowerplasma_mat - xlowerwall_mat) + + np.square(ylowerplasma_mat - ylowerwall_mat) + ) + upper_gap = uppergap_mat.min(axis=(0, 1)) + lower_gap = lowergap_mat.min(axis=(0, 1)) + + # Interpolate to the requested timebase + upper_gap = interp1(efittime, upper_gap, params.times) + lower_gap = interp1(efittime, lower_gap, params.times) + + return {"upper_gap": upper_gap, "lower_gap": lower_gap} From cca6419d867bc921f7e820469d7f25a6a2eb0a76 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sun, 30 Mar 2025 05:58:16 +0800 Subject: [PATCH 12/28] Split get_efit_gaps into public and hidden methods. Add pupper_gap and plower_gap to public method outputs --- disruption_py/machine/east/physics.py | 94 ++++++++++++++++----------- 1 file changed, 57 insertions(+), 37 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index dcd01894a..183b22bbb 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -1603,48 +1603,16 @@ def get_h98(params: PhysicsMethodParams): return {"h98": h98_y2} @staticmethod - @physics_method( - columns=["upper_gap", "lower_gap"], - tokamak=Tokamak.EAST, - ) - def get_efit_gaps(params: PhysicsMethodParams): + def _get_efit_gaps(params: PhysicsMethodParams, tree: str = "_efit_tree"): """ - This script calculates upper and lower gaps from EFIT information - on plasma boundary and first wall geometry. This is needed because - the EAST version of EFIT does not output the upper and lower gaps - directly. - - Parameters - ---------- - params : PhysicsMethodParams - Parameters containing MDS connection and shot information - - Returns - ------- - dict - A dictionary containing the following keys: - - 'upper_gap' : array - Upper gap [m]. - - 'lower_gap' : array - Lower gap [m]. - - References - ------- - https://github.com/MIT-PSFC/disruption-py/blob/matlab/EAST/get_EFIT_gaps.m - - Original Authors - ---------------- - Robert Granetz - Jiaxiang Zhu - - Last major update: 2014/11/22 by William Wei + Hidden method to calculate the EFIT and P-EFIT gaps """ upper_gap = [np.nan] lower_gap = [np.nan] # Get plasma boundary data data, efittime = params.mds_conn.get_data_with_dims( - r"\top.results.geqdsk:bdry", tree_name="_efit_tree" + r"\top.results.geqdsk:bdry", tree_name=tree ) # Convert the order of indices to MATLAB order # MATLAB (0, 1, 2) -> Python (2, 1, 0) @@ -1654,10 +1622,10 @@ def get_efit_gaps(params: PhysicsMethodParams): # Get first wall geometry data xfirstwall = params.mds_conn.get_data( - r"\top.results.geqdsk:xlim", tree_name="_efit_tree" + r"\top.results.geqdsk:xlim", tree_name=tree ) yfirstwall = params.mds_conn.get_data( - r"\top.results.geqdsk:ylim", tree_name="_efit_tree" + r"\top.results.geqdsk:ylim", tree_name=tree ) seed = np.ones((len(xcoords), 1)) xfirstwall = np.reshape(xfirstwall, (-1, 1)) @@ -1716,3 +1684,55 @@ def get_efit_gaps(params: PhysicsMethodParams): lower_gap = interp1(efittime, lower_gap, params.times) return {"upper_gap": upper_gap, "lower_gap": lower_gap} + + @staticmethod + @physics_method( + columns=["upper_gap", "lower_gap", "pupper_gap", "plower_gap"], + tokamak=Tokamak.EAST, + ) + def get_efit_gaps(params: PhysicsMethodParams): + """ + This script calculates upper and lower gaps from the EFIT and P-EFIT + information on plasma boundary and first wall geometry. This is needed + because the EAST version of EFIT does not output the upper and lower + gaps directly. + + Parameters + ---------- + params : PhysicsMethodParams + Parameters containing MDS connection and shot information + + Returns + ------- + dict + A dictionary containing the following keys: + - 'upper_gap' : array + Upper gap [m]. + - 'lower_gap' : array + Lower gap [m]. + - 'pupper_gap' : array + Upper gap computed from the P-EFIT tree [m]. + - 'plower_gap' : array + Lower gap computed from the P-EFIT tree [m]. + + References + ------- + https://github.com/MIT-PSFC/disruption-py/blob/matlab/EAST/get_EFIT_gaps.m + https://github.com/MIT-PSFC/disruption-py/blob/matlab/EAST/get_PEFIT_gaps.m + + Original Authors + ---------------- + Robert Granetz + Jiaxiang Zhu + + Last major update: 2015/03/29 by William Wei + """ + efit_gaps = EastPhysicsMethods._get_efit_gaps(params=params, tree="_efit_tree") + pefit_gaps = EastPhysicsMethods._get_efit_gaps(params=params, tree="pefit_east") + + return { + "upper_gap": efit_gaps["upper_gap"], + "lower_gap": efit_gaps["lower_gap"], + "pupper_gap": pefit_gaps["upper_gap"], + "plower_gap": pefit_gaps["lower_gap"], + } From 165abffdcd97d6303c142a05b703b5ca127232de Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sun, 30 Mar 2025 06:04:27 +0800 Subject: [PATCH 13/28] Remove efit signals that are not in SQL table from test --- disruption_py/machine/east/config.toml | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 5da585f33..0b787c745 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -12,15 +12,9 @@ port = 3306 [east.tests] expected_failure_columns = [ - "aminor", - "area", "beta_n", "beta_p", "beta_p_rt", - "chisq", - "dbetap_dt", - "dli_dt", - "dwmhd_dt", "greenwald_fraction", "ip_error_rt", "kappa", @@ -41,17 +35,11 @@ expected_failure_columns = [ "wmhd_rt", ] test_columns = [ - "aminor", - "area", "beta_n", "beta_p", "beta_p_rt", "btor", - "chisq", - "dbetap_dt", "dipprog_dt", - "dli_dt", - "dwmhd_dt", "greenwald_fraction", "ip", "ip_error", From 9dd9fbcc98479a8b7d52901209b424a3d92fffac Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Sun, 30 Mar 2025 06:16:20 +0800 Subject: [PATCH 14/28] Remove unnecessary newlines --- disruption_py/machine/east/physics.py | 35 ++++++--------------------- 1 file changed, 7 insertions(+), 28 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 183b22bbb..53f74c7c8 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -216,10 +216,7 @@ def get_ip_parameters(params: PhysicsMethodParams): } @staticmethod - @physics_method( - columns=["v_loop"], - tokamak=Tokamak.EAST, - ) + @physics_method(columns=["v_loop"], tokamak=Tokamak.EAST) def get_v_loop(params: PhysicsMethodParams): """ This routine gets the loop voltage from the EAST tree. The signal in the @@ -372,8 +369,7 @@ def get_z_error(params: PhysicsMethodParams): @staticmethod @physics_method( - columns=["n_e", "greenwald_fraction", "dn_dt"], - tokamak=Tokamak.EAST, + columns=["n_e", "greenwald_fraction", "dn_dt"], tokamak=Tokamak.EAST ) def get_density_parameters(params: PhysicsMethodParams): r""" @@ -713,10 +709,7 @@ def get_heating_power(nodes, tree): return output @staticmethod - @physics_method( - columns=["p_oh"], - tokamak=Tokamak.EAST, - ) + @physics_method(columns=["p_oh"], tokamak=Tokamak.EAST) def get_p_ohm(params: PhysicsMethodParams): r""" This script calculates the ohmic power, p_ohm. We use the following @@ -1170,14 +1163,7 @@ def get_pcs_parameters(params: PhysicsMethodParams): return output @staticmethod - @physics_method( - columns=[ - "p_rad_rt", - "p_lh_rt", - "p_nbi_rt", - ], - tokamak=Tokamak.EAST, - ) + @physics_method(columns=["p_rad_rt", "p_lh_rt", "p_nbi_rt"], tokamak=Tokamak.EAST) def get_pcs_power(params: PhysicsMethodParams): """ This function gets the real time power signals that are actually used @@ -1320,10 +1306,7 @@ def get_pcs_power(params: PhysicsMethodParams): } @staticmethod - @physics_method( - columns=["prad_peaking"], - tokamak=Tokamak.EAST, - ) + @physics_method(columns=["prad_peaking"], tokamak=Tokamak.EAST) def get_prad_peaking(params: PhysicsMethodParams): """ This routine calculates the peaking factor of the profiles of radiated @@ -1406,8 +1389,7 @@ def get_prad_peaking(params: PhysicsMethodParams): @staticmethod @physics_method( - columns=["mirnov_std", "mirnov_std_normalized"], - tokamak=Tokamak.EAST, + columns=["mirnov_std", "mirnov_std_normalized"], tokamak=Tokamak.EAST ) def get_mirnov_std(params: PhysicsMethodParams): """ @@ -1565,10 +1547,7 @@ def get_n1rms_n2rms(params: PhysicsMethodParams): } @staticmethod - @physics_method( - columns=["h98"], - tokamak=Tokamak.EAST, - ) + @physics_method(columns=["h98"], tokamak=Tokamak.EAST) def get_h98(params: PhysicsMethodParams): """ Get the H98y2 energy confinement time parameter. From 4a11cd5454e3635b49a385696b07480b672c8581 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Mon, 31 Mar 2025 23:47:39 +0800 Subject: [PATCH 15/28] minor corrections --- disruption_py/machine/east/physics.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 53f74c7c8..a47de563a 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -554,15 +554,15 @@ def get_radiated_power(params: PhysicsMethodParams): ) def get_power(params: PhysicsMethodParams): """ - This function gets the input heating powers -- ohmic (p_OHM),electron - cyclotron resonance heating (p_ECRH), neutral beam injection system (p_NBI) - ion cyclotron (p_ICRF), and lower hybrid (p_LH) -- as well as the radiated - output power (p_RAD). If any of the auxiliary heating powers are not - available (there was no ICRF or LH), then this function returns an array - of zeros for them. The ohmic heating power is obtained by calling a - separate function, get_P_ohm.m. If either the ohmic heating power or the - radiated power is unavailable, then arrays of NaN (Not-a-Number) are - returned for them, and for the radiated power fraction. + This function gets the axuillary heating powers -- electron cyclotron + resonance heating (p_ECRH), neutral beam injection system (p_NBI) + ion cyclotron (p_ICRF), and lower hybrid (p_LH). If any of the auxiliary + heating powers are not available (there was no ICRF or LH), then this + function returns an array of zeros for them. + + This function also calculates the total input power, the radiated input + fraction, and the radiated loss function by calling get_p_ohm and + get_radiated_power. Parameters ---------- @@ -573,8 +573,6 @@ def get_power(params: PhysicsMethodParams): ------- dict A dictionary containing the following keys: - - 'p_rad' : array - Radiated power [W]. - 'p_ecrh' : array Electron cyclotron resonance heating power [W]. - 'p_lh' : array @@ -602,10 +600,8 @@ def get_power(params: PhysicsMethodParams): Last major update: 11/20/24 by William Wei """ - p_rad = [np.nan] p_ecrh = [np.nan] p_lh = [np.nan] - p_ohm = [np.nan] p_icrf = [np.nan] p_nbi = [np.nan] rad_input_frac = [np.nan] @@ -764,6 +760,7 @@ def get_p_ohm(params: PhysicsMethodParams): ) # [A] sign_ip = np.sign(sum(ip)) dipdt = np.gradient(ip, ip_time) + # TODO: switch to causal boxcar smoothing dipdt_smoothed = smooth(dipdt, 11) # Use 11-point boxcar smoothing # Interpolate fetched signals to the requested timebase @@ -775,7 +772,6 @@ def get_p_ohm(params: PhysicsMethodParams): ) # Calculate p_ohm - # TODO: check DIV/0 error # TODO: replace 1.85 with fetched r0 signal inductance = 4e-7 * np.pi * 1.85 * li / 2 # For EAST use R0 = 1.85 m v_inductive = -inductance * dipdt_smoothed @@ -1596,8 +1592,7 @@ def _get_efit_gaps(params: PhysicsMethodParams, tree: str = "_efit_tree"): # Convert the order of indices to MATLAB order # MATLAB (0, 1, 2) -> Python (2, 1, 0) data = np.transpose(data, [2, 1, 0]) - xcoords = np.reshape(data[0, :, :], (-1, len(efittime))) - ycoords = np.reshape(data[1, :, :], (-1, len(efittime))) + xcoords, ycoords = data[0, :, :], data[1, :, :] # Get first wall geometry data xfirstwall = params.mds_conn.get_data( From 4be3a4007ec6952e895baafc8e6ea6ada88509c5 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 2 Apr 2025 02:35:06 +0800 Subject: [PATCH 16/28] Simplify unpacking --- disruption_py/machine/east/physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index a47de563a..d08802fc0 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -1592,7 +1592,7 @@ def _get_efit_gaps(params: PhysicsMethodParams, tree: str = "_efit_tree"): # Convert the order of indices to MATLAB order # MATLAB (0, 1, 2) -> Python (2, 1, 0) data = np.transpose(data, [2, 1, 0]) - xcoords, ycoords = data[0, :, :], data[1, :, :] + xcoords, ycoords = data # Get first wall geometry data xfirstwall = params.mds_conn.get_data( From 49b06ae2dbe504eafa2a4de2c63db852a24bc236 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 2 Apr 2025 02:36:03 +0800 Subject: [PATCH 17/28] fix typo --- disruption_py/machine/east/physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index d08802fc0..acb6bd0b7 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -554,7 +554,7 @@ def get_radiated_power(params: PhysicsMethodParams): ) def get_power(params: PhysicsMethodParams): """ - This function gets the axuillary heating powers -- electron cyclotron + This function gets the auxiliary heating powers -- electron cyclotron resonance heating (p_ECRH), neutral beam injection system (p_NBI) ion cyclotron (p_ICRF), and lower hybrid (p_LH). If any of the auxiliary heating powers are not available (there was no ICRF or LH), then this From 2efa36dab6ad7470316cbf54e69bb28b15975e55 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 2 Apr 2025 02:39:58 +0800 Subject: [PATCH 18/28] Add fallback option for get_v_loop --- disruption_py/machine/east/physics.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index acb6bd0b7..bcbcf0ac6 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -218,14 +218,17 @@ def get_ip_parameters(params: PhysicsMethodParams): @staticmethod @physics_method(columns=["v_loop"], tokamak=Tokamak.EAST) def get_v_loop(params: PhysicsMethodParams): - """ - This routine gets the loop voltage from the EAST tree. The signal in the + r""" + This routine gets the loop voltage (\vp1_s) from the EAST tree. The signal in the tree is derived by taking the time derivative of a flux loop near the inboard midplane. Two possible signals are available, one digitised at a high rate (50 kHz), and the other sub-sampled down to 1 kHz. This routine reads in the 1 kHz signal. It linearly interpolates the loop voltage signal onto the specified timebase. + If \vp1_s isn't available, the method will fall back to using \pcvloop from the + pcs_east tree instead. + Parameters ---------- params : PhysicsMethodParams @@ -250,10 +253,17 @@ def get_v_loop(params: PhysicsMethodParams): # Get "\vp1_s" signal from the EAST tree. (This signal is a sub-sampled # version of "vp1".) - # TODO: Fallback to \pcvloop when \vp1_s isn't available (e.g. shot 81548)? - v_loop, v_loop_time = params.mds_conn.get_data_with_dims( - r"\vp1_s", tree_name="east" - ) + try: + v_loop, v_loop_time = params.mds_conn.get_data_with_dims( + r"\vp1_s", tree_name="east" + ) + except mdsExceptions.TreeException: + params.logger.verbose( + r"v_loop: Failed to get \vp1_s data. Use \pcvloop from pcs_east instead." + ) + v_loop, v_loop_time = params.mds_conn.get_data_with_dims( + r"\pcvloop", tree_name="pcs_east" + ) # [V] # Interpolate the signal onto the requested timebase v_loop = interp1(v_loop_time, v_loop, params.times) From ba4c56492a800985becd2e2fef811a96883cf33b Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 2 Apr 2025 02:41:12 +0800 Subject: [PATCH 19/28] Change get_raw_axuv_data to a private method --- disruption_py/machine/east/physics.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index bcbcf0ac6..ad1a6d874 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -448,7 +448,7 @@ def get_density_parameters(params: PhysicsMethodParams): return {"n_e": ne, "greenwald_fraction": greenwald_fraction, "dn_dt": dn_dt} @staticmethod - def get_raw_axuv_data(params: PhysicsMethodParams): + def _get_raw_axuv_data(params: PhysicsMethodParams): """ Get the raw (uncalibrated) data from the AXUV arrays for calculating the total radiated power and prad_peaking. The AXUV diagnostic contains @@ -458,8 +458,6 @@ def get_raw_axuv_data(params: PhysicsMethodParams): get_prad_peaking.m are slightly different in where the calibration factors are applied. This difference isn't explained in the scripts so I keep these functions different. - - TODO: should this function be put in EastUtilMethods? """ # Get XUV data (xuvtime,) = params.mds_conn.get_dims(r"\pxuv1", tree_name="east_1") # [s] @@ -509,7 +507,7 @@ def get_radiated_power(params: PhysicsMethodParams): """ # Get the raw AXUV data try: - xuv, xuvtime = EastPhysicsMethods.get_raw_axuv_data(params) # [W], [s] + xuv, xuvtime = EastPhysicsMethods._get_raw_axuv_data(params) # [W], [s] except mdsExceptions.MdsException: return {"p_rad": [np.nan]} @@ -1350,7 +1348,7 @@ def get_prad_peaking(params: PhysicsMethodParams): Last major update: 2014/11/22 by William Wei """ - xuv, xuvtime = EastPhysicsMethods.get_raw_axuv_data(params) + xuv, xuvtime = EastPhysicsMethods._get_raw_axuv_data(params) # Get calibration factors calib_factors = EastUtilMethods.get_axuv_calib_factors() From 3fca4c008235e40777766828cf0f6b703f2637cd Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 2 Apr 2025 02:48:26 +0800 Subject: [PATCH 20/28] Remove some TODOs & NOTEs --- disruption_py/machine/east/physics.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index ad1a6d874..b93faa2e6 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -471,7 +471,6 @@ def _get_raw_axuv_data(params: PhysicsMethodParams): for iarray in range(4): for ichan in range(16): ichord = 16 * iarray + ichan - # TODO: confirm the actual node of each chord signal = params.mds_conn.get_data( r"\pxuv" + str(ichord + 1), tree_name="east_1" ) @@ -502,8 +501,6 @@ def get_radiated_power(params: PhysicsMethodParams): References ------- - Prad_bulk_xuv2014_2016.m (Not in repo) - - TODO: Update docstring """ # Get the raw AXUV data try: @@ -522,7 +519,7 @@ def get_radiated_power(params: PhysicsMethodParams): * calib_factors["Fac4"] ) # Correction for bad channels (from Duan Yanmin's program) - # NOTE: why not [35] = ([34]+[36])/2 when [35] is the bad channel? + # TODO: why not [35] = ([34]+[36])/2 when [35] is the bad channel? xuv[:, 11] = 0.5 * (xuv[:, 10] + xuv[:, 12]) xuv[:, 35] = 0.5 * (xuv[:, 34] + xuv[:, 35]) # Apply geometric calibrations and Fac5 @@ -542,7 +539,7 @@ def get_radiated_power(params: PhysicsMethodParams): + list(range(32, 48)) + list(range(53, 59)) ) - p_rad = np.sum(xuv[:, core_indices], axis=1) # TODO: check nans + p_rad = np.nansum(xuv[:, core_indices], axis=1) # Interpolate to requested time base p_rad = interp1(xuvtime, p_rad, params.times, kind="linear", fill_value=0) return {"p_rad": p_rad} @@ -696,10 +693,9 @@ def get_heating_power(nodes, tree): # Calculate p_input, rad_input_frac, and rad_loss_frac p_input = p_ohm + p_lh + p_icrf + p_ecrh + p_nbi - rad_input_frac = ( - p_rad / p_input - ) # TODO: div/0 error? Use with np.errorstate(...) - rad_loss_frac = p_rad / (p_input - dwmhd_dt) + with np.errstate(divide="ignore", invalid="ignore"): + rad_input_frac = p_rad / p_input + rad_loss_frac = p_rad / (p_input - dwmhd_dt) output = { "p_ecrh": p_ecrh, @@ -904,7 +900,7 @@ def get_n_equal_1_data(params: PhysicsMethodParams): rmp_n_equal_1 = amplitude[:, 1] rmp_n_equal_1_phase = phase[:, 1] # Interpolate onto the requested timebase - # TODO: interpolate phase may cause unexpected error + # TODO: figure out how to interpolate phase rmp_n_equal_1 = interp1(saddletime, rmp_n_equal_1, params.times) rmp_n_equal_1_phase = interp1(saddletime, rmp_n_equal_1_phase, params.times) @@ -1144,7 +1140,6 @@ def get_pcs_parameters(params: PhysicsMethodParams): timearray, signal, params.times, kind="linear", bounds_error=0 ) output[name] = signal - # TODO: Specify error type except mdsExceptions.MdsException: output[name] = [np.nan] @@ -1160,7 +1155,6 @@ def get_pcs_parameters(params: PhysicsMethodParams): q95_rt_time, q95_rt, params.times, kind="linear", bounds_error=0 ) output["q95_rt"] = q95_rt - # TODO: Specify error type except mdsExceptions.MdsException: output["q95_rt"] = [np.nan] @@ -1354,7 +1348,7 @@ def get_prad_peaking(params: PhysicsMethodParams): calib_factors = EastUtilMethods.get_axuv_calib_factors() # Apply calibration factors for i in range(64): - # NOTE: Original script has Del_r(ichan) which should be a mistake + # Original script has Del_r(ichan) which should be a mistake xuv[:, i] *= ( calib_factors["Fac1"][i % 16] * calib_factors["Fac2"][i // 16][i % 16] @@ -1367,7 +1361,7 @@ def get_prad_peaking(params: PhysicsMethodParams): # * calib_factors["Fac5"] ) # Correction for bad channels (from Duan Yanmin's program) - # NOTE: get_radiated_power does this first and then apply geometric correction. + # get_radiated_power does these 2 lines first and then apply geometric correction. xuv[:, 11] = 0.5 * (xuv[:, 10] + xuv[:, 12]) # xuv[:, 35] = 0.5 * (xuv[:, 34] + xuv[:, 35]) @@ -1522,7 +1516,6 @@ def get_n1rms_n2rms(params: PhysicsMethodParams): btor = EastPhysicsMethods.get_btor(params)["btor"] # Compute the n=1 and n=2 mode signals from the Mirnov array signals - # TODO: Verify the output structure of scipy.fft.fft; MATLAB one has index 3 (so 0, 1, 2?) mir_fft_output = scipy.fft.fft(mir, axis=1) amplitude = abs(mir_fft_output) / len(mir_fft_output[0]) amplitude[:, 1:] *= 2 # TODO: Why? From f98aa7a02db18241530a6ab01954570d6100f6b4 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Wed, 2 Apr 2025 02:56:12 +0800 Subject: [PATCH 21/28] Change TreeException to MdsException in get_v_loop --- disruption_py/machine/east/physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index b93faa2e6..1dccb8ee6 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -257,7 +257,7 @@ def get_v_loop(params: PhysicsMethodParams): v_loop, v_loop_time = params.mds_conn.get_data_with_dims( r"\vp1_s", tree_name="east" ) - except mdsExceptions.TreeException: + except mdsExceptions.MdsException: params.logger.verbose( r"v_loop: Failed to get \vp1_s data. Use \pcvloop from pcs_east instead." ) From c0c79af8da3a4343ae96bf551bce174d84ed5242 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Mon, 14 Apr 2025 22:32:33 +0800 Subject: [PATCH 22/28] Simplify rad_input_frac & rad_loss_frac computations --- disruption_py/machine/east/physics.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 1dccb8ee6..e1a045e56 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -693,9 +693,15 @@ def get_heating_power(nodes, tree): # Calculate p_input, rad_input_frac, and rad_loss_frac p_input = p_ohm + p_lh + p_icrf + p_ecrh + p_nbi - with np.errstate(divide="ignore", invalid="ignore"): - rad_input_frac = p_rad / p_input - rad_loss_frac = p_rad / (p_input - dwmhd_dt) + p_loss = p_input - dwmhd_dt + + rad_input_frac = np.full(len(params.times), np.nan) + (valid_indices,) = np.where((p_input != 0) & ~np.isnan(p_input)) + rad_input_frac[valid_indices] = p_rad[valid_indices] / p_input[valid_indices] + + rad_loss_frac = np.full(len(params.times), np.nan) + (valid_indices,) = np.where((p_loss != 0) & ~np.isnan(p_loss)) + rad_loss_frac[valid_indices] = p_rad[valid_indices] / p_loss[valid_indices] output = { "p_ecrh": p_ecrh, From e14558ab5004fea9a5349dd4dd9be2ac02d41faa Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Mon, 14 Apr 2025 22:46:58 +0800 Subject: [PATCH 23/28] Clean up comments --- disruption_py/machine/east/physics.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index e1a045e56..a85d363b8 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -882,11 +882,9 @@ def get_n_equal_1_data(params: PhysicsMethodParams): saddle[:, i] = params.mds_conn.get_data(node, tree_name="east") except mdsExceptions.MdsException: saddle[:, i] = 0 - # \sad_no not operational in 2015 sad_lo = params.mds_conn.get_data(r"\sad_lo", tree_name="east") sad_lm = params.mds_conn.get_data(r"\sad_lm", tree_name="east") saddle[:, 7] = sad_lo - sad_lm - # End of get_rmp_and_saddle_signals # Calculate RMP n=1 Fourier component amplitude and phase (on the timebase # of the saddle signals) From 06eadc2bfd8e21d15b2cf32ec5d7317e6059a60e Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Mon, 14 Apr 2025 23:17:06 +0800 Subject: [PATCH 24/28] Fix bug introduced by github --- disruption_py/machine/east/config.toml | 2 +- disruption_py/machine/east/physics.py | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 0b787c745..35aa65279 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -49,6 +49,7 @@ test_columns = [ "kappa_area", "li", "li_rt", + "lower_gap", "mirnov_std", "mirnov_std_normalized", "n_e", @@ -83,7 +84,6 @@ test_columns = [ "rmp_n_equal_1_phase", "time_until_disrupt", "upper_gap", - "lower_gap", "v_loop", "wmhd", "wmhd_rt", diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index c43028a80..eb114cc41 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -477,7 +477,9 @@ def _get_raw_axuv_data(params: PhysicsMethodParams): # Subtract baseline signal = signal - np.mean(signal[:100]) # TODO: change this to causal smoothing - xuv[:, ichord] = smooth(signal, smoothing_window) * 1e3 # [kW] -> [W] + xuv[:, ichord] = ( + matlab_smooth(signal, smoothing_window) * 1e3 + ) # [kW] -> [W] return xuv, xuvtime @staticmethod From 7b2f50b599337beac54ad4add0ea4a4a467d58bb Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 6 May 2025 01:52:02 +0800 Subject: [PATCH 25/28] fail-proof get_power --- disruption_py/machine/east/physics.py | 125 ++++++++++++++++++-------- 1 file changed, 90 insertions(+), 35 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index eb114cc41..b2336fb34 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -571,6 +571,20 @@ def get_power(params: PhysicsMethodParams): fraction, and the radiated loss function by calling get_p_ohm and get_radiated_power. + Notes + ---------- + For the moment we assume that the MDSplus trees always have the corresponding + tree nodes for each of the powers, even if the heating source or diagnostic + was not turned on in that shot. If the method is not able to access a particular + node, then we assume the tree is broken and skips computing that power, rather than + considering that heating source or diagnostic was turned off during that shot and assign + an array of 0 to that power. This was not the case in several shots on CMOD, and we need + to verify this on EAST in the future. + + For now, if we encounter an MdsException error for any of the powers, we return + that power with nans and skip calculating the radiated fractions and p_input + if applicable. + Parameters ---------- params : PhysicsMethodParams @@ -649,17 +663,27 @@ def get_heating_power(nodes, tree): r"\plhi1*1e3", # LHW 2.45 GHz data r"\plhi2*1e3", # LHW 4.66 GHz data ] - p_lh = get_heating_power(lh_nodes, "east_1") # [W] + try: + p_lh = get_heating_power(lh_nodes, "east_1") # [W] + except mdsExceptions.MdsException: + params.logger.warning("Failed to get LH heating power") + p_lh = [np.nan] # Get ECRH power - p_ecrh, ecrh_time = params.mds_conn.get_data_with_dims( - r"\pecrh1i*1e3", tree_name="analysis" - ) # [W], [s] - (baseline_indices,) = np.where(ecrh_time < 0) - if len(baseline_indices) > 0: - p_ecrh_baseline = np.mean(p_ecrh[baseline_indices]) - p_ecrh -= p_ecrh_baseline - p_ecrh = interp1(ecrh_time, p_ecrh, params.times, kind="linear", fill_value=0) + try: + p_ecrh, ecrh_time = params.mds_conn.get_data_with_dims( + r"\pecrh1i*1e3", tree_name="analysis" + ) # [W], [s] + (baseline_indices,) = np.where(ecrh_time < 0) + if len(baseline_indices) > 0: + p_ecrh_baseline = np.mean(p_ecrh[baseline_indices]) + p_ecrh -= p_ecrh_baseline + p_ecrh = interp1( + ecrh_time, p_ecrh, params.times, kind="linear", fill_value=0 + ) + except mdsExceptions.MdsException: + params.logger.warning("Failed to get ECRH heating power") + p_ecrh = [np.nan] # Get ICRF power # p_ICRF = (p_icrfii - p_icrfir) + (p_icrfbi - p_icrfbr) @@ -669,7 +693,11 @@ def get_heating_power(nodes, tree): r"\picrfbi*1e3", r"\picrfbr*-1e3", ] - p_icrf = get_heating_power(icrf_nodes, "icrf_east") # [W] + try: + p_icrf = get_heating_power(icrf_nodes, "icrf_east") # [W] + except mdsExceptions.MdsException: + params.logger.warning("Failed to get ICRF heating power") + p_icrf = [np.nan] # Get NBI power nbi_nodes = [ @@ -678,7 +706,11 @@ def get_heating_power(nodes, tree): r"\pnbi2lsource*1e3", r"\pnbi2rsource*1e3", ] - p_nbi = get_heating_power(nbi_nodes, "nbi_east") # [W] + try: + p_nbi = get_heating_power(nbi_nodes, "nbi_east") # [W] + except mdsExceptions.MdsException: + params.logger.warning("Failed to get NBI heating power") + p_nbi = [np.nan] # Get ohmic power p_ohm = EastPhysicsMethods.get_p_ohm(params)["p_oh"] # [W] @@ -686,24 +718,43 @@ def get_heating_power(nodes, tree): # Get radiated power p_rad = EastPhysicsMethods.get_radiated_power(params)["p_rad"] # [W] - # Get Wmhd and calculate dWmhd_dt - wmhd, efittime = params.mds_conn.get_data_with_dims( - r"\efit_aeqdsk:wplasm", tree_name="_efit_tree" - ) # [W], [s] - dwmhd_dt = np.gradient(wmhd, efittime) - dwmhd_dt = interp1(efittime, dwmhd_dt, params.times) + # Calculate p_input + if not ( + np.isnan(p_ohm).all() + and np.isnan(p_lh).all() + and np.isnan(p_icrf).all() + and np.isnan(p_ecrh).all() + and np.isnan(p_nbi).all() + ): + p_input = p_ohm + p_lh + p_icrf + p_ecrh + p_nbi + else: + p_input = [np.nan] - # Calculate p_input, rad_input_frac, and rad_loss_frac - p_input = p_ohm + p_lh + p_icrf + p_ecrh + p_nbi - p_loss = p_input - dwmhd_dt + # Get Wmhd, calculate dWmhd_dt, and calculate p_loss + try: + wmhd, efittime = params.mds_conn.get_data_with_dims( + r"\efit_aeqdsk:wplasm", tree_name="_efit_tree" + ) # [W], [s] + dwmhd_dt = np.gradient(wmhd, efittime) + dwmhd_dt = interp1(efittime, dwmhd_dt, params.times) + if not np.isnan(p_input).all(): + p_loss = p_input - dwmhd_dt + else: + p_loss = [np.nan] + except mdsExceptions.MdsException: + params.logger.warning("Failed to get proper signals to compute p_loss") + p_loss = [np.nan] rad_input_frac = np.full(len(params.times), np.nan) - (valid_indices,) = np.where((p_input != 0) & ~np.isnan(p_input)) - rad_input_frac[valid_indices] = p_rad[valid_indices] / p_input[valid_indices] - rad_loss_frac = np.full(len(params.times), np.nan) - (valid_indices,) = np.where((p_loss != 0) & ~np.isnan(p_loss)) - rad_loss_frac[valid_indices] = p_rad[valid_indices] / p_loss[valid_indices] + if ~(np.isnan(p_rad).all() and np.isnan(p_input).all()): + (valid_indices,) = np.where((p_input != 0) & ~np.isnan(p_input)) + rad_input_frac[valid_indices] = ( + p_rad[valid_indices] / p_input[valid_indices] + ) + if ~(np.isnan(p_rad).all() and np.isnan(p_loss).all()): + (valid_indices,) = np.where((p_loss != 0) & ~np.isnan(p_loss)) + rad_loss_frac[valid_indices] = p_rad[valid_indices] / p_loss[valid_indices] output = { "p_ecrh": p_ecrh, @@ -760,16 +811,20 @@ def get_p_ohm(params: PhysicsMethodParams): Last major update: 2014/11/21 by William Wei """ # Get raw signals - vloop, vloop_time = params.mds_conn.get_data_with_dims( - r"\pcvloop", tree_name="pcs_east" - ) # [V] - li, li_time = params.mds_conn.get_data_with_dims( - r"\efit_aeqdsk:li", tree_name="_efit_tree" - ) # [H] - # Fetch raw ip signal to calculate dip_dt and apply smoothing - ip, ip_time = params.mds_conn.get_data_with_dims( - r"\pcrl01", tree_name="pcs_east" - ) # [A] + try: + vloop, vloop_time = params.mds_conn.get_data_with_dims( + r"\pcvloop", tree_name="pcs_east" + ) # [V] + li, li_time = params.mds_conn.get_data_with_dims( + r"\efit_aeqdsk:li", tree_name="_efit_tree" + ) # [H] + # Fetch raw ip signal to calculate dip_dt and apply smoothing + ip, ip_time = params.mds_conn.get_data_with_dims( + r"\pcrl01", tree_name="pcs_east" + ) # [A] + except mdsExceptions.MdsException: + params.logger.warning("Failed to get necessary signals to compute p_ohm") + return {"p_oh": [np.nan]} sign_ip = np.sign(sum(ip)) dipdt = np.gradient(ip, ip_time) # TODO: switch to causal boxcar smoothing From 0af1c92f8ea2285db20b759336eb73bad2fdb1dd Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 6 May 2025 02:57:59 +0800 Subject: [PATCH 26/28] Fix logics, add additional logger --- disruption_py/machine/east/physics.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index b2336fb34..864b7292d 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -508,6 +508,7 @@ def get_radiated_power(params: PhysicsMethodParams): try: xuv, xuvtime = EastPhysicsMethods._get_raw_axuv_data(params) # [W], [s] except mdsExceptions.MdsException: + params.logger.warning("Failed to get raw AXUS data") return {"p_rad": [np.nan]} # Get calibration factors @@ -719,12 +720,12 @@ def get_heating_power(nodes, tree): p_rad = EastPhysicsMethods.get_radiated_power(params)["p_rad"] # [W] # Calculate p_input - if not ( + if ~( np.isnan(p_ohm).all() - and np.isnan(p_lh).all() - and np.isnan(p_icrf).all() - and np.isnan(p_ecrh).all() - and np.isnan(p_nbi).all() + or np.isnan(p_lh).all() + or np.isnan(p_icrf).all() + or np.isnan(p_ecrh).all() + or np.isnan(p_nbi).all() ): p_input = p_ohm + p_lh + p_icrf + p_ecrh + p_nbi else: @@ -747,12 +748,12 @@ def get_heating_power(nodes, tree): rad_input_frac = np.full(len(params.times), np.nan) rad_loss_frac = np.full(len(params.times), np.nan) - if ~(np.isnan(p_rad).all() and np.isnan(p_input).all()): + if ~(np.isnan(p_rad).all() or np.isnan(p_input).all()): (valid_indices,) = np.where((p_input != 0) & ~np.isnan(p_input)) rad_input_frac[valid_indices] = ( p_rad[valid_indices] / p_input[valid_indices] ) - if ~(np.isnan(p_rad).all() and np.isnan(p_loss).all()): + if ~(np.isnan(p_rad).all() or np.isnan(p_loss).all()): (valid_indices,) = np.where((p_loss != 0) & ~np.isnan(p_loss)) rad_loss_frac[valid_indices] = p_rad[valid_indices] / p_loss[valid_indices] From e86d6cbc60e09604add80693085e5fb0e65ede85 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 6 May 2025 03:02:15 +0800 Subject: [PATCH 27/28] Chamge wmhd source from wplasm (all 0) to wmhd --- disruption_py/machine/east/physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/disruption_py/machine/east/physics.py b/disruption_py/machine/east/physics.py index 864b7292d..b8717bae8 100644 --- a/disruption_py/machine/east/physics.py +++ b/disruption_py/machine/east/physics.py @@ -734,7 +734,7 @@ def get_heating_power(nodes, tree): # Get Wmhd, calculate dWmhd_dt, and calculate p_loss try: wmhd, efittime = params.mds_conn.get_data_with_dims( - r"\efit_aeqdsk:wplasm", tree_name="_efit_tree" + r"\efit_aeqdsk:wmhd", tree_name="_efit_tree" ) # [W], [s] dwmhd_dt = np.gradient(wmhd, efittime) dwmhd_dt = interp1(efittime, dwmhd_dt, params.times) From ce8014381cf8394037e0607426477723c7ad96b5 Mon Sep 17 00:00:00 2001 From: yumouwei <46117079+yumouwei@users.noreply.github.com> Date: Tue, 6 May 2025 03:04:00 +0800 Subject: [PATCH 28/28] Add XFAILs due to east/east_1 tree issues --- disruption_py/machine/east/config.toml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/disruption_py/machine/east/config.toml b/disruption_py/machine/east/config.toml index 35aa65279..3f3ce7ac4 100644 --- a/disruption_py/machine/east/config.toml +++ b/disruption_py/machine/east/config.toml @@ -23,14 +23,22 @@ expected_failure_columns = [ "li_rt", "mirnov_std", "mirnov_std_normalized", + "n_equal_1_mode", # east/east_1 tree issues + "n_equal_1_normalized", # east/east_1 tree issues + "n_equal_1_phase", # east/east_1 tree issues "p_input", + "p_lh", # east/east_1 tree issues "p_oh", + "p_rad", # east/east_1 tree issues "prad_peaking", "q0", "q95", "qstar", "rad_input_frac", "rad_loss_frac", + "rmp_n_equal_1", # east/east_1 tree issues + "rmp_n_equal_1_phase", # east/east_1 tree issues + "v_loop", # east/east_1 tree issues "wmhd", "wmhd_rt", ]