 |
 |
The mission of this script is to process data from Mirnov coils, determine plasma position and thus help with plasma stabilization. Mirnov coils are type of magnetic diagnostics and are used to the plasma position determination on GOLEM. Four Mirnov coils are placed inside the tokamak 93 mm from the center of the chamber. The effective area of each coil is $A_{eff}=3.8\cdot10^{-3}$ m. Coils $mc_9$ and $mc_1$ are used to determination the horizontal plasma position and coils $mc_5$ and $mc_{13}$ to determination vertical plasma position.
|
|
import os
import glob
#remove files from the previous discharge
def remove_old(name):
if os.path.exists(name):
os.remove(name)
names = ['analysis.html','icon-fig.png', 'plasma_position.png']
for name in names:
remove_old(name)
if os.path.exists('Results/'):
file=glob.glob('Results/*')
for f in file:
os.remove(f)
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
#execute time
from timeit import default_timer as timer
start_execute = timer()
BDdata_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/{identifier}.csv"
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/{identifier}.csv" #NI Standart (DAS) new address
shot_no = 45781 # to be replaced by the actual discharge number
vacuum_shot = 45744 # to be replaced by the discharge command line paramater "vacuum_shot"
ds = np.DataSource(destpath='')
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, data_type='lvm'):
file = open_remote(shot_no, identifier, url)
return pd.read_csv(file, names=['Time', identifier],
index_col = 'Time').squeeze()
The coils measure voltage induced by changes in poloidal magnetic field. In order to obtain measured quantity (poloidal magnetic field), it is necessary to integrate the measured voltage and multiply by constant: $B(t)=\frac{1}{A_{eff}}\int_{0}^{t} U_{sig}(\tau)d\tau$
Ideally axis of the coils is perpendicular to the toroidal magnetic field, but in fact they are slightly deflected and measure toroidal magnetic field too. To determine the position of the plasma, we must eliminate this additional signal. For this we use vacuum discharge with the same parameters.
def elimination (shot_no, identifier, vacuum_shot):
#load data
mc = read_signal(shot_no, identifier, data_URL)
mc = mc.replace([np.inf, -np.inf, np.nan], value = 0)
konst = 1/(3.8e-03)
mc_vacuum = read_signal(vacuum_shot, identifier, data_URL)
mc_vacuum -= mc_vacuum.loc[:0.9e-3].mean() #remove offset
mc_vacuum = mc_vacuum.replace([np.inf, -np.inf], 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_start = float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/t_plasma_start').text)
plasma_end = float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/t_plasma_end').text)
print ('Plasma start =', round(plasma_start, 3), 'ms')
print ('Plasma end =', round(plasma_end, 3), 'ms')
Plasma start = 2.758 ms Plasma end = 26.648 ms
To calculate displacement of plasma column, we can use approximation of the straight conductor. If we measure the poloidal field on the two opposite sides of the column, its displacement can be expressed as: $\Delta=\frac{B_{\Theta=0}-B_{\Theta=\pi}}{B_{\Theta=0}+B_{\Theta=\pi}}\cdot b$
def horpol(shot_no,vacuum_shot):
mc1 = elimination(shot_no, 'U_mc1', vacuum_shot)
mc9 = elimination (shot_no, 'U_mc9', vacuum_shot)
b = 93
r = ((mc1-mc9)/(mc1+mc9))*b
r = r.replace([np.nan], value = 0)
return r.loc[plasma_start:plasma_end]
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)
<matplotlib.lines.Line2D at 0x7fb544366040>
def vertpol(shot_no, vacuum_shot):
mc5 = elimination(shot_no, 'U_mc5', vacuum_shot)
mc13 = elimination (shot_no, 'U_mc13', vacuum_shot)
b = 93
z = ((mc5-mc13)/(mc5+mc13))*b
z = z.replace([np.nan], value = 0)
return z.loc[plasma_start:plasma_end]
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)
<matplotlib.lines.Line2D at 0x7fb4f83ef340>
def plasma_radius(shot_no, vacuum_shot):
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:plasma_end]
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))
[(0.0, 85.0), (2.758, 26.648), Text(0.5, 0, 'Time [ms]'), Text(0, 0.5, '$a$ [mm]'), Text(0.5, 1.0, 'Plasma column radius #45781')]
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.757 ms end = 26.459 ms
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
| r | z | a | |
|---|---|---|---|
| Time | |||
| 2.758 | 8.537578 | 12.296587 | 70.030154 |
| 2.760 | 8.603504 | 12.408869 | 69.900320 |
| 2.762 | 8.750328 | 12.505458 | 69.737146 |
| 2.764 | 8.824954 | 12.718657 | 69.519559 |
| 2.766 | 8.759815 | 12.964624 | 69.353408 |
| ... | ... | ... | ... |
| 26.450 | -40.416906 | 73.414848 | 1.195071 |
| 26.450 | -40.416906 | 73.414848 | 1.195071 |
| 26.450 | -40.416906 | 73.414848 | 1.195071 |
| 26.450 | -40.416906 | 73.414848 | 1.195071 |
| 26.450 | -40.416906 | 73.414848 | 1.195071 |
11848 rows × 3 columns
savedata = 'plasma_position.csv'
df_processed.to_csv(savedata)
Markdown("[Plasma position data - $r, z, a$ ](./{})".format(savedata))
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
df_processed.describe()
| r | z | a | |
|---|---|---|---|
| count | 11848.000000 | 11848.000000 | 11848.000000 |
| mean | -9.520174 | 18.321851 | 57.975079 |
| std | 19.141883 | 10.681102 | 13.283431 |
| min | -61.368936 | 1.298354 | 0.069794 |
| 25% | -26.456829 | 8.915683 | 53.824164 |
| 50% | 0.760000 | 17.012791 | 61.789282 |
| 75% | 7.709776 | 23.937170 | 65.473336 |
| max | 18.234037 | 74.582332 | 78.755400 |
z_end=round(z_cut.iat[-25],2)
z_start=round(z_cut.iat[25],2)
r_end=round(r_cut.iat[-25],2)
r_start=round(r_cut.iat[25],2)
name=('z_start', 'z_end', 'r_start', 'r_end','plasma_start', 'plasma_end')
values=(z_start, z_end, r_start, r_end, start, end)
data={'Values':[z_start, z_end, r_start, r_end, start, end]}
results=pd.DataFrame(data, index=name)
dictionary='Results/'
j=0
for i in values :
with open(dictionary+name[j], 'w') as f:
f.write(str(i))
j+=1
results
| Values | |
|---|---|
| z_start | 9.970 |
| z_end | 63.740 |
| r_start | 17.220 |
| r_end | -42.530 |
| plasma_start | 2.757 |
| plasma_end | 26.459 |
# TODO: missing data on current in the stabilization coils
# responsestab = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/PowSup4InnerStab"%shot_no)
# stab = responsestab.text
# if not '-1' in stab or 'Not Found' in stab and exStab==False:
# fig, axs = plt.subplots(4, 1, sharex=True, dpi=200)
# dataNI = (read_signal(shot_no, 'Rog', data_URL, 'lvm')) #data from NI
# dataNI = dataNI.set_index('Time')
# dataNI = dataNI.set_index(dataNI.index*1e3)
# Irog = dataNI['RogQuadr']
# # Irog=abs(Irog)
# Irog*=1/5e-3
# r.plot(grid=True, ax=axs[0])
# z.plot(grid=True, ax=axs[1])
# a.plot(grid=True, ax=axs[2])
# Irog.plot(grid=True, ax=axs[3])
# axs[3].set(xlim=(start,end), ylabel = 'I [A]')
# 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))
# else:
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')
tms = np.arange(np.round(start),np.round(end)-1,1)
circle1 = plt.Circle((0, 0), 100, color='black', fill = False, lw = 5)
circle2 = plt.Circle((0, 0), 85, color='black', fill = False, lw = 3)
fig, ax = plt.subplots(figsize = (10,10))
ax.axvline(x=0, color = 'black', alpha = 0.5)
ax.axhline(y=0, color = 'black', alpha = 0.5)
ax.set_xlim(-100,100)
ax.set_ylim(-100,100)
ax.add_patch(circle1)
ax.add_patch(circle2)
plt.annotate('vessel',xy=(27-100, 170-100), xytext=(2-100, 190-100),
arrowprops=dict(facecolor='black', shrink=0.05), fontsize = 12)
plt.annotate('limiter',xy=(175-100, 140-100), xytext=(170-100, 190-100),
arrowprops=dict(facecolor='black', shrink=0.05), fontsize = 12)
color_idx = np.linspace(0, 1, tms.size)
a = 0
for i in tms:
mask_pos = ((df_processed.index>i) & (df_processed.index<i+1))
a_mean = df_processed['a'][mask_pos].mean()
r_mean = df_processed['r'][mask_pos].mean()
z_mean = df_processed['z'][mask_pos].mean()
circle3 = plt.Circle((r_mean, z_mean), a_mean, fill = False,
color=plt.cm.jet(color_idx[a]),label = str(i) +'-' + str(i+1) + ' ms')
ax.add_patch(circle3)
a+=1
plt.grid(True)
plt.title('Estimated plasma position during discharge', fontweight = 'bold', fontsize = 20)
leg=plt.legend(bbox_to_anchor = (1,1))
for text in leg.get_texts():
plt.setp(text, fontsize = 16)
ax.set_xlabel('r [mm]',fontweight = 'bold', fontsize = 15)
ax.set_ylabel('z [mm]', fontweight = 'bold', fontsize = 15)
plt.xticks( fontweight = 'bold', fontsize = 14)
plt.yticks( fontweight = 'bold', fontsize = 14)
plt.savefig('plasma_position.png', bbox_inches='tight')
# TODO: missing data on current in the stabilization coils
# if not '-1' in stab and not 'Not Found' in stab:
# dataNI = (read_signal(shot_no, 'Rog', data_URL, 'lvm')) #data from NI
# dataNI = dataNI.set_index('Time')
# dataNI=dataNI.set_index(dataNI.index*1e3)
# Irog = dataNI['RogQuadr']
# Irog*=1/5e-3
# Irog=Irog.loc[start:end]
# data = pd.concat([df_processed, Irog], axis='columns')
# data.to_csv('plasma_position.csv')
# data['RogQuadr'].plot()