🚧 Implement new suite of bootstrap current scalings

documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md

@@ -7,7 +7,7 @@ Some more info can be found [here](https://wiki.fusion.ciemat.es/wiki/Bootstrap_
 
 The fraction of the plasma current provided by the bootstrap effect
 can be either input into the code directly, or calculated using one of five
-methods, as summarised here. Note that methods `i_bootstrap_current = 1-3 & 5` do not take into account the 
+methods, as summarised here. Note that methods `i_bootstrap_current = 1-3 & 5-6` do not take into account the
 existence of pedestals, whereas the Sauter et al. scaling 
 (`i_bootstrap_current = 4`) allows general profiles to be used. 
 
@@ -460,7 +460,7 @@ The $-1$ subscript in this case refers to the value of the variable in the previ
 
 ### Sakai Scaling | `bootstrap_fraction_sakai()`
 
-Is selected by setting `i_bootstrap_current = 5`[^5]
+Is selected by setting `i_bootstrap_current = 5`[^6]
 
 $$
 f_{\text{BS}} = 10^{0.951 \epsilon - 0.948} \cdot \beta_p^{1.226 \epsilon + 1.584} \cdot l_i^{-0.184\epsilon - 0.282} \cdot \left(\frac{q_{95}}{q_0}\right)^{-0.042 \epsilon - 0.02} \\
@@ -469,6 +469,102 @@ $$
 
 The model includes the toroidal diamagnetic current in the calculation due to the dataset, so `i_diamagnetic_current = 0` can only be used with it
 
+-------------------
+
+### ARIES Scaling | `bootstrap_fraction_aries()`
+
+Is selected by setting `i_bootstrap_current = 6`[^8]
+
+The source reference[^8] does not provide any info about the derivation of the formula. It is only stated like that shown below.
+
+$$
+a_1 = 1.10-1.165l_{\text{i}}+0.47l_{\text{i}}^2
+$$
+
+$$
+b_1 = 0.806-0.885l_{\text{i}}+0.297l_{\text{i}}^2
+$$
+
+$$
+C_{\text{BS}} = a_1+b_1\left(\frac{n(0)}{\langle n \rangle}\right)
+$$
+
+$$
+f_{\text{BS}} = C_{\text{BS}} \sqrt{\epsilon}\beta_{\text{p}}
+$$
+
+---------------------
+
+### Andrade Scaling | `bootstrap_fraction_andrade()`
+
+Is selected by setting `i_bootstrap_current = 7`[^9]
+
+$$
+C_{\text{BS}} = 0.2340 \pm 0.0007
+$$
+
+$$
+c_p = \frac{p_0}{\langle p \rangle}
+$$
+
+$$
+f_{\text{BS}} = C_{\text{BS}} \sqrt{\epsilon}\beta_{\text{p}}c_p^{0.8}
+$$
+
+---------------------
+
+### Hoang Scaling | `bootstrap_fraction_hoang()`
+
+Is selected by setting `i_bootstrap_current = 8`[^10]
+
+$$
+C_{\text{BS}} = \sqrt{\frac{\alpha_{\text{p}}}{\alpha_{\text{j}}}}
+$$
+
+$$
+f_{\text{BS}} = 0.4C_{\text{BS}} \sqrt{\epsilon}\beta_{\text{p}}^{0.9}
+$$
+
+---------------------
+
+### Wong Scaling | `bootstrap_fraction_wong()`
+
+Is selected by setting `i_bootstrap_current = 9`[^11]
+
+$$
+C_{\text{BS}} = 0.773+0.019\kappa
+$$
+
+$$
+f_{\text{peak}} = \left(\int_0^1 \left(1-\rho^2 \right)^{\alpha_{\text{T}}} \left(1-\rho^2 \right)^{\alpha_{\text{n}}} \ \ \mathrm{d\rho}\right)^{-1}
+$$
+
+The integral above is set by the definite solution below
+
+$$
+\frac{\operatorname{B}\left(\frac{1}{2},{\alpha}_{n} + {\alpha}_{T} + 1\right)}{2}
+$$
+
+Assuming that $\alpha_{\text{n}} + \alpha_{\text{T}} > 0,   \alpha_{\text{n}} + \alpha_{\text{T}} + 1 > 0$
+
+$$
+f_{\text{BS}} = C_{\text{BS}} f_{\text{peak}}^{0.25}\beta_{\text{p}}\sqrt{\epsilon}
+$$
+
+---------------------
+
+### Gi Scaling | `bootstrap_fraction_gi()`
+
+Is selected by setting `i_bootstrap_current = 10`[^12]
+
+$$
+C_{\text{BS}} = 0.474 \epsilon^{-0.1} \alpha_{\text{p}}^{0.974} \alpha_{\text{T}}^{-0.416} Z_{\text{eff}}^{0.178} \left(\frac{q_{95}}{q_0}\right)^{-0.133}
+$$
+
+$$
+f_{\text{BS}} = C_{\text{BS}} \beta_{\text{p}}\sqrt{\epsilon}
+$$
+
 ---------------------
 
 ## Setting of maximum desirable bootstrap current fraction
@@ -500,5 +596,8 @@ Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN
 [^5]: O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: “Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime” [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052  
 [^6]: Ryosuke Sakai, Takaaki Fujita, Atsushi Okamoto, Derivation of bootstrap current fraction scaling formula for 0-D system code analysis, Fusion Engineering and Design, Volume 149, 2019, 111322, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2019.111322.
 [^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 
-
-
+[^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015.
+[^9]: M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, doi: https://doi.org/10.1088/0741-3335/50/6/065001.
+[^10]: G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, Jan. 1997, doi: https://doi.org/10.1063/1.53414.
+[^11]: C.-P. Wong, J. C. Wesley, R. D. Stambaugh, and E. T. Cheng, “Toroidal reactor designs as a function of aspect ratio and elongation,” vol. 42, no. 5, pp. 547–556, May 2002, doi: https://doi.org/10.1088/0029-5515/42/5/307.
+[^12]: K. Gi, M. Nakamura, Kenji Tobita, and Y. Ono, “Bootstrap current fraction scaling for a tokamak reactor design study,” Fusion Engineering and Design, vol. 89, no. 11, pp. 2709–2715, Aug. 2014, doi: https://doi.org/10.1016/j.fusengdes.2014.07.009.

documentation/proc-pages/physics-models/profiles/plasma_profiles.md

@@ -407,10 +407,17 @@ $$
 With the coefficients used to turn the temperature from $\text{keV}$ back to Joules.
 
 Since we multiply the density and temperature profiles together to get the pressure we can find the pressure profile factor by adding the profile factors for temperature and density.
+
 $$
 \alpha_p = \alpha_n + \alpha_T
 $$
 
+The volume averaged pressure can then be set if we assume the pressure also has a parabolic profile. Using the standard relation used for both density and temeprature we can set the volume averaged pressure as:
+
+$$
+\langle p \rangle =   \frac{p_0}{\alpha_p+1}
+$$
+
 ???+ note "Pressure profile factor"
 
     The calculation of $\alpha_p$ is only valid assuming a parabolic profile case. The calculatio of $p_0$ is still true as the core values are calculated independantly for each profile type. 

process/io/plot_proc.py

@@ -1876,14 +1876,14 @@ def plot_tf_wp(axis, mfile_data, scan: int) -> None:
                 )
             )
 
-        plt.minorticks_on()
-        plt.xlim(0.0, r_tf_inboard_out * 1.1)
-        plt.ylim((y14[-1] * 1.25), (-y14[-1] * 1.25))
+        axis.minorticks_on()
+        axis.set_xlim(0.0, r_tf_inboard_out * 1.1)
+        axis.set_ylim((y14[-1] * 1.25), (-y14[-1] * 1.25))
 
-        plt.title("Top-down view of inboard TF coil at midplane")
-        plt.xlabel("Radial distance [m]")
-        plt.ylabel("Toroidal distance [m]")
-        plt.legend(bbox_to_anchor=(0.0, -0.25), loc="upper left")
+        axis.set_title("Top-down view of inboard TF coil at midplane")
+        axis.set_xlabel("Radial distance [m]")
+        axis.set_ylabel("Toroidal distance [m]")
+        axis.legend(bbox_to_anchor=(0.0, -0.25), loc="upper left")
 
 
 def plot_tf_turn(axis, mfile_data, scan: int) -> None:
@@ -1974,8 +1974,8 @@ def plot_tf_turn(axis, mfile_data, scan: int) -> None:
                 edgecolor="black",
             ),
         )
