Tokamak GOLEM Basic diagnostics¶
Prerequisities: function definitions¶
Load libraries
In [1]:
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy import constants, integrate, signal, interpolate
#import sqlalchemy # high-level library for SQL in Python
import pandas as pd
For interactive web figures
In [2]:
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
For conditional rich-text boxes
In [3]:
from IPython.display import Markdown
Define global constants.
In [44]:
base_URL = "http://golem.fjfi.cvut.cz/utils/data/{shot_no}/{identifier}" #remote access
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/StandardDAS/{identifier}.csv" # TODO workaround
destination=''
#os.makedirs(destination, exist_ok=True );
#shot_no = os.environ.get('SHOT_NO', 0)
shot_no = 32183
The DataSource downloads and caches data (by full URL) in a temporary directory and makes them accessible as files.
In [45]:
ds = np.DataSource(destpath='/tmp')
In [46]:
def print_and_save(PhysQuant, Value, Format ):
print(PhysQuant+" = %.5f" % Value)
with open(destination+PhysQuant, 'w') as f:
f.write("%.3f" % Value)
#engine = sqlalchemy.create_engine('postgresql://golem:rabijosille@192.168.2.116/golem_database')
#engine.execute("""UPDATE shots SET pressure=50 WHERE shot_no IN(SELECT max(shot_no) FROM shots)""")
def open_remote(shot_no, identifier, url_template=base_URL):
return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))
def read_value(shot_no, identifier):
"""Return the value for given shot as a number if possible"""
value = open_remote(shot_no, identifier, base_URL).read()
return pd.to_numeric(value, errors='ignore')
def read_signal(shot_no, identifier):
file = open_remote(shot_no, identifier, data_URL)
return pd.read_csv(file, names=['Time',identifier],
index_col='Time', squeeze=True) # squeeze makes simple 1-column signals a Series
$U_l$ management¶
Check the data availability¶
In [47]:
loop_voltage = read_signal(shot_no, 'LoopVoltageCoil_raw')
ax = loop_voltage.plot(grid=True)
ax.set(xlabel="Time [s]", ylabel="$U_l$ [V]", title="Loop voltage $U_l$ #{}".format(shot_no));
Getting $U_l$ mean¶
In [16]:
loop_voltage_mean = loop_voltage.abs().mean()
print_and_save("UloopMean", loop_voltage_mean, "%.3f")
Getting $U_l$ max¶
The 99% quantile is used to neglect large peaks. This method also ignores NaNs.
In [17]:
loop_voltage_max = loop_voltage.quantile(0.99)
print_and_save("UloopMax", loop_voltage_max, "%.3f")
$B_t$ calculation¶
Check the data availability¶
It is as magnetic measurement, so the raw data only give $\frac{dB_t}{dt}$
In [18]:
dBt = read_signal(shot_no,'BtCoil_raw')
ax = dBt.plot(grid=True)
ax.set(xlabel="Time [s]", ylabel="$dU_{B_t}/dt$ [V]", title="BtCoil_raw signal #{}".format(shot_no));
Integration (it is a magnetic diagnostic) & calibration¶
In [19]:
K_BtCoil = read_value(shot_no, 'K_BtCoil') # Get BtCoil calibration factor
print('BtCoil calibration factor K_BtCoil={} T/(Vs)'.format(K_BtCoil))
In [20]:
Bt = pd.Series(integrate.cumtrapz(dBt, x=dBt.index, initial=0) * K_BtCoil,
index=dBt.index, name='Bt')
ax = Bt.plot(grid=True)
ax.set(xlabel="Time [s]", ylabel="$B_t$ [T]", title="Toroidal magnetic field $B_t$ #{}".format(shot_no));
Getting $B_t$ mean¶
In [21]:
Bt_mean = Bt.abs().mean()
print_and_save("BtMean", Bt_mean, "%.3f")
Getting $B_t$ max¶
In [22]:
Bt_max = Bt.quantile(0.99)
print_and_save("BtMax", Bt_max, "%.3f")
Plasma + Chamber current $I_{p+ch}$ calculation¶
Check the data availability¶
Because it is a magnetic measurement, the raw data only gives $\frac{dI_{p+ch}}{dt}$
In [23]:
dIpch = read_signal(shot_no, 'RogowskiCoil_raw')
dIpch -= dIpch.loc[:0.9e-3].mean() # subtract offset
ax = dIpch.plot(grid=True)
ax.set(xlabel="Time [s]", ylabel="$dU_{I_{p+ch}}/dt$ [V]", title="RogowskiCoil_raw signal #{}".format(shot_no));
Integration (it is a magnetic diagnostics) & calibration¶
In [24]:
K_RogowskiCoil = read_value(shot_no, 'K_RogowskiCoil') # Get RogowskiCoil calibration factor
print('RogowskiCoil calibration factor K_RogowskiCoil={} A/(Vs)'.format(K_RogowskiCoil))
In [25]:
Ipch = pd.Series(integrate.cumtrapz(dIpch, x=dIpch.index, initial=0) * K_RogowskiCoil * (-1),
index=dIpch.index, name='Ipch')
ax = Ipch.plot(grid=True)
ax.set(xlabel="Time [s]", ylabel="$I_{p+ch}$ [A]", title="Total (plasma+chamber) current $I_{{p+ch}}$ #{}".format(shot_no));
Plasma current $I_{p}$ calculation¶
In [26]:
R_chamber = read_value(shot_no, 'R_chamber') # Get Chamber resistivity
print('Chamber resistivity R_chamber={} Ohm'.format(R_chamber))
In [27]:
L_chamber = 1.2e-6 # H
The chamber current $I_{ch}$ satisfies the equation (neglecting the mutual inductance with the plasma) $$U_l = R_{ch} I_{ch} + L_{ch} \frac{d I_{ch}}{dt}$$ Therefore, the following initial value problem must be solved to take into account the chamber inductance properly $$\frac{d I_{ch}}{dt} = \frac{1}{L_{ch}}\left( U_l - R_{ch} I_{ch}\right), \quad I_{ch}(t=0)=0$$
In [28]:
U_l_func = interpolate.interp1d(loop_voltage.index, loop_voltage) # 1D interpolator
def dIch_dt(t, Ich):
return (U_l_func(t) - R_chamber * Ich) / L_chamber
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')
In [29]:
Ip_naive = Ipch - loop_voltage/R_chamber # creates a new Series
Ip = Ipch - Ich
Ip.name = 'Ip'
Ip_naive.plot(grid=True, label='naive $I_{ch}=U_l/R_{ch}$')
ax = Ip.plot(grid=True, label=r'using $U_l = R_{ch} I_{ch} - L_{ch} \frac{d I_{ch}}{dt}$')
ax.legend()
ax.set(xlabel="Time [s]", ylabel="$I_{p}$ [A]", title="Plasma current $I_{{p}}$ #{}".format(shot_no));
Getting $I_p$ mean¶
In [30]:
Ip_mean = Ip.abs().mean()
print_and_save("IpMean", Ip_mean, "%.3f")
Getting $I_p$ max¶
In [31]:
Ip_max = Ip.quantile(0.99)
print_and_save("IpMax", Ip_max, "%.3f")
Plasma detect¶
The plasma detection is based on finding the time points of the greatest peaks in the first (smooth) derivative of $I_p$ which corresponds to the breaks in slope of $I_p$, i.e. the rise or fall to 0, respectively. The smooth derivative is found using the Savitzky-Golay filter.
In [32]:
Ip_detect = Ip.loc[0.0025:] # select time > 2.5 ms to avoid thyristors
dt = (Ip_detect.index[-1] - Ip_detect.index[0]) / (Ip_detect.index.size) # estimated sampling step
print('Sampling step {dt:.3g} s, sampling frequency {fs:.3g} Hz'.format(dt=dt, fs=1/dt))
window_length = int(0.5e-3/dt) # window fit over 1 ms
if window_length % 2 == 0: # must not be odd
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 # [MA/s]
threshold = 0.05
plasma_start = dIp[dIp > dIp.max()*threshold].index[0] # plasma actually starts before the sharpest rise
plasma_end = dIp.idxmin() # due to inductance smoothing the greatest fall is the best predictor
is_plasma = int(plasma_start < plasma_end and dIp.abs().max() > 0.5)
if not is_plasma:
plasma_start = plasma_end = np.nan
print_and_save("Plasma", is_plasma, "%.3f")
print_and_save("PlasmaStart", plasma_start, "%.3f")
print_and_save("PlasmaEnd", plasma_end, "%.3f")
print_and_save("PlasmaTimeLength", plasma_end - plasma_start, "%.3f")
In [33]:
fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)
dIp.plot(grid=True, ax=axs[0])
axs[0].set(xlabel='', ylabel='$dI_p/dt$ [MA/s]', title='Plasma detection from $d I_p/dt$ in #{}'.format(shot_no))
Ip.plot(grid=True, ax=axs[1])
axs[1].set(ylabel="$I_{p}$ [A]")
loop_voltage.plot(grid=True, ax=axs[2])
axs[2].set(xlabel="Time [s]", ylabel="$U_{l}$ [V]")
for ax in axs:
for t in (plasma_start, plasma_end):
ax.axvline(t, color='k', linestyle='--', linewidth=2)
Overview graphs and parameters¶
In [34]:
df_processed = pd.concat(
[loop_voltage.rename('Ul'), Bt, Ip*1e-3], axis='columns')
df_processed.index = df_processed.index * 1e3 # to ms
df_processed.head()
Out[34]:
In [35]:
plasma_lines = hv.VLine(plasma_start*1e3) * hv.VLine(plasma_end*1e3)
layout = hv.Layout([df_processed[v].hvplot.line(
xlabel='', ylabel=l, legend=False, title='', grid=True).relabel(group=v) * plasma_lines
for (v, l) in [('Ul', 'Uₗ [V]'), ('Bt', 'Bₜ [T]'), ('Ip', 'Iₚ [kA]')] ])
layout.cols(1).opts(hv.opts.Curve(width=600, height=200),
hv.opts.VLine(color='black', alpha=0.7, line_dash='dashed'),
hv.opts.Curve('Ip', xlabel='time [ms]'))
In [36]:
fname = 'basig_diagnostics_processed.csv'
In [37]:
df_processed.to_csv(fname)
In [38]:
Markdown("[Time series in graph in CSV format](./{})".format(fname))
Out[38]:
In [39]:
units = ['V', 'T', 'kA']
In [40]:
plasma_start_ms = plasma_start *1e3
plasma_end_ms = plasma_end *1e3
In [41]:
if is_plasma:
heading = Markdown("### Basic diagnostics values during plasma\n\n"
"plasma lifetime of {:.3g} ms, from {:.3g} ms to {:.3g} ms".format(
plasma_end_ms-plasma_start_ms, plasma_start_ms, plasma_end_ms))
else:
heading = Markdown("### No plasma detected (vacuum discharge)")
heading
Out[41]:
In [42]:
df_during_plasma = df_processed.loc[plasma_start_ms:plasma_end_ms]
df_overview = pd.DataFrame()
for n in ('mean', 'min', 'max'):
df_overview[n] = df_during_plasma.aggregate(n)
df_overview['start'] = df_during_plasma.iloc[0]
df_overview['end'] = df_during_plasma.iloc[-1]
In [43]:
df_overview.insert(0,'units', units)
df_overview.T
Out[43]:
End of overview row
TODO¶
- filter out thyristor noise in coils during start-up
- handle polarities (should only some saved scalars use
abs()?)
In [ ]: