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 = 33456 # 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)
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]:
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]:
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]:
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')
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]:
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]:
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')