-        plt.xlim(-turn_width * 0.05, turn_width * 1.05)
-        plt.ylim(-turn_width * 0.05, turn_width * 1.05)
+        axis.set_xlim(-turn_width * 0.05, turn_width * 1.05)
+        axis.set_ylim(-turn_width * 0.05, turn_width * 1.05)
 
     # Non square turns
     elif integer_turns == 1:
@@ -2026,14 +2026,14 @@ def plot_tf_turn(axis, mfile_data, scan: int) -> None:
             ),
         )
 
-        plt.xlim(-turn_width * 0.05, turn_width * 1.05)
-        plt.ylim(-turn_height * 0.05, turn_height * 1.05)
+        axis.set_xlim(-turn_width * 0.05, turn_width * 1.05)
+        axis.set_ylim(-turn_height * 0.05, turn_height * 1.05)
 
-    plt.minorticks_on()
-    plt.title("WP Turn Structure")
-    plt.xlabel("X [mm]")
-    plt.ylabel("Y [mm]")
-    plt.legend(bbox_to_anchor=(0.0, -0.25), loc="upper left")
+    axis.minorticks_on()
+    axis.set_title("WP Turn Structure")
+    axis.set_xlabel("X [mm]")
+    axis.set_ylabel("Y [mm]")
+    axis.legend(loc="upper right", bbox_to_anchor=(1.0, -0.25))
 
 
 def plot_pf_coils(axis, mfile_data, scan, colour_scheme):
