diff --git a/documentation/source/physics-models/plasma_current/plasma_current.md b/documentation/source/physics-models/plasma_current/plasma_current.md index fdd994fd9..9a635a43e 100644 --- a/documentation/source/physics-models/plasma_current/plasma_current.md +++ b/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,12 +489,12 @@ $$ ----------------------------- -### 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}}$. @@ -503,7 +502,7 @@ The effect of toroidicity and of shaping of the plasma cross section is to modif $$ -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. @@ -511,7 +510,7 @@ It is not possible to derive a general analytic expression for $F$, so it has be 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}} diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 17168438b..88444f218 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -62,6 +62,7 @@ ) from process.models.physics.current_drive import ElectronBernstein, ElectronCyclotron from process.models.physics.impurity_radiation import read_impurity_file +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 @@ -3023,17 +3024,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, @@ -3049,8 +3050,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", @@ -13406,23 +13407,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, @@ -13432,7 +13434,7 @@ 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, @@ -13440,19 +13442,19 @@ def main_plot( ) # 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) @@ -13460,52 +13462,52 @@ def main_plot( # 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_pipe_bends(figs[29], m_file, scan) + ax_blanket = figs[29].add_subplot(122, aspect="equal") plot_blanket(ax_blanket, m_file, scan, radial_build, colour_scheme) plot_firstwall(ax_blanket, m_file, scan, radial_build, colour_scheme) ax_blanket.set_xlabel("Radial position [m]") @@ -13548,13 +13550,13 @@ def main_plot( ) 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): @@ -13631,7 +13633,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( diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 72e4fd564..1cee819fe 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -961,6 +961,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]""" @@ -1611,6 +1636,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, \ @@ -1893,6 +1927,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 diff --git a/process/main.py b/process/main.py index a2b243d10..9a8f6fc2f 100644 --- a/process/main.py +++ b/process/main.py @@ -102,6 +102,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 @@ -706,6 +707,7 @@ def __init__(self): self.plasma_bootstrap_current = PlasmaBootstrapCurrent( plasma_profile=self.plasma_profile ) + self.plasma_current = PlasmaCurrent() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, @@ -714,6 +716,7 @@ def __init__(self): plasma_density_limit=self.plasma_density_limit, plasma_exhaust=self.plasma_exhaust, plasma_bootstrap_current=self.plasma_bootstrap_current, + plasma_current=self.plasma_current, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 84c686704..accecbcb9 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -31,6 +31,52 @@ logger = logging.getLogger(__name__) +def calculate_cylindrical_safety_factor( + rmajor: float, + rminor: float, + plasma_current: float, + b_plasma_toroidal_on_axis: float, + kappa95: float, + triang95: float, +) -> float: + """Calculate the cylindrical safety factor from the IPDG89 guidelines. + + Parameters + ---------- + rmajor : float + Major radius of the tokamak in meters. + rminor : float + Minor radius of the tokamak in meters. + plasma_current : float + Plasma current in amperes. + b_plasma_toroidal_on_axis : float + Toroidal magnetic field on axis in tesla. + kappa95 : float + Elongation at 95% of the plasma boundary. + triang95 : float + Triangularity at 95% of the plasma boundary. + Returns + ------- + float + Cylindrical safety factor (dimensionless). + Notes + ----- + The cylindrical safety factor is calculated following the IPDG89 guidelines. + The formula accounts for plasma elongation and triangularity effects on the + safety factor through the kappa95 and triang95 parameters. + + """ + + # Calculate cyclindrical safety factor from IPDG89 + return ( + ((2 * np.pi) / constants.RMU0) + * rminor**2 + / (rmajor * plasma_current / b_plasma_toroidal_on_axis) + * 0.5 + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + ) + + @nb.jit(nopython=True, cache=True) def rether( alphan, @@ -84,490 +130,6 @@ def rether( return conie * (temp_plasma_ion_vol_avg_kev - te) / (te**1.5) -# ----------------------------------------------------- -# Plasma Current & Poloidal Field Calculations -# ----------------------------------------------------- - - -@nb.jit(nopython=True, cache=True) -def _plascar_bpol( - aspect: float, eps: float, kappa: float, delta: float -) -> tuple[float, float, float, float]: - """Calculate the poloidal field coefficients for determining the plasma current - and poloidal field. - - - This internal function calculates the poloidal field coefficients, - which is used to calculate the poloidal field and the plasma current. - - Parameters - ---------- - aspect : - plasma aspect ratio - eps : - inverse aspect ratio - kappa : - plasma elongation - delta : - plasma triangularity - - Returns - ------- - : - coefficients ff1, ff2, d1, d2 - - References - ---------- - - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). - 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), - 1729-1738. https://doi.org/10.13182/FST92-A29971 - - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - - """ - # Original coding, only suitable for TARTs [STAR Code] - - c1 = (kappa**2 / (1.0 + delta)) + delta - c2 = (kappa**2 / (1.0 - delta)) - delta - - d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 - d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 - - c1_aspect = ((c1 * eps) - 1.0) if aspect < c1 else (1.0 - (c1 * eps)) - - y1 = np.sqrt(c1_aspect / (1.0 + eps)) * ((1.0 + delta) / kappa) - y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * ((1.0 - delta) / kappa) - - h2 = (1.0 + (c2 - 1.0) * (eps / 2.0)) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) - f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * ((c2 * eps) + 1.0)) - g = (eps * kappa) / (1.0 - (eps * delta)) - ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) - - h1 = (1.0 + (1.0 - c1) * (eps / 2.0)) / np.sqrt((1.0 + eps) * c1_aspect) - f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - - if aspect < c1: - ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) - else: - ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) - - return ff1, ff2, d1, d2 - - -def calculate_poloidal_field( - i_plasma_current: int, - ip: float, - q95: float, - aspect: float, - eps: float, - b_plasma_toroidal_on_axis: float, - kappa: float, - delta: float, - perim: float, - rmu0: float, -) -> float: - """Function to calculate poloidal field from the plasma current - - This function calculates the poloidal field from the plasma current in Tesla, - using a simple calculation using Ampere's law for conventional - tokamaks, or for TARTs, a scaling from Peng, Galambos and - Shipe (1992). - - Parameters - ---------- - i_plasma_current : - current scaling model to use - ip : - plasma current (A) - q95 : - 95% flux surface safety factor - aspect : - plasma aspect ratio - eps : - inverse aspect ratio - b_plasma_toroidal_on_axis : - toroidal field on axis (T) - kappa : - plasma elongation - delta : - plasma triangularity - perim : - plasma perimeter (m) - rmu0 : - vacuum permeability (H/m) - - Returns - ------- - : - poloidal field in Tesla - - - References - ---------- - - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). - 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), - 1729-1738. https://doi.org/10.13182/FST92-A29971 - - """ - # Use Ampere's law using the plasma poloidal cross-section - if i_plasma_current != 2: - return rmu0 * ip / perim - # Use the relation from Peng, Galambos and Shipe (1992) [STAR code] otherwise - ff1, ff2, _, _ = _plascar_bpol(aspect, eps, kappa, delta) - - # Transform q95 to qbar - qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 - - return b_plasma_toroidal_on_axis * (ff1 + ff2) / (2.0 * np.pi * qbar) - - -def calculate_current_coefficient_peng( - eps: float, len_plasma_poloidal: float, rminor: float -) -> float: - """Calculate the plasma current scaling coefficient for the Peng scaling from the STAR code. - - Parameters - ---------- - eps : float - Plasma inverse aspect ratio. - len_plasma_poloidal : float - Plasma poloidal perimeter length [m]. - rminor : float - Plasma minor radius [m]. - - Returns - ------- - float - The plasma current scaling coefficient. - - """ - - return ( - (1.22 - 0.68 * eps) - / ((1.0 - eps * eps) ** 2) - * (len_plasma_poloidal / (2.0 * np.pi * rminor)) ** 2 - ) - - -def calculate_plasma_current_peng( - q95: float, - aspect: float, - eps: float, - rminor: float, - b_plasma_toroidal_on_axis: float, - kappa: float, - delta: float, -) -> float: - """Function to calculate plasma current (Peng scaling from the STAR code) - - This function calculates the plasma current in MA, - using a scaling from Peng, Galambos and Shipe (1992). - It is primarily used for Tight Aspect Ratio Tokamaks and is - selected via i_plasma_current=2. - - Parameters - ---------- - q95 : - 95% flux surface safety factor - aspect : - plasma aspect ratio - eps : - inverse aspect ratio - rminor : - plasma minor radius (m) - b_plasma_toroidal_on_axis : - toroidal field on axis (T) - kappa : - plasma elongation - delta : - plasma triangularity - - - Returns - ------- - : - plasma current in MA - - - References - ---------- - - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). - 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), - 1729-1738. https://doi.org/10.13182/FST92-A29971 - - """ - - # Transform q95 to qbar - qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 - - ff1, ff2, d1, d2 = _plascar_bpol(aspect, eps, kappa, delta) - - e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) - e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) - - return ( - rminor - * b_plasma_toroidal_on_axis - / qbar - * 5.0 - * kappa - / (2.0 * np.pi**2) - * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) - * (ff1 + ff2) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_ipdg89( - eps: float, kappa95: float, triang95: float -) -> float: - """Calculate the fq coefficient from the IPDG89 guidlines used in the plasma current scaling. - - This function calculates the fq coefficient used in the IPDG89 plasma current scaling, - based on the given plasma parameters. - - Parameters - ---------- - eps : - plasma inverse aspect ratio - kappa95 : - plasma elongation 95% - triang95 : - plasma triangularity 95% - - Returns - ------- - : - the fq plasma current coefficient - - - References - ---------- - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - return ( - 0.5 - * (1.17 - 0.65 * eps) - / ((1.0 - eps * eps) ** 2) - * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_todd( - eps: float, kappa95: float, triang95: float, model: int -) -> float: - """Calculate the fq coefficient used in the two Todd plasma current scalings. - - This function calculates the fq coefficient based on the given plasma parameters for the two Todd scalings. - - Parameters - ---------- - eps : - plasma inverse aspect ratio - kappa95 : - plasma elongation 95% - triang95 : - plasma triangularity 95% - model: - - Returns - ------- - : - the fq plasma current coefficient - - References - ---------- - - D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - # Calculate the Todd scaling based on the model - base_scaling = ( - (1.0 + 2.0 * eps**2) - * ((1.0 + kappa95**2) / 0.5) - * (1.24 - 0.54 * kappa95 + 0.3 * (kappa95**2 + triang95**2) + 0.125 * triang95) - ) - if model == 1: - return base_scaling - if model == 2: - return base_scaling * (1.0 + (abs(kappa95 - 1.2)) ** 3) - raise ProcessValueError(f"model = {model} is an invalid option") - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_hastie( - alphaj: float, - alphap: float, - b_plasma_toroidal_on_axis: float, - delta95: float, - eps: float, - kappa95: float, - pres_plasma_on_axis: float, - rmu0: float, -) -> float: - """Routine to calculate the f_q coefficient for the Connor-Hastie model used for scaling the plasma current. - - This routine calculates the f_q coefficient used for scaling the plasma current, - using the Connor-Hastie scaling - - - Parameters - ---------- - alphaj : - the current profile index - alphap : - the pressure profile index - b_plasma_toroidal_on_axis : - the toroidal field on axis (T) - delta95 : - the plasma triangularity 95% - eps : - the inverse aspect ratio - kappa95 : - the plasma elongation 95% - pres_plasma_on_axis : - the central plasma pressure (Pa) - rmu0 : - the vacuum permeability (H/m) - - Returns - ------- - : - the F coefficient - - - Reference: - - J.W.Connor and R.J.Hastie, Culham Lab Report CLM-M106 (1985). - https://scientific-publications.ukaea.uk/wp-content/uploads/CLM-M106-1.pdf - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - # Exponent in Connor-Hastie current profile - lamda = alphaj - - # Exponent in Connor-Hastie pressure profile - nu = alphap - - # Central plasma beta - beta0 = 2.0 * rmu0 * pres_plasma_on_axis / (b_plasma_toroidal_on_axis**2) - - # Plasma internal inductance - lamp1 = 1.0 + lamda - li = lamp1 / lamda * (lamp1 / lamda * np.log(lamp1) - 1.0) - - # T/r in AEA FUS 172 - kap1 = kappa95 + 1.0 - tr = kappa95 * delta95 / kap1**2 - - # E/r in AEA FUS 172 - er = (kappa95 - 1.0) / kap1 - - # T primed in AEA FUS 172 - tprime = 2.0 * tr * lamp1 / (1.0 + 0.5 * lamda) - - # E primed in AEA FUS 172 - eprime = er * lamp1 / (1.0 + lamda / 3.0) - - # Delta primed in AEA FUS 172 - deltap = (0.5 * kap1 * eps * 0.5 * li) + (beta0 / (0.5 * kap1 * eps)) * lamp1**2 / ( - 1.0 + nu - ) - - # Delta/R0 in AEA FUS 172 - deltar = beta0 / 6.0 * (1.0 + 5.0 * lamda / 6.0 + 0.25 * lamda**2) + ( - 0.5 * kap1 * eps - ) ** 2 * 0.125 * (1.0 - (lamda**2) / 3.0) - - # F coefficient - return (0.5 * kap1) ** 2 * ( - 1.0 - + eps**2 * (0.5 * kap1) ** 2 - + 0.5 * deltap**2 - + 2.0 * deltar - + 0.5 * (eprime**2 + er**2) - + 0.5 * (tprime**2 + 4.0 * tr**2) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_sauter( - eps: float, - kappa: float, - triang: float, -) -> float: - """Routine to calculate the f_q coefficient for the Sauter model used for scaling the plasma current. - - Parameters - ---------- - eps : - inverse aspect ratio - kappa : - plasma elongation at the separatrix - triang : - plasma triangularity at the separatrix - - Returns - ------- - : - the fq coefficient - - Reference: - - O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, - Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, - ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2016.04.033. - - """ - w07 = 1.0 # zero squareness - can be modified later if required - - return ( - (4.1e6 / 5.0e6) - * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) - * (1.0 + 0.09 * triang + 0.16 * triang**2) - * (1.0 + 0.45 * triang * eps) - / (1.0 - 0.74 * eps) - * (1.0 + 0.55 * (w07 - 1.0)) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_fiesta( - eps: float, kappa: float, triang: float -) -> float: - """Calculate the fq coefficient used in the FIESTA plasma current scaling. - - This function calculates the fq coefficient based on the given plasma parameters for the FIESTA scaling. - - Parameters - ---------- - eps : - plasma inverse aspect ratio - kappa : - plasma elongation at the separatrix - triang : - plasma triangularity at the separatrix - - Returns - ------- - : - the fq plasma current coefficient - - - References - ---------- - - S.Muldrew et.al,“PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, - Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. - - """ - return 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060 - - # ----------------------------------------------------- # Diamagnetic and Pfirsch-Schlüter Current Calculations # ----------------------------------------------------- @@ -641,6 +203,7 @@ def __init__( plasma_density_limit, plasma_exhaust, plasma_bootstrap_current: PlasmaBootstrapCurrent, + plasma_current, ): self.outfile = constants.NOUT self.mfile = constants.MFILE @@ -651,6 +214,7 @@ def __init__( self.density_limit = plasma_density_limit self.exhaust = plasma_exhaust self.plasma_bootstrap_current = plasma_bootstrap_current + self.current = plasma_current def physics(self): """Routine to calculate tokamak plasma physics information @@ -700,11 +264,7 @@ def physics(self): ) # Calculate plasma current - ( - physics_variables.b_plasma_poloidal_average, - physics_variables.qstar, - physics_variables.plasma_current, - ) = self.calculate_plasma_current( + physics_variables.plasma_current = self.current.calculate_plasma_current( physics_variables.alphaj, physics_variables.alphap, physics_variables.b_plasma_toroidal_on_axis, @@ -721,6 +281,30 @@ def physics(self): physics_variables.triang95, ) + physics_variables.qstar = calculate_cylindrical_safety_factor( + rmajor=physics_variables.rmajor, + rminor=physics_variables.rminor, + plasma_current=physics_variables.plasma_current, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + kappa95=physics_variables.kappa95, + triang95=physics_variables.triang95, + ) + + # Calculate the poloidal field generated by the plasma current + physics_variables.b_plasma_poloidal_average = ( + self.current.calculate_poloidal_field( + i_plasma_current=physics_variables.i_plasma_current, + ip=physics_variables.plasma_current, + q95=physics_variables.q95, + aspect=physics_variables.aspect, + eps=physics_variables.eps, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + kappa=physics_variables.kappa, + delta=physics_variables.triang, + perim=physics_variables.len_plasma_poloidal, + ) + ) + # ----------------------------------------------------- # Plasma Current Profile # ----------------------------------------------------- @@ -2026,211 +1610,6 @@ def plasma_ohmic_heating( return pden_plasma_ohmic_mw, p_plasma_ohmic_mw, f_res_plasma_neo, res_plasma - @staticmethod - def calculate_plasma_current( - alphaj: float, - alphap: float, - b_plasma_toroidal_on_axis: float, - eps: float, - i_plasma_current: int, - kappa: float, - kappa95: float, - pres_plasma_on_axis: float, - len_plasma_poloidal: float, - q95: float, - rmajor: float, - rminor: float, - triang: float, - triang95: float, - ) -> tuple[float, float, float, float, float]: - """Calculate the plasma current. - - This routine calculates the plasma current based on the edge safety factor q95. - It will also make the current profile parameters consistent with the q-profile if required. - - Parameters - ---------- - alphaj : float - Current profile index. - alphap : float - Pressure profile index. - b_plasma_toroidal_on_axis : float - Toroidal field on axis (T). - eps : float - Inverse aspect ratio. - i_plasma_current : int - Current scaling model to use. - 1 = Peng analytic fit - 2 = Peng divertor scaling (TART,STAR) - 3 = Simple ITER scaling - 4 = IPDG89 scaling - 5 = Todd empirical scaling I - 6 = Todd empirical scaling II - 7 = Connor-Hastie model - 8 = Sauter scaling (allowing negative triangularity) - 9 = FIESTA ST scaling - kappa : float - Plasma elongation. - kappa95 : float - Plasma elongation at 95% surface. - pres_plasma_on_axis : float - Central plasma pressure (Pa). - len_plasma_poloidal : float - Plasma perimeter length (m). - q95 : float - Plasma safety factor at 95% flux (= q-bar for i_plasma_current=2). - ind_plasma_internal_norm : float - Plasma normalised internal inductance. - rmajor : float - Major radius (m). - rminor : float - Minor radius (m). - triang : float - Plasma triangularity. - triang95 : float - Plasma triangularity at 95% surface. - - Returns - ------- - Tuple[float, float, float,] - Tuple containing b_plasma_poloidal_average, qstar, plasma_current, - - Raises - ------ - ValueError - If invalid value for i_plasma_current is provided. - ValueError - If invalid value for i_plasma_current is provided. - ValueError - If invalid value for i_plasma_current is provided. - - References - ---------- - - J D Galambos, STAR Code - Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document - - J D Galambos, STAR Code - Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document - - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). - 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), - 1729-1738. https://doi.org/10.13182/FST92-A29971 - - ITER Physics Design Guidelines - 1989 [IPDG89], N. A. Uckan et al, ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 - - M. Kovari et al, 2014, "PROCESS" - A systems code for fusion power plants - Part 1: Physics - - M. Kovari et al, 2014, "PROCESS" - A systems code for fusion power plants - Part 1: Physics - - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO - - M. Kovari et al, 2014, "PROCESS" - A systems code for fusion power plants - Part 1: Physics - - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO - - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant - - M. Kovari et al, 2014, "PROCESS" - A systems code for fusion power plants - Part 1: Physics - - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO - - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant - - Sauter, Geometric formulas for systems codes..., FED 2016 - - """ - # Aspect ratio - aspect_ratio = 1.0 / eps - - # Only the Sauter scaling (i_plasma_current=8) is suitable for negative triangularity: - if i_plasma_current != 8 and triang < 0.0: - raise ProcessValueError( - f"Triangularity is negative without i_plasma_current = 8 selected: {triang=}, {i_plasma_current=}" - ) - - # Calculate the function Fq that scales the edge q from the - # circular cross-section cylindrical case - - # Peng analytical fit - if i_plasma_current == 1: - fq = calculate_current_coefficient_peng(eps, len_plasma_poloidal, rminor) - - # Peng scaling for double null divertor; TARTs [STAR Code] - elif i_plasma_current == 2: - plasma_current = 1.0e6 * calculate_plasma_current_peng( - q95, aspect_ratio, eps, rminor, b_plasma_toroidal_on_axis, kappa, triang - ) - - # Simple ITER scaling (simply the cylindrical case) - elif i_plasma_current == 3: - fq = 1.0 - - # ITER formula (IPDG89) - elif i_plasma_current == 4: - fq = calculate_current_coefficient_ipdg89(eps, kappa95, triang95) - - # Todd empirical scalings - # D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 - elif i_plasma_current in [5, 6]: - fq = calculate_current_coefficient_todd(eps, kappa95, triang95, model=1) - - if i_plasma_current == 6: - fq = calculate_current_coefficient_todd(eps, kappa95, triang95, model=2) - - # Connor-Hastie asymptotically-correct expression - elif i_plasma_current == 7: - fq = calculate_current_coefficient_hastie( - alphaj, - alphap, - b_plasma_toroidal_on_axis, - triang95, - eps, - kappa95, - pres_plasma_on_axis, - constants.RMU0, - ) - - # Sauter scaling allowing negative triangularity [FED May 2016] - # https://doi.org/10.1016/j.fusengdes.2016.04.033. - elif i_plasma_current == 8: - # Assumes zero squareness, note takes kappa, delta at separatrix not _95 - fq = calculate_current_coefficient_sauter(eps, kappa, triang) - - # FIESTA ST scaling - # https://doi.org/10.1016/j.fusengdes.2020.111530. - elif i_plasma_current == 9: - fq = calculate_current_coefficient_fiesta(eps, kappa, triang) - else: - raise ProcessValueError(f"Invalid value {i_plasma_current=}") - - # Main plasma current calculation using the fq value from the different settings - if i_plasma_current != 2: - plasma_current = ( - (2.0 * np.pi / constants.RMU0) - * rminor**2 - / (rmajor * q95) - * fq - * b_plasma_toroidal_on_axis - ) - # i_plasma_current == 2 case covered above - - # Calculate cyclindrical safety factor from IPDG89 - qstar = ( - 5.0e6 - * rminor**2 - / (rmajor * plasma_current / b_plasma_toroidal_on_axis) - * 0.5 - * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) - ) - - # Calculate the poloidal field generated by the plasma current - b_plasma_poloidal_average = calculate_poloidal_field( - i_plasma_current, - plasma_current, - q95, - aspect_ratio, - eps, - b_plasma_toroidal_on_axis, - kappa, - triang, - len_plasma_poloidal, - constants.RMU0, - ) - - return b_plasma_poloidal_average, qstar, plasma_current - def outtim(self): po.oheadr(self.outfile, "Times") @@ -2599,10 +1978,12 @@ def outplas(self): "OP ", ) + self.current.output_plasma_current_models() + if physics_variables.i_alphaj == 1: po.ovarrf( self.outfile, - "Current density profile factor", + "\nCurrent density profile factor", "(alphaj)", physics_variables.alphaj, "OP ", diff --git a/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py new file mode 100644 index 000000000..3a6973ede --- /dev/null +++ b/process/models/physics/plasma_current.py @@ -0,0 +1,928 @@ +import logging +from enum import IntEnum + +import matplotlib.pyplot as plt +import numpy as np + +import process.core.io.mfile as mf +from process.core import constants +from process.core import process_output as po +from process.core.exceptions import ProcessValueError +from process.data_structure import ( + physics_variables, +) + +logger = logging.getLogger(__name__) + + +class PlasmaCurrentModel(IntEnum): + PENG_ANALYTIC_FIT = (1, "Peng analytic fit") + PENG_DIVERTOR_SCALING = (2, "Peng divertor scaling") + ITER_SCALING = (3, "Simple ITER scaling (cylindrical case)") + IPDG89_SCALING = (4, "IPDG89 scaling") + TODD_EMPIRICAL_SCALING_I = (5, "Todd empirical scaling I") + TODD_EMPIRICAL_SCALING_II = (6, "Todd empirical scaling II") + CONNOR_HASTIE_MODEL = (7, "Connor-Hastie model") + SAUTER_SCALING = (8, "Sauter scaling") + FIESTA_ST_SCALING = (9, "FIESTA ST scaling") + + def __new__(cls, value, full_name): + obj = int.__new__(cls, value) + obj._value_ = value + obj.full_name = full_name + return obj + + +class PlasmaCurrent: + """Class to hold plasma current calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + def calculate_plasma_current( + self, + alphaj: float, + alphap: float, + b_plasma_toroidal_on_axis: float, + eps: float, + i_plasma_current: int, + kappa: float, + kappa95: float, + pres_plasma_on_axis: float, + len_plasma_poloidal: float, + q95: float, + rmajor: float, + rminor: float, + triang: float, + triang95: float, + ) -> tuple[float, float, float, float, float]: + """Calculate the plasma current. + + Args: + alphaj (float): Current profile index. + alphap (float): Pressure profile index. + b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). + eps (float): Inverse aspect ratio. + i_plasma_current (int): Current scaling model to use. + 1 = Peng analytic fit + 2 = Peng divertor scaling (TART,STAR) + 3 = Simple ITER scaling + 4 = IPDG89 scaling + 5 = Todd empirical scaling I + 6 = Todd empirical scaling II + 7 = Connor-Hastie model + 8 = Sauter scaling (allowing negative triangularity) + 9 = FIESTA ST scaling + kappa (float): Plasma elongation. + kappa95 (float): Plasma elongation at 95% surface. + pres_plasma_on_axis (float): Central plasma pressure (Pa). + len_plasma_poloidal (float): Plasma perimeter length (m). + q95 (float): Plasma safety factor at 95% flux (= q-bar for i_plasma_current=2). + rmajor (float): Major radius (m). + rminor (float): Minor radius (m). + triang (float): Plasma triangularity. + triang95 (float): Plasma triangularity at 95% surface. + + Returns: + Tuple[float, float, float,]: Tuple containing b_plasma_poloidal_average, qstar, plasma_current, + + Raises: + ValueError: If invalid value for i_plasma_current is provided. + + Notes: + This routine calculates the plasma current based on the edge safety factor q95. + It will also make the current profile parameters consistent with the q-profile if required. + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729-1738. https://doi.org/10.13182/FST92-A29971 + - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - Part 1: Physics + - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO + - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant + - Sauter, Geometric formulas for systems codes..., FED 2016 + """ + # Aspect ratio + aspect_ratio = 1.0 / eps + + # Only the Sauter scaling (i_plasma_current=8) is suitable for negative triangularity: + if i_plasma_current != 8 and triang < 0.0: + raise ProcessValueError( + f"Triangularity is negative without i_plasma_current = 8 selected: {triang=}, {i_plasma_current=}" + ) + try: + model = PlasmaCurrentModel(int(physics_variables.i_plasma_current)) + # Calculate the function Fq that scales the edge q from the + # circular cross-section cylindrical case + + # Peng analytical fit + if model == PlasmaCurrentModel.PENG_ANALYTIC_FIT: + fq = self.calculate_current_coefficient_peng( + eps, len_plasma_poloidal, rminor + ) + + # Peng scaling for double null divertor; TARTs [STAR Code] + elif model == PlasmaCurrentModel.PENG_DIVERTOR_SCALING: + plasma_current = 1.0e6 * self.calculate_plasma_current_peng( + q95, + aspect_ratio, + eps, + rminor, + b_plasma_toroidal_on_axis, + kappa, + triang, + ) + + # Simple ITER scaling (simply the cylindrical case) + elif model == PlasmaCurrentModel.ITER_SCALING: + fq = 1.0 + + # ITER formula (IPDG89) + elif model == PlasmaCurrentModel.IPDG89_SCALING: + fq = self.calculate_current_coefficient_ipdg89(eps, kappa95, triang95) + + # Todd empirical scalings + # D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + elif model in [ + PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_I, + PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II, + ]: + fq = self.calculate_current_coefficient_todd( + eps, kappa95, triang95, model=1 + ) + + if model == PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II: + fq = self.calculate_current_coefficient_todd( + eps, kappa95, triang95, model=2 + ) + + # Connor-Hastie asymptotically-correct expression + elif model == PlasmaCurrentModel.CONNOR_HASTIE_MODEL: + fq = self.calculate_current_coefficient_hastie( + alphaj, + alphap, + b_plasma_toroidal_on_axis, + triang95, + eps, + kappa95, + pres_plasma_on_axis, + constants.RMU0, + ) + + # Sauter scaling allowing negative triangularity [FED May 2016] + # https://doi.org/10.1016/j.fusengdes.2016.04.033. + elif model == PlasmaCurrentModel.SAUTER_SCALING: + # Assumes zero squareness, note takes kappa, delta at separatrix not _95 + fq = self.calculate_current_coefficient_sauter(eps, kappa, triang) + + # FIESTA ST scaling + # https://doi.org/10.1016/j.fusengdes.2020.111530. + elif model == PlasmaCurrentModel.FIESTA_ST_SCALING: + fq = self.calculate_current_coefficient_fiesta(eps, kappa, triang) + + except ValueError: + raise ProcessValueError( + "Illegal value of i_plasma_current", + i_plasma_current=physics_variables.i_plasma_current, + ) from None + + # Main plasma current calculation using the fq value from the different settings + if i_plasma_current != 2: + plasma_current = ( + self.calculate_cyclindrical_plasma_current( + rminor=rminor, + rmajor=rmajor, + q95=q95, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + ) + * fq + ) + + return plasma_current + + def calculate_all_plasma_current_models( + self, + alphaj: float, + alphap: float, + b_plasma_toroidal_on_axis: float, + eps: float, + kappa: float, + kappa95: float, + pres_plasma_on_axis: float, + len_plasma_poloidal: float, + q95: float, + rmajor: float, + rminor: float, + triang: float, + triang95: float, + ) -> dict[PlasmaCurrentModel, float]: + """Calculate the plasma current for all models. + + This function calculates the plasma current for all models and returns a dictionary of the results. + + Parameters + ---------- + alphaj : + current profile index + alphap : + pressure profile index + b_plasma_toroidal_on_axis : + toroidal field on axis (T) + eps : + inverse aspect ratio + kappa : + plasma elongation + kappa95 : + plasma elongation at 95% surface + pres_plasma_on_axis : + central plasma pressure (Pa) + len_plasma_poloidal : + plasma perimeter length (m) + q95 : + plasma safety factor at 95% flux + rmajor : + major radius (m) + rminor : + minor radius (m) + triang : + plasma triangularity + triang95 : + plasma triangularity at 95% surface + Returns + ------- + dict[PlasmaCurrentModel, float] + Dictionary containing the plasma current for each model. + """ + results = {} + for model in PlasmaCurrentModel: + physics_variables.i_plasma_current = model.value + results[model] = self.calculate_plasma_current( + alphaj=alphaj, + alphap=alphap, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + eps=eps, + i_plasma_current=model.value, + kappa=kappa, + kappa95=kappa95, + pres_plasma_on_axis=pres_plasma_on_axis, + len_plasma_poloidal=len_plasma_poloidal, + q95=q95, + rmajor=rmajor, + rminor=rminor, + triang=triang, + triang95=triang95, + ) + return results + + def output_plasma_current_models(self) -> None: + """Output the plasma current for all models. + + This function outputs the plasma current for all models to the output file. + """ + plasma_currents = self.calculate_all_plasma_current_models( + alphaj=physics_variables.alphaj, + alphap=physics_variables.alphap, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + eps=physics_variables.eps, + kappa=physics_variables.kappa, + kappa95=physics_variables.kappa95, + pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis, + len_plasma_poloidal=physics_variables.len_plasma_poloidal, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + rminor=physics_variables.rminor, + triang=physics_variables.triang, + triang95=physics_variables.triang95, + ) + + physics_variables.c_plasma_peng_analytic = plasma_currents.get( + PlasmaCurrentModel.PENG_ANALYTIC_FIT + ) + physics_variables.c_plasma_peng_double_null = plasma_currents.get( + PlasmaCurrentModel.PENG_DIVERTOR_SCALING + ) + physics_variables.c_plasma_cyclindrical = plasma_currents.get( + PlasmaCurrentModel.ITER_SCALING + ) + physics_variables.c_plasma_ipdg89 = plasma_currents.get( + PlasmaCurrentModel.IPDG89_SCALING + ) + physics_variables.c_plasma_todd_empirical_i = plasma_currents.get( + PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_I + ) + physics_variables.c_plasma_todd_empirical_ii = plasma_currents.get( + PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II + ) + physics_variables.c_plasma_connor_hastie = plasma_currents.get( + PlasmaCurrentModel.CONNOR_HASTIE_MODEL + ) + physics_variables.c_plasma_sauter = plasma_currents.get( + PlasmaCurrentModel.SAUTER_SCALING + ) + physics_variables.c_plasma_fiesta_st = plasma_currents.get( + PlasmaCurrentModel.FIESTA_ST_SCALING + ) + + po.osubhd(self.outfile, "Plasma Currents using different models :") + + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.PENG_ANALYTIC_FIT.full_name}", + "(c_plasma_peng_analytic)", + physics_variables.c_plasma_peng_analytic, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.PENG_DIVERTOR_SCALING.full_name}", + "(c_plasma_peng_double_null)", + physics_variables.c_plasma_peng_double_null, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.ITER_SCALING.full_name}", + "(c_plasma_cyclindrical)", + physics_variables.c_plasma_cyclindrical, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.IPDG89_SCALING.full_name}", + "(c_plasma_ipdg89)", + physics_variables.c_plasma_ipdg89, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_I.full_name}", + "(c_plasma_todd_empirical_i)", + physics_variables.c_plasma_todd_empirical_i, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II.full_name}", + "(c_plasma_todd_empirical_ii)", + physics_variables.c_plasma_todd_empirical_ii, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.CONNOR_HASTIE_MODEL.full_name}", + "(c_plasma_connor_hastie)", + physics_variables.c_plasma_connor_hastie, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.SAUTER_SCALING.full_name}", + "(c_plasma_sauter)", + physics_variables.c_plasma_sauter, + "OP ", + ) + po.ovarre( + self.outfile, + f"{PlasmaCurrentModel.FIESTA_ST_SCALING.full_name}", + "(c_plasma_fiesta_st)", + physics_variables.c_plasma_fiesta_st, + "OP ", + ) + + @staticmethod + def plot_plasma_current_comparison(axis: plt.Axes, mfile: mf.MFile, scan: int): + """Function to plot a scatter box plot of different plasma current comparisons. + + Parameters + ---------- + axis : + Axis object to plot to. + mfile : + MFILE data object. + scan : + Scan number to use. + """ + c_plasma_peng_analytic = mfile.get("c_plasma_peng_analytic", scan=scan) + c_plasma_peng_double_null = mfile.get("c_plasma_peng_double_null", scan=scan) + c_plasma_cyclindrical = mfile.get("c_plasma_cyclindrical", scan=scan) + c_plasma_ipdg89 = mfile.get("c_plasma_ipdg89", scan=scan) + c_plasma_todd_empirical_i = mfile.get("c_plasma_todd_empirical_i", scan=scan) + c_plasma_todd_empirical_ii = mfile.get("c_plasma_todd_empirical_ii", scan=scan) + c_plasma_connor_hastie = mfile.get("c_plasma_connor_hastie", scan=scan) + c_plasma_sauter = mfile.get("c_plasma_sauter", scan=scan) + c_plasma_fiesta_st = mfile.get("c_plasma_fiesta_st", scan=scan) + + # Data for the box plot + data = { + f"{PlasmaCurrentModel.PENG_ANALYTIC_FIT.full_name}": c_plasma_peng_analytic, + f"{PlasmaCurrentModel.PENG_DIVERTOR_SCALING.full_name}": c_plasma_peng_double_null, + f"{PlasmaCurrentModel.ITER_SCALING.full_name}": c_plasma_cyclindrical, + f"{PlasmaCurrentModel.IPDG89_SCALING.full_name}": c_plasma_ipdg89, + f"{PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_I.full_name}": c_plasma_todd_empirical_i, + f"{PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II.full_name}": c_plasma_todd_empirical_ii, + f"{PlasmaCurrentModel.CONNOR_HASTIE_MODEL.full_name}": c_plasma_connor_hastie, + f"{PlasmaCurrentModel.SAUTER_SCALING.full_name}": c_plasma_sauter, + f"{PlasmaCurrentModel.FIESTA_ST_SCALING.full_name}": c_plasma_fiesta_st, + } + + # 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=(-0.9, 1)) + + # Calculate average, standard deviation, and median + data_values = list(data.values()) + avg_density_limit = np.mean(data_values) + std_density_limit = np.std(data_values) + median_density_limit = np.median(data_values) + + # Plot average, standard deviation, and median as text + axis.text( + -0.45, + 0.15, + rf"Average: {avg_density_limit * 1e-6:.4f}", + transform=axis.transAxes, + fontsize=9, + ) + axis.text( + -0.45, + 0.1, + rf"Standard Dev: {std_density_limit * 1e-6:.4f}", + transform=axis.transAxes, + fontsize=9, + ) + axis.text( + -0.45, + 0.05, + rf"Median: {median_density_limit * 1e-6:.4f}", + transform=axis.transAxes, + fontsize=9, + ) + + axis.set_title("Plasma Current Comparison") + axis.set_ylabel(r"Plasma Current [MA]") + axis.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x * 1e-6:.1f}")) + axis.set_xlim([0.5, 1.5]) + axis.set_xticks([]) + axis.set_xticklabels([]) + axis.set_facecolor("#f0f0f0") + + @staticmethod + def calculate_cyclindrical_plasma_current( + rminor: float, rmajor: float, q95: float, b_plasma_toroidal_on_axis: float + ) -> float: + """Calculate the plasma current for a cylindrical plasma. + + Parameters + ---------- + rminor : + plasma minor radius (m) + rmajor : + plasma major radius (m) + q95 : + plasma safety factor at 95% flux + b_plasma_toroidal_on_axis : + toroidal field on axis (T) + + Returns + ------- + float + plasma current (A) + + """ + + return ( + (2.0 * np.pi / constants.RMU0) + * rminor**2 + / (rmajor * q95) + * b_plasma_toroidal_on_axis + ) + + @staticmethod + def _plascar_bpol( + aspect: float, eps: float, kappa: float, delta: float + ) -> tuple[float, float, float, float]: + """Calculate the poloidal field coefficients for determining the plasma current + and poloidal field. + + + This internal function calculates the poloidal field coefficients, + which is used to calculate the poloidal field and the plasma current. + + Parameters + ---------- + aspect : + plasma aspect ratio + eps : + inverse aspect ratio + kappa : + plasma elongation + delta : + plasma triangularity + + Returns + ------- + : + coefficients ff1, ff2, d1, d2 + + References + ---------- + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729-1738. https://doi.org/10.13182/FST92-A29971 + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + + """ + # Original coding, only suitable for TARTs [STAR Code] + + c1 = (kappa**2 / (1.0 + delta)) + delta + c2 = (kappa**2 / (1.0 - delta)) - delta + + d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 + d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 + + c1_aspect = ((c1 * eps) - 1.0) if aspect < c1 else (1.0 - (c1 * eps)) + + y1 = np.sqrt(c1_aspect / (1.0 + eps)) * ((1.0 + delta) / kappa) + y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * ((1.0 - delta) / kappa) + + h2 = (1.0 + (c2 - 1.0) * (eps / 2.0)) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) + f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * ((c2 * eps) + 1.0)) + g = (eps * kappa) / (1.0 - (eps * delta)) + ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) + + h1 = (1.0 + (1.0 - c1) * (eps / 2.0)) / np.sqrt((1.0 + eps) * c1_aspect) + f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) + + if aspect < c1: + ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) + else: + ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) + + return ff1, ff2, d1, d2 + + def calculate_poloidal_field( + self, + i_plasma_current: int, + ip: float, + q95: float, + aspect: float, + eps: float, + b_plasma_toroidal_on_axis: float, + kappa: float, + delta: float, + perim: float, + ) -> float: + """Function to calculate poloidal field from the plasma current + + This function calculates the poloidal field from the plasma current in Tesla, + using a simple calculation using Ampere's law for conventional + tokamaks, or for TARTs, a scaling from Peng, Galambos and + Shipe (1992). + + Parameters + ---------- + i_plasma_current : + current scaling model to use + ip : + plasma current (A) + q95 : + 95% flux surface safety factor + aspect : + plasma aspect ratio + eps : + inverse aspect ratio + b_plasma_toroidal_on_axis : + toroidal field on axis (T) + kappa : + plasma elongation + delta : + plasma triangularity + perim : + plasma perimeter (m) + + Returns + ------- + : + poloidal field in Tesla + + + References + ---------- + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729-1738. https://doi.org/10.13182/FST92-A29971 + + """ + # Use Ampere's law using the plasma poloidal cross-section + if i_plasma_current != 2: + return constants.RMU0 * ip / perim + # Use the relation from Peng, Galambos and Shipe (1992) [STAR code] otherwise + ff1, ff2, _, _ = self._plascar_bpol(aspect, eps, kappa, delta) + + # Transform q95 to qbar + qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 + + return b_plasma_toroidal_on_axis * (ff1 + ff2) / (2.0 * np.pi * qbar) + + @staticmethod + def calculate_current_coefficient_peng( + eps: float, len_plasma_poloidal: float, rminor: float + ) -> float: + """ + Calculate the plasma current scaling coefficient for the Peng scaling from the STAR code. + + :param eps: Plasma inverse aspect ratio. + :type eps: float + :param len_plasma_poloidal: Plasma poloidal perimeter length [m]. + :type len_plasma_poloidal: float + :param rminor: Plasma minor radius [m]. + :type rminor: float + + :return: The plasma current scaling coefficient. + :rtype: float + + :references: None + """ + + return ( + (1.22 - 0.68 * eps) + / ((1.0 - eps * eps) ** 2) + * (len_plasma_poloidal / (2.0 * np.pi * rminor)) ** 2 + ) + + def calculate_plasma_current_peng( + self, + q95: float, + aspect: float, + eps: float, + rminor: float, + b_plasma_toroidal_on_axis: float, + kappa: float, + delta: float, + ) -> float: + """ + Function to calculate plasma current (Peng scaling from the STAR code) + + Parameters: + - q95: float, 95% flux surface safety factor + - aspect: float, plasma aspect ratio + - eps: float, inverse aspect ratio + - rminor: float, plasma minor radius (m) + - b_plasma_toroidal_on_axis: float, toroidal field on axis (T) + - kappa: float, plasma elongation + - delta: float, plasma triangularity + + Returns: + - float, plasma current in MA + + This function calculates the plasma current in MA, + using a scaling from Peng, Galambos and Shipe (1992). + It is primarily used for Tight Aspect Ratio Tokamaks and is + selected via i_plasma_current=2. + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729-1738. https://doi.org/10.13182/FST92-A29971 + """ + + # Transform q95 to qbar + qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 + + ff1, ff2, d1, d2 = self._plascar_bpol(aspect, eps, kappa, delta) + + e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) + e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) + + return ( + rminor + * b_plasma_toroidal_on_axis + / qbar + * 5.0 + * kappa + / (2.0 * np.pi**2) + * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) + * (ff1 + ff2) + ) + + @staticmethod + def calculate_current_coefficient_ipdg89( + eps: float, kappa95: float, triang95: float + ) -> float: + """ + Calculate the fq coefficient from the IPDG89 guidlines used in the plasma current scaling. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa95: float, plasma elongation 95% + - triang95: float, plasma triangularity 95% + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient used in the IPDG89 plasma current scaling, + based on the given plasma parameters. + + References: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 0.5 + * (1.17 - 0.65 * eps) + / ((1.0 - eps * eps) ** 2) + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + ) + + @staticmethod + def calculate_current_coefficient_todd( + eps: float, kappa95: float, triang95: float, model: int + ) -> float: + """ + Calculate the fq coefficient used in the two Todd plasma current scalings. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa95: float, plasma elongation 95% + - triang95: float, plasma triangularity 95% + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient based on the given plasma parameters for the two Todd scalings. + + References: + - D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + # Calculate the Todd scaling based on the model + base_scaling = ( + (1.0 + 2.0 * eps**2) + * ((1.0 + kappa95**2) / 2) + * ( + 1.24 + - 0.54 * kappa95 + + 0.3 * (kappa95**2 + triang95**2) + + 0.125 * triang95 + ) + ) + if model == 1: + return base_scaling + if model == 2: + return base_scaling * (1.0 + (abs(kappa95 - 1.2)) ** 3) + raise ProcessValueError(f"model = {model} is an invalid option") + + @staticmethod + def calculate_current_coefficient_hastie( + alphaj: float, + alphap: float, + b_plasma_toroidal_on_axis: float, + delta95: float, + eps: float, + kappa95: float, + pres_plasma_on_axis: float, + rmu0: float, + ) -> float: + """ + Routine to calculate the f_q coefficient for the Connor-Hastie model used for scaling the plasma current. + + Parameters: + - alphaj: float, the current profile index + - alphap: float, the pressure profile index + - b_plasma_toroidal_on_axis: float, the toroidal field on axis (T) + - delta95: float, the plasma triangularity 95% + - eps: float, the inverse aspect ratio + - kappa95: float, the plasma elongation 95% + - pres_plasma_on_axis: float, the central plasma pressure (Pa) + - rmu0: float, the vacuum permeability (H/m) + + Returns: + - float, the F coefficient + + This routine calculates the f_q coefficient used for scaling the plasma current, + using the Connor-Hastie scaling + + Reference: + - J.W.Connor and R.J.Hastie, Culham Lab Report CLM-M106 (1985). + https://scientific-publications.ukaea.uk/wp-content/uploads/CLM-M106-1.pdf + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + # Exponent in Connor-Hastie current profile + lamda = alphaj + + # Exponent in Connor-Hastie pressure profile + nu = alphap + + # Central plasma beta + beta0 = 2.0 * rmu0 * pres_plasma_on_axis / (b_plasma_toroidal_on_axis**2) + + # Plasma internal inductance + lamp1 = 1.0 + lamda + li = lamp1 / lamda * (lamp1 / lamda * np.log(lamp1) - 1.0) + + # T/r in AEA FUS 172 + kap1 = kappa95 + 1.0 + tr = kappa95 * delta95 / kap1**2 + + # E/r in AEA FUS 172 + er = (kappa95 - 1.0) / kap1 + + # T primed in AEA FUS 172 + tprime = 2.0 * tr * lamp1 / (1.0 + 0.5 * lamda) + + # E primed in AEA FUS 172 + eprime = er * lamp1 / (1.0 + lamda / 3.0) + + # Delta primed in AEA FUS 172 + deltap = (0.5 * kap1 * eps * 0.5 * li) + ( + beta0 / (0.5 * kap1 * eps) + ) * lamp1**2 / (1.0 + nu) + + # Delta/R0 in AEA FUS 172 + deltar = beta0 / 6.0 * (1.0 + 5.0 * lamda / 6.0 + 0.25 * lamda**2) + ( + 0.5 * kap1 * eps + ) ** 2 * 0.125 * (1.0 - (lamda**2) / 3.0) + + # F coefficient + return (0.5 * kap1) ** 2 * ( + 1.0 + + eps**2 * (0.5 * kap1) ** 2 + + 0.5 * deltap**2 + + 2.0 * deltar + + 0.5 * (eprime**2 + er**2) + + 0.5 * (tprime**2 + 4.0 * tr**2) + ) + + @staticmethod + def calculate_current_coefficient_sauter( + eps: float, + kappa: float, + triang: float, + ) -> float: + """ + Routine to calculate the f_q coefficient for the Sauter model used for scaling the plasma current. + + Parameters: + - eps: float, inverse aspect ratio + - kappa: float, plasma elongation at the separatrix + - triang: float, plasma triangularity at the separatrix + + Returns: + - float, the fq coefficient + + Reference: + - O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, + Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, + ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2016.04.033. + """ + w07 = 1.0 # zero squareness - can be modified later if required + + return ( + (4.1e6 / 5.0e6) + * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) + * (1.0 + 0.09 * triang + 0.16 * triang**2) + * (1.0 + 0.45 * triang * eps) + / (1.0 - 0.74 * eps) + * (1.0 + 0.55 * (w07 - 1.0)) + ) + + @staticmethod + def calculate_current_coefficient_fiesta( + eps: float, kappa: float, triang: float + ) -> float: + """ + Calculate the fq coefficient used in the FIESTA plasma current scaling. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa: float, plasma elongation at the separatrix + - triang: float, plasma triangularity at the separatrix + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient based on the given plasma parameters for the FIESTA scaling. + + References: + - S.Muldrew et.al,“PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, + Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. + """ + return 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060 diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 52f5799ef..0c3a778ed 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -28,15 +28,13 @@ Physics, PlasmaBeta, PlasmaInductance, - calculate_current_coefficient_hastie, - calculate_plasma_current_peng, - calculate_poloidal_field, diamagnetic_fraction_hender, diamagnetic_fraction_scene, ps_fraction_scene, res_diff_time, rether, ) +from process.models.physics.plasma_current import PlasmaCurrent from process.models.physics.plasma_profiles import PlasmaProfile @@ -62,6 +60,7 @@ def physics(): PlasmaDensityLimit(), PlasmaExhaust(), PlasmaBootstrapCurrent(plasma_profile=PlasmaProfile()), + PlasmaCurrent(), ) @@ -1153,10 +1152,6 @@ class PlasmaCurrentParam(NamedTuple): expected_normalised_total_beta: Any = None - expected_bp: Any = None - - expected_qstar: Any = None - expected_plasma_current: Any = None @@ -1181,8 +1176,6 @@ class PlasmaCurrentParam(NamedTuple): triang=0.5, triang95=0.33333333333333331, expected_normalised_total_beta=2.4784688886891844, - expected_bp=0.96008591022564971, - expected_qstar=2.900802902105021, expected_plasma_current=18398455.678867526, ), PlasmaCurrentParam( @@ -1203,8 +1196,6 @@ class PlasmaCurrentParam(NamedTuple): triang=0.5, triang95=0.33333333333333331, expected_normalised_total_beta=2.4784688886891844, - expected_bp=0.96008591022564971, - expected_qstar=2.900802902105021, expected_plasma_current=18398455.678867526, ), ), @@ -1224,11 +1215,11 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): monkeypatch.setattr(physics_variables, "beta_total_vol_avg", plasmacurrentparam.beta) - b_plasma_poloidal_average, qstar, plasma_current = physics.calculate_plasma_current( + plasma_current = physics.current.calculate_plasma_current( i_plasma_current=plasmacurrentparam.i_plasma_current, alphaj=plasmacurrentparam.alphaj, alphap=plasmacurrentparam.alphap, - b_plasma_toroidal_on_axis=plasmacurrentparam.b_plasma_toroidal_on_axis, + b_plasma_toroidal_on_axis=(plasmacurrentparam.b_plasma_toroidal_on_axis), eps=plasmacurrentparam.eps, kappa=plasmacurrentparam.kappa, kappa95=plasmacurrentparam.kappa95, @@ -1241,10 +1232,6 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): triang95=plasmacurrentparam.triang95, ) - assert b_plasma_poloidal_average == pytest.approx(plasmacurrentparam.expected_bp) - - assert qstar == pytest.approx(plasmacurrentparam.expected_qstar) - assert plasma_current == pytest.approx(plasmacurrentparam.expected_plasma_current) @@ -1277,8 +1264,10 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): ), ), ) -def test_calculate_plasma_current_peng(arguments, expected): - assert calculate_plasma_current_peng(**arguments) == pytest.approx(expected) +def test_calculate_plasma_current_peng(arguments, expected, physics): + assert physics.current.calculate_plasma_current_peng(**arguments) == pytest.approx( + expected + ) @pytest.mark.parametrize( @@ -1295,7 +1284,6 @@ def test_calculate_plasma_current_peng(arguments, expected): "kappa": 1.85, "delta": 0.5, "perim": 24, - "rmu0": constants.RMU0, }, 3.4401302177092803, ), @@ -1310,7 +1298,6 @@ def test_calculate_plasma_current_peng(arguments, expected): "kappa": 1.85, "delta": 0.5, "perim": 24, - "rmu0": constants.RMU0, }, 2.9310284627233205, ), @@ -1325,14 +1312,15 @@ def test_calculate_plasma_current_peng(arguments, expected): "kappa": 1.85, "delta": 0.5, "perim": 24, - "rmu0": constants.RMU0, }, 0.8377580413333333, ), ), ) -def test_calculate_poloidal_field(arguments, expected): - assert calculate_poloidal_field(**arguments) == pytest.approx(expected) +def test_calculate_poloidal_field(arguments, expected, physics): + assert physics.current.calculate_poloidal_field(**arguments) == pytest.approx( + expected + ) def test_calculate_beta_limit(): @@ -1341,8 +1329,8 @@ def test_calculate_beta_limit(): ) == pytest.approx(0.0297619) -def test_conhas(): - assert calculate_current_coefficient_hastie( +def test_conhas(physics): + assert physics.current.calculate_current_coefficient_hastie( 5, 5, 12, 0.5, 0.33, 1.85, 2e3, constants.RMU0 ) == pytest.approx(2.518876726889116) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 6bec26a1f..43ea5a0d1 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -37,6 +37,7 @@ PlasmaBeta, PlasmaInductance, ) +from process.models.physics.plasma_current import PlasmaCurrent from process.models.physics.plasma_profiles import PlasmaProfile from process.models.power import Power from process.models.stellarator.build import st_build @@ -92,6 +93,7 @@ def stellarator(): PlasmaDensityLimit(), PlasmaExhaust(), PlasmaBootstrapCurrent(plasma_profile=PlasmaProfile()), + PlasmaCurrent(), ), Neoclassics(), plasma_beta=PlasmaBeta(),