Create plasma current class

documentation/source/physics-models/plasma_current/plasma_current.md

@@ -1,4 +1,4 @@
-# Plasma Current 
+# Plasma Current | `PlasmaCurrent`
 
 ## Overview
 
@@ -92,8 +92,7 @@ Instead of $q_a$, $q_{95}$ is used as in plasma configurations with divertors th
 
 ## Plasma Current Calculation | `calculate_plasma_current()`
 
-This function calculates the plasma current shaping factor ($f_q$), plasma current ($I_{\text{p}}$), qstar ($q^*$) and then poloidal field and the profile settings for $\texttt{alphaj}$ ($\alpha_J$) and $\texttt{ind_plasma_internal_norm}$ ($l_{\text{i}}$). This is done in 5 separate steps which are shown in the following numbered sections.
-
+This function calculates the plasma current shaping factor ($f_q$) and then plasma current ($I_{\text{p}}$). 
 
 $$\begin{aligned}
 I_{\text{p}} = f_q \frac{2\pi}{\mu_0}  \frac{a^2 B_{\text{T}}}{R \ q_{95}}
@@ -490,28 +489,28 @@ $$
 
 -----------------------------
 
-### 2. Calculate the cylindrical safety factor
+## Calculate the cylindrical safety factor | `calculate_cylindrical_safety_factor()`
 
 In reactor design the total plasma current $I_{\text{p}}$, is limited by considerations of tearing mode stability, and major disruptions. These considerations place a limit on the safety factor at the plasma boundary, $q_{\psi}$. In a cylindrical plasma model there is a simple relation between this value, $q_{\text{cyl}}$, and the plasma current, given by
 
 $$
-q_{\psi} = q_{\text{cyl}} \equiv \frac{5a^2B_0}{R_0I_p}
+q_{\psi} = q_{\text{cyl}} \equiv \frac{2\pi}{\mu_0}\frac{a^2B_0}{R_0I_p}
 $$
 
 where the minor and major radii are expressed in metres, the toroidal magnetic field $B_0$, in Tesla, and the plasma current $I_{\text{p}}$, in megaamps. Thus a stability requirement of the form  $q_{\psi} > q_{\text{crit}}$ (usually considered to be ~ 2.0) places a limit on $I_{\text{p}}$.
 The effect of toroidicity and of shaping of the plasma cross section is to modify the relation above between $q_{\psi}$ and the plasma current, $I_{\text{p}}$. In general this relation may be written in the form:
 
 
 $$
-q_{\psi} =  \frac{5a^2B_0}{R_0I_p}F\left(\epsilon,\kappa,\delta; \text{profiles}\right)
+q_{\psi} =  \frac{2\pi}{\mu_0}\frac{a^2B_0}{R_0I_p}F\left(\epsilon,\kappa,\delta; \text{profiles}\right)
 $$
 
 It is not possible to derive a general analytic expression for $F$, so it has been customary to employ an empirical, or semi-empirical expression involving $\epsilon$,$\kappa$ and $\delta$, chosen to fit data taken from numerical solutions of the Grad Shafranov equation.
 
 The cylindrical equivalent $q$ definition used currently is the one given below from IPDG89[^5]
 
 $$
-q_* = \frac{5 \times 10^6a^2B_T}{RI_{\text{p}}}\frac{(1+\kappa_{95}^2(1+2\delta_{95}^2-1.2\delta_{95}^3))}{2}
+q_* = \frac{2\pi}{\mu_0}\frac{a^2B_T}{RI_{\text{p}}}\frac{(1+\kappa_{95}^2(1+2\delta_{95}^2-1.2\delta_{95}^3))}{2}
 $$
 
 
@@ -524,9 +523,9 @@ $$
 --------------
 
 
-### 3. Plasma Current Poloidal Field
+## Plasma Current Poloidal Field | `calculate_poloidal_field()`
 
-For calculating the poloidal magnetic field created due to the presence of the plasma current, [Ampere's law](https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law) can be used. In this case the poloidal field is simply returned as:
+For calculating the poloidal magnetic field created due to the presence of the plasma current, [Ampere's law](https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law) can be used. In this case the plasma edge average poloidal field is simply returned as:
 
 $$
 B_{\text{p}} = \frac{\mu_0 I_{\text{p}}}{\texttt{len_plasma_poloidal}}

process/core/io/plot_proc.py

@@ -66,6 +66,7 @@
 from process.models.physics.current_drive import ElectronBernstein, ElectronCyclotron
 from process.models.physics.impurity_radiation import read_impurity_file
 from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel
+from process.models.physics.plasma_current import PlasmaCurrent, PlasmaCurrentModel
 from process.models.tfcoil.superconducting import SUPERCONDUCTING_TF_TYPES
 
 mpl.rcParams["figure.max_open_warning"] = 40
@@ -3027,17 +3028,17 @@ def plot_main_plasma_information(
 
     # Add plasma current information
     textstr_currents = (
-        f"          $\\mathbf{{Plasma\\ currents:}}$\n\n"
-        f"          Plasma current {mfile.get('plasma_current_ma', scan=scan):.4f} MA\n"
-        f"            - Bootstrap fraction {mfile.get('f_c_plasma_bootstrap', scan=scan):.4f}\n"
-        f"            - Diamagnetic fraction {mfile.get('f_c_plasma_diamagnetic', scan=scan):.4f}\n"
-        f"            - Pfirsch-Schlüter fraction {mfile.get('f_c_plasma_pfirsch_schluter', scan=scan):.4f}\n"
-        f"            - Auxiliary fraction {mfile.get('f_c_plasma_auxiliary', scan=scan):.4f}\n"
-        f"            - Inductive fraction {mfile.get('f_c_plasma_inductive', scan=scan):.4f}"
+        f"$\\mathbf{{Plasma\\ currents:}}$\n\n"
+        f"Plasma current ({PlasmaCurrentModel(int(mfile.get('i_plasma_current', scan=scan))).full_name}): {mfile.get('plasma_current_ma', scan=scan):.4f} MA\n"
+        f"  - Bootstrap fraction {mfile.get('f_c_plasma_bootstrap', scan=scan):.4f}\n"
+        f"  - Diamagnetic fraction {mfile.get('f_c_plasma_diamagnetic', scan=scan):.4f}\n"
+        f"  - Pfirsch-Schlüter fraction {mfile.get('f_c_plasma_pfirsch_schluter', scan=scan):.4f}\n"
+        f"  - Auxiliary fraction {mfile.get('f_c_plasma_auxiliary', scan=scan):.4f}\n"
+        f"  - Inductive fraction {mfile.get('f_c_plasma_inductive', scan=scan):.4f}"
     )
 
     axis.text(
-        0.765,
+        0.72,
         0.975,
         textstr_currents,
         fontsize=9,
@@ -3053,8 +3054,8 @@ def plot_main_plasma_information(
 
     # Add plasma current label
     axis.text(
-        0.78,
         0.925,
+        0.9,
         "$I_{\\text{p}} $",
         fontsize=23,
         verticalalignment="top",
@@ -13554,23 +13555,24 @@ def main_plot(
     plot_ebw_ecrh_coupling_graph(figs[12].add_subplot(111), m_file, scan)
 
     plot_bootstrap_comparison(figs[13].add_subplot(221), m_file, scan)
-    plot_h_threshold_comparison(figs[13].add_subplot(224), m_file, scan)
+    PlasmaCurrent.plot_plasma_current_comparison(figs[13].add_subplot(224), m_file, scan)
+    plot_h_threshold_comparison(figs[14].add_subplot(224), m_file, scan)
     plot_density_limit_comparison(figs[14].add_subplot(221), m_file, scan)
-    plot_confinement_time_comparison(figs[14].add_subplot(224), m_file, scan)
+    plot_confinement_time_comparison(figs[15].add_subplot(224), m_file, scan)
 
-    plot_debye_length_profile(figs[15].add_subplot(232), m_file, scan)
-    plot_velocity_profile(figs[15].add_subplot(233), m_file, scan)
-    plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan)
+    plot_debye_length_profile(figs[16].add_subplot(232), m_file, scan)
+    plot_velocity_profile(figs[16].add_subplot(233), m_file, scan)
+    plot_plasma_coloumb_logarithms(figs[16].add_subplot(231), m_file, scan)
 
-    plot_electron_frequency_profile(figs[15].add_subplot(212), m_file, scan)
+    plot_electron_frequency_profile(figs[16].add_subplot(212), m_file, scan)
 
-    plot_ion_frequency_profile(figs[16].add_subplot(311), m_file, scan)
+    plot_ion_frequency_profile(figs[17].add_subplot(311), m_file, scan)
 
-    plot_larmor_radius_profile(figs[16].add_subplot(313), m_file, scan)
+    plot_larmor_radius_profile(figs[17].add_subplot(313), m_file, scan)
 
     # Plot poloidal cross-section
     poloidal_cross_section(
-        figs[17].add_subplot(121, aspect="equal"),
+        figs[18].add_subplot(121, aspect="equal"),
         m_file,
         scan,
         demo_ranges,
@@ -13580,90 +13582,90 @@ def main_plot(
 
     # Plot toroidal cross-section
     toroidal_cross_section(
-        figs[17].add_subplot(122, aspect="equal"),
+        figs[18].add_subplot(122, aspect="equal"),
         m_file,
         scan,
         demo_ranges,
         colour_scheme,
     )
 
     # Plot color key
-    ax17 = figs[17].add_subplot(222)
+    ax17 = figs[18].add_subplot(222)
     ax17.set_position([0.5, 0.5, 0.5, 0.5])
     color_key(ax17, m_file, scan, colour_scheme)
 
-    ax18 = figs[18].add_subplot(211)
+    ax18 = figs[19].add_subplot(211)
     ax18.set_position([0.1, 0.33, 0.8, 0.6])
     plot_radial_build(ax18, m_file, colour_scheme)
 
     # Make each axes smaller vertically to leave room for the legend
-    ax185 = figs[19].add_subplot(211)
+    ax185 = figs[20].add_subplot(211)
     ax185.set_position([0.1, 0.61, 0.8, 0.32])
 
-    ax18b = figs[19].add_subplot(212)
+    ax18b = figs[20].add_subplot(212)
     ax18b.set_position([0.1, 0.13, 0.8, 0.32])
     plot_upper_vertical_build(ax185, m_file, colour_scheme)
     plot_lower_vertical_build(ax18b, m_file, colour_scheme)
 
     # Can only plot WP and turn structure if superconducting coil at the moment
     if m_file.get("i_tf_sup", scan=scan) == 1:
         # TF coil with WP
-        ax19 = figs[20].add_subplot(221, aspect="equal")
+        ax19 = figs[21].add_subplot(221, aspect="equal")
         ax19.set_position([
             0.025,
             0.45,
             0.5,
             0.5,
         ])  # Half height, a bit wider, top left
-        plot_superconducting_tf_wp(ax19, m_file, scan, figs[20])
+        plot_superconducting_tf_wp(ax19, m_file, scan, figs[21])
 
         # TF coil turn structure
-        ax20 = figs[21].add_subplot(325, aspect="equal")
+        ax20 = figs[22].add_subplot(325, aspect="equal")
         ax20.set_position([0.025, 0.5, 0.4, 0.4])
-        plot_tf_cable_in_conduit_turn(ax20, figs[21], m_file, scan)
-        plot_205 = figs[21].add_subplot(223, aspect="equal")
+        plot_tf_cable_in_conduit_turn(ax20, figs[22], m_file, scan)
+        plot_205 = figs[22].add_subplot(223, aspect="equal")
         plot_205.set_position([0.075, 0.1, 0.3, 0.3])
-        plot_cable_in_conduit_cable(plot_205, figs[21], m_file, scan)
+        plot_cable_in_conduit_cable(plot_205, figs[22], m_file, scan)
     else:
-        ax19 = figs[20].add_subplot(211, aspect="equal")
+        ax19 = figs[21].add_subplot(211, aspect="equal")
         ax19.set_position([0.06, 0.55, 0.675, 0.4])
-        plot_resistive_tf_wp(ax19, m_file, scan, figs[20])
-        plot_resistive_tf_info(ax19, m_file, scan, figs[20])
+        plot_resistive_tf_wp(ax19, m_file, scan, figs[21])
+        plot_resistive_tf_info(ax19, m_file, scan, figs[21])
     plot_tf_coil_structure(
-        figs[22].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme
+        figs[23].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme
     )
 
-    plot_plasma_outboard_toroidal_ripple_map(figs[23], m_file, scan)
+    plot_plasma_outboard_toroidal_ripple_map(figs[24], m_file, scan)
 
-    plot_tf_stress(figs[24].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file)
+    plot_tf_stress(figs[25].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file)
 
-    plot_current_profiles_over_time(figs[25].add_subplot(111), m_file, scan)
+    plot_current_profiles_over_time(figs[26].add_subplot(111), m_file, scan)
 
     plot_cs_coil_structure(
-        figs[26].add_subplot(121, aspect="equal"), figs[26], m_file, scan
+        figs[27].add_subplot(121, aspect="equal"), figs[27], m_file, scan
     )
     plot_cs_turn_structure(
-        figs[26].add_subplot(224, aspect="equal"), figs[26], m_file, scan
+        figs[27].add_subplot(224, aspect="equal"), figs[27], m_file, scan
     )
 
     plot_first_wall_top_down_cross_section(
-        figs[27].add_subplot(221, aspect="equal"), m_file, scan
+        figs[28].add_subplot(221, aspect="equal"), m_file, scan
     )
-    plot_first_wall_poloidal_cross_section(figs[27].add_subplot(122), m_file, scan)
-    plot_fw_90_deg_pipe_bend(figs[27].add_subplot(337), m_file, scan)
+    plot_first_wall_poloidal_cross_section(figs[28].add_subplot(122), m_file, scan)
+    plot_fw_90_deg_pipe_bend(figs[28].add_subplot(337), m_file, scan)
 
-    plot_blkt_pipe_bends(figs[28], m_file, scan)
-    ax_blanket = figs[28].add_subplot(122, aspect="equal")
-    plot_blkt_structure(ax_blanket, figs[28], m_file, scan, radial_build, colour_scheme)
+    plot_blkt_pipe_bends(figs[29], m_file, scan)
+    ax_blanket = figs[29].add_subplot(122, aspect="equal")
+    plot_blkt_structure(ax_blanket, figs[29], m_file, scan, radial_build, colour_scheme)
 
     plot_main_power_flow(
-        figs[29].add_subplot(111, aspect="equal"), m_file, scan, figs[29]
+        figs[30].add_subplot(111, aspect="equal"), m_file, scan, figs[30]
     )
 
-    ax24 = figs[30].add_subplot(111)
+    ax24 = figs[31].add_subplot(111)
     # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height
     ax24.set_position([0.08, 0.35, 0.84, 0.57])
-    plot_system_power_profiles_over_time(ax24, m_file, scan, figs[30])
+    plot_system_power_profiles_over_time(ax24, m_file, scan, figs[31])
 
 
 def create_thickness_builds(m_file, scan: int):
@@ -13740,7 +13742,7 @@ def main(args=None):
 
     # create main plot
     # Increase range when adding new page
-    pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(31)]
+    pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(32)]
 
     # run main_plot
     main_plot(

process/data_structure/physics_variables.py

@@ -906,6 +906,31 @@
 plasma_current: float = None
 """plasma current (A)"""
 
+c_plasma_peng_analytic: float = None
+"""Peng analytic plasma current (A)"""
+
+c_plasma_peng_double_null: float = None
+"""Peng double null divertor plasma current (A)"""
+
+c_plasma_cyclindrical: float = None
+"""Cylindrical plasma current (A)"""
+
+c_plasma_ipdg89: float = None
+"""ITER IPDG89 plasma current (A)"""
+
+c_plasma_todd_empirical_i: float = None
+"""Todd empirical plasma current I (A)"""
+
+c_plasma_todd_empirical_ii: float = None
+"""Todd empirical plasma current II (A)"""
+c_plasma_connor_hastie: float = None
+"""Connor-Hastie plasma current (A)"""
+
+c_plasma_sauter: float = None
+"""Sauter plasma current (A)"""
+
+c_plasma_fiesta_st: float = None
+"""FIESTA ST plasma current (A)"""
 
 p_plasma_neutron_mw: float = None
 """Neutron fusion power from just the plasma [MW]"""
@@ -1556,6 +1581,15 @@ def init_physics_variables():
         pflux_fw_rad_mw, \
         pden_ion_electron_equilibration_mw, \
         plasma_current, \
+        c_plasma_peng_analytic, \
+        c_plasma_peng_double_null, \
+        c_plasma_cyclindrical, \
+        c_plasma_ipdg89, \
+        c_plasma_todd_empirical_i, \
+        c_plasma_todd_empirical_ii, \
+        c_plasma_connor_hastie, \
+        c_plasma_sauter, \
+        c_plasma_fiesta_st, \
         p_plasma_neutron_mw, \
         p_neutron_total_mw, \
         pden_neutron_total_mw, \
@@ -1838,6 +1872,15 @@ def init_physics_variables():
     pflux_fw_rad_mw = 0.0
     pden_ion_electron_equilibration_mw = 0.0
     plasma_current = 0.0
+    c_plasma_peng_analytic = 0.0
+    c_plasma_peng_double_null = 0.0
+    c_plasma_cyclindrical = 0.0
+    c_plasma_ipdg89 = 0.0
+    c_plasma_todd_empirical_i = 0.0
+    c_plasma_todd_empirical_ii = 0.0
+    c_plasma_connor_hastie = 0.0
+    c_plasma_sauter = 0.0
+    c_plasma_fiesta_st = 0.0
     p_plasma_neutron_mw = 0.0
     p_neutron_total_mw = 0.0
     pden_neutron_total_mw = 0.0

process/main.py

@@ -104,6 +104,7 @@
     PlasmaBeta,
     PlasmaInductance,
 )
+from process.models.physics.plasma_current import PlasmaCurrent
 from process.models.physics.plasma_geometry import PlasmaGeom
 from process.models.physics.plasma_profiles import PlasmaProfile
 from process.models.power import Power
@@ -710,6 +711,7 @@ def __init__(self):
         )
         self.plasma_confinement = PlasmaConfinementTime()
         self.plasma_transition = PlasmaConfinementTransition()
+        self.plasma_current = PlasmaCurrent()
         self.physics = Physics(
             plasma_profile=self.plasma_profile,
             current_drive=self.current_drive,
@@ -720,6 +722,7 @@ def __init__(self):
             plasma_bootstrap_current=self.plasma_bootstrap_current,
             plasma_confinement=self.plasma_confinement,
             plasma_transition=self.plasma_transition,
+            plasma_current=self.plasma_current,
         )
         self.physics_detailed = DetailedPhysics(
             plasma_profile=self.plasma_profile,