Plasma column position

Adapted version of InnerStabilisation.ipynb by Daniela Kropáčková, Honza Stöckel and Vojtěch Svoboda.

This notebook loads Mirnov coil data and calculates the plasma column position from them.

Import libraries

In [1]:
import os

#delete files from the previous discharge
if os.path.exists('icon-fig.png'):
    os.remove('icon-fig.png')

import numpy as np
import matplotlib.pyplot as plt
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
from scipy import integrate, signal, interpolate
import pandas as pd
import requests
from IPython.display import Markdown

Set discharge number

To calculate the plasma position, $B_t$ contribution must be subtracted from the Mirnov coil signal. This contribution can be acquired by two ways: from a dedicated vacuum shot with the same toroidal magnetic field, or it can be estimated from the standard diagnostics. If a suitable vacuum shot is available, its number is stored in the variable vacuum_shot. If it isn't available, the vacuum_shot variable shall be set to False and the $B_t$ contribution shall be estimated.

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}/Diagnostics/BasicDiagnostics/{identifier}.csv" #BD = basic diagnostic
# data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/022Daniela0InnerStabilization/{identifier}.csv"  
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/DAS_raw_data_dir/Nidatap_6133.lvm" #NI Standart (DAS)

from urllib.request import urlopen

def get_file(url, shot, silent=False):
    URL = 'http://golem.fjfi.cvut.cz/shots/%i/%s' % (shot, url)
    if not silent:
        print(URL)
    f = urlopen(URL)
    try:
        return np.loadtxt(f)
    except ValueError: # the data couldn't be converted to a row of numbers
        print('The data could not be read! Using default values.')
        return np.array([35022])

shot_no = 35806
#shot_no = 35021
vacuum_shot = int(get_file('/Diagnostics/LimiterMirnovCoils/Parameters/vacuum_shot', 0))
#vacuum_shot = 35022 #number of the vacuum shot or 'False'
print('Vacuum shot: %i' % vacuum_shot)

ds = np.DataSource(destpath='')  
http://golem.fjfi.cvut.cz/shots/0//Diagnostics/LimiterMirnovCoils/Parameters/vacuum_shot
Vacuum shot: 35785

Data access

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, data_type='csv'): 
    file = open_remote(shot_no, identifier, url)
    if data_type == 'lvm':
        channels = ['Time', 'mc1', 'mc5', 'mc9', 'mc13', '5', 'RogQuadr', '7', '8', '9', '10', '11', '12']
        return pd.read_table(file, sep = '\t', names=channels)
    else:
        return pd.read_csv(file, names=['Time', identifier],
                     index_col = 'Time', squeeze=True)

def get_basic_diag_data(shot):
    url = ( 'http://golem.fjfi.cvut.cz/shots/{}/Diagnostics/BasicDiagnostics/'
            'basig_diagnostics_processed.csv'.format(shot) )
    print(url)
    return pd.read_csv(url, header=0, index_col='Time')

Mirnov coil signal processing: integration and $B_t$ elimination

Load the data of a given Mirnov coil, integrate it and subtract the $B_t$ contribution from them.

In [4]:
def elimination (shot_no, identifier, vacuum_shot = False):
    #load data 
    mc = (read_signal(shot_no, identifier, data_URL, 'lvm'))
    mc = mc.set_index('Time')
    mc = mc.replace([np.inf, -np.inf, np.nan], value = 0)
    mc = mc[identifier]
    
    konst = 1/(3.8e-03)

    if vacuum_shot == False: 
        signal_start = mc.index[0]
        length = len(mc)
        basic_diag = get_basic_diag_data(shot_no)
        Bt = basic_diag['Bt']
        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, data_URL, 'lvm'))
        mc_vacuum = mc_vacuum.set_index('Time')
        mc_vacuum = mc_vacuum[identifier]
        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

Find plasma start and end

In [5]:
plasma_start = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_start', shot_no))
plasma_end = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_end', shot_no))

print ('Plasma start =', round(plasma_start, 3), 'ms')
print ('Plasma end =', round(plasma_end, 3), 'ms')
http://golem.fjfi.cvut.cz/shots/35806/Diagnostics/BasicDiagnostics/Results/t_plasma_start
http://golem.fjfi.cvut.cz/shots/35806/Diagnostics/BasicDiagnostics/Results/t_plasma_end
Plasma start = 6.613 ms
Plasma end = 20.655 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: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)
Out[6]:
<matplotlib.lines.Line2D at 0x7f8bb887bed0>

Vertical plasma position $\Delta z$ calculation

In [7]:
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: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)
Out[7]:
<matplotlib.lines.Line2D at 0x7f8bb4090c50>

Plasma column radius $a$ calculation

In [8]:
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: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))
Out[8]:
[(0, 85),
 Text(0, 0.5, '$a$ [mm]'),
 (6.613, 20.655),
 Text(0.5, 0, 'Time [ms]'),
 Text(0.5, 1.0, 'Plasma column radius #35806')]

Save the results

In [9]:
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 = 6.613 ms
end = 20.647 ms
In [10]:
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[10]:
r z a
Time
6.614 17.250590 15.382756 61.886973
6.615 17.243683 15.406650 61.876214
6.616 17.238403 15.408256 61.879081
6.617 17.251596 15.369543 61.895013
6.618 17.262596 15.327243 61.914905
... ... ... ...
20.643 -68.205629 26.613487 11.786029
20.644 -68.900116 25.822943 11.419769
20.645 -72.573759 23.502680 8.715490
20.646 -75.952597 23.671396 5.444158
20.647 -76.863654 24.236286 4.405839

14034 rows × 3 columns

In [11]:
savedata = 'plasma_position.csv'
df_processed.to_csv(savedata)

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

Plot the results

In [12]:
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[12]:
In [13]:
responsestab = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/PowSup4InnerStab"%shot_no)
stab = responsestab.text

if '-1' in stab or 'Not Found' in stab:
    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))

else:
    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))


plt.savefig('icon-fig')