Add C-Mod Thermal Quench Labeler

[hwietfeldt]

Overview

I created a C-Mod physics method get_thermal_quench_onset_time() that labels the onset of the thermal quench (TQ) for a given shot (NaN for non-disruptive shots) using a vertical SXR array. The TQ time is returned as a column with identical values for each time slice. I used SXR instead of ECE because ECE is less reliable (especially for shots at low Bt) and because the SXR arrays can track the plasma core during a hot VDE. I used a vertical array because the horizontal array wasn't working for most 2005 shots. For 100 manually labeled flattop TQ onsets, this method labeled the onset within 1 ms of the manual label in 90% of cases. For the remaining 10%, the automated TQ label is usually still a more accurate disruption time metric compared to the current t_disrupt based on the maximum of |dIp/dt|. This method has also worked well in my UFO labeling workflow. I have not tested this method during ramp-up or ramp-down

Steps of Method

1. Read Ip and SXR chords from MDSplus
2. Subtract background from each SXR chord
3. Discard bad chords dominated by noise by removing chords with low autocorrelation.
4. Filter the SXR chords using a non-causal Butterworth low pass filter.
   This helps to smooth transient, positive SXR spikes during the TQ, which I attribute to recombination radiation. There is a tradeoff when choosing the cutoff frequency between maintaining decent temporal resolution and removing the recombination spikes. I found a cutoff frequency of 1 kHz works best.
5. Select the SXR chord viewing the core by taking the maximum of SXR across all chords for each time slice.
6. Search for the onset of the CQ based on Ip and t_disrupt.
   This enables us to search for the TQ in a small time window prior to the CQ onset. We need a small search window to avoid labeling sawtooth crashes as the TQ while enabling us to label the TQ for shots with long CQs (see shots 1050311013, 1050802017).
7. Find the midpoint of the TQ by finding the minimum of dSXR/dt within a 5 ms window prior to the CQ onset.
8. Find the TQ onset within a 0.5 ms window prior to the TQ midpoint.
   Note that almost all TQs have a duration less than 1 ms.

Case that works well

We correctly label the TQ onset, which is around 7 ms prior to $t_{disrupt}$. The $T_{e0}$ signal is junk because ECE was not working on this shot, which motivates the use of SXR for TQ labeling on C-Mod.

Results from Testing

I tested this method on 100 disruptive shots for which I manually labeled the TQ onset and end using ECE, SXR, H alpha, and other diagnostics. The plot below compared the error in the automated TQ onset vs the error we currently get when we use t_disrupt as the onset of the disruption compared to my manually labeled TQ onsets. From the y-values, we see that the automated TQ labeling routine labeled the TQ onset within 0.5 ms of the actual onset for most shots. As guaranteed by the algorithm, no automated TQ labels are later than t_disrupt. There are a several labels that are early or late by more than a ms. In these cases, there may be some ambiguity in the TQ onset (multiple partial thermal collapses), or there were strong recombination spikes. The search windows and filtering cutoff frequency have tradeoffs, so I was not able to eliminate these outliers. Nevertheless, I believe that the TQ labels provide a much more accurate metric of the onset time of the disruption for many of our workflows compared to the current t_disrupt.

Most early label

In this shot, the TQ time is automatically labeled around 0.511 s, which is at the onset of a sawtooth or small thermal collapse. I manually labeled the TQ onset around 0.514 s, which is the main thermal quench according to ECE (second panel) and the Ip spike (first panel). Again, I don't think that is 3 ms early label for the TQ is too bad. Using t_disrupt may be worse for some workflows as it is clearly after the disruption onset. If using t_disrupt, the user risks leaking disruption information into their disruption prediction model.

Most late label

Here a very strong SXR recombination spike coupled with a long delay between the Ip spike and when Ip crashes causes a late label. The TQ onset is not much better than t_disrupt, but it's also not worse.

Timing

I timed various steps of the TQ labeler. It can be slow for disruptive shots from 2012-2016. Note that for 2005 shots, the SXR was digitized at 5 kHz, not 250 kHz, so the algorithm is faster. Reading all 38 chords from MDSplus is a bottleneck. I tried using MDSplus's getMany() function to read the chords in one batch call; however, that was slower. For non-disruptive shots, the algorithm is fast because it immediately returns NaNs if t_disrupt is 0.

Timing averages for 2005 shots
Reading MDSplus: 0.20789 s
Calculating autocorrelation: 0.0264536 s
Butterworth Filter: 0.00573944 s
Finding TQ: 0.000120203 s
Total: 0.242485 s

Timing averages for 2012-2016 shots
Reading MDSplus: 3.87387 s
Calculating autocorrelation: 0.0920277 s
Butterworth Filter: 0.194109 s
Finding TQ: 0.00131993 s
Total: 4.16501 s

[hwietfeldt]

@yumouwei @AlexSaperstein

Once we are satisfied with the method and testing, I'll clean up the method and turn the draft PR into a normal PR.

[gtrevisan]

I think we only need to discuss about the SQL-based disruption time vs our python-based current quench time, and then we can merge!

[gtrevisan]