@@ -2898,10 +2898,102 @@ def plot_current_drive_info(axis, mfile_data, scan):
     plot_info(axis, data, mfile_data, scan)
 
 
+def plot_bootstrap_comparison(axis, mfile_data, scan):
+    """Function to plot a scatter box plot of bootstrap current fractions.
+
+    Arguments:
+        axis --> axis object to plot to
+        mfile_data --> MFILE data object
+        scan --> scan number to use
+    """
+
+    boot_ipdg = mfile_data.data["bscf_iter89"].get_scan(scan)
+    boot_sauter = mfile_data.data["bscf_sauter"].get_scan(scan)
+    boot_nenins = mfile_data.data["bscf_nevins"].get_scan(scan)
+    boot_wilson = mfile_data.data["bscf_wilson"].get_scan(scan)
+    boot_sakai = mfile_data.data["bscf_sakai"].get_scan(scan)
+    boot_aries = mfile_data.data["bscf_aries"].get_scan(scan)
+    boot_andrade = mfile_data.data["bscf_andrade"].get_scan(scan)
+    boot_hoang = mfile_data.data["bscf_hoang"].get_scan(scan)
+    boot_wong = mfile_data.data["bscf_wong"].get_scan(scan)
+    boot_gi = mfile_data.data["bscf_gi"].get_scan(scan)
+
+    # Data for the box plot
+    data = [
+        boot_ipdg,
+        boot_sauter,
+        boot_nenins,
+        boot_wilson,
+        boot_sakai,
+        boot_aries,
+        boot_andrade,
+        boot_hoang,
+        boot_wong,
+        boot_gi,
+    ]
+    labels = [
+        "IPDG",
+        "Sauter",
+        "Nevins",
+        "Wilson",
+        "Sakai",
+        "ARIES",
+        "Andrade",
+        "Hoang",
+        "Wong",
+        "Gi",
+    ]
+
+    x = np.ones(len(data))
+
+    # Create the violin plot
+    axis.violinplot(data, showextrema=False)
+
+    # Create the box plot
+    axis.boxplot(data, showfliers=True, showmeans=True, meanline=True, widths=0.3)
+
+    # Scatter plot for each data point
+    colors = plt.cm.plasma(np.linspace(0, 1, len(data)))
+    for i, value in enumerate(data):
+        axis.scatter(x[i], value, color=colors[i], label=labels[i], alpha=1.0)
+    axis.legend(loc="upper left", bbox_to_anchor=(1, 1))
+
+    # Calculate average, standard deviation, and median
+    avg_bootstrap = np.mean(data)
+    std_bootstrap = np.std(data)
+    median_bootstrap = np.median(data)
+
+    # Plot average, standard deviation, and median as text
+    axis.text(
+        0.6, 0.9, f"Average: {avg_bootstrap:.4f}", transform=axis.transAxes, fontsize=9
+    )
+    axis.text(
+        0.6,
+        0.85,
+        f"Standard Dev: {std_bootstrap:.4f}",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+    axis.text(
+        0.6,
+        0.8,
+        f"Median: {median_bootstrap:.4f}",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+
+    axis.set_title("Bootstrap Current Fraction Comparison")
+    axis.set_ylabel("Bootstrap Current Fraction")
+    axis.set_xlim([0.5, 1.5])
+    axis.set_xticks([])
+    axis.set_xticklabels([])
+
+
 def main_plot(
     fig1,
     fig2,
     fig3,
+    fig4,
     m_file_data,
     scan,
     imp="../data/lz_non_corona_14_elements/",
@@ -2988,7 +3080,7 @@ def main_plot(
     plot_current_drive_info(plot_6, m_file_data, scan)
     fig1.subplots_adjust(wspace=0.25, hspace=0.25)
 
-    # Can only plot WP and turn sturcutre if superconducting coil at the moment
+    # Can only plot WP and turn structure if superconducting coil at the moment
     if m_file_data.data["i_tf_sup"].get_scan(scan) == 1:
         # TF coil with WP
         plot_7 = fig3.add_subplot(321)
@@ -2998,6 +3090,9 @@ def main_plot(
         plot_8 = fig3.add_subplot(322, aspect="equal")
         plot_tf_turn(plot_8, m_file_data, scan)
 
+    plot_9 = fig4.add_subplot(221)
+    plot_bootstrap_comparison(plot_9, m_file_data, scan)
+
 
 def main(args=None):
     # TODO The use of globals here isn't ideal, but is required to get main()
@@ -3252,12 +3347,14 @@ def main(args=None):
     page1 = plt.figure(figsize=(12, 9), dpi=80)
     page2 = plt.figure(figsize=(12, 9), dpi=80)
     page3 = plt.figure(figsize=(12, 9), dpi=80)
+    page4 = plt.figure(figsize=(12, 9), dpi=80)
 
     # run main_plot
     main_plot(
         page1,
         page2,
         page3,
+        page4,
         m_file,
         scan=scan,
         demo_ranges=demo_ranges,
@@ -3269,6 +3366,7 @@ def main(args=None):
         pdf.savefig(page1)
         pdf.savefig(page2)
         pdf.savefig(page3)
+        pdf.savefig(page4)
 
     # show fig if option used
     if args.show:
@@ -3277,6 +3375,7 @@ def main(args=None):
     plt.close(page1)
     plt.close(page2)
     plt.close(page3)
+    plt.close(page4)
 
 
 if __name__ == "__main__":