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 = 35790
#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/35790/Diagnostics/BasicDiagnostics/Results/t_plasma_start
http://golem.fjfi.cvut.cz/shots/35790/Diagnostics/BasicDiagnostics/Results/t_plasma_end
Plasma start = 4.641 ms
Plasma end = 19.037 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 0x7f61152e2650>

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

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]'),
 (4.641, 19.037),
 Text(0.5, 0, 'Time [ms]'),
 Text(0.5, 1.0, 'Plasma column radius #35790')]

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 = 4.641 ms
end = 19.035 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
4.642 -8.876519 55.113327 29.176426
4.643 -8.775680 55.225235 29.081853
4.644 -8.594790 55.311568 29.024648
4.645 -8.356172 55.312369 29.059998
4.646 -8.093418 55.192683 29.217067
... ... ... ...
19.031 -39.553957 10.248606 44.139880
19.032 -39.659334 9.937809 44.114516
19.033 -39.564361 10.893909 43.963237
19.034 -36.653731 10.259948 46.937387
19.035 -35.966709 6.559604 48.440014

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