Implement generic current time computation method

[yumouwei]

Implemented changes

• Added GenericPhysicsMethods.get_current_quench_time(...) and GenericUtilMethod.get_end_of_current(...) based on Bob's original scripts in IDL & MATLAB (see implement disruption time computation #223).
• The computed current_quench_time is returned as a column with the same value for each shot.

Testing

• So far, I've run spot checks on the testing shots in CMOD & DIII-D, and a few DIII-D shots in the 201xxx range which are recently used in training the DPRF model for the AI Task Force experiment. All shots except cmod 1150805014 (see below) agrees with the existing t_disrupt from the cached SQL database.

TODOs

• Figure out how to return the current_quench_time in the dataset/dataframe.
• Implement unit test for this method.
• Apply this method to MAST & HBT-EP and determine the proper threshold values.
  • Complete implementation for MAST & HBT-EP
  • Determine MAST thresholds
  • Determine HBT-EP thresholds

References

• implement disruption time computation #223

Note on cmod 1150805014

cmod 1150805014 is considered a non-disruptive shot in the SQL database. However, the parameters of this shot do satisfy all of the 6 testing criteria. The shot does look like a proper ramp down disruption with sufficient pre-disruption ip, so I believe it should be relabeled as disruptive.

[yumouwei]

Adapting the method to MAST & HBT-EP

MAST

I copied in the thresholds from CMOD and tested them with the testing shots. Right now the method is able to locate the current quench event in the disruptive shots, but it occasionally marks a non-disruptive shot as disruptive (see mast 29695 below). I tried adjusting some of the thresholds but it didn't help.

I think the problem is that right now this method has some hard-coded assumptions about the machine's disruption & rampdown timescales (e.g. L219-222 where it checks for the pre-ip-spike maximum current based on the computed shot duration - 50 ms, or L232-237 where it searches for the candidate t_disrupt within +-50ms of the shot duration). So far these assumptions have worked out fine across CMOD, DIII-D, & EAST, but they should be generalized and become machine-specific so that we can adapt this method to other machines with either faster or slower timescales. I believe this should be done in the next revamp of this method when we plan to create a disruption dataset of MAST.

disruption-py/disruption_py/machine/generic/physics.py

Lines 219 to 222 in 08cea0a

	# Find the maximum plasma current excluding the current spike
	ip_upright = ip * polarity
	(time_indices,) = np.where((t_ip > 0) & (t_ip < duration - 0.050))
	ip_max = max(ip_upright[time_indices]) * polarity

disruption-py/disruption_py/machine/generic/physics.py

Lines 232 to 237 in 08cea0a

	# Compute dI/dt during the latter part of the discharge
	(time_indices,) = np.where((t_ip > duration - 0.05) & (t_ip < duration + 0.05))
	didt_upright = np.diff(ip_upright[time_indices]) / np.diff(t_ip[time_indices])
	indx = np.argmin(didt_upright)
	candidate_max_didt = didt_upright[indx] * polarity
	candidate_t_disrupt = t_ip[time_indices[indx]]

HBT-EP

For HBT-EP, the problem with the timescale assumptions is more severe. An HBT shot usually lasts between 5~10 ms. It always ends with a disruption, and the disruption precursor usually occurs within the last 1 ms prior to the current quench. These shorter timescales break every assumption about the timescales of current quench events in the method, which makes the method in its current form incompatible with HBT regardless of the threshold values. Similar to the suggestions for MAST, I believe this should be done in a subsequent revamp after a more throughout study of the disruption dynamics in HBT-EP.

[yumouwei]

Comparison of current_quench_time vs t_disrupt

• The following figures compare the current_quench_time computed using this new physics method with t_disrupt previously computed by Bob Granetz using all shots from the disruption_warning database.

1. "Confusion Matrix"

2. Difference between current_quench_time vs t_disrupt

• These figures include shots that have both current_quench_time and t_disrupt (i.e. lower right category in the "confusion matrix")

Analysis

• The get_current_quench_time physics method is mostly robust for D3D shots. This is expected as the method was translated mainly based on the D3D version of test_for_disruption.m and end_of_current_d3d.m.
  • I also expect it to work well for EAST (& KSTAR once we incorporate it) since their versions of the MATLAB scripts are nearly identical to the D3D version.
• For C-Mod, right now there are a lot of "false positives" (shots that don't have t_disrupt but have current_quench_time). I spot checked a couple of these shots and many of them are what I'd consider as "disruptive" based on the descriptions in implement disruption time computation #223 similar to cmod 1150805014 shown above. One exception I found so far is cmod 1120502001 shown below.
• Upon further investigation, I found out that there are actually a couple differences between test_disrupt.pro (for CMOD) and the other MATLAB scripts：
  • C-Mod uses different search windows to get the intermediate variables like ip_max, ip_final and candidate_di_dt. This is the exact same timescale problem I mentioned earlier about MAST and HBT-EP.
  • C-Mod uses a different method to find the pre-CQ Ip (candidate_ip0) by backtracking from candidate_t_disrupt until dI/dt crosses a threshold value, as opposed to computing the mean Ip between a window based on the end of current.

Based on these analysis I'd expect this method to take more time to complete. Alternatively we can merge it "as-is", keep everything that uses t_disrupt intact, and implement the revamps before we decouple with SQL DBs.

Additional figure

• cmod 1120502001 (should be labeled as non-disruptive)

• cmod 1140522013 (Ip0 < 100kA so should be considered as non-disruptive)