this is using the disruption time as read from the disruptions SQL table.
is this what we want?
I'd imagine we want to rely on our python-controlled current quench time, now.

[gtrevisan]

same logic here, to be discussed.
ideally you'd like to see the different between the DB disruption time and our python CQ time for a handful of your shots.

[gtrevisan]

and again -- you are using the SQL-based disruption_time extensively, yet its calculation is not within our control. we might want to discuss this.

[hwietfeldt]

I changed the scipy import to just one import scipy with full qualifiers for each package at each function call. In addition to changes in get_thermal_quench_time(), there was also a change in one line of _get_te_profile_params_ece(), which used scipy.constants. I performed a diff compare on DisruptionPy outputs for both functions on 1000 C-Mod shots before and after the changes to the scipy import. Results were identical

[gtrevisan]

after discussions with Henry, this is our to-do list:

• move from old-school disruption_time to our new current_quench_time
• evaluate cyclic import
• evaluate timing

[hwietfeldt]

Switch from SQL disruption time to Generic CQ time

I tested the differences in thermal quench time labels using the generic current_quench_time vs the SQL disruption_time on all C-Mod shots that reached flattop from 2012-2016 (~8000 shots). Below are plots similar to #545:

More shots with thermal quench

Using the generic CQ method instead of the SQL disruption_time, more shots are labeled to have a TQ. This is because the generic CQ method labels more shots to have CQs than the routines used to build the SQL table containing disruption_time (see #545). A shot should be labeled to have a TQ if and only if it's labeled to have a CQ, so this behavior is desired.

Differences in TQ time values are basically unchanged

The above plot is a histogram in the differences in the output thermal_quench_time when using the generic current_quench_time vs the SQL disruption_time for all 2840 shots in which both choices result in a labeled TQ. The differences are $\ll 1 \mu$s, which is much less than the error based on comparison to manual labels (see above posts), so these differences are negligible. This finding is expected because according to #545, differences in labeled CQ times are fractions of a ms, and the TQ time labeler should not be sensitive to differences in the CQ of < 1 ms.

[hwietfeldt]

Speed-up for Thermal Quench Labeling

The computation of the thermal quench time is limited by the read-time for the soft X-ray (SXR) chords.
The main issue is that for shots in the years 2012-2016, the SXR chords were digitized to 250 kHz,
resulting in very long arrays.

Summary of speed-ups

Prior to the speed-ups, the read time per disrupted shot (2012-2016) was 3.7 s on the MFE workstations.
Note that for shots without a current quench (about half of shots), get_thermal_quench_time immediately returns NaN.

• Read the time base of only one chord
• Only read 21 chords instead of 38 (excluding the outer 17 views did not result in meaningful changes)
• Subtract SXR chord background within a TDI expression
• Read only snippets of SXR chords (only minor speed-up of 0.06 s / shot)
• vectorized autocorrelation calculation (minor speed-up; not relevant for reading signals)

These improvements result in a 0.7 s read time per disrupted shot (2012-2016) on the MFE workstations.
The first three bullets each significantly contribute to the speed-up.
The output thermal_quench_time is not meaningully different with these speed-ups for 50 shots tested.

TDI Expression for reading chord data

The data for each chord is now read with a TDI expression that grabs only a snippet of the data and subtracts the background signal. This was found to be the fastest way to get background-corrected chord data.

sig = f"{array_path}:CHORD_{idx_first_chord+i+1:02}"
tdi_expr = f"""
    _s = data({sig});
    _y = _s[{j_chord_start} : {j_chord_end-1}]
        - mean(_s[{j_bgrnd_start} : {j_t0-1}]);
    _y
    """
chord = params.get_data(tdi_expr, tree_name="xtomo")

Changes to output thermal quenches

For about 1400 flattop disruptions in the C-Mod disruption warning database, I calculated the differences in the output thermal_quench_time compared to those prior to the speed-up related commits. The vast majority of shots had no changes. There were five cases of shots in which the magnitude of the differences was over 1 ms. On these shots, the SNR of the SXR chords was low (so the algorithm struggles anyways), and minor diffs in the autocorrelation calculation on smaller snippets propagated into the output. Importantly, there were no shots for which the difference was larger than the error range based on manual testing (see plot in previous comments).

Attempted speed-ups that did not work

Resampling in TDI expression

I attempted to resample the signals with the following TDI expression:

tdi_expr = (
    f"data(resample(make_signal("
    f"data({sig})[{i_chord_start}:{i_chord_end-1}], *, "
    f"dim_of({sig})[{i_chord_start}:{i_chord_end-1}]), {t_wndw_start}, {t_wndw_end}, {0.001}))"
)
chord = params.get_data(tdi_expr, tree_name="xtomo")

This resulted in 1.25 s read time per disruptive shot compared to 0.8 s / shot in the final implementation.
The likely explanation is that the resampling function still needs to decompress the entire SXR signal before resampling and addtionally needs to read the time base.

Batch call of SXR chords with fall-back

I attempted to get all chords in one TDI expression containing an array of data() calls; however, this did not speed-up the reads.

Removing background subtraction

Background subtraction is necessary when filtering out SXR chords dominated by white noise (see C-Mod shot 1050311013).