Inner Stabilisation stuff

Basic info

co-authors: DanielaK, HonzaS, VojtechS

Description ...

Main result: Other results ...

Procedure (This notebook to download)

In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt

from scipy import integrate, signal, interpolate
import pandas as pd

import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
import requests

from IPython.display import Markdown
In [2]:
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/LimiterMirnovCoils/{identifier}.csv"  #Mirnov coils and quadrupole
BDdata_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/StandardDAS/{identifier}.csv" #BD = basic diagnostic


shot_no = 33445 # to be replaced by the actual discharge number
#shot_no = 32607 # Test High performance shot
# shot_no = 32660 # Test Low performance shot
# shot_no = 32947
vacuum_shot = 33440  # to be replaced by the discharge command line paramater "vacuum_shot"
# vacuum_shot = 32929 #number of the vacuum shot or 'False'


ds = np.DataSource(destpath='') #/tmp
In [3]:
def open_remote(shot_no, identifier, url_template=data_URL):
    return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))

def read_signal(shot_no, identifier, url = data_URL): 
    file = open_remote(shot_no, identifier, url)
    return pd.read_csv(file, names=['Time', identifier],
                     index_col = 'Time', squeeze=True)

Data integration and $B_t$ elimination

In [4]:
def elimination (shot_no, identifier, vacuum_shot = False):
    #load data 
    mc = (read_signal(shot_no, identifier))
    mc = mc.replace([np.inf, -np.inf, np.nan], value = 0)
    
    konst = 1/(3.8e-03)
    
       
    if vacuum_shot == False: 
        signal_start = mc.index[0]
        length = len(mc)
        Bt = read_signal(shot_no, 'BtCoil_integrated', BDdata_URL).loc[signal_start:signal_start+length*1e-06]
        if len(Bt)>len(mc):
            Bt = Bt.iloc[:length]            
        if len(Bt)<len(mc):
            mc = mc.iloc[:len(Bt)]
        
        if identifier == 'mc1':
            k=300
        elif identifier == 'mc5':
            k= 14
        elif identifier == 'mc9':
            k = 31
        elif identifier == 'mc13':
            k = -100 
        
        mc_vacuum = Bt/k
    else:
        mc_vacuum = (read_signal(vacuum_shot, identifier))
        mc_vacuum -= mc_vacuum.loc[:0.9e-3].mean()    #remove offset
        mc_vacuum = mc_vacuum.replace([np.inf, -np.inf, np.nan], value = 0)
        mc_vacuum = pd.Series(integrate.cumtrapz(mc_vacuum, x=mc_vacuum.index, initial=0) * konst,
                    index=mc_vacuum.index*1000, name= identifier)    #integration

    mc -= mc.loc[:0.9e-3].mean()  #remove offset
       
    mc = pd.Series(integrate.cumtrapz(mc, x=mc.index, initial=0) * konst,
                    index=mc.index*1000, name= identifier)    #integration
    
    #Bt elimination
    mc_vacuum = np.array(mc_vacuum) 
    mc_elim = mc - mc_vacuum
    
    return mc_elim

Plasma life time

In [5]:
loop_voltage = read_signal(shot_no, 'LoopVoltageCoil_raw', BDdata_URL)

dIpch = read_signal(shot_no, 'RogowskiCoil_raw', BDdata_URL)

dIpch -= dIpch.loc[:0.9e-3].mean()

Ipch = pd.Series(integrate.cumtrapz(dIpch, x=dIpch.index, initial=0) * (-5.3*1e06),
                index=dIpch.index, name='Ipch')

U_l_func = interpolate.interp1d(loop_voltage.index, loop_voltage)  
def dIch_dt(t, Ich):
    return (U_l_func(t) - 0.0097 * Ich) / (1.2e-6/2)
t_span = loop_voltage.index[[0, -1]]
solution = integrate.solve_ivp(dIch_dt, t_span, [0], t_eval=loop_voltage.index, )
Ich = pd.Series(solution.y[0], index=loop_voltage.index, name='Ich')
Ip = Ipch - Ich
Ip.name = 'Ip'

Ip_detect = Ip.loc[0.0025:]

dt = (Ip_detect.index[-1] - Ip_detect.index[0]) / (Ip_detect.index.size) 

window_length = int(0.5e-3/dt)  
if window_length % 2 == 0:  
    window_length += 1
dIp = pd.Series(signal.savgol_filter(Ip_detect, window_length, polyorder=3, deriv=1, delta=dt),
                name='dIp', index=Ip_detect.index) / 1e6 

threshold = 0.05

CD = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/CD_orientation" % shot_no)
CD_orientation = CD.text

if "ACW" in CD_orientation:
    plasma_start = dIp[dIp < dIp.min()*threshold].index[0]*1e3 
    plasma_end = dIp.idxmax()*1e3 
else: 
    plasma_start = dIp[dIp > dIp.max()*threshold].index[0]*1e3 
    plasma_end = dIp.idxmin()*1e3     


print ('Plasma start =', round(plasma_start, 3), 'ms')
print ('Plasma end =', round(plasma_end, 3), 'ms')
# print (CD_orientation)

Plasma start = 2.5 ms
Plasma end = 12.484 ms

