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 = 33448 # 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 = 13.048 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 0x7f9b01543b50>

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 0x7f9b31ec0710>

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.50036784, 13.0483678),
 Text(0.5, 0, 'Time [ms]'),
 Text(0.5, 1.0, 'Plasma column radius #33448')]
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.55 ms
end = 13.047 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.550558 64.048671 54.484945 0.911704
2.551558 63.210823 50.239757 4.255705
2.552558 62.411990 46.162823 7.370993
2.553558 61.582456 42.461409 10.197793
2.554558 60.709687 39.300787 12.679754
... ... ... ...
13.041560 -36.484591 46.457986 25.928265
13.042560 -36.555268 46.494981 25.855508
13.043560 -36.573893 46.500962 25.839294
13.044560 -36.590232 46.535338 25.802173
13.045560 -36.613679 46.603307 25.734245

10496 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')