{ "cells": [ { "cell_type": "markdown", "metadata": { "golem": {} }, "source": [ "# Tokamak GOLEM Basic diagnostics\n", "\n", "
\n", "\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Procedure (This notebook to download)
\n", "bash wrapper, Error log" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisities: function definitions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load libraries" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import os\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy import constants, integrate, signal, interpolate\n", "import sqlalchemy # high-level library for SQL in Python\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For interactive web figures" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": {} }, "outputs": [], "source": [ "import holoviews as hv\n", "hv.extension('bokeh')\n", "import hvplot.pandas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For conditional rich-text boxes" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.display import Markdown" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define global constants." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "base_URL = \"http://golem.fjfi.cvut.cz/utils/data/{shot_no}/{identifier}\" #remote access\n", "data_URL = \"http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/StandardDAS/{identifier}.csv\" # TODO workaround\n", "parameters_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Production/Parameters/{identifier}'\n", "\n", "destination=''\n", "#os.makedirs(destination, exist_ok=True );\n", "# try to get thot number form SHOT_NO envirnoment variable, otherwise use the specified one\n", "shot_no = os.environ.get('SHOT_NO', 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `DataSource` downloads and caches data (by full URL) in a temporary directory and makes them accessible as files." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds = np.DataSource(destpath='/tmp')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def print_and_save(phys_quant, value, format_str='%.3f'):\n", " print(phys_quant+\" = %.5f\" % value)\n", " with open(destination+phys_quant, 'w') as f:\n", " f.write(format_str % value)\n", " update_db_current_shot(phys_quant,value) \n", " \n", "def update_db_current_shot(field_name, value):\n", " try:\n", " engine = sqlalchemy.create_engine('postgresql://golem:rabijosille@192.168.2.116/golem_database')\n", " except:\n", " return\n", " engine.execute(f\"\"\"UPDATE shots SET \"{field_name}\"={value} WHERE shot_no IN(SELECT max(shot_no) FROM shots)\"\"\") \n", " \n", " \n", "def open_remote(shot_no, identifier, url_template=base_URL):\n", " return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))\n", "\n", "def read_value(shot_no, identifier):\n", " \"\"\"Return the value for given shot as a number if possible\"\"\"\n", " value = open_remote(shot_no, identifier, base_URL).read()\n", " return pd.to_numeric(value, errors='ignore')\n", "\n", "def read_parameter(shot_no, identifier):\n", " return open_remote(shot_no, identifier, parameters_URL).read().strip()\n", " \n", "def read_signal(shot_no, identifier): \n", " file = open_remote(shot_no, identifier, data_URL)\n", " return pd.read_csv(file, names=['Time',identifier],\n", " index_col='Time', squeeze=True) # squeeze makes simple 1-column signals a Series" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def correct_inf(signal):\n", " \"\"\"Inteprolate Inf values\"\"\"\n", " signal = signal.replace([-np.inf, np.inf], np.nan).interpolate()\n", " return signal" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# TODO hardcoded until properly recorded in database\n", "t_Bt = 1e-3\n", "t_CD = 2e-3\n", "offset_sl = slice(None, min(t_Bt, t_CD) - 1e-4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $U_l$ management" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check the data availability" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "loop_voltage = read_signal(shot_no, 'LoopVoltageCoil_raw')\n", "polarity_CD = read_parameter(shot_no, 'CD_orientation')\n", "if polarity_CD != 'CW': # TODO hardcoded for now!\n", " loop_voltage *= -1 # make positive\n", "loop_voltage = correct_inf(loop_voltage)\n", "loop_voltage.loc[:t_CD] = 0\n", "ax = loop_voltage.plot(grid=True)\n", "ax.set(xlabel=\"Time [s]\", ylabel=\"$U_l$ [V]\", title=\"Loop voltage $U_l$ #{}\".format(shot_no));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $B_t$ calculation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check the data availability\n", "It is as magnetic measurement, so the raw data only give $\\frac{dB_t}{dt}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dBt = read_signal(shot_no,'BtCoil_raw')\n", "polarity_Bt = read_parameter(shot_no, 'Bt_orientation')\n", "if polarity_Bt != 'CW': # TODO hardcoded for now!\n", " dBt *= -1 # make positive\n", "dBt = correct_inf(dBt)\n", "dBt -= dBt.loc[offset_sl].mean()\n", "ax = dBt.plot(grid=True)\n", "ax.set(xlabel=\"Time [s]\", ylabel=\"$dU_{B_t}/dt$ [V]\", title=\"BtCoil_raw signal #{}\".format(shot_no));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integration (it is a magnetic diagnostic) & calibration" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "K_BtCoil = read_value(shot_no, 'K_BtCoil') # Get BtCoil calibration factor\n", "print('BtCoil calibration factor K_BtCoil={} T/(Vs)'.format(K_BtCoil))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Bt = pd.Series(integrate.cumtrapz(dBt, x=dBt.index, initial=0) * K_BtCoil, \n", " index=dBt.index, name='Bt')\n", "ax = Bt.plot(grid=True)\n", "ax.set(xlabel=\"Time [s]\", ylabel=\"$B_t$ [T]\", title=\"Toroidal magnetic field $B_t$ #{}\".format(shot_no));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plasma + Chamber current $I_{p+ch}$ calculation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check the data availability\n", "\n", "Because it is a magnetic measurement, the raw data only gives $\\frac{dI_{p+ch}}{dt}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "polarity_CD" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dIpch = read_signal(shot_no, 'RogowskiCoil_raw')\n", "if polarity_CD == 'CW': # TODO hardcoded for now!\n", " dIpch *= -1 # make positive\n", "dIpch = correct_inf(dIpch)\n", "dIpch -= dIpch.loc[offset_sl].mean() # subtract offset\n", "dIpch.loc[:t_CD] = 0\n", "ax = dIpch.plot(grid=True)\n", "ax.set(xlabel=\"Time [s]\", ylabel=\"$dU_{I_{p+ch}}/dt$ [V]\", title=\"RogowskiCoil_raw signal #{}\".format(shot_no));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integration (it is a magnetic diagnostics) & calibration" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "K_RogowskiCoil = read_value(shot_no, 'K_RogowskiCoil') # Get RogowskiCoil calibration factor\n", "print('RogowskiCoil calibration factor K_RogowskiCoil={} A/(Vs)'.format(K_RogowskiCoil))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Ipch = pd.Series(integrate.cumtrapz(dIpch, x=dIpch.index, initial=0) * K_RogowskiCoil,\n", " index=dIpch.index, name='Ipch')\n", "ax = Ipch.plot(grid=True)\n", "ax.set(xlabel=\"Time [s]\", ylabel=\"$I_{p+ch}$ [A]\", title=\"Total (plasma+chamber) current $I_{{p+ch}}$ #{}\".format(shot_no));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plasma current $I_{p}$ calculation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "R_chamber = read_value(shot_no, 'R_chamber') # Get Chamber resistivity\n", "print('Chamber resistivity R_chamber={} Ohm'.format(R_chamber))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L_chamber = 1.2e-6/2 # H" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The chamber current $I_{ch}$ satisfies the equation (neglecting the mutual inductance with the plasma)\n", "$$U_l = R_{ch} I_{ch} + L_{ch} \\frac{d I_{ch}}{dt}$$\n", "Therefore, the following initial value problem must be solved to take into account the chamber inductance properly\n", "$$\\frac{d I_{ch}}{dt} = \\frac{1}{L_{ch}}\\left( U_l - R_{ch} I_{ch}\\right), \\quad I_{ch}(t=0)=0$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "U_l_func = interpolate.interp1d(loop_voltage.index, loop_voltage) # 1D interpolator\n", "def dIch_dt(t, Ich):\n", " return (U_l_func(t) - R_chamber * Ich) / L_chamber\n", "t_span = loop_voltage.index[[0, -1]]\n", "solution = integrate.solve_ivp(dIch_dt, t_span, [0], t_eval=loop_voltage.index, )\n", "Ich = pd.Series(solution.y[0], index=loop_voltage.index, name='Ich')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Ip_naive = Ipch - loop_voltage/R_chamber # creates a new Series\n", "Ip = Ipch - Ich\n", "Ip.name = 'Ip'\n", "Ip_naive.plot(grid=True, label='naive $I_{ch}=U_l/R_{ch}$')\n", "ax = Ip.plot(grid=True, label=r'using $U_l = R_{ch} I_{ch} - L_{ch} \\frac{d I_{ch}}{dt}$')\n", "ax.legend()\n", "ax.set(xlabel=\"Time [s]\", ylabel=\"$I_{p}$ [A]\", title=\"Plasma current $I_{{p}}$ #{}\".format(shot_no));" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for I in [Ich.rename('$I_{ch}$'), Ipch.rename('$I_{ch}+I_p$'), Ip.rename('$I_p$')]:\n", " ax = I.plot()\n", "ax.legend()\n", "ax.set(xlabel='Time [s]', ylabel='$I$ [A]', title='estimated plasma and chamber current and measured total')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plasma detect\n", "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.\n", "The smooth derivative is found using the Savitzky-Golay filter." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Ip_detect = Ip.loc[0.0025:] # select time > 2.5 ms to avoid thyristors\n", "dt = (Ip_detect.index[-1] - Ip_detect.index[0]) / (Ip_detect.index.size) # estimated sampling step\n", "print('Sampling step {dt:.3g} s, sampling frequency {fs:.3g} Hz'.format(dt=dt, fs=1/dt))\n", "window_length = int(0.5e-3/dt) # window fit over 1 ms\n", "if window_length % 2 == 0: # must not be odd\n", " window_length += 1\n", "dIp = pd.Series(signal.savgol_filter(Ip_detect, window_length, polyorder=3, deriv=1, delta=dt),\n", " name='dIp', index=Ip_detect.index) / 1e6 # [MA/s]\n", "\n", "threshold = 0.05\n", "plasma_start = dIp[dIp > dIp.max()*threshold].index[0] # plasma actually starts before the sharpest rise\n", "plasma_end = dIp.idxmin() # due to inductance smoothing the greatest fall is the best predictor\n", "is_plasma = int(plasma_start < plasma_end and dIp.abs().max() > 0.5)\n", "print_and_save(\"b_discharge\", is_plasma, \"%.3f\")\n", "if not is_plasma:\n", " plasma_start = plasma_end = np.nan\n", " plasma_lifetime = 0\n", " print_and_save(\"t_plasma_start\", -1, \"%.3f\")\n", " print_and_save(\"t_plasma_end\", -1, \"%.3f\")\n", " print_and_save(\"t_discharge_duration\", -1, \"%.3f\")\n", "else:\n", " plasma_lifetime = plasma_end - plasma_start\n", " print_and_save(\"t_plasma_start\", plasma_start*1000, \"%.3f\")\n", " print_and_save(\"t_plasma_end\", plasma_end*1000, \"%.3f\")\n", " print_and_save(\"t_discharge_duration\", plasma_lifetime*1000, \"%.3f\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)\n", "dIp.plot(grid=True, ax=axs[0])\n", "axs[0].set(xlabel='', ylabel='$dI_p/dt$ [MA/s]', title='Plasma detection from $d I_p/dt$ in #{}'.format(shot_no))\n", "Ip.plot(grid=True, ax=axs[1])\n", "axs[1].set(ylabel=\"$I_{p}$ [A]\")\n", "loop_voltage.plot(grid=True, ax=axs[2])\n", "axs[2].set(xlabel=\"Time [s]\", ylabel=\"$U_{l}$ [V]\")\n", "\n", "for ax in axs:\n", " for t in (plasma_start, plasma_end):\n", " ax.axvline(t, color='k', linestyle='--', linewidth=2)\n", " \n", "plt.savefig('icon-fig') " ] }, { "cell_type": "markdown", "metadata": { "golem": { "row": "start" } }, "source": [ "## Overview graphs and parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "col": "col-lg-8", "row": "start" } }, "outputs": [], "source": [ "df_processed = pd.concat(\n", " [loop_voltage.rename('U_loop'), Bt, Ip*1e-3], axis='columns')\n", "df_processed.index = df_processed.index * 1e3 # to ms\n", "df_processed.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "col": "col-12" } }, "outputs": [], "source": [ "if is_plasma:\n", " plasma_lines = hv.VLine(plasma_start*1e3) * hv.VLine(plasma_end*1e3)\n", "else:\n", " plasma_lines = hv.Curve([])\n", "layout = hv.Layout([df_processed[v].hvplot.line(\n", " xlabel='', ylabel=l, legend=False, title='', grid=True, group_label=v) * plasma_lines\n", " for (v, l) in [('U_loop', 'Uₗ [V]'), ('Bt', 'Bₜ [T]'), ('Ip', 'Iₚ [kA]')] ])\n", "\n", "plot=layout.cols(1).opts(hv.opts.Curve(width=600, height=200),\n", " hv.opts.VLine(color='black', alpha=0.7, line_dash='dashed'),\n", " hv.opts.Curve('Ip', xlabel='time [ms]'))\n", "hvplot.save(plot, 'homepage_figure.html')\n", "plot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fname = 'basig_diagnostics_processed.csv'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df_processed.to_csv(fname)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "col": "col-12" } }, "outputs": [], "source": [ "Markdown(\"[Time series in graph in CSV format](./{})\".format(fname))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "row": "end" } }, "outputs": [], "source": [ "units = ['V', 'T', 'kA']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "col": "col-lg-4", "row": "start" } }, "outputs": [], "source": [ "plasma_start_ms = plasma_start *1e3\n", "plasma_end_ms = plasma_end *1e3\n", "plasma_lifetime_ms = plasma_lifetime * 1e3" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "col": "col-12" } }, "outputs": [], "source": [ "if is_plasma:\n", " heading = Markdown(\"### Basic diagnostics values during plasma\\n\\n\"\n", " \"plasma lifetime of {:.3g} ms, from {:.3g} ms to {:.3g} ms\".format(\n", " plasma_lifetime_ms, plasma_start_ms, plasma_end_ms))\n", "else:\n", " heading = Markdown(\"### No plasma detected (vacuum discharge)\")\n", "heading" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if is_plasma:\n", " plasma_sl = slice(plasma_start_ms, plasma_end_ms)\n", "else:\n", " plasma_sl = slice() # TODO really use whole discharge ?\n", "df_during_plasma = df_processed.loc[plasma_sl]\n", "df_overview = df_during_plasma.quantile([0.01, 0.5, 0.99]) # use quantiles to skip peaks\n", "df_overview.index = ['min', 'mean', 'max'] # actually quantiles, but easier to understand\n", "if is_plasma:\n", " df_overview.loc['start'] = df_during_plasma.iloc[0]\n", " df_overview.loc['end'] = df_during_plasma.iloc[-1]\n", "else:\n", " df_overview.loc['start'] = df_overview.loc['end'] = np.nan\n", "df_overview.loc['units'] = units\n", "# make units row first\n", "df_overview = df_overview.iloc[np.roll(np.arange(df_overview.shape[0]), 1)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "golem": { "col": "col-12" } }, "outputs": [], "source": [ "df_overview" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for agg in ('mean', 'max'):\n", " for quantity, value in df_overview.loc[agg].iteritems():\n", " print_and_save(quantity+'_'+agg, value)" ] }, { "cell_type": "markdown", "metadata": { "golem": { "row": "end" } }, "source": [ "End of overview row" ] } ], "metadata": { "celltoolbar": "Edit Metadata", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" } }, "nbformat": 4, "nbformat_minor": 4 }