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)

[gtrevisan]

I ran the new disruption-errors script for C-MOD and got:

disruption-errors.sh plasmas -m get_current_quench_time

      4  get_current_quench_time: UnboundLocalError("cannot access local variable 'ip' where it is not associated with a value")
     39  get_current_quench_time: ValueError('max() iterable argument is empty')

[ ERROR ] #1120918022 | get_current_quench_time: UnboundLocalError("cannot access local variable 'ip' where it is not associated with a value")
[ ERROR ] #1050505004 | get_current_quench_time: ValueError('max() iterable argument is empty')

1050505004 1120918022

[gtrevisan]

I ran the new disruption-errors script for DIII-D and got:

disruption-errors.sh plasmas -m get_current_quench_time

     47  get_current_quench_time: ValueError('attempt to get argmin of an empty sequence')
    443  Failed retrieval! KeyError('param')
   8545  get_current_quench_time: ValueError('max() iterable argument is empty')

[CRITICAL] #70581 | Failed retrieval! KeyError('param')
[ ERROR ] #78257 | get_current_quench_time: ValueError('attempt to get argmin of an empty sequence')
[ ERROR ] #70415 | get_current_quench_time: ValueError('max() iterable argument is empty')

70415 70581 78257

(actually, I ran with 1 cpu for roughly 65k out of 75k shots, then it was grinding almost to a halt for some MDSplus reason.)

[yumouwei]

Analysis of the 5 identified problematic shots

Shot	Identified cause	Implemented change
cmod 1050505004	get_end_of_current calculates duration=0.0914 which does not satisfy the bona fide plasma condition by Bob, then returns duration=0.	Add check for duration==0 (i.e. not bona fide plasma) in get_current_quench_time
cmod 1120918022	ip = array of nan (no ip signal)	Remove the try-except block around get_data_with_dims and let the framework catch MDSplus errors
d3d 70415	ip was not digitized to the end of the discharge. get_end_of_current marks this by setting duration = -duration	Check for negative duration and raise CalculationError
d3d 70581	Typo in ValueError	Corrected typo, change it to CalculationError
d3d 78257	ip_upright[time_indices] has only 1 element, thus didt_upright = [].	Check for len(didt_upright)==0 as well as other outputs of np.where, raise CalculationError accordingly.

[zapatace]

I reviewed 300 random shots from MAST. It works reasonably well.

I found some corner cases in which the method fails returning all NaNs. 18 cases out of 300 so 6%.

Few examples:

The list of pathological shots is 16553, 19792, 23730, 25287, 25459, 25735, 26606, 26612, 26617, 26842, 28219, 28232, 29077, 29136, 29146, 29645, 30011, 30278

I'll continue investigating...