diff --git a/documentation/proc-pages/images/plot_proc_3.PNG b/documentation/proc-pages/images/plot_proc_3.PNG index e026084ca2..d50faea82f 100644 Binary files a/documentation/proc-pages/images/plot_proc_3.PNG and b/documentation/proc-pages/images/plot_proc_3.PNG differ diff --git a/documentation/proc-pages/images/plot_proc_4.PNG b/documentation/proc-pages/images/plot_proc_4.PNG new file mode 100644 index 0000000000..0504ef94d7 Binary files /dev/null and b/documentation/proc-pages/images/plot_proc_4.PNG differ diff --git a/documentation/proc-pages/io/utilities.md b/documentation/proc-pages/io/utilities.md index ddc559b860..1a5e8b9caa 100644 --- a/documentation/proc-pages/io/utilities.md +++ b/documentation/proc-pages/io/utilities.md @@ -69,7 +69,7 @@ A `.csv` file will be saved to the directory of the input file. > `process/io/plot_proc.py` -A utility to produce a three-page PDF summary of the output from PROCESS, including the major parameters, poloidal and toroidal cross-sections, temperature and density profiles and TF coil layout and turn strucuture. +A utility to produce a four-page PDF summary of the output from PROCESS, including the major parameters, poloidal and toroidal cross-sections, temperature and density profiles and TF coil layout and turn strucuture. ### Usage diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index e3c4472200..13ac886a03 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -6,10 +6,8 @@ Some more info can be found [here](https://wiki.fusion.ciemat.es/wiki/Bootstrap_ ## Selection 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 -existence of pedestals, whereas the Sauter et al. scaling -(`i_bootstrap_current = 4`) allows general profiles to be used. +can be either input into the code directly, or calculated using one of eleven +methods, as summarized below: -------------- @@ -460,7 +458,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 +467,230 @@ $$ 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 + +Is selected by setting `i_bootstrap_current = 7`[^9] + +Based off of 350 plasma profiles from Experimento Tokamak Esferico (ETE) spherical tokamak. +Profiles were taken to be Gaussian shaped functions. + +The range of parameters from the discharges in ETE can be found in the table below: + +| Parameter | Value | Parameter | Value | +|-----------------------|-----------|-----------------------|-----------| +| Aspect ratio, $A$ | 1.5 | Electron temperature profile index, $\alpha_{\text{T}_\text{e}}$ | 0.02 | +| Major radius, $R_0$ | 0.3 $\text{[m]}$ | Ion temperature profile index, $\alpha_{\text{T}_\text{i}}$ | 2 | +| Core plasma pressure, $p(0)$ | 15 $\text{[kPa]}$ | Elongation, $\kappa(a)$ | 2 | +| Core electron/ion temperature, $T_{\text{e,i}}(0)$ | 1 $\text{[keV]}$ | Triangularity, $\delta$ | 0.3 | +| Edge electron/ion temperature $T_{\text{e,i}}(a)$ | 0.1 $\text{[keV]}$ | Plasma Current, $I_{\text{p}}$ | 200 $\text{[kA]}$ | +| Pressure profile index, $\alpha_{\text{p}}$ | 3 | Toroidal field on-axis, $B_0$ | 0.4 $\text{[T]}$ | +| Plasma total beta, $\beta$ | 4-10% | Plasma effective charge, $Z_{\text{eff}}$ | 1 | + +Errors mostly up to the order of 10% are obtained when both expressions are compared with the equilibrium estimates for the bootstrap current in ETE. + +#### Scaling 1 + +!!! note "Applicability of 1st scaling" + + Andrade et.al[^9] actually presents two scalings, the first is: + + $$ + \frac{I_{\text{BS}}}{I_\text{p}}5C_{\text{bs}}c_{\text{p}}^{\lambda}\frac{\beta_{\text{N}}q_{\text{cyl}}}{\epsilon^{1/2}l_i^{\gamma}}\frac{R_0}{R_{\text{m}}} + $$ + + $R_{\text{m}}$ in this case is the major radius of the magnetic axis which `PROCESS` does not have as it does not currently calculate the [Shafranov shift](https://wiki.fusion.ciemat.es/wiki/Shafranov_shift). If this is implemented then the first Andrade scaling can be properly implemented. + +#### Scaling 2 | `bootstrap_fraction_andrade()` + +This form of the Andrade scaling in terms of $\beta_{\text{p}}$ found using a least-square fit to 347 of the equilibria points. + +$$ +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] + +This scaling is based off of 170 discharges from TFTR with the neoclassical bootstrap current being calculated by the TRANSP plasma analysis code. +The plasma parameters of the discharges can be seen in the table below: + +| Parameter | Value | +|-----------------------|-----------| +| Plasma Current, $I_{\text{p}}$ | 0.6 - 2.7 $\text{[MA]}$ | +| Edge safety factor, $q(a)$ | 2.8 - 11.0 | +| Injected NBI power, $P_{\text{NBI}}$ | 2.0 - 35.0 $\text{[MW]}$ | +| Injected ICRH power, $P_{\text{ICRH}}$ | 1.5 - 6.0 $\text{[MW]}$ | +| Toroidal field on-axis, $B_{\text{T}}$ | 1.9 - 5.7 $\text{[T]}$ | +| Core electron density, $n_{\text{e,0}}$ | 0.2 - 1.2 $[10^{20} \text{m}^{-3}]$ | + +A wide variety of discharge regimes are included, such as: L-mode supershots, discharges +with reversed shear (RS) and enhanced reversed shear (ERS), and discharges with increased-$l_i$. +Discharges with both monotonic $q$ profiles and with reversed shear are included in the dataset. + +For an example ERS discharge the bootstrap current is driven by when the thermal particles surpasses 1 MA ($f_{\text{boot}}$ = 63%). Some ERS discharges in TFTR achieved $f_{\text{boot}}$ greater than 100% transiently. + + +!!! note "Change of profile index definition" + + Hoang et.al uses a different definition for the profile indexes such that + $\alpha_{\text{p}}$ is defined as the ratio of the central and the volume-averaged values, and the peakedness of the density of the total plasma current (defined as ratio of the central value and $I_{\text{p}}$), $\alpha_{\text{J}}$ + + Assuming that the pressure and current profile is parabolic we can represent these ratios as $\frac{p_0}{\langle p \rangle}= \alpha_{\text{p}}+1$ + + **This could lead to large changes in the value depending on interpretation of the profile index** + +$$ +C_{\text{BS}} = \sqrt{\frac{\alpha_{\text{p}}+1}{\alpha_{\text{j}}+1}} +$$ + +$$ +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] + +This scaling data is based off of equilibria from Miller et.al.[^12]. +The equilibria from Miller et.al. are in the range of $A$ = 1.2 - 3 that are stable to infinite $n$ ballooning and low $n$ kink modes at a bootstrap fraction of 99% for $\kappa$ = 2, 2.5, 3.0. The results were parameterized as a function of aspect ratio and elongation. + +The parametric dependency of $\beta_{\text{p}}$ and $\beta_{\text{T}}$ are based on the fitting of the DIII-D high equivalent DT yield results. + +$$ +\beta_{\text{N}}=\frac{\left(3.09+\frac{3.35}{A}+\frac{3.87}{A^{0.5}}\right)\left(\frac{\kappa}{3}\right)^{0.5}}{f_{\text{peak}}^{0.5}} +$$ + +$$ +\beta_{\text{T}}=\frac{25}{\beta_p}\left(\frac{1+\kappa^2}{2}\right)\left(\frac{\beta_N}{100}\right)^2 +$$ + +Here $\beta_{\text{p}}$ is given by + +$$ +\beta_p=f_{b s} \frac{\sqrt{A}}{C_{b s} f_{\text{peak}}^{0.25}} +$$ + +$$ +C_{\text{BS}} = 0.773+0.019\kappa +$$ + +Parabolic profiles should be used for best results as the pressure peaking value is calculated as the product of a parabolic temperature and density profile. + +$$ +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's + +This scaling is found by solving the Hirshman-Sigmar bootstrap current model using the matrix inversion method to create bootstrap current scalings with variables given explicitly in the TPC systems code[^13]. +A 8800 point database for the bootstrap current fraction using the bootstrap current density calculation module in the ACCOME code is used, using the variable ranges in the table below: + +The fitting of the variable exponents is done using the least squares method with a $R^2$ value of > 0.98 for [scaling one](#scaling-1--bootstrap_fraction_gi_i) and > 0.96 for [scaling two](#scaling-2--bootstrap_fraction_gi_ii) compared to the ACCOME data. + + +| Parameter | Range | Points | +|----------------------------|----------------|--------| +| Major radius, $R$ | 5.0 $[\mathrm{m}]$ | 1 | +| Aspect ratio, $A$ | 1.3, 1.5, 1.7, 2.0, 2.2, 2.5, 3.0, 3.5, 4.0, 5.0 | 10 | +| Elongation, $\kappa$ | $\sim$ 2 | 1 | +| Triangularity, $\delta$ | $\sim$ 0.3 | 1 | +| Density profile index, $a_{\text{n}}$ | 0.1-0.8 | 8 | +| Temperature profile index, $a_{\text{T}}$ | 1.0-3.0 | 11 | +| Effective charge, $Z_{\text{eff}}$ | 1.2-3.0 | 10 | + +The plasma parameters for each point in the aspect ratio scan can be seen in the table below: + +| Aspect ratio A | 1.3 | 1.5 | 1.7 | 2.0 | 2.2 | 2.5 | 3.0 | 3.5 | 4.0 | 5.0 | +|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:| +| Electron density at axis, $n_{\text{e0}}\left[10^{20} \mathrm{~m}^{-3}\right]$ | 1.0 | 1.0 | 1.5 | 1.5 | 1.5 | 2.0 | 2.0 | 2.0 | 2.0 | 1.0 | +| Electron temperature at axis, $T_{\text{e0}}[\mathrm{keV}]$ | 40 | 20 | 30 | 30 | 30 | 40 | 20 | 40 | 20 | 30 | +| Plasma current, $I_p$ $[\mathrm{MA}]$ | 20 | 15 | 20 | 15 | 15 | 15 | 10 | 15 | 10 | 5 | +| Toroidal magnetic field at axis, $B_{\text{T}}$ $[\mathrm{T}]$ | 3.0 | 2.0 | 3.0 | 2.0 | 4.0 | 2.0 | 2.0 | 6.0 | 2.0 | 5.0 | +| Poloidal beta, $\beta_{\text{p}}$ | 0.9-2.6 | 0.6-1.8 | 0.6-1.8 | 0.8-1.9 | 0.7-2.0 | 0.9-2.7 | 0.7-2.2 | 0.5-1.4 | 0.4-1.2 | 0.8-2.3 | + + + +#### Scaling 1 | `bootstrap_fraction_gi_I()` + +Is selected by setting `i_bootstrap_current = 10` + +Scaling 1 has better accuracy than Scaling 2. However, Scaling 1 overestimated the $f_{\text{BS}}$ value for reversed shear equilibrium. Although Scaling 2 does not have an internal current profile term, it can predict the $f_{\text{BS}}$ values to a certain extent for the high-$f_{\text{BS}}$ equilibria which are expected for next fusion devices. + +$$ +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} +$$ + +#### Scaling 2 | `bootstrap_fraction_gi_II()` + +Is selected by setting `i_bootstrap_current = 11` + +This scaling has the $q$ profile dependance removed to obtain a scaling formula with much more flexible variables than that by a single profile factor for internal current profile. + +$$ +C_{\text{BS}} = 0.382 \epsilon^{-0.242} \alpha_{\text{p}}^{0.974} \alpha_{\text{T}}^{-0.416} Z_{\text{eff}}^{0.178} +$$ + +$$ +f_{\text{BS}} = C_{\text{BS}} \beta_{\text{p}}\sqrt{\epsilon} +$$ + --------------------- ## Setting of maximum desirable bootstrap current fraction @@ -500,5 +722,9 @@ 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]: Miller, R L, "Stable bootstrap-current driven equilibria for low aspect ratio tokamaks". Switzerland: N. p., 1996. Web.https://fusion.gat.com/pubs-ext/MISCONF96/A22433.pdf +[^13]: 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. \ No newline at end of file diff --git a/documentation/proc-pages/physics-models/profiles/plasma_profiles.md b/documentation/proc-pages/physics-models/profiles/plasma_profiles.md index 7c26247223..e211f3d223 100644 --- a/documentation/proc-pages/physics-models/profiles/plasma_profiles.md +++ b/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. diff --git a/documentation/proc-pages/usage/plotting.md b/documentation/proc-pages/usage/plotting.md index 79074c509f..703ff7ae7e 100644 --- a/documentation/proc-pages/usage/plotting.md +++ b/documentation/proc-pages/usage/plotting.md @@ -29,6 +29,11 @@ process -i path/to/IN.DAT --plot --mfile path/to/MFILE.DAT
Figure 3: plot_proc TF coil and turn structure page
+
+![plot_proc_4](../images/plot_proc_4.PNG){ width="100%"} +
Figure 4: plot_proc bootstrap current model comparison page
+
+

Scan files

`plot_scans` is a tool to show the change in variables as a scan variable is varied. diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 822202d38c..77ffaca6aa 100755 --- a/process/io/plot_proc.py +++ b/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): @@ -2404,7 +2404,7 @@ def plot_physics_info(axis, mfile_data, scan): data = [ ("fusion_power", "Fusion power", "MW"), ("bigq", "$Q_{p}$", ""), - ("plasma_current_MA", "$I_p$", "MA"), + ("plasma_current_ma", "$I_p$", "MA"), ("bt", "Vacuum $B_T$ at $R_0$", "T"), ("q95", r"$q_{\mathrm{95}}$", ""), ("normalised_thermal_beta", r"$\beta_N$, thermal", "% m T MA$^{-1}$"), @@ -2898,10 +2898,92 @@ 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_I = mfile_data.data["bscf_gi_i"].get_scan(scan) + boot_gi_II = mfile_data.data["bscf_gi_ii"].get_scan(scan) + + # Data for the box plot + data = { + "IPDG": boot_ipdg, + "Sauter": boot_sauter, + "Nevins": boot_nenins, + "Wilson": boot_wilson, + "Sakai": boot_sakai, + "ARIES": boot_aries, + "Andrade": boot_andrade, + "Hoang": boot_hoang, + "Wong": boot_wong, + "Gi-I": boot_gi_I, + "Gi-II": boot_gi_II, + } + # Create the violin plot + axis.violinplot(data.values(), showextrema=False) + + # Create the box plot + axis.boxplot( + data.values(), 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.values()))) + for index, (key, value) in enumerate(data.items()): + axis.scatter(1, value, color=colors[index], label=key, alpha=1.0) + axis.legend(loc="upper left", bbox_to_anchor=(1, 1)) + + # Calculate average, standard deviation, and median + data_values = list(data.values()) + avg_bootstrap = np.mean(data_values) + std_bootstrap = np.std(data_values) + median_bootstrap = np.median(data_values) + + # Plot average, standard deviation, and median as text + axis.text( + 0.65, 0.9, f"Average: {avg_bootstrap:.4f}", transform=axis.transAxes, fontsize=9 + ) + axis.text( + 0.65, + 0.85, + f"Standard Dev: {std_bootstrap:.4f}", + transform=axis.transAxes, + fontsize=9, + ) + axis.text( + 0.65, + 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 +3070,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 +3080,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 +3337,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 +3356,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 +3365,7 @@ def main(args=None): plt.close(page1) plt.close(page2) plt.close(page3) + plt.close(page4) if __name__ == "__main__": diff --git a/process/physics.py b/process/physics.py index 9afa7f19fe..7bedecbf45 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1,6 +1,7 @@ import math from typing import Tuple import numpy as np +import scipy import numba as nb from scipy.optimize import root_scalar from process.utilities.f2py_string_patch import f2py_compatible_to_string @@ -1726,6 +1727,69 @@ def physics(self): ) ) + current_drive_variables.bscf_aries = ( + current_drive_variables.cboot + * self.bootstrap_fraction_aries( + betap=physics_variables.betap, + rli=physics_variables.rli, + core_density=physics_variables.ne0, + average_density=physics_variables.dene, + inverse_aspect=physics_variables.eps, + ) + ) + + current_drive_variables.bscf_andrade = ( + current_drive_variables.cboot + * self.bootstrap_fraction_andrade( + betap=physics_variables.betap, + core_pressure=physics_variables.p0, + average_pressure=physics_variables.vol_avg_pressure, + inverse_aspect=physics_variables.eps, + ) + ) + current_drive_variables.bscf_hoang = ( + current_drive_variables.cboot + * self.bootstrap_fraction_hoang( + betap=physics_variables.betap, + pressure_index=physics_variables.alphap, + current_index=physics_variables.alphaj, + inverse_aspect=physics_variables.eps, + ) + ) + current_drive_variables.bscf_wong = ( + current_drive_variables.cboot + * self.bootstrap_fraction_wong( + betap=physics_variables.betap, + density_index=physics_variables.alphan, + temperature_index=physics_variables.alphat, + inverse_aspect=physics_variables.eps, + elongation=physics_variables.kappa, + ) + ) + current_drive_variables.bscf_gi_I = ( + current_drive_variables.cboot + * self.bootstrap_fraction_gi_I( + betap=physics_variables.betap, + pressure_index=physics_variables.alphap, + temperature_index=physics_variables.alphat, + inverse_aspect=physics_variables.eps, + effective_charge=physics_variables.zeff, + q95=physics_variables.q95, + q0=physics_variables.q0, + ) + ) + + current_drive_variables.bscf_gi_II = ( + current_drive_variables.cboot + * self.bootstrap_fraction_gi_II( + betap=physics_variables.betap, + pressure_index=physics_variables.alphap, + temperature_index=physics_variables.alphat, + inverse_aspect=physics_variables.eps, + effective_charge=physics_variables.zeff, + ) + ) + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: current_drive_variables.bootstrap_current_fraction = abs( current_drive_variables.bootstrap_current_fraction_max @@ -1756,6 +1820,30 @@ def physics(self): current_drive_variables.bootstrap_current_fraction = ( current_drive_variables.bscf_sakai ) + elif physics_variables.i_bootstrap_current == 6: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_aries + ) + elif physics_variables.i_bootstrap_current == 7: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_andrade + ) + elif physics_variables.i_bootstrap_current == 8: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_hoang + ) + elif physics_variables.i_bootstrap_current == 9: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_wong + ) + elif physics_variables.i_bootstrap_current == 10: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_gi_I + ) + elif physics_variables.i_bootstrap_current == 11: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_gi_II + ) else: error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) @@ -3640,6 +3728,20 @@ def outplas(self): physics_variables.dnla, "OP ", ) + po.ovarre( + self.outfile, + "Plasma pressure on axis (Pa)", + "(p0)", + physics_variables.p0, + "OP ", + ) + po.ovarre( + self.outfile, + "Volume averaged plasma pressure (Pa)", + "(vol_avg_pressure)", + physics_variables.vol_avg_pressure, + "OP ", + ) if stellarator_variables.istell == 0: po.ovarre( @@ -3891,6 +3993,12 @@ def outplas(self): "(tbeta)", physics_variables.tbeta, ) + po.ovarrf( + self.outfile, + "Pressure profile index", + "(alphap)", + physics_variables.alphap, + ) if stellarator_variables.istell == 0: po.osubhd(self.outfile, "Density Limit using different models :") @@ -5024,7 +5132,6 @@ def outplas(self): current_drive_variables.bscf_wilson, "OP ", ) - po.ovarrf( self.outfile, "Bootstrap fraction (Sakai)", @@ -5032,6 +5139,49 @@ def outplas(self): current_drive_variables.bscf_sakai, "OP ", ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (ARIES)", + "(bscf_aries)", + current_drive_variables.bscf_aries, + "OP ", + ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Andrade)", + "(bscf_andrade)", + current_drive_variables.bscf_andrade, + "OP ", + ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Hoang)", + "(bscf_hoang)", + current_drive_variables.bscf_hoang, + "OP ", + ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Wong)", + "(bscf_wong)", + current_drive_variables.bscf_wong, + "OP ", + ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Gi I)", + "(bscf_gi_I)", + current_drive_variables.bscf_gi_I, + "OP ", + ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Gi II)", + "(bscf_gi_II)", + current_drive_variables.bscf_gi_II, + "OP ", + ) + po.ovarrf( self.outfile, "Diamagnetic fraction (Hender)", @@ -5088,6 +5238,36 @@ def outplas(self): self.outfile, " (Sakai et al bootstrap current fraction model used)", ) + elif physics_variables.i_bootstrap_current == 6: + po.ocmmnt( + self.outfile, + " (ARIES bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 7: + po.ocmmnt( + self.outfile, + " (Andrade et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 8: + po.ocmmnt( + self.outfile, + " (Hoang et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 9: + po.ocmmnt( + self.outfile, + " (Wong et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 10: + po.ocmmnt( + self.outfile, + " (Gi-I et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 11: + po.ocmmnt( + self.outfile, + " (Gi-II et al bootstrap current fraction model used)", + ) if physics_variables.i_diamagnetic_current == 0: po.ocmmnt( @@ -5723,6 +5903,7 @@ def bootstrap_fraction_sakai( alphan (float): Density profile index alphat (float): Temperature profile index eps (float): Inverse aspect ratio. + rli (float): Internal Inductance Returns: float: The calculated bootstrap fraction. @@ -5752,6 +5933,276 @@ def bootstrap_fraction_sakai( * alphat ** (0.502 * eps - 0.273) ) + @staticmethod + def bootstrap_fraction_aries( + betap: float, + rli: float, + core_density: float, + average_density: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the ARIES formula. + + Parameters: + betap (float): Plasma poloidal beta. + rli (float): Plasma normalized internal inductance. + core_density (float): Core plasma density. + average_density (float): Average plasma density. + inverse_aspect (float): Inverse aspect ratio. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + - The source reference does not provide any info about the derivation of the formula. It is only stated + + References: + - 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. + + """ + # Using the standard variable naming from the ARIES paper + a_1 = 1.10 - 1.165 * rli + 0.47 * rli**2 + b_1 = 0.806 - 0.885 * rli + 0.297 * rli**2 + + c_bs = a_1 + b_1 * (core_density / average_density) + + return c_bs * np.sqrt(inverse_aspect) * betap + + @staticmethod + def bootstrap_fraction_andrade( + betap: float, + core_pressure: float, + average_pressure: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Andrade et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + core_pressure (float): Core plasma pressure. + average_pressure (float): Average plasma pressure. + inverse_aspect (float): Inverse aspect ratio. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + - Based off 350 plasma profiles from Experimento Tokamak Esferico (ETE) spherical tokamak + - A = 1.5, R_0 = 0.3m, I_p = 200kA, B_0=0.4T, beta = 4-10%. Profiles taken as Gaussian shaped functions. + - Errors mostly up to the order of 10% are obtained when both expressions are compared with the equilibrium estimates for the + bootstrap current in ETE + + References: + - 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. + + """ + + # Using the standard variable naming from the Andrade et.al. paper + c_p = core_pressure / average_pressure + + # Error +- 0.0007 + c_bs = 0.2340 + + return c_bs * np.sqrt(inverse_aspect) * betap * c_p**0.8 + + @staticmethod + def bootstrap_fraction_hoang( + betap: float, + pressure_index: float, + current_index: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Hoang et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + pressure_index (float): Pressure profile index. + current_index (float): Current profile index. + inverse_aspect (float): Inverse aspect ratio. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + - Based off of TFTR data calculated using the TRANSP plasma analysis code + - 170 discharges which was assembled to study the tritium influx and transport in discharges with D-only neutral beam + injection (NBI) + - Contains L-mode, supershots, reversed shear, enhanced reversed shear and increased li discharges + - Discharges with monotonic flux profiles with reversed shear are also included + - Is applied to circular cross-section plasmas + + References: + - 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. + """ + + # Using the standard variable naming from the Hoang et.al. paper + # Hoang et.al uses a different definition for the profile indexes such that + # alpha_p is defined as the ratio of the central and the volume-averaged values, and the peakedness of the density of the total plasma current + # (defined as ratio of the central value and I_p), alpha_j$ + + # We assume the pressure and current profile is parabolic and use the (profile_index +1) term in lieu + # The term represents the ratio of the the core to volume averaged value + + # This could lead to large changes in the value depending on interpretation of the profile index + + c_bs = np.sqrt((pressure_index + 1) / (current_index + 1)) + + return 0.4 * np.sqrt(inverse_aspect) * betap**0.9 * c_bs + + @staticmethod + def bootstrap_fraction_wong( + betap: float, + density_index: float, + temperature_index: float, + inverse_aspect: float, + elongation: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Wong et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + density_index (float): Density profile index. + temperature_index (float): Temperature profile index. + inverse_aspect (float): Inverse aspect ratio. + elongation (float): Plasma elongation. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + - Data is based off of equilibria from Miller et al. + - A: 1.2 - 3.0 and stable to n ballooning and low n kink modes at a bootstrap fraction of 99% for kappa = 2, 2.5 and 3 + - The results were parameterized as a function of aspect ratio and elongation + - The parametric dependency of beta_p and beta_T are based on fitting of the DIII-D high equivalent DT yield results + - Parabolic profiles should be used for best results as the pressure peaking value is calculated as the product of a parabolic + temperature and density profile + + References: + - 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. + + - Miller, R L, "Stable bootstrap-current driven equilibria for low aspect ratio tokamaks". + Switzerland: N. p., 1996. Web.https://fusion.gat.com/pubs-ext/MISCONF96/A22433.pdf + """ + # Using the standard variable naming from the Wong et.al. paper + f_peak = 2.0 / scipy.special.beta(0.5, density_index + temperature_index + 1) + + c_bs = 0.773 + 0.019 * elongation + + return c_bs * f_peak**0.25 * betap * np.sqrt(inverse_aspect) + + @staticmethod + def bootstrap_fraction_gi_I( + betap: float, + pressure_index: float, + temperature_index: float, + inverse_aspect: float, + effective_charge: float, + q95: float, + q0: float, + ) -> float: + """ + Calculate the bootstrap fraction using the first scaling from the Gi et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + pressure_index (float): Pressure profile index. + temperature_index (float): Temperature profile index. + inverse_aspect (float): Inverse aspect ratio. + effective_charge (float): Plasma effective charge. + q95 (float): Safety factor at 95% of the plasma radius. + q0 (float): Safety factor at the magnetic axis. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + - Scaling found by solving the Hirshman-Sigmar bootstrap modelusing the matrix inversion method + - Method was done to put the scaling into parameters compatible with the TPC systems code + - Uses the ACCOME code to create bootstrap current fractions without using the itrative calculations of the + curent drive and equilibrium models in the scan + - R = 5.0 m, A = 1.3 - 5.0, kappa = 2, traing = 0.3, alpha_n = 0.1 - 0.8, alpha_t = 1.0 - 3.0, Z_eff = 1.2 - 3.0 + - Uses parabolic plasma profiles only. + - Scaling 1 has better accuracy than Scaling 2. However, Scaling 1 overestimated the f_BS value for reversed shear + equilibrium. + + References: + - 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. + """ + + # Using the standard variable naming from the Gi et.al. paper + + c_bs = ( + 0.474 + * inverse_aspect**-0.1 + * pressure_index**0.974 + * temperature_index**-0.416 + * effective_charge**0.178 + * (q95 / q0) ** -0.133 + ) + + return c_bs * np.sqrt(inverse_aspect) * betap + + @staticmethod + def bootstrap_fraction_gi_II( + betap: float, + pressure_index: float, + temperature_index: float, + inverse_aspect: float, + effective_charge: float, + ) -> float: + """ + Calculate the bootstrap fraction using the second scaling from the Gi et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + pressure_index (float): Pressure profile index. + temperature_index (float): Temperature profile index. + inverse_aspect (float): Inverse aspect ratio. + effective_charge (float): Plasma effective charge. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + - Scaling found by solving the Hirshman-Sigmar bootstrap modelusing the matrix inversion method + - Method was done to put the scaling into parameters compatible with the TPC systems code + - Uses the ACCOME code to create bootstrap current fractions without using the itrative calculations of the + curent drive and equilibrium models in the scan + - R = 5.0 m, A = 1.3 - 5.0, kappa = 2, traing = 0.3, alpha_n = 0.1 - 0.8, alpha_t = 1.0 - 3.0, Z_eff = 1.2 - 3.0 + - Uses parabolic plasma profiles only. + - This scaling has the q profile dependance removed to obtain a scaling formula with much more flexible variables than + that by a single profile factor for internal current profile. + + References: + - 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. + """ + + # Using the standard variable naming from the Gi et.al. paper + + c_bs = ( + 0.382 + * inverse_aspect**-0.242 + * pressure_index**0.974 + * temperature_index**-0.416 + * effective_charge**0.178 + ) + + return c_bs * np.sqrt(inverse_aspect) * betap + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index a24ec353fa..0322576da0 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -264,6 +264,12 @@ def calculate_profile_factors(self) -> None: physics_variables.alphap = physics_variables.alphan + physics_variables.alphat + # Shall assume that the pressure profile is parabolic. Can find volume average from + # profile index and core value the same as for density and temperature + physics_variables.vol_avg_pressure = physics_variables.p0 / ( + physics_variables.alphap + 1 + ) + @staticmethod def calculate_parabolic_profile_factors() -> None: """ diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index f313216951..5818527958 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -45,8 +45,26 @@ module current_drive_variables real(dp) :: bscf_sakai !! Bootstrap current fraction, Sakai et al model + real(dp) :: bscf_aries + !! Bootstrap current fraction, ARIES model + + real(dp) :: bscf_andrade + !! Bootstrap current fraction, Andrade et al model + + real(dp) :: bscf_hoang + !! Bootstrap current fraction, Hoang et al model + + real(dp) :: bscf_wong + !! Bootstrap current fraction, Wong et al model + + real(dp) :: bscf_gi_I + !! Bootstrap current fraction, first Gi et al model + + real(dp) :: bscf_gi_II + !! Bootstrap current fraction, second Gi et al model + real(dp) :: cboot - !! bootstrap current fraction multiplier (`i_bootstrap_current=1`) + !! bootstrap current fraction multiplier real(dp) :: beam_current !! neutral beam current (A) @@ -243,6 +261,13 @@ subroutine init_current_drive_variables bscf_nevins = 0.0D0 bscf_sauter = 0.0D0 bscf_wilson = 0.0D0 + bscf_sakai = 0.0D0 + bscf_aries = 0.0D0 + bscf_andrade = 0.0D0 + bscf_hoang = 0.0D0 + bscf_wong = 0.0D0 + bscf_gi_I = 0.0D0 + bscf_gi_II = 0.0D0 cboot = 1.0D0 beam_current = 0.0D0 diacf_hender = 0.0D0 diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 35cdb77bf7..5c4051825b 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -621,7 +621,7 @@ subroutine parse_input_file(in_file,out_file,show_changes) call parse_real_variable('taumax', taumax, 0.1D0, 100.0D0, & 'Maximum allowed energy confinement time (s)') case ('i_bootstrap_current') - call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 5, & + call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 11, & 'Switch for bootstrap scaling') case ('iculbl') call parse_int_variable('iculbl', iculbl, 0, 3, & diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index 3d2c298872..69dbdf09f8 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -253,6 +253,12 @@ module physics_variables !! - =3 for Wilson et al numerical scaling !! - =4 for Sauter et al scaling !! - =5 for Sakai et al scaling + !! - =6 for ARIES scaling + !! - =7 for Andrade et al scaling + !! - =8 for Hoang et al scaling + !! - =9 for Wong et al scaling + !! - =10 for Gi-I et al scaling + !! - =11 for Gi-II et al scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) @@ -553,6 +559,9 @@ module physics_variables real(dp) :: p0 !! central total plasma pressure (Pa) + real(dp) :: vol_avg_pressure + !! Volume averaged plasma pressure (Pa) + real(dp) :: f_dd_branching_trit !! branching ratio for DD -> T diff --git a/tests/integration/ref_dicts.json b/tests/integration/ref_dicts.json index 948132a67f..f01dc0c31a 100644 --- a/tests/integration/ref_dicts.json +++ b/tests/integration/ref_dicts.json @@ -1163,6 +1163,13 @@ "bscf_nevins": 0.0, "bscf_sauter": 0.0, "bscf_wilson": 0.0, + "bscf_sakai": 0.0, + "bscf_aries": 0.0, + "bscf_andrade": 0.0, + "bscf_hoang": 0.0, + "bscf_wong": 0.0, + "bscf_gi_I": 0.0, + "bscf_gi_II": 0.0, "bootstrap_current_fraction_max": 0.9, "bt": 5.68, "btot": 0.0, @@ -8947,6 +8954,13 @@ "bscf_nevins": "bootstrap current fraction, Nevins et al model", "bscf_sauter": "bootstrap current fraction, Sauter et al model", "bscf_wilson": "bootstrap current fraction, Wilson et al model", + "bscf_sakai": "bootstrap current fraction, Sakai model", + "bscf_aries": "bootstrap current fraction, ARIES model", + "bscf_andrade": "bootstrap current fraction, Andrade model", + "bscf_hoang": "bootstrap current fraction, Hoang model", + "bscf_wong": "bootstrap current fraction, Wong model", + "bscf_gi_I": "bootstrap current fraction, GI model I", + "bscf_gi_II": "bootstrap current fraction, GI model II", "bootstrap_current_fraction_max": "maximum fraction of plasma current from bootstrap; if `bootstrap_current_fraction_max < 0`,\n bootstrap fraction = abs(bootstrap_current_fraction_max)", "bt": "toroidal field on axis (T) (`iteration variable 2`)", "btot": "total toroidal + poloidal field (T)", @@ -17772,6 +17786,13 @@ "bscf_nevins", "bscf_sauter", "bscf_wilson", + "bscf_sakai", + "bscf_aries", + "bscf_andrade", + "bscf_hoang", + "bscf_wong", + "bscf_gi_I", + "bscf_gi_II", "cboot", "beam_current", "diacf_hender", diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 81e6935475..e94a44fdd5 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -566,6 +566,294 @@ def test_bootstrap_fraction_sakai(bootstrapfractionsakaiparam, monkeypatch, phys assert bfs == pytest.approx(bootstrapfractionsakaiparam.expected_bfs) +class BootstrapFractionAriesParam(NamedTuple): + betap: Any = None + + rli: Any = None + + core_density: Any = None + + average_density: Any = None + + inverse_aspect: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractionariesparam", + ( + BootstrapFractionAriesParam( + betap=1.2708883332338736, + rli=1.4279108047138775, + core_density=1.0695994460047332e20, + average_density=8.1317358967210131e19, + inverse_aspect=1 / 3, + expected_bfs=4.3237405809568441e-01, + ), + ), +) +def test_bootstrap_fraction_aries(bootstrapfractionariesparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_aries. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_aries( + betap=bootstrapfractionariesparam.betap, + rli=bootstrapfractionariesparam.rli, + core_density=bootstrapfractionariesparam.core_density, + average_density=bootstrapfractionariesparam.average_density, + inverse_aspect=bootstrapfractionariesparam.inverse_aspect, + ) + + assert bfs == pytest.approx(bootstrapfractionariesparam.expected_bfs) + + +class BootstrapFractionAndradeParam(NamedTuple): + betap: Any = None + + core_pressure: Any = None + + average_pressure: Any = None + + inverse_aspect: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractionandradeparam", + ( + BootstrapFractionAndradeParam( + betap=1.2708883332338736, + core_pressure=8.3049163275475602e05, + average_pressure=2.4072221239268288e05, + inverse_aspect=1 / 3, + expected_bfs=4.6240007834873120e-01, + ), + ), +) +def test_bootstrap_fraction_andrade(bootstrapfractionandradeparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_andrade. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_andrade( + betap=bootstrapfractionandradeparam.betap, + core_pressure=bootstrapfractionandradeparam.core_pressure, + average_pressure=bootstrapfractionandradeparam.average_pressure, + inverse_aspect=bootstrapfractionandradeparam.inverse_aspect, + ) + + assert bfs == pytest.approx(bootstrapfractionandradeparam.expected_bfs) + + +class BootstrapFractionHoangParam(NamedTuple): + betap: Any = None + + pressure_index: Any = None + + current_index: Any = None + + inverse_aspect: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractionhoangparam", + ( + BootstrapFractionHoangParam( + betap=1.2708883332338736, + pressure_index=2.4500000000000002e00, + current_index=2.8314361644755763e00, + inverse_aspect=1 / 3, + expected_bfs=0.27190974213794156, + ), + ), +) +def test_bootstrap_fraction_hoang(bootstrapfractionhoangparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_hoang. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_hoang( + betap=bootstrapfractionhoangparam.betap, + pressure_index=bootstrapfractionhoangparam.pressure_index, + current_index=bootstrapfractionhoangparam.current_index, + inverse_aspect=bootstrapfractionhoangparam.inverse_aspect, + ) + + assert bfs == pytest.approx(bootstrapfractionhoangparam.expected_bfs) + + +class BootstrapFractionWongParam(NamedTuple): + betap: Any = None + + density_index: Any = None + + temperature_index: Any = None + + inverse_aspect: Any = None + + elongation: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractionwongparam", + ( + BootstrapFractionWongParam( + betap=1.2708883332338736, + density_index=1.0000000000000000e00, + temperature_index=1.4500000000000000e00, + inverse_aspect=1 / 3, + elongation=1.8500000000000001e00, + expected_bfs=7.0706527916080808e-01, + ), + ), +) +def test_bootstrap_fraction_wong(bootstrapfractionwongparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_wong. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_wong( + betap=bootstrapfractionwongparam.betap, + density_index=bootstrapfractionwongparam.density_index, + temperature_index=bootstrapfractionwongparam.temperature_index, + inverse_aspect=bootstrapfractionwongparam.inverse_aspect, + elongation=bootstrapfractionwongparam.elongation, + ) + + assert bfs == pytest.approx(bootstrapfractionwongparam.expected_bfs) + + +class BootstrapFractionGiIParam(NamedTuple): + betap: Any = None + + pressure_index: Any = None + + temperature_index: Any = None + + inverse_aspect: Any = None + + effective_charge: Any = None + + q95: Any = None + + q0: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractiongiiparam", + ( + BootstrapFractionGiIParam( + betap=1.2708883332338736, + pressure_index=2.4500000000000002e00, + temperature_index=1.4500000000000000e00, + inverse_aspect=1 / 3, + effective_charge=2.5368733516769737e00, + q95=3.4656394133756647e00, + q0=1.0, + expected_bfs=7.9639753138719782e-01, + ), + ), +) +def test_bootstrap_fraction_gi_I(bootstrapfractiongiiparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_gi. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_gi_I( + betap=bootstrapfractiongiiparam.betap, + pressure_index=bootstrapfractiongiiparam.pressure_index, + temperature_index=bootstrapfractiongiiparam.temperature_index, + inverse_aspect=bootstrapfractiongiiparam.inverse_aspect, + effective_charge=bootstrapfractiongiiparam.effective_charge, + q95=bootstrapfractiongiiparam.q95, + q0=bootstrapfractiongiiparam.q0, + ) + + assert bfs == pytest.approx(bootstrapfractiongiiparam.expected_bfs) + + +class BootstrapFractionGiIIParam(NamedTuple): + betap: Any = None + + pressure_index: Any = None + + temperature_index: Any = None + + inverse_aspect: Any = None + + effective_charge: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractiongiiiparam", + ( + BootstrapFractionGiIIParam( + betap=1.2708883332338736, + pressure_index=2.4500000000000002e00, + temperature_index=1.4500000000000000e00, + inverse_aspect=1 / 3, + effective_charge=2.5368733516769737e00, + expected_bfs=8.8502865710180589e-01, + ), + ), +) +def test_bootstrap_fraction_gi_II(bootstrapfractiongiiiparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_gi. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_gi_II( + betap=bootstrapfractiongiiiparam.betap, + pressure_index=bootstrapfractiongiiiparam.pressure_index, + temperature_index=bootstrapfractiongiiiparam.temperature_index, + inverse_aspect=bootstrapfractiongiiiparam.inverse_aspect, + effective_charge=bootstrapfractiongiiiparam.effective_charge, + ) + + assert bfs == pytest.approx(bootstrapfractiongiiiparam.expected_bfs) + + class PlasmaCurrentParam(NamedTuple): normalised_total_beta: Any = None