Horizontal plasma position $\Delta r$ calculation

In [6]:
def horpol(shot_no, vacuum_shot=False):
    mc1 = elimination(shot_no, 'mc1', vacuum_shot)
    mc9 = elimination (shot_no, 'mc9', vacuum_shot)
    
    b = 93
    
    r = ((mc1-mc9)/(mc1+mc9))*b
    r = r.replace([np.nan], value = 0)
    
#     return r.loc[plasma_start:]
    return r.loc[plasma_start:plasma_end]
#     return r
In [7]:
r = horpol(shot_no, vacuum_shot)
ax = r.plot()
ax.set(ylim=(-85,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\Delta$r [mm]', title = 'Horizontal plasma position #{}'.format(shot_no))
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
Out[7]:
<matplotlib.lines.Line2D at 0x7fd0f93e1710>

Vertical plasma position $\Delta z$ calculation

In [8]:
def vertpol(shot_no, vacuum_shot = False):
    mc5 = elimination(shot_no, 'mc5', vacuum_shot)
    mc13 = elimination (shot_no, 'mc13', vacuum_shot)
    
    b = 93
    
    z = ((mc5-mc13)/(mc5+mc13))*b
    z = z.replace([np.nan], value = 0)
#     return z.loc[plasma_start:]
    return z.loc[plasma_start:plasma_end]
#     return z
In [9]:
z = vertpol (shot_no, vacuum_shot)
ax = z.plot()
ax.set(ylim=(-85, 85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\Delta$z [mm]', title = 'Vertical plasma position #{}'.format(shot_no))
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
Out[9]:
<matplotlib.lines.Line2D at 0x7fd1240c6090>

Plasma column radius $a$ calculation

In [10]:
def plasma_radius(shot_no, vacuum_shot=False):
    r = horpol(shot_no, vacuum_shot) 
    z = vertpol(shot_no, vacuum_shot) 
    
    a0 = 85
    a = a0 - np.sqrt((r**2)+(z**2)) 
    a = a.replace([np.nan], value = 0)
#     return a.loc[plasma_start:]
    return a.loc[plasma_start:plasma_end]
#     return a
In [11]:
a = plasma_radius(shot_no,vacuum_shot)
ax = a.plot()
ax.set(ylim=(0,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]', title = 'Plasma column radius #{}'.format(shot_no))
Out[11]:
[(0, 85),
 Text(0, 0.5, '$a$ [mm]'),
 (2.5002904000000004, 12.4842904),
 Text(0.5, 0, 'Time [ms]'),
 Text(0.5, 1.0, 'Plasma column radius #33445')]
In [12]:
plasma_time = []
t = 0
for i in a:
    if i>85 or i <0:
        a = a.replace(i, value = 0)
    else:

        plasma_time.append(a.index[t])

    t+=1
start = plasma_time[0]-1e-03 
end = plasma_time[-1]-1e-03 
print('start =', round(start, 3), 'ms')
print('end =', round(end, 3), 'ms')

start = 2.5 ms
end = 12.483 ms

Graphs

In [13]:
r_cut = r.loc[start:end]
a_cut = a.loc[start:end]
z_cut = z.loc[start:end]
df_processed = pd.concat(
    [r_cut.rename('r'), z_cut.rename('z'), a_cut.rename('a')], axis= 'columns')
df_processed
Out[13]:
r z a
Time
2.500952 24.288464 -4.303127 60.333294
2.501952 24.309423 -3.897117 60.380180
2.502952 24.391292 -3.537700 60.353491
2.503952 24.534909 -3.238185 60.252322
2.504952 24.708397 -2.983324 60.112149
... ... ... ...
12.478950 -45.025625 42.152406 23.322353
12.479950 -44.862510 42.194620 23.412413
12.480950 -44.772532 42.232435 23.451985
12.481950 -44.822731 42.276133 23.385483
12.482950 -45.053277 42.360885 23.159541

9983 rows × 3 columns

In [14]:
savedata = 'plasma_position_%i.csv' %shot_no 
df_processed.to_csv(savedata)

Data to download

In [15]:
Markdown("[Plasma position data - r, z, a ](./{})".format(savedata))
Out[15]:

Plasma position data - r, z, a

In [16]:
hline = hv.HLine(0)
hline.opts(
    color='k', 
    line_dash='dashed',
    alpha = 0.4,
    line_width=1.0)

layout = hv.Layout([df_processed[v].hvplot.line(
    xlabel='', ylabel=l,ylim=(-85,85), xlim=(start,end),legend=False, title='', grid=True, group_label=v)
                    for (v, l) in [('r', ' r [mm]'), ('z', 'z [mm]'), ('a', 'a [mm]')] ])*hline

plot=layout.cols(1).opts(hv.opts.Curve(width=600, height=200),  
                    hv.opts.Curve('a', xlabel='time [ms]'))
plot
Out[16]:
In [17]:
fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)
r.plot(grid=True, ax=axs[0])
z.plot(grid=True, ax=axs[1])
a.plot(grid=True, ax=axs[2])
axs[2].set(ylim=(0,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]')
axs[1].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\Delta$z [mm]')
axs[0].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\Delta$r [mm]', title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))


plt.savefig('icon-fig')