diff --git a/process/main.py b/process/main.py index 07b0ae868..db5c837f6 100644 --- a/process/main.py +++ b/process/main.py @@ -93,7 +93,7 @@ ) from process.log import logging_model_handler, show_errors from process.pfcoil import PFCoil -from process.physics import DetailedPhysics, Physics, PlasmaBeta, PlasmaInductance +from process.physics import DetailedPhysics, Physics, PlasmaBeta, PlasmaInductance, PlasmaBootstrapCurrent from process.plasma_geometry import PlasmaGeom from process.plasma_profiles import PlasmaProfile from process.power import Power @@ -684,11 +684,13 @@ def __init__(self): ) self.plasma_beta = PlasmaBeta() self.plasma_inductance = PlasmaInductance() + self.plasma_bootstrap = PlasmaBootstrapCurrent() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, plasma_beta=self.plasma_beta, plasma_inductance=self.plasma_inductance, + plasma_bootstrap=self.plasma_bootstrap, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, @@ -707,6 +709,7 @@ def __init__(self): neoclassics=self.neoclassics, plasma_beta=self.plasma_beta, plasma_inductance=self.plasma_inductance, + plasma_bootstrap=self.plasma_bootstrap, ) self.dcll = DCLL(fw=self.fw) diff --git a/process/physics.py b/process/physics.py index de5922c93..360766167 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1387,13 +1387,14 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance): + def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance, plasma_bootstrap): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive self.beta = plasma_beta self.inductance = plasma_inductance + self.bootstrap = plasma_bootstrap def physics(self): """ @@ -1684,7 +1685,7 @@ def physics(self): # Calculate bootstrap current fraction using various models current_drive_variables.f_c_plasma_bootstrap_iter89 = ( current_drive_variables.cboot - * self.bootstrap_fraction_iter89( + * self.bootstrap.bootstrap_fraction_iter89( physics_variables.aspect, physics_variables.beta_total_vol_avg, physics_variables.b_plasma_total, @@ -1698,7 +1699,7 @@ def physics(self): current_drive_variables.f_c_plasma_bootstrap_nevins = ( current_drive_variables.cboot - * self.bootstrap_fraction_nevins( + * self.bootstrap.bootstrap_fraction_nevins( physics_variables.alphan, physics_variables.alphat, physics_variables.beta_toroidal_vol_avg, @@ -1717,7 +1718,7 @@ def physics(self): # Wilson scaling uses thermal poloidal beta, not total current_drive_variables.f_c_plasma_bootstrap_wilson = ( current_drive_variables.cboot - * self.bootstrap_fraction_wilson( + * self.bootstrap.bootstrap_fraction_wilson( physics_variables.alphaj, physics_variables.alphap, physics_variables.alphat, @@ -1732,14 +1733,14 @@ def physics(self): ( current_drive_variables.f_c_plasma_bootstrap_sauter, physics_variables.j_plasma_bootstrap_sauter_profile, - ) = self.bootstrap_fraction_sauter(self.plasma_profile) + ) = self.bootstrap.bootstrap_fraction_sauter(self.plasma_profile) current_drive_variables.f_c_plasma_bootstrap_sauter *= ( current_drive_variables.cboot ) current_drive_variables.f_c_plasma_bootstrap_sakai = ( current_drive_variables.cboot - * self.bootstrap_fraction_sakai( + * self.bootstrap.bootstrap_fraction_sakai( beta_poloidal=physics_variables.beta_poloidal_vol_avg, q95=physics_variables.q95, q0=physics_variables.q0, @@ -1752,7 +1753,7 @@ def physics(self): current_drive_variables.f_c_plasma_bootstrap_aries = ( current_drive_variables.cboot - * self.bootstrap_fraction_aries( + * self.bootstrap.bootstrap_fraction_aries( beta_poloidal=physics_variables.beta_poloidal_vol_avg, ind_plasma_internal_norm=physics_variables.ind_plasma_internal_norm, core_density=physics_variables.nd_plasma_electron_on_axis, @@ -1763,7 +1764,7 @@ def physics(self): current_drive_variables.f_c_plasma_bootstrap_andrade = ( current_drive_variables.cboot - * self.bootstrap_fraction_andrade( + * self.bootstrap.bootstrap_fraction_andrade( beta_poloidal=physics_variables.beta_poloidal_vol_avg, core_pressure=physics_variables.pres_plasma_thermal_on_axis, average_pressure=physics_variables.pres_plasma_thermal_vol_avg, @@ -1772,7 +1773,7 @@ def physics(self): ) current_drive_variables.f_c_plasma_bootstrap_hoang = ( current_drive_variables.cboot - * self.bootstrap_fraction_hoang( + * self.bootstrap.bootstrap_fraction_hoang( beta_poloidal=physics_variables.beta_poloidal_vol_avg, pressure_index=physics_variables.alphap, current_index=physics_variables.alphaj, @@ -1781,7 +1782,7 @@ def physics(self): ) current_drive_variables.f_c_plasma_bootstrap_wong = ( current_drive_variables.cboot - * self.bootstrap_fraction_wong( + * self.bootstrap.bootstrap_fraction_wong( beta_poloidal=physics_variables.beta_poloidal_vol_avg, density_index=physics_variables.alphan, temperature_index=physics_variables.alphat, @@ -1791,7 +1792,7 @@ def physics(self): ) current_drive_variables.bscf_gi_i = ( current_drive_variables.cboot - * self.bootstrap_fraction_gi_I( + * self.bootstrap.bootstrap_fraction_gi_I( beta_poloidal=physics_variables.beta_poloidal_vol_avg, pressure_index=physics_variables.alphap, temperature_index=physics_variables.alphat, @@ -1804,7 +1805,7 @@ def physics(self): current_drive_variables.bscf_gi_ii = ( current_drive_variables.cboot - * self.bootstrap_fraction_gi_II( + * self.bootstrap.bootstrap_fraction_gi_II( beta_poloidal=physics_variables.beta_poloidal_vol_avg, pressure_index=physics_variables.alphap, temperature_index=physics_variables.alphat, @@ -1814,7 +1815,7 @@ def physics(self): ) current_drive_variables.f_c_plasma_bootstrap_sugiyama_l = ( current_drive_variables.cboot - * self.bootstrap_fraction_sugiyama_l_mode( + * self.bootstrap.bootstrap_fraction_sugiyama_l_mode( eps=physics_variables.eps, beta_poloidal=physics_variables.beta_poloidal_vol_avg, alphan=physics_variables.alphan, @@ -1826,7 +1827,7 @@ def physics(self): ) current_drive_variables.f_c_plasma_bootstrap_sugiyama_h = ( current_drive_variables.cboot - * self.bootstrap_fraction_sugiyama_h_mode( + * self.bootstrap.bootstrap_fraction_sugiyama_h_mode( eps=physics_variables.eps, beta_poloidal=physics_variables.beta_poloidal_vol_avg, alphan=physics_variables.alphan, @@ -1842,31 +1843,19 @@ def physics(self): ) ) - bootstrap_map = { - 0: current_drive_variables.f_c_plasma_bootstrap, - 1: current_drive_variables.f_c_plasma_bootstrap_iter89, - 2: current_drive_variables.f_c_plasma_bootstrap_nevins, - 3: current_drive_variables.f_c_plasma_bootstrap_wilson, - 4: current_drive_variables.f_c_plasma_bootstrap_sauter, - 5: current_drive_variables.f_c_plasma_bootstrap_sakai, - 6: current_drive_variables.f_c_plasma_bootstrap_aries, - 7: current_drive_variables.f_c_plasma_bootstrap_andrade, - 8: current_drive_variables.f_c_plasma_bootstrap_hoang, - 9: current_drive_variables.f_c_plasma_bootstrap_wong, - 10: current_drive_variables.bscf_gi_i, - 11: current_drive_variables.bscf_gi_ii, - 12: current_drive_variables.f_c_plasma_bootstrap_sugiyama_l, - 13: current_drive_variables.f_c_plasma_bootstrap_sugiyama_h, - } - if int(physics_variables.i_bootstrap_current) in bootstrap_map: - current_drive_variables.f_c_plasma_bootstrap = bootstrap_map[ + # Calculate beta_norm_max based on i_beta_norm_max + try: + model = BootstrapCurrentFractionModel( int(physics_variables.i_bootstrap_current) - ] - else: + ) + current_drive_variables.f_c_plasma_bootstrap = ( + self.bootstrap.get_bootstrap_current_fraction_value(model) + ) + except ValueError: raise ProcessValueError( "Illegal value of i_bootstrap_current", i_bootstrap_current=physics_variables.i_bootstrap_current, - ) + ) from None physics_variables.err242 = 0 if ( @@ -5522,254 +5511,7 @@ def outplas(self): self.inductance.output_volt_second_information() - po.ovarrf( - self.outfile, - "Bootstrap current fraction multiplier", - "(cboot)", - current_drive_variables.cboot, - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (ITER 1989)", - "(f_c_plasma_bootstrap_iter89)", - current_drive_variables.f_c_plasma_bootstrap_iter89, - "OP ", - ) - - po.ovarrf( - self.outfile, - "Bootstrap fraction (Sauter et al)", - "(f_c_plasma_bootstrap_sauter)", - current_drive_variables.f_c_plasma_bootstrap_sauter, - "OP ", - ) - for point in range(len(physics_variables.j_plasma_bootstrap_sauter_profile)): - po.ovarrf( - self.mfile, - f"Sauter et al bootstrap current density profile at point {point}", - f"(j_plasma_bootstrap_sauter_profile{point})", - physics_variables.j_plasma_bootstrap_sauter_profile[point], - "OP ", - ) - - po.ovarrf( - self.outfile, - "Bootstrap fraction (Nevins et al)", - "(f_c_plasma_bootstrap_nevins)", - current_drive_variables.f_c_plasma_bootstrap_nevins, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (Wilson)", - "(f_c_plasma_bootstrap_wilson)", - current_drive_variables.f_c_plasma_bootstrap_wilson, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (Sakai)", - "(f_c_plasma_bootstrap_sakai)", - current_drive_variables.f_c_plasma_bootstrap_sakai, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (ARIES)", - "(f_c_plasma_bootstrap_aries)", - current_drive_variables.f_c_plasma_bootstrap_aries, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (Andrade)", - "(f_c_plasma_bootstrap_andrade)", - current_drive_variables.f_c_plasma_bootstrap_andrade, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (Hoang)", - "(f_c_plasma_bootstrap_hoang)", - current_drive_variables.f_c_plasma_bootstrap_hoang, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (Wong)", - "(f_c_plasma_bootstrap_wong)", - current_drive_variables.f_c_plasma_bootstrap_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, - "Bootstrap fraction (Sugiyama L-mode)", - "(f_c_plasma_bootstrap_sugiyama_l)", - current_drive_variables.f_c_plasma_bootstrap_sugiyama_l, - "OP ", - ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (Sugiyama H-mode)", - "(f_c_plasma_bootstrap_sugiyama_h)", - current_drive_variables.f_c_plasma_bootstrap_sugiyama_h, - "OP ", - ) - - po.ovarrf( - self.outfile, - "Diamagnetic fraction (Hender)", - "(f_c_plasma_diamagnetic_hender)", - current_drive_variables.f_c_plasma_diamagnetic_hender, - "OP ", - ) - po.ovarrf( - self.outfile, - "Diamagnetic fraction (SCENE)", - "(f_c_plasma_diamagnetic_scene)", - current_drive_variables.f_c_plasma_diamagnetic_scene, - "OP ", - ) - po.ovarrf( - self.outfile, - "Pfirsch-Schlueter fraction (SCENE)", - "(f_c_plasma_pfirsch_schluter_scene)", - current_drive_variables.f_c_plasma_pfirsch_schluter_scene, - "OP ", - ) - # Error to catch if bootstap fraction limit has been enforced - if physics_variables.err242 == 1: - logger.error("Bootstrap fraction upper limit enforced") - - # Error to catch if self-driven current fraction limit has been enforced - if physics_variables.err243 == 1: - logger.error( - "Predicted plasma driven current is more than upper limit on non-inductive fraction" - ) - - if physics_variables.i_bootstrap_current == 0: - po.ocmmnt( - self.outfile, " (User-specified bootstrap current fraction used)" - ) - elif physics_variables.i_bootstrap_current == 1: - po.ocmmnt( - self.outfile, " (ITER 1989 bootstrap current fraction model used)" - ) - elif physics_variables.i_bootstrap_current == 2: - po.ocmmnt( - self.outfile, - " (Nevins et al bootstrap current fraction model used)", - ) - elif physics_variables.i_bootstrap_current == 3: - po.ocmmnt( - self.outfile, " (Wilson bootstrap current fraction model used)" - ) - elif physics_variables.i_bootstrap_current == 4: - po.ocmmnt( - self.outfile, - " (Sauter et al bootstrap current fraction model used)", - ) - elif physics_variables.i_bootstrap_current == 5: - po.ocmmnt( - 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( - self.outfile, " (Diamagnetic current fraction not calculated)" - ) - # Error to show if diamagnetic current is above 1% but not used - if current_drive_variables.f_c_plasma_diamagnetic_scene > 0.01e0: - logger.error( - "Diamagnetic fraction is more than 1%, but not calculated. " - "Consider using i_diamagnetic_current=2 and i_pfirsch_schluter_current=1" - ) - - elif physics_variables.i_diamagnetic_current == 1: - po.ocmmnt( - self.outfile, " (Hender diamagnetic current fraction scaling used)" - ) - elif physics_variables.i_diamagnetic_current == 2: - po.ocmmnt( - self.outfile, " (SCENE diamagnetic current fraction scaling used)" - ) - - if physics_variables.i_pfirsch_schluter_current == 0: - po.ocmmnt( - self.outfile, " Pfirsch-Schluter current fraction not calculated" - ) - elif physics_variables.i_pfirsch_schluter_current == 1: - po.ocmmnt( - self.outfile, - " (SCENE Pfirsch-Schluter current fraction scaling used)", - ) - - po.ovarrf( - self.outfile, - "Bootstrap fraction (enforced)", - "(f_c_plasma_bootstrap.)", - current_drive_variables.f_c_plasma_bootstrap, - "OP ", - ) - po.ovarrf( - self.outfile, - "Diamagnetic fraction (enforced)", - "(f_c_plasma_diamagnetic.)", - current_drive_variables.f_c_plasma_diamagnetic, - "OP ", - ) - po.ovarrf( - self.outfile, - "Pfirsch-Schlueter fraction (enforced)", - "(f_c_plasma_pfirsch_schluter.)", - current_drive_variables.f_c_plasma_pfirsch_schluter, - "OP ", - ) + self.bootstrap.output_bootstrap_current_information() po.osubhd(self.outfile, "Fuelling :") po.ovarre( @@ -5927,3344 +5669,3637 @@ def output_confinement_comparison(self, istell: int) -> None: po.oblnkl(self.outfile) po.ostars(self.outfile, 110) - @staticmethod - def bootstrap_fraction_iter89( - aspect: float, - beta: float, - b_plasma_toroidal_on_axis: float, - plasma_current: float, - q95: float, - q0: float, - rmajor: float, - vol_plasma: float, - ) -> float: + def find_other_h_factors(self, i_confinement_time: int) -> float: """ - Calculate the bootstrap-driven fraction of the plasma current. + Function to find H-factor for the equivalent confinement time in other scalings. Args: - aspect (float): Plasma aspect ratio. - beta (float): Plasma total beta. - b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). - plasma_current (float): Plasma current (A). - q95 (float): Safety factor at 95% surface. - q0 (float): Central safety factor. - rmajor (float): Plasma major radius (m). - vol_plasma (float): Plasma volume (m3). + i_confinement_time (int): Index of the confinement time scaling to use. Returns: - float: The bootstrap-driven fraction of the plasma current. - - This function performs the original ITER calculation of the plasma current bootstrap fraction. - - Reference: - - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, - - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + float: The calculated H-factor. """ - # Calculate the bootstrap current coefficient - c_bs = 1.32 - 0.235 * (q95 / q0) + 0.0185 * (q95 / q0) ** 2 + def fhz(hfact: float) -> float: + """ + Function used to find power balance. - # Calculate the average minor radius - average_a = np.sqrt(vol_plasma / (2 * np.pi**2 * rmajor)) + Args: + hfact (float): H-factor to be used in the calculation. - b_pa = (plasma_current / 1e6) / (5 * average_a) + Returns: + float: The difference between the calculated power and the required power for balance. + """ + ( + ptrez, + ptriz, + _, + _, + _, + _, + ) = self.calculate_confinement_time( + physics_variables.m_fuel_amu, + physics_variables.p_alpha_total_mw, + physics_variables.aspect, + physics_variables.b_plasma_toroidal_on_axis, + physics_variables.nd_plasma_ions_total_vol_avg, + physics_variables.nd_plasma_electrons_vol_avg, + physics_variables.nd_plasma_electron_line, + physics_variables.eps, + hfact, + i_confinement_time, + physics_variables.i_plasma_ignited, + physics_variables.kappa, + physics_variables.kappa95, + physics_variables.p_non_alpha_charged_mw, + current_drive_variables.p_hcd_injected_total_mw, + physics_variables.plasma_current, + physics_variables.pden_plasma_core_rad_mw, + physics_variables.rmajor, + physics_variables.rminor, + physics_variables.temp_plasma_electron_density_weighted_kev, + physics_variables.temp_plasma_ion_density_weighted_kev, + physics_variables.q95, + physics_variables.qstar, + physics_variables.vol_plasma, + physics_variables.n_charge_plasma_effective_vol_avg, + ) - # Calculate the poloidal beta for bootstrap current - betapbs = beta * (b_plasma_toroidal_on_axis / b_pa) ** 2 + # At power balance, fhz is zero. + fhz_value = ( + ptrez + + ptriz + - physics_variables.f_p_alpha_plasma_deposited + * physics_variables.pden_alpha_total_mw + - physics_variables.pden_non_alpha_charged_mw + - physics_variables.pden_plasma_ohmic_mw + ) - # Ensure betapbs is positive - if betapbs <= 0.0: - return 0.0 + # Take into account whether injected power is included in tau_e calculation (i.e. whether device is ignited) + if physics_variables.i_plasma_ignited == 0: + fhz_value -= ( + current_drive_variables.p_hcd_injected_total_mw + / physics_variables.vol_plasma + ) - # Calculate and return the bootstrap current fraction - return c_bs * (betapbs / np.sqrt(aspect)) ** 1.3 + # Include the radiation power if requested + if physics_variables.i_rad_loss == 0: + fhz_value += physics_variables.pden_plasma_rad_mw + elif physics_variables.i_rad_loss == 1: + fhz_value += physics_variables.pden_plasma_core_rad_mw - @staticmethod - def bootstrap_fraction_wilson( - alphaj: float, - alphap: float, - alphat: float, - betpth: float, - q0: float, - q95: float, + return fhz_value + + return root_scalar(fhz, bracket=(0.01, 150), xtol=0.001).root + + @staticmethod + def calculate_confinement_time( + m_fuel_amu: float, + p_alpha_total_mw: float, + aspect: float, + b_plasma_toroidal_on_axis: float, + nd_plasma_ions_total_vol_avg: float, + nd_plasma_electrons_vol_avg: float, + nd_plasma_electron_line: float, + eps: float, + hfact: float, + i_confinement_time: int, + i_plasma_ignited: int, + kappa: float, + kappa95: float, + p_non_alpha_charged_mw: float, + p_hcd_injected_total_mw: float, + plasma_current: float, + pden_plasma_core_rad_mw: float, rmajor: float, rminor: float, - ) -> float: + temp_plasma_electron_density_weighted_kev: float, + temp_plasma_ion_density_weighted_kev: float, + q95: float, + qstar: float, + vol_plasma: float, + zeff: float, + ) -> tuple[float, float, float, float, float, float, float]: """ - Bootstrap current fraction from Wilson et al scaling - - Args: - alphaj (float): Current profile index. - alphap (float): Pressure profile index. - alphat (float): Temperature profile index. - betpth (float): Thermal component of poloidal beta. - q0 (float): Safety factor on axis. - q95 (float): Edge safety factor. - rmajor (float): Major radius (m). - rminor (float): Minor radius (m). - - Returns: - float: The bootstrap current fraction. + Calculate the confinement times and the transport power loss terms. - This function calculates the bootstrap current fraction using the numerically fitted algorithm written by Howard Wilson. + :param m_fuel_amu: Average mass of fuel (amu) + :param p_alpha_total_mw: Alpha particle power (MW) + :param aspect: Aspect ratio + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T) + :param nd_plasma_ions_total_vol_avg: Total ion density (/m3) + :param nd_plasma_electrons_vol_avg: Volume averaged electron density (/m3) + :param nd_plasma_electron_line: Line-averaged electron density (/m3) + :param eps: Inverse aspect ratio + :param hfact: H factor on energy confinement scalings + :param i_confinement_time: Switch for energy confinement scaling to use + :param i_plasma_ignited: Switch for ignited calculation + :param kappa: Plasma elongation + :param kappa95: Plasma elongation at 95% surface + :param p_non_alpha_charged_mw: Non-alpha charged particle fusion power (MW) + :param p_hcd_injected_total_mw: Auxiliary power to ions and electrons (MW) + :param plasma_current: Plasma current (A) + :param pden_plasma_core_rad_mw: Total core radiation power (MW/m3) + :param q95: Edge safety factor (tokamaks), or rotational transform iotabar (stellarators) + :param qstar: Equivalent cylindrical edge safety factor + :param rmajor: Plasma major radius (m) + :param rminor: Plasma minor radius (m) + :param temp_plasma_electron_density_weighted_kev: Density weighted average electron temperature (keV) + :param temp_plasma_ion_density_weighted_kev: Density weighted average ion temperature (keV) + :param vol_plasma: Plasma volume (m3) + :param a_plasma_poloidal: Plasma cross-sectional area (m2) + :param zeff: Plasma effective charge - Reference: - - AEA FUS 172: Physics Assessment for the European Reactor Study, 1989 - - H. R. Wilson, Nuclear Fusion 32 (1992) 257 + :return: Tuple containing: + - pden_electron_transport_loss_mw (float): Electron transport power (MW/m3) + - pden_ion_transport_loss_mw (float): Ion transport power (MW/m3) + - t_electron_energy_confinement (float): Electron energy confinement time (s) + - t_ion_energy_confinement (float): Ion energy confinement time (s) + - t_energy_confinement (float): Global energy confinement time (s) + - p_plasma_loss_mw (float): Heating power (MW) assumed in calculation """ - term1 = np.log(0.5) - term2 = np.log(q0 / q95) - - # Re-arranging of parabolic profile to be equal to (r/a)^2 where the profile value is half of the core value - - termp = 1.0 - 0.5 ** (1.0 / alphap) - termt = 1.0 - 0.5 ** (1.0 / alphat) - termj = 1.0 - 0.5 ** (1.0 / alphaj) + # ======================================================================== - # Assuming a parabolic safety factor profile of the form q = q0 + (q95 - q0) * (r/a)^2 - # Substitute (r/a)^2 term from temperature,pressure and current profiles into q profile when values is 50% of core value - # Take natural log of q profile over q95 and q0 to get the profile index + # Calculate heating power (MW) + p_plasma_loss_mw = ( + physics_variables.f_p_alpha_plasma_deposited * p_alpha_total_mw + + p_non_alpha_charged_mw + + physics_variables.p_plasma_ohmic_mw + ) - alfpnw = term1 / np.log(np.log((q0 + (q95 - q0) * termp) / q95) / term2) - alftnw = term1 / np.log(np.log((q0 + (q95 - q0) * termt) / q95) / term2) - aj = term1 / np.log(np.log((q0 + (q95 - q0) * termj) / q95) / term2) + # If the device is not ignited, add the injected auxiliary power + if i_plasma_ignited == 0: + p_plasma_loss_mw = p_plasma_loss_mw + p_hcd_injected_total_mw - # Crude check for NaN errors or other illegal values. - if np.isnan(aj) or np.isnan(alfpnw) or np.isnan(alftnw) or aj < 0: - raise ProcessValueError( - "Illegal profile value found", aj=aj, alfpnw=alfpnw, alftnw=alftnw + # Include the radiation as a loss term if requested + if physics_variables.i_rad_loss == 0: + p_plasma_loss_mw = ( + p_plasma_loss_mw - physics_variables.pden_plasma_rad_mw * vol_plasma ) + elif physics_variables.i_rad_loss == 1: + p_plasma_loss_mw = ( + p_plasma_loss_mw - pden_plasma_core_rad_mw * vol_plasma + ) # shouldn't this be vol_core instead of vol_plasma? + # else do not adjust p_plasma_loss_mw for radiation - # Ratio of ionic charge to electron charge + # Ensure heating power is positive (shouldn't be necessary) + p_plasma_loss_mw = max(p_plasma_loss_mw, 1.0e-3) - z = 1.0 + # ======================================================================== - # Inverse aspect ratio: r2 = maximum plasma radius, r1 = minimum - # This is the definition used in the paper - r2 = rmajor + rminor - r1 = rmajor - rminor - eps1 = (r2 - r1) / (r2 + r1) + # Line averaged electron density in scaled units + dnla20 = nd_plasma_electron_line * 1.0e-20 + dnla19 = nd_plasma_electron_line * 1.0e-19 - # Coefficients fitted using least squares techniques + # Volume averaged electron density in units of 10**20 m**-3 + n20 = nd_plasma_electrons_vol_avg / 1.0e20 - # Square root of current profile index term - saj = np.sqrt(aj) + # Plasma current in MA + pcur = plasma_current / 1.0e6 - a = np.array([ - 1.41 * (1.0 - 0.28 * saj) * (1.0 + 0.12 / z), - 0.36 * (1.0 - 0.59 * saj) * (1.0 + 0.8 / z), - -0.27 * (1.0 - 0.47 * saj) * (1.0 + 3.0 / z), - 0.0053 * (1.0 + 5.0 / z), - -0.93 * (1.0 - 0.34 * saj) * (1.0 + 0.15 / z), - -0.26 * (1.0 - 0.57 * saj) * (1.0 - 0.27 * z), - 0.064 * (1.0 - 0.6 * aj + 0.15 * aj * aj) * (1.0 + 7.6 / z), - -0.0011 * (1.0 + 9.0 / z), - -0.33 * (1.0 - aj + 0.33 * aj * aj), - -0.26 * (1.0 - 0.87 / saj - 0.16 * aj), - -0.14 * (1.0 - 1.14 / saj - 0.45 * saj), - -0.0069, - ]) + # Separatrix kappa defined with plasma volume for IPB scalings + # Updated version of kappa used by the IPB98 scalings correction in: - seps1 = np.sqrt(eps1) + # None Otto Kardaun, N. K. Thomsen, and None Alexander Chudnovskiy, + # “Corrections to a sequence of papers in Nuclear Fusion,” Nuclear Fusion, + # vol. 48, no. 9, pp. 099801099801, Aug. 2008, + # doi: https://doi.org/10.1088/0029-5515/48/9/099801. - b = np.array([ - 1.0, - alfpnw, - alftnw, - alfpnw * alftnw, - seps1, - alfpnw * seps1, - alftnw * seps1, - alfpnw * alftnw * seps1, - eps1, - alfpnw * eps1, - alftnw * eps1, - alfpnw * alftnw * eps1, - ]) + physics_variables.kappa_ipb = (vol_plasma / (2.0 * np.pi * rmajor)) / ( + np.pi * rminor**2 + ) - # Empirical bootstrap current fraction - return seps1 * betpth * (a * b).sum() + # Electron energy confinement times - @staticmethod - def bootstrap_fraction_nevins( - alphan: float, - alphat: float, - beta_toroidal: float, - b_plasma_toroidal_on_axis: float, - nd_plasma_electrons_vol_avg: float, - plasma_current: float, - q95: float, - q0: float, - rmajor: float, - rminor: float, - te: float, - zeff: float, - ) -> float: - """ - Calculate the bootstrap current fraction from Nevins et al scaling. + # ======================================================================== - Args: - alphan (float): Density profile index. - alphat (float): Temperature profile index. - beta_toroidal (float): Toroidal plasma beta. - b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). - nd_plasma_electrons_vol_avg (float): Electron density (/m3). - plasma_current (float): Plasma current (A). - q0 (float): Central safety factor. - q95 (float): Safety factor at 95% surface. - rmajor (float): Plasma major radius (m). - rminor (float): Plasma minor radius (m). - te (float): Volume averaged plasma temperature (keV). - zeff (float): Plasma effective charge. + # User defined confinement time + if i_confinement_time == 0: # t_electron_energy_confinement is an input + t_electron_confinement = physics_variables.tauee_in - Returns: - float: The bootstrap current fraction. + # ======================================================================== - This function calculates the bootstrap current fraction using the Nevins et al method. + # Nec-Alcator(NA) OH scaling + if i_confinement_time == 1: + t_electron_confinement = confinement.neo_alcator_confinement_time( + n20, rminor, rmajor, qstar + ) - Reference: See appendix of: - Keii Gi, Makoto Nakamura, Kenji Tobita, Yasushi Ono, - Bootstrap current fraction scaling for a tokamak reactor design study, - Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, - ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2014.07.009. + # ======================================================================== - Nevins, W. M. "Summary report: ITER specialists' meeting on heating and current drive." - ITER-TN-PH-8-4, June 1988. 1988. + # "Mirnov"-like scaling (H-mode) + elif i_confinement_time == 2: # Mirnov scaling (H-mode) + t_electron_confinement = confinement.mirnov_confinement_time( + rminor, kappa95, pcur + ) - """ - # Calculate peak electron beta at plasma centre, this is not the form used in the paper - # The paper assumes parabolic profiles for calculating core values with the profile indexes. - # We instead use the directly calculated electron density and temperature values at the core. - # So that it is compatible with all profiles + # ======================================================================== - betae0 = ( - physics_variables.nd_plasma_electron_on_axis - * physics_variables.temp_plasma_electron_on_axis_kev - * 1.0e3 - * constants.ELECTRON_CHARGE - / (b_plasma_toroidal_on_axis**2 / (2.0 * constants.RMU0)) - ) + # Merezhkin-Mukhovatov (MM) OH/L-mode scaling + elif i_confinement_time == 3: + t_electron_confinement = confinement.merezhkin_muhkovatov_confinement_time( + rmajor, + rminor, + kappa95, + qstar, + dnla20, + m_fuel_amu, + temp_plasma_electron_density_weighted_kev, + ) - # Call integration routine using definite integral routine from scipy + # ======================================================================== - ainteg, _ = integrate.quad( - lambda y: _nevins_integral( - y, - nd_plasma_electrons_vol_avg, - te, - b_plasma_toroidal_on_axis, - rminor, + # Shimomura (S) optimized H-mode scaling + elif i_confinement_time == 4: + t_electron_confinement = confinement.shimomura_confinement_time( + rmajor, rminor, b_plasma_toroidal_on_axis, kappa95, m_fuel_amu + ) + + # ======================================================================== + + # Kaye-Goldston scaling (L-mode) + elif i_confinement_time == 5: + t_electron_confinement = confinement.kaye_goldston_confinement_time( + pcur, rmajor, - zeff, - alphat, - alphan, - q0, - q95, - beta_toroidal, - ), - 0, # Lower bound - 1.0, # Upper bound - ) + rminor, + kappa, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ) - # Calculate bootstrap current and fraction + # ======================================================================== - aibs = 2.5 * betae0 * rmajor * b_plasma_toroidal_on_axis * q95 * ainteg - return 1.0e6 * aibs / plasma_current + # ITER Power scaling - ITER 89-P (L-mode) + elif i_confinement_time == 6: + t_electron_confinement = confinement.iter_89p_confinement_time( + pcur, + rmajor, + rminor, + kappa, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ) - @staticmethod - def bootstrap_fraction_sauter(plasma_profile: float) -> float: - """Calculate the bootstrap current fraction from the Sauter et al scaling. + # ======================================================================== - Args: - plasma_profile (PlasmaProfile): The plasma profile object. + # ITER Offset linear scaling - ITER 89-O (L-mode) + elif i_confinement_time == 7: + t_electron_confinement = confinement.iter_89_0_confinement_time( + pcur, + rmajor, + rminor, + kappa, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ) + # ======================================================================== - Returns: - float: The bootstrap current fraction. + # Rebut-Lallia offset linear scaling (L-mode) + elif i_confinement_time == 8: + t_electron_confinement = confinement.rebut_lallia_confinement_time( + rminor, + rmajor, + kappa, + m_fuel_amu, + pcur, + zeff, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + ) - This function calculates the bootstrap current fraction using the Sauter, Angioni, and Lin-Liu scaling. + # ======================================================================== - Reference: - - O. Sauter, C. Angioni, Y. R. Lin-Liu; - Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. - Phys. Plasmas 1 July 1999; 6 (7): 2834-2839. https://doi.org/10.1063/1.873240 - - 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 + # Goldston scaling (L-mode) + elif i_confinement_time == 9: # Goldston scaling (L-mode) + t_electron_confinement = confinement.goldston_confinement_time( + pcur, rmajor, rminor, kappa95, m_fuel_amu, p_plasma_loss_mw + ) - Note: - The code was supplied by Emiliano Fable, IPP Garching (private communication). - """ + # ======================================================================== - # Radial points from 0 to 1 seperated by 1/profile_size - roa = plasma_profile.neprofile.profile_x + # T-10 scaling (L-mode) + elif i_confinement_time == 10: + t_electron_confinement = confinement.t10_confinement_time( + dnla20, + rmajor, + qstar, + b_plasma_toroidal_on_axis, + rminor, + kappa95, + p_plasma_loss_mw, + zeff, + pcur, + ) - # Local circularised minor radius - rho = np.sqrt(physics_variables.a_plasma_poloidal / np.pi) * roa + # ======================================================================== - # Square root of local aspect ratio - sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) + # JAERI / Odajima-Shimomura L-mode scaling + elif i_confinement_time == 11: # JAERI scaling + t_electron_confinement = confinement.jaeri_confinement_time( + kappa95, + rminor, + m_fuel_amu, + n20, + pcur, + b_plasma_toroidal_on_axis, + rmajor, + qstar, + zeff, + p_plasma_loss_mw, + ) - # Calculate electron and ion density profiles - ne = plasma_profile.neprofile.profile_y * 1e-19 - ni = ( - physics_variables.nd_plasma_ions_total_vol_avg - / physics_variables.nd_plasma_electrons_vol_avg - ) * ne + # ======================================================================== - # Calculate electron and ion temperature profiles - tempe = plasma_profile.teprofile.profile_y - tempi = ( - physics_variables.temp_plasma_ion_vol_avg_kev - / physics_variables.temp_plasma_electron_vol_avg_kev - ) * tempe + # Kaye "big" L-mode scaling (based only on big tokamak data) + elif i_confinement_time == 12: + t_electron_confinement = confinement.kaye_big_confinement_time( + rmajor, + rminor, + b_plasma_toroidal_on_axis, + kappa95, + pcur, + n20, + m_fuel_amu, + p_plasma_loss_mw, + ) - # Flat Zeff profile assumed - # Return tempi like array object filled with zeff - zeff = np.full_like(tempi, physics_variables.n_charge_plasma_effective_vol_avg) + # ======================================================================== - # inverse_q = 1/safety factor - # Parabolic q profile assumed - inverse_q = 1 / ( - physics_variables.q0 - + (physics_variables.q95 - physics_variables.q0) * roa**2 - ) - # Create new array of average mass of fuel portion of ions - amain = np.full_like(inverse_q, physics_variables.m_ions_total_amu) + # ITER H90-P H-mode scaling + elif i_confinement_time == 13: + t_electron_confinement = confinement.iter_h90_p_confinement_time( + pcur, + rmajor, + rminor, + kappa, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ) - # Create new array of average main ion charge - zmain = np.full_like(inverse_q, 1.0 + physics_variables.f_plasma_fuel_helium3) + # ======================================================================== - # Calculate total bootstrap current (MA) by summing along profiles - # Looping from 2 because _calculate_l31_coefficient() etc should return 0 @ j == 1 - radial_elements = np.arange(2, plasma_profile.profile_size) + # Minimum of ITER 89-P and ITER 89-O + elif i_confinement_time == 14: + t_electron_confinement = min( + confinement.iter_89p_confinement_time( + pcur, + rmajor, + rminor, + kappa, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ), + confinement.iter_89_0_confinement_time( + pcur, + rmajor, + rminor, + kappa, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ), + ) - # Change in localised minor radius to be used as delta term in derivative - drho = rho[radial_elements] - rho[radial_elements - 1] + # ======================================================================== - # Area of annulus, assuming circular plasma cross-section - da = 2 * np.pi * rho[radial_elements - 1] * drho # area of annulus + # Riedel scaling (L-mode) + elif i_confinement_time == 15: + t_electron_confinement = confinement.riedel_l_confinement_time( + pcur, + rmajor, + rminor, + kappa95, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + ) - # Create the partial derivatives using numpy gradient (central differences) - dlogte_drho = np.gradient(np.log(tempe), rho)[radial_elements - 1] - dlogti_drho = np.gradient(np.log(tempi), rho)[radial_elements - 1] - dlogne_drho = np.gradient(np.log(ne), rho)[radial_elements - 1] + # ======================================================================== - jboot = ( - 0.5 - * ( - _calculate_l31_coefficient( - radial_elements, - plasma_profile.profile_size, - physics_variables.rmajor, - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.triang, - ne, - ni, - tempe, - tempi, - inverse_q, - rho, - zeff, - sqeps, - ) - * dlogne_drho - + _calculate_l31_32_coefficient( - radial_elements, - plasma_profile.profile_size, - physics_variables.rmajor, - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.triang, - ne, - ni, - tempe, - tempi, - inverse_q, - rho, - zeff, - sqeps, - ) - * dlogte_drho - + _calculate_l34_alpha_31_coefficient( - radial_elements, - plasma_profile.profile_size, - physics_variables.rmajor, - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.triang, - inverse_q, - sqeps, - tempi, - tempe, - amain, - zmain, - ni, - ne, - rho, - zeff, - ) - * dlogti_drho - ) - * 1.0e6 - * ( - -physics_variables.b_plasma_toroidal_on_axis - / (0.2 * np.pi * physics_variables.rmajor) - * rho[radial_elements - 1] - * inverse_q[radial_elements - 1] + # Christiansen et al scaling (L-mode) + elif i_confinement_time == 16: + t_electron_confinement = confinement.christiansen_confinement_time( + pcur, + rmajor, + rminor, + kappa95, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + m_fuel_amu, ) - ) # A/m2 - - return (np.sum(da * jboot, axis=0) / physics_variables.plasma_current), jboot - - @staticmethod - def bootstrap_fraction_sakai( - beta_poloidal: float, - q95: float, - q0: float, - alphan: float, - alphat: float, - eps: float, - ind_plasma_internal_norm: float, - ) -> float: - """ - Calculate the bootstrap fraction using the Sakai formula. - - Parameters: - beta_poloidal (float): Plasma poloidal beta. - q95 (float): Safety factor at 95% of the plasma radius. - q0 (float): Safety factor at the magnetic axis. - alphan (float): Density profile index - alphat (float): Temperature profile index - eps (float): Inverse aspect ratio. - ind_plasma_internal_norm (float): Plasma normalised internal inductance - Returns: - float: The calculated bootstrap fraction. + # ======================================================================== - Notes: - The profile assumed for the alphan and alpat indexes is only a parabolic profile without a pedestal (L-mode). - The Root Mean Squared Error for the fitting database of this formula was 0.025 - Concentrating on the positive shear plasmas using the ACCOME code equilibria with the fully non-inductively driven - conditions with neutral beam (NB) injection only are calculated. - The electron temperature and the ion temperature were assumed to be equal - This can be used for all apsect ratios. - The diamagnetic fraction is included in this formula. + # Lackner-Gottardi scaling (L-mode) + elif i_confinement_time == 17: + t_electron_confinement = confinement.lackner_gottardi_confinement_time( + pcur, + rmajor, + rminor, + kappa95, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + ) - References: - 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. - """ - # Sakai states that the ACCOME dataset used has the toridal diamagnetic current included in the bootstrap current - # So the diamganetic current should not be calculated with this. i_diamagnetic_current = 0 - return ( - 10 ** (0.951 * eps - 0.948) - * beta_poloidal ** (1.226 * eps + 1.584) - * ind_plasma_internal_norm ** (-0.184 * eps - 0.282) - * (q95 / q0) ** (-0.042 * eps - 0.02) - * alphan ** (0.13 * eps + 0.05) - * alphat ** (0.502 * eps - 0.273) - ) + # ======================================================================== - @staticmethod - def bootstrap_fraction_aries( - beta_poloidal: float, - ind_plasma_internal_norm: float, - core_density: float, - average_density: float, - inverse_aspect: float, - ) -> float: - """ - Calculate the bootstrap fraction using the ARIES formula. + # Neo-Kaye scaling (L-mode) + elif i_confinement_time == 18: + t_electron_confinement = confinement.neo_kaye_confinement_time( + pcur, + rmajor, + rminor, + kappa95, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + ) - Parameters: - beta_poloidal (float): Plasma poloidal beta. - ind_plasma_internal_norm (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. + # Riedel scaling (H-mode) + elif i_confinement_time == 19: + t_electron_confinement = confinement.riedel_h_confinement_time( + pcur, + rmajor, + rminor, + kappa95, + dnla20, + b_plasma_toroidal_on_axis, + m_fuel_amu, + p_plasma_loss_mw, + ) - 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. + # Amended version of ITER H90-P law + elif i_confinement_time == 20: + t_electron_confinement = confinement.iter_h90_p_amended_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + m_fuel_amu, + rmajor, + p_plasma_loss_mw, + kappa, + ) - """ - # Using the standard variable naming from the ARIES paper - a_1 = ( - 1.10 - 1.165 * ind_plasma_internal_norm + 0.47 * ind_plasma_internal_norm**2 - ) - b_1 = ( - 0.806 - - 0.885 * ind_plasma_internal_norm - + 0.297 * ind_plasma_internal_norm**2 - ) + # ========================================================================== - c_bs = a_1 + b_1 * (core_density / average_density) + # Sudo et al. scaling (stellarators/heliotron) + elif i_confinement_time == 21: + t_electron_confinement = confinement.sudo_et_al_confinement_time( + rmajor, + rminor, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + ) - return c_bs * np.sqrt(inverse_aspect) * beta_poloidal + # ========================================================================== - @staticmethod - def bootstrap_fraction_andrade( - beta_poloidal: float, - core_pressure: float, - average_pressure: float, - inverse_aspect: float, - ) -> float: - """ - Calculate the bootstrap fraction using the Andrade et al formula. + # Gyro-reduced Bohm scaling + elif i_confinement_time == 22: + t_electron_confinement = confinement.gyro_reduced_bohm_confinement_time( + b_plasma_toroidal_on_axis, + dnla20, + p_plasma_loss_mw, + rminor, + rmajor, + ) - Parameters: - beta_poloidal (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. + # Lackner-Gottardi stellarator scaling + elif i_confinement_time == 23: + t_electron_confinement = ( + confinement.lackner_gottardi_stellarator_confinement_time( + rmajor, + rminor, + dnla20, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + q95, + ) + ) - 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. + # ITER_93 ELM-free H-mode scaling + elif i_confinement_time == 24: + t_electron_confinement = confinement.iter_93h_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + m_fuel_amu, + rmajor, + dnla20, + aspect, + kappa, + ) - """ + # ========================================================================== + # Scaling removed + elif i_confinement_time == 25: + raise ProcessValueError("Scaling removed") + # ========================================================================== - # Using the standard variable naming from the Andrade et.al. paper - c_p = core_pressure / average_pressure + # ELM-free: ITERH-97P + elif i_confinement_time == 26: + t_electron_confinement = confinement.iter_h97p_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + dnla19, + rmajor, + aspect, + kappa, + m_fuel_amu, + ) - # Error +- 0.0007 - c_bs = 0.2340 + # ========================================================================== - return c_bs * np.sqrt(inverse_aspect) * beta_poloidal * c_p**0.8 + # ELMy: ITERH-97P(y) + elif i_confinement_time == 27: + t_electron_confinement = confinement.iter_h97p_elmy_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + dnla19, + rmajor, + aspect, + kappa, + m_fuel_amu, + ) - @staticmethod - def bootstrap_fraction_hoang( - beta_poloidal: float, - pressure_index: float, - current_index: float, - inverse_aspect: float, - ) -> float: - """ - Calculate the bootstrap fraction using the Hoang et al formula. + # ========================================================================== - Parameters: - beta_poloidal (float): Plasma poloidal beta. - pressure_index (float): Pressure profile index. - current_index (float): Current profile index. - inverse_aspect (float): Inverse aspect ratio. + # ITER-96P (= ITER-97L) L-mode scaling + elif i_confinement_time == 28: + t_electron_confinement = confinement.iter_96p_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + kappa95, + rmajor, + aspect, + dnla19, + m_fuel_amu, + p_plasma_loss_mw, + ) - 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 + # Valovic modified ELMy-H mode scaling + # WARNING: No reference found for this scaling. This may not be its real name + elif i_confinement_time == 29: + t_electron_confinement = confinement.valovic_elmy_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + m_fuel_amu, + rmajor, + rminor, + kappa, + p_plasma_loss_mw, + ) - 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$ + # Kaye PPPL Workshop April 1998 L-mode scaling + # WARNING: No reference found for this scaling. This may not be its real name + elif i_confinement_time == 30: + t_electron_confinement = confinement.kaye_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + kappa, + rmajor, + aspect, + dnla19, + m_fuel_amu, + p_plasma_loss_mw, + ) - # 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 + # ITERH-PB98P(y), ELMy H-mode scaling + # WARNING: No reference found for this scaling. This may not be its real name + elif i_confinement_time == 31: + t_electron_confinement = confinement.iter_pb98py_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) - c_bs = np.sqrt((pressure_index + 1) / (current_index + 1)) + # ========================================================================== - return 0.4 * np.sqrt(inverse_aspect) * beta_poloidal**0.9 * c_bs + # IPB98(y), ELMy H-mode scaling + elif i_confinement_time == 32: + t_electron_confinement = confinement.iter_ipb98y_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + kappa, + aspect, + m_fuel_amu, + ) - @staticmethod - def bootstrap_fraction_wong( - beta_poloidal: float, - density_index: float, - temperature_index: float, - inverse_aspect: float, - elongation: float, - ) -> float: - """ - Calculate the bootstrap fraction using the Wong et al formula. + # ========================================================================== - Parameters: - beta_poloidal (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. + # IPB98(y,1), ELMy H-mode scaling + elif i_confinement_time == 33: + t_electron_confinement = confinement.iter_ipb98y1_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) - 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 + # IPB98(y,2), ELMy H-mode scaling + elif i_confinement_time == 34: + t_electron_confinement = confinement.iter_ipb98y2_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) - 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) + # IPB98(y,3), ELMy H-mode scaling + elif i_confinement_time == 35: + t_electron_confinement = confinement.iter_ipb98y3_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) - c_bs = 0.773 + 0.019 * elongation + # ========================================================================== - return c_bs * f_peak**0.25 * beta_poloidal * np.sqrt(inverse_aspect) + # IPB98(y,4), ELMy H-mode scaling + elif i_confinement_time == 36: + t_electron_confinement = confinement.iter_ipb98y4_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) - @staticmethod - def bootstrap_fraction_gi_I( # noqa: N802 - beta_poloidal: 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: - beta_poloidal (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. + # ISS95 stellarator scaling + elif i_confinement_time == 37: + # dummy argument q95 is actual argument iotabar for stellarators + iotabar = q95 + t_electron_confinement = confinement.iss95_stellarator_confinement_time( + rminor, + rmajor, + dnla19, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + iotabar, + ) - 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. + # ISS04 stellarator scaling + elif i_confinement_time == 38: + # dummy argument q95 is actual argument iotabar for stellarators + iotabar = q95 + t_electron_confinement = confinement.iss04_stellarator_confinement_time( + rminor, + rmajor, + dnla19, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + iotabar, + ) - 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 + # DS03 beta-independent H-mode scaling + elif i_confinement_time == 39: + t_electron_confinement = confinement.ds03_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + kappa95, + aspect, + m_fuel_amu, + ) - 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) * beta_poloidal + # Murari "Non-power law" scaling + elif i_confinement_time == 40: + t_electron_confinement = confinement.murari_confinement_time( + pcur, + rmajor, + physics_variables.kappa_ipb, + dnla19, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + ) - @staticmethod - def bootstrap_fraction_gi_II( # noqa: N802 - beta_poloidal: 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: - beta_poloidal (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. + # Petty08, beta independent dimensionless scaling + elif i_confinement_time == 41: + t_electron_confinement = confinement.petty08_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + ) - 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. + # Lang high density relevant confinement scaling + elif i_confinement_time == 42: + t_electron_confinement = confinement.lang_high_density_confinement_time( + plasma_current, + b_plasma_toroidal_on_axis, + nd_plasma_electron_line, + p_plasma_loss_mw, + rmajor, + rminor, + q95, + qstar, + aspect, + m_fuel_amu, + physics_variables.kappa_ipb, + ) - 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 + # Hubbard 2017 I-mode confinement time scaling - nominal + elif i_confinement_time == 43: + t_electron_confinement = confinement.hubbard_nominal_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla20, + p_plasma_loss_mw, + ) - 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) * beta_poloidal + # Hubbard 2017 I-mode confinement time scaling - lower + elif i_confinement_time == 44: + t_electron_confinement = confinement.hubbard_lower_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla20, + p_plasma_loss_mw, + ) - @staticmethod - def bootstrap_fraction_sugiyama_l_mode( - eps: float, - beta_poloidal: float, - alphan: float, - alphat: float, - zeff: float, - q95: float, - q0: float, - ) -> float: - """ - Calculate the bootstrap fraction using the L-mode scaling from the Sugiyama et al formula. + # ========================================================================== - :param eps: Inverse aspect ratio. - :type eps: float - :param beta_poloidal: Plasma poloidal beta. - :type beta_poloidal: float - :param alphan: Density profile index. - :type alphan: float - :param alphat: Temperature profile index. - :type alphat: float - :param zeff: Plasma effective charge. - :type zeff: float - :param q95: Safety factor at 95% of the plasma radius. - :type q95: float - :param q0: Safety factor at the magnetic axis. - :type q0: float + # Hubbard 2017 I-mode confinement time scaling - upper + elif i_confinement_time == 45: + t_electron_confinement = confinement.hubbard_upper_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla20, + p_plasma_loss_mw, + ) - :returns: The calculated bootstrap fraction. - :rtype: float + # ========================================================================== - :notes: - - This scaling is derived for L-mode plasmas. - - Ion and electron temperature are the same - - Z_eff has a uniform profile, with only fully stripped carbon impurity + # Menard NSTX, ELMy H-mode scaling + elif i_confinement_time == 46: + t_electron_confinement = confinement.menard_nstx_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) - :references: - - S. Sugiyama, T. Goto, H. Utoh, and Y. Sakamoto, “Improvement of core plasma power and - current balance models for tokamak systems code considering H-mode plasma profiles,” - Fusion Engineering and Design, vol. 216, p. 115022, Jul. 2025, doi: - https://doi.org/10.1016/j.fusengdes.2025.115022. - """ + # ========================================================================== - return ( - 0.740 - * eps**0.418 - * beta_poloidal**0.904 - * alphan**0.06 - * alphat**-0.138 - * zeff**0.230 - * (q95 / q0) ** -0.142 - ) + # Menard NSTX-Petty08 Hybrid + elif i_confinement_time == 47: + t_electron_confinement = ( + confinement.menard_nstx_petty08_hybrid_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.kappa_ipb, + aspect, + m_fuel_amu, + ) + ) - @staticmethod - def bootstrap_fraction_sugiyama_h_mode( - eps: float, - beta_poloidal: float, - alphan: float, - alphat: float, - tbeta: float, - zeff: float, - q95: float, - q0: float, - radius_plasma_pedestal_density_norm: float, - nd_plasma_pedestal_electron: float, - n_greenwald: float, - temp_plasma_pedestal_kev: float, - ) -> float: - """ - Calculate the bootstrap fraction using the H-mode scaling from the Sugiyama et al formula. + # ========================================================================== - :param eps: Inverse aspect ratio. - :type eps: float - :param beta_poloidal: Plasma poloidal beta. - :type beta_poloidal: float - :param alphan: Density profile index. - :type alphan: float - :param alphat: Temperature profile index. - :type alphat: float - :param tbeta: Second temperature profile index. - :type tbeta: float - :param zeff: Plasma effective charge. - :type zeff: float - :param q95: Safety factor at 95% of the plasma radius. - :type q95: float - :param q0: Safety factor at the magnetic axis. - :type q0: float - :param radius_plasma_pedestal_density_norm: Normalised plasma radius of density pedestal. - :type radius_plasma_pedestal_density_norm: float - :param nd_plasma_pedestal_electron: Electron number density at the pedestal [m^-3]. - :type nd_plasma_pedestal_electron: float - :param n_greenwald: Greenwald density limit [m^-3]. - :type n_greenwald: float - :param temp_plasma_pedestal_kev: Electron temperature at the pedestal [keV]. - :type temp_plasma_pedestal_kev: float + # NSTX gyro-Bohm (Buxton) + elif i_confinement_time == 48: + t_electron_confinement = confinement.nstx_gyro_bohm_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + rmajor, + dnla20, + ) - :returns: The calculated bootstrap fraction. - :rtype: float + # ========================================================================== - :notes: - - This scaling is derived for H-mode plasmas. - - The temperature and density pedestal positions are the same - - Separatrix temperature and density are zero - - Ion and electron temperature are the same - - Z_eff has a uniform profile, with only fully stripped carbon impurity + # ITPA20 H-mode scaling + elif i_confinement_time == 49: + t_electron_confinement = confinement.itpa20_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + dnla19, + p_plasma_loss_mw, + rmajor, + physics_variables.triang, + physics_variables.kappa_ipb, + eps, + physics_variables.m_ions_total_amu, + ) - :references: - - S. Sugiyama, T. Goto, H. Utoh, and Y. Sakamoto, “Improvement of core plasma power and - current balance models for tokamak systems code considering H-mode plasma profiles,” - Fusion Engineering and Design, vol. 216, p. 115022, Jul. 2025, doi: - https://doi.org/10.1016/j.fusengdes.2025.115022. - """ + # ========================================================================== - return ( - 0.789 - * eps**0.606 - * beta_poloidal**0.960 - * alphan**0.0319 - * alphat**0.00822 - * tbeta**-0.0783 - * zeff**0.241 - * (q95 / q0) ** -0.103 - * radius_plasma_pedestal_density_norm**0.367 - * (nd_plasma_pedestal_electron / n_greenwald) ** -0.174 - * temp_plasma_pedestal_kev**0.0552 - ) + # ITPA20-IL confinement time scaling + elif i_confinement_time == 50: + t_electron_confinement = confinement.itpa20_il_confinement_time( + pcur, + b_plasma_toroidal_on_axis, + p_plasma_loss_mw, + dnla19, + physics_variables.m_ions_total_amu, + rmajor, + physics_variables.triang, + physics_variables.kappa_ipb, + ) - def find_other_h_factors(self, i_confinement_time: int) -> float: - """ - Function to find H-factor for the equivalent confinement time in other scalings. + # ========================================================================== - Args: - i_confinement_time (int): Index of the confinement time scaling to use. + else: + raise ProcessValueError( + "Illegal value for i_confinement_time", + i_confinement_time=i_confinement_time, + ) - Returns: - float: The calculated H-factor. - """ + # Apply H-factor correction to chosen scaling + t_electron_energy_confinement = hfact * t_electron_confinement - def fhz(hfact: float) -> float: - """ - Function used to find power balance. + # Ion energy confinement time + t_ion_energy_confinement = t_electron_energy_confinement - Args: - hfact (float): H-factor to be used in the calculation. + # Calculation of the transport power loss terms + # Transport losses in Watts/m3 are 3/2 * n.e.T / tau , with T in eV + # (here, temp_plasma_ion_density_weighted_kev and temp_plasma_electron_density_weighted_kev are in keV, and pden_electron_transport_loss_mw and pden_ion_transport_loss_mw are in MW/m3) - Returns: - float: The difference between the calculated power and the required power for balance. - """ - ( - ptrez, - ptriz, - _, - _, - _, - _, - ) = self.calculate_confinement_time( - physics_variables.m_fuel_amu, - physics_variables.p_alpha_total_mw, - physics_variables.aspect, - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.nd_plasma_ions_total_vol_avg, - physics_variables.nd_plasma_electrons_vol_avg, - physics_variables.nd_plasma_electron_line, - physics_variables.eps, - hfact, - i_confinement_time, - physics_variables.i_plasma_ignited, - physics_variables.kappa, - physics_variables.kappa95, - physics_variables.p_non_alpha_charged_mw, - current_drive_variables.p_hcd_injected_total_mw, - physics_variables.plasma_current, - physics_variables.pden_plasma_core_rad_mw, - physics_variables.rmajor, - physics_variables.rminor, - physics_variables.temp_plasma_electron_density_weighted_kev, - physics_variables.temp_plasma_ion_density_weighted_kev, - physics_variables.q95, - physics_variables.qstar, - physics_variables.vol_plasma, - physics_variables.n_charge_plasma_effective_vol_avg, - ) + # The transport losses is just the electron and ion thermal energies divided by the confinement time. + pden_ion_transport_loss_mw = ( + (3 / 2) + * (constants.ELECTRON_CHARGE / 1e3) + * nd_plasma_ions_total_vol_avg + * temp_plasma_ion_density_weighted_kev + / t_ion_energy_confinement + ) + pden_electron_transport_loss_mw = ( + (3 / 2) + * (constants.ELECTRON_CHARGE / 1e3) + * nd_plasma_electrons_vol_avg + * temp_plasma_electron_density_weighted_kev + / t_electron_energy_confinement + ) - # At power balance, fhz is zero. - fhz_value = ( - ptrez - + ptriz - - physics_variables.f_p_alpha_plasma_deposited - * physics_variables.pden_alpha_total_mw - - physics_variables.pden_non_alpha_charged_mw - - physics_variables.pden_plasma_ohmic_mw - ) + ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * ( + temp_plasma_ion_density_weighted_kev + / temp_plasma_electron_density_weighted_kev + ) - # Take into account whether injected power is included in tau_e calculation (i.e. whether device is ignited) - if physics_variables.i_plasma_ignited == 0: - fhz_value -= ( - current_drive_variables.p_hcd_injected_total_mw - / physics_variables.vol_plasma - ) + # Global energy confinement time - # Include the radiation power if requested - if physics_variables.i_rad_loss == 0: - fhz_value += physics_variables.pden_plasma_rad_mw - elif physics_variables.i_rad_loss == 1: - fhz_value += physics_variables.pden_plasma_core_rad_mw + t_energy_confinement = (ratio + 1.0e0) / ( + ratio / t_ion_energy_confinement + 1.0e0 / t_electron_energy_confinement + ) - return fhz_value + # For comparison directly calculate the confinement time from the stored energy calculated + # from the total plasma beta and the loss power used above. + physics_variables.t_energy_confinement_beta = ( + physics_variables.e_plasma_beta / 1e6 + ) / p_plasma_loss_mw - return root_scalar(fhz, bracket=(0.01, 150), xtol=0.001).root + return ( + pden_electron_transport_loss_mw, + pden_ion_transport_loss_mw, + t_electron_energy_confinement, + t_ion_energy_confinement, + t_energy_confinement, + p_plasma_loss_mw, + ) @staticmethod - def calculate_confinement_time( + def calculate_plasma_masses( m_fuel_amu: float, - p_alpha_total_mw: float, - aspect: float, - b_plasma_toroidal_on_axis: float, + m_ions_total_amu: float, nd_plasma_ions_total_vol_avg: float, - nd_plasma_electrons_vol_avg: float, - nd_plasma_electron_line: float, - eps: float, - hfact: float, - i_confinement_time: int, - i_plasma_ignited: int, - kappa: float, - kappa95: float, - p_non_alpha_charged_mw: float, - p_hcd_injected_total_mw: float, - plasma_current: float, - pden_plasma_core_rad_mw: float, - rmajor: float, - rminor: float, - temp_plasma_electron_density_weighted_kev: float, - temp_plasma_ion_density_weighted_kev: float, - q95: float, - qstar: float, + nd_plasma_fuel_ions_vol_avg: float, + nd_plasma_alphas_vol_avg: float, vol_plasma: float, - zeff: float, - ) -> tuple[float, float, float, float, float, float, float]: + nd_plasma_electrons_vol_avg: float, + ) -> tuple[float, float, float, float, float]: """ - Calculate the confinement times and the transport power loss terms. + Calculate the plasma masses. - :param m_fuel_amu: Average mass of fuel (amu) - :param p_alpha_total_mw: Alpha particle power (MW) - :param aspect: Aspect ratio - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T) - :param nd_plasma_ions_total_vol_avg: Total ion density (/m3) - :param nd_plasma_electrons_vol_avg: Volume averaged electron density (/m3) - :param nd_plasma_electron_line: Line-averaged electron density (/m3) - :param eps: Inverse aspect ratio - :param hfact: H factor on energy confinement scalings - :param i_confinement_time: Switch for energy confinement scaling to use - :param i_plasma_ignited: Switch for ignited calculation - :param kappa: Plasma elongation - :param kappa95: Plasma elongation at 95% surface - :param p_non_alpha_charged_mw: Non-alpha charged particle fusion power (MW) - :param p_hcd_injected_total_mw: Auxiliary power to ions and electrons (MW) - :param plasma_current: Plasma current (A) - :param pden_plasma_core_rad_mw: Total core radiation power (MW/m3) - :param q95: Edge safety factor (tokamaks), or rotational transform iotabar (stellarators) - :param qstar: Equivalent cylindrical edge safety factor - :param rmajor: Plasma major radius (m) - :param rminor: Plasma minor radius (m) - :param temp_plasma_electron_density_weighted_kev: Density weighted average electron temperature (keV) - :param temp_plasma_ion_density_weighted_kev: Density weighted average ion temperature (keV) - :param vol_plasma: Plasma volume (m3) - :param a_plasma_poloidal: Plasma cross-sectional area (m2) - :param zeff: Plasma effective charge + :param m_fuel_amu: Average mass of fuel (amu). + :type m_fuel_amu: float + :param m_ions_total_amu: Average mass of all ions (amu). + :type m_ions_total_amu: float + :param nd_plasma_ions_total_vol_avg: Total ion density (/m3). + :type nd_plasma_ions_total_vol_avg: float + :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (/m3). + :type nd_plasma_fuel_ions_vol_avg: float + :param nd_plasma_alphas_vol_avg: Alpha ash density (/m3). + :type nd_plasma_alphas_vol_avg: float + :param vol_plasma: Plasma volume (m3). + :type vol_plasma: float + :param nd_plasma_electrons_vol_avg: Volume averaged electron density (/m3). + :type nd_plasma_electrons_vol_avg: float - :return: Tuple containing: - - pden_electron_transport_loss_mw (float): Electron transport power (MW/m3) - - pden_ion_transport_loss_mw (float): Ion transport power (MW/m3) - - t_electron_energy_confinement (float): Electron energy confinement time (s) - - t_ion_energy_confinement (float): Ion energy confinement time (s) - - t_energy_confinement (float): Global energy confinement time (s) - - p_plasma_loss_mw (float): Heating power (MW) assumed in calculation + :returns: A tuple containing: + :rtype: tuple[float, float, float, float, float] """ - # ======================================================================== + # Calculate mass of fuel ions + m_plasma_fuel_ions = (m_fuel_amu * constants.ATOMIC_MASS_UNIT) * ( + nd_plasma_fuel_ions_vol_avg * vol_plasma + ) - # Calculate heating power (MW) - p_plasma_loss_mw = ( - physics_variables.f_p_alpha_plasma_deposited * p_alpha_total_mw - + p_non_alpha_charged_mw - + physics_variables.p_plasma_ohmic_mw + m_plasma_ions_total = (m_ions_total_amu * constants.ATOMIC_MASS_UNIT) * ( + nd_plasma_ions_total_vol_avg * vol_plasma ) - # If the device is not ignited, add the injected auxiliary power - if i_plasma_ignited == 0: - p_plasma_loss_mw = p_plasma_loss_mw + p_hcd_injected_total_mw + m_plasma_alpha = (nd_plasma_alphas_vol_avg * vol_plasma) * constants.ALPHA_MASS - # Include the radiation as a loss term if requested - if physics_variables.i_rad_loss == 0: - p_plasma_loss_mw = ( - p_plasma_loss_mw - physics_variables.pden_plasma_rad_mw * vol_plasma - ) - elif physics_variables.i_rad_loss == 1: - p_plasma_loss_mw = ( - p_plasma_loss_mw - pden_plasma_core_rad_mw * vol_plasma - ) # shouldn't this be vol_core instead of vol_plasma? - # else do not adjust p_plasma_loss_mw for radiation + m_plasma_electron = constants.ELECTRON_MASS * ( + nd_plasma_electrons_vol_avg * vol_plasma + ) - # Ensure heating power is positive (shouldn't be necessary) - p_plasma_loss_mw = max(p_plasma_loss_mw, 1.0e-3) + m_plasma = m_plasma_electron + m_plasma_ions_total - # ======================================================================== + return ( + m_plasma_fuel_ions, + m_plasma_ions_total, + m_plasma_alpha, + m_plasma_electron, + m_plasma, + ) - # Line averaged electron density in scaled units - dnla20 = nd_plasma_electron_line * 1.0e-20 - dnla19 = nd_plasma_electron_line * 1.0e-19 - # Volume averaged electron density in units of 10**20 m**-3 - n20 = nd_plasma_electrons_vol_avg / 1.0e20 +def res_diff_time(rmajor, res_plasma, kappa95): + """Calculates resistive diffusion time - # Plasma current in MA - pcur = plasma_current / 1.0e6 + Author: James Morris (UKAEA) - # Separatrix kappa defined with plasma volume for IPB scalings - # Updated version of kappa used by the IPB98 scalings correction in: + :param rmajor: plasma major radius (m) + :param res_plasma: plasma resistivity (Ohms) + :param kappa95: plasma elongation at 95% flux surface + """ - # None Otto Kardaun, N. K. Thomsen, and None Alexander Chudnovskiy, - # “Corrections to a sequence of papers in Nuclear Fusion,” Nuclear Fusion, - # vol. 48, no. 9, pp. 099801099801, Aug. 2008, - # doi: https://doi.org/10.1088/0029-5515/48/9/099801. + return 2 * constants.RMU0 * rmajor / (res_plasma * kappa95) - physics_variables.kappa_ipb = (vol_plasma / (2.0 * np.pi * rmajor)) / ( - np.pi * rminor**2 - ) - # Electron energy confinement times +def l_h_threshold_power( + nd_plasma_electron_line: float, + b_plasma_toroidal_on_axis: float, + rmajor: float, + rminor: float, + kappa: float, + a_plasma_surface: float, + m_ions_total_amu: float, + aspect: float, + plasma_current: float, +) -> list[float]: + """ + L-mode to H-mode power threshold calculation. - # ======================================================================== + :param nd_plasma_electron_line: Line-averaged electron density (/m3) + :type nd_plasma_electron_line: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T) + :type b_plasma_toroidal_on_axis: float + :param rmajor: Plasma major radius (m) + :type rmajor: float + :param rminor: Plasma minor radius (m) + :type rminor: float + :param kappa: Plasma elongation + :type kappa: float + :param a_plasma_surface: Plasma surface area (m2) + :type a_plasma_surface: float + :param m_ions_total_amu: Average mass of all ions (amu) + :type m_ions_total_amu: float + :param aspect: Aspect ratio + :type aspect: float + :param plasma_current: Plasma current (A) + :type plasma_current: float - # User defined confinement time - if i_confinement_time == 0: # t_electron_energy_confinement is an input - t_electron_confinement = physics_variables.tauee_in + :returns: Array of power thresholds + :rtype: list[float] + """ - # ======================================================================== + dnla20 = 1e-20 * nd_plasma_electron_line - # Nec-Alcator(NA) OH scaling - if i_confinement_time == 1: - t_electron_confinement = confinement.neo_alcator_confinement_time( - n20, rminor, rmajor, qstar - ) + # ======================================================================== - # ======================================================================== + # ITER-1996 H-mode power threshold database + # Fit to 1996 H-mode power threshold database: nominal - # "Mirnov"-like scaling (H-mode) - elif i_confinement_time == 2: # Mirnov scaling (H-mode) - t_electron_confinement = confinement.mirnov_confinement_time( - rminor, kappa95, pcur - ) + # i_l_h_threshold = 1 + iterdd = transition.calculate_iter1996_nominal( + dnla20, b_plasma_toroidal_on_axis, rmajor + ) - # ======================================================================== + # Fit to 1996 H-mode power threshold database: upper bound + # i_l_h_threshold = 2 + iterdd_ub = transition.calculate_iter1996_upper( + dnla20, b_plasma_toroidal_on_axis, rmajor + ) - # Merezhkin-Mukhovatov (MM) OH/L-mode scaling - elif i_confinement_time == 3: - t_electron_confinement = confinement.merezhkin_muhkovatov_confinement_time( - rmajor, - rminor, - kappa95, - qstar, - dnla20, - m_fuel_amu, - temp_plasma_electron_density_weighted_kev, - ) + # Fit to 1996 H-mode power threshold database: lower bound + # i_l_h_threshold = 3 + iterdd_lb = transition.calculate_iter1996_lower( + dnla20, b_plasma_toroidal_on_axis, rmajor + ) - # ======================================================================== + # ======================================================================== - # Shimomura (S) optimized H-mode scaling - elif i_confinement_time == 4: - t_electron_confinement = confinement.shimomura_confinement_time( - rmajor, rminor, b_plasma_toroidal_on_axis, kappa95, m_fuel_amu - ) + # Snipes 1997 ITER H-mode power threshold - # ======================================================================== + # i_l_h_threshold = 4 + snipes_1997 = transition.calculate_snipes1997_iter( + dnla20, b_plasma_toroidal_on_axis, rmajor + ) - # Kaye-Goldston scaling (L-mode) - elif i_confinement_time == 5: - t_electron_confinement = confinement.kaye_goldston_confinement_time( - pcur, - rmajor, - rminor, - kappa, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ) + # i_l_h_threshold = 5 + snipes_1997_kappa = transition.calculate_snipes1997_kappa( + dnla20, b_plasma_toroidal_on_axis, rmajor, kappa + ) - # ======================================================================== + # ======================================================================== - # ITER Power scaling - ITER 89-P (L-mode) - elif i_confinement_time == 6: - t_electron_confinement = confinement.iter_89p_confinement_time( - pcur, - rmajor, - rminor, - kappa, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ) + # Martin et al (2008) for recent ITER scaling, with mass correction + # and 95% confidence limits - # ======================================================================== + # i_l_h_threshold = 6 + martin_nominal = transition.calculate_martin08_nominal( + dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu + ) - # ITER Offset linear scaling - ITER 89-O (L-mode) - elif i_confinement_time == 7: - t_electron_confinement = confinement.iter_89_0_confinement_time( - pcur, - rmajor, - rminor, - kappa, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ) - # ======================================================================== - - # Rebut-Lallia offset linear scaling (L-mode) - elif i_confinement_time == 8: - t_electron_confinement = confinement.rebut_lallia_confinement_time( - rminor, - rmajor, - kappa, - m_fuel_amu, - pcur, - zeff, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - ) - - # ======================================================================== - - # Goldston scaling (L-mode) - elif i_confinement_time == 9: # Goldston scaling (L-mode) - t_electron_confinement = confinement.goldston_confinement_time( - pcur, rmajor, rminor, kappa95, m_fuel_amu, p_plasma_loss_mw - ) + # i_l_h_threshold = 7 + martin_ub = transition.calculate_martin08_upper( + dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu + ) - # ======================================================================== + # i_l_h_threshold = 8 + martin_lb = transition.calculate_martin08_lower( + dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu + ) - # T-10 scaling (L-mode) - elif i_confinement_time == 10: - t_electron_confinement = confinement.t10_confinement_time( - dnla20, - rmajor, - qstar, - b_plasma_toroidal_on_axis, - rminor, - kappa95, - p_plasma_loss_mw, - zeff, - pcur, - ) + # ======================================================================== - # ======================================================================== + # Snipes et al (2000) scaling with mass correction + # Nominal, upper and lower - # JAERI / Odajima-Shimomura L-mode scaling - elif i_confinement_time == 11: # JAERI scaling - t_electron_confinement = confinement.jaeri_confinement_time( - kappa95, - rminor, - m_fuel_amu, - n20, - pcur, - b_plasma_toroidal_on_axis, - rmajor, - qstar, - zeff, - p_plasma_loss_mw, - ) + # i_l_h_threshold = 9 + snipes_2000 = transition.calculate_snipes2000_nominal( + dnla20, b_plasma_toroidal_on_axis, rmajor, rminor, m_ions_total_amu + ) - # ======================================================================== + # i_l_h_threshold = 10 + snipes_2000_ub = transition.calculate_snipes2000_upper( + dnla20, b_plasma_toroidal_on_axis, rmajor, rminor, m_ions_total_amu + ) - # Kaye "big" L-mode scaling (based only on big tokamak data) - elif i_confinement_time == 12: - t_electron_confinement = confinement.kaye_big_confinement_time( - rmajor, - rminor, - b_plasma_toroidal_on_axis, - kappa95, - pcur, - n20, - m_fuel_amu, - p_plasma_loss_mw, - ) + # i_l_h_threshold = 11 + snipes_2000_lb = transition.calculate_snipes2000_lower( + dnla20, b_plasma_toroidal_on_axis, rmajor, rminor, m_ions_total_amu + ) - # ======================================================================== + # ======================================================================== - # ITER H90-P H-mode scaling - elif i_confinement_time == 13: - t_electron_confinement = confinement.iter_h90_p_confinement_time( - pcur, - rmajor, - rminor, - kappa, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ) + # Snipes et al (2000) scaling (closed divertor) with mass correction + # Nominal, upper and lower - # ======================================================================== + # i_l_h_threshold = 12 + snipes_2000_cd = transition.calculate_snipes2000_closed_divertor_nominal( + dnla20, b_plasma_toroidal_on_axis, rmajor, m_ions_total_amu + ) - # Minimum of ITER 89-P and ITER 89-O - elif i_confinement_time == 14: - t_electron_confinement = min( - confinement.iter_89p_confinement_time( - pcur, - rmajor, - rminor, - kappa, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ), - confinement.iter_89_0_confinement_time( - pcur, - rmajor, - rminor, - kappa, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ), - ) + # i_l_h_threshold = 13 + snipes_2000_cd_ub = transition.calculate_snipes2000_closed_divertor_upper( + dnla20, b_plasma_toroidal_on_axis, rmajor, m_ions_total_amu + ) - # ======================================================================== + # i_l_h_threshold = 14 + snipes_2000_cd_lb = transition.calculate_snipes2000_closed_divertor_lower( + dnla20, b_plasma_toroidal_on_axis, rmajor, m_ions_total_amu + ) - # Riedel scaling (L-mode) - elif i_confinement_time == 15: - t_electron_confinement = confinement.riedel_l_confinement_time( - pcur, - rmajor, - rminor, - kappa95, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - ) + # ======================================================================== - # ======================================================================== + # Hubbard et al. 2012 L-I threshold scaling - # Christiansen et al scaling (L-mode) - elif i_confinement_time == 16: - t_electron_confinement = confinement.christiansen_confinement_time( - pcur, - rmajor, - rminor, - kappa95, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - m_fuel_amu, - ) + # i_l_h_threshold = 15 + hubbard_2012 = transition.calculate_hubbard2012_nominal(plasma_current, dnla20) - # ======================================================================== + # i_l_h_threshold = 16 + hubbard_2012_lb = transition.calculate_hubbard2012_lower(plasma_current, dnla20) - # Lackner-Gottardi scaling (L-mode) - elif i_confinement_time == 17: - t_electron_confinement = confinement.lackner_gottardi_confinement_time( - pcur, - rmajor, - rminor, - kappa95, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - ) + # i_l_h_threshold = 17 + hubbard_2012_ub = transition.calculate_hubbard2012_upper(plasma_current, dnla20) - # ======================================================================== + # ======================================================================== - # Neo-Kaye scaling (L-mode) - elif i_confinement_time == 18: - t_electron_confinement = confinement.neo_kaye_confinement_time( - pcur, - rmajor, - rminor, - kappa95, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - ) + # Hubbard et al. 2017 L-I threshold scaling - # ======== ================================================================ + # i_l_h_threshold = 18 + hubbard_2017 = transition.calculate_hubbard2017( + dnla20, a_plasma_surface, b_plasma_toroidal_on_axis + ) - # Riedel scaling (H-mode) - elif i_confinement_time == 19: - t_electron_confinement = confinement.riedel_h_confinement_time( - pcur, - rmajor, - rminor, - kappa95, - dnla20, - b_plasma_toroidal_on_axis, - m_fuel_amu, - p_plasma_loss_mw, - ) + # ======================================================================== - # ======================================================================== + # Aspect ratio corrected Martin et al (2008) - # Amended version of ITER H90-P law - elif i_confinement_time == 20: - t_electron_confinement = confinement.iter_h90_p_amended_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - m_fuel_amu, - rmajor, - p_plasma_loss_mw, - kappa, - ) + # i_l_h_threshold = 19 + martin_nominal_aspect = transition.calculate_martin08_aspect_nominal( + dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu, aspect + ) - # ========================================================================== + # i_l_h_threshold = 20 + martin_ub_aspect = transition.calculate_martin08_aspect_upper( + dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu, aspect + ) - # Sudo et al. scaling (stellarators/heliotron) - elif i_confinement_time == 21: - t_electron_confinement = confinement.sudo_et_al_confinement_time( - rmajor, - rminor, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - ) + # i_l_h_threshold = 21 + martin_lb_aspect = transition.calculate_martin08_aspect_lower( + dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu, aspect + ) - # ========================================================================== + # ======================================================================== - # Gyro-reduced Bohm scaling - elif i_confinement_time == 22: - t_electron_confinement = confinement.gyro_reduced_bohm_confinement_time( - b_plasma_toroidal_on_axis, - dnla20, - p_plasma_loss_mw, - rminor, - rmajor, - ) + return [ + iterdd, + iterdd_ub, + iterdd_lb, + snipes_1997, + snipes_1997_kappa, + martin_nominal, + martin_ub, + martin_lb, + snipes_2000, + snipes_2000_ub, + snipes_2000_lb, + snipes_2000_cd, + snipes_2000_cd_ub, + snipes_2000_cd_lb, + hubbard_2012, + hubbard_2012_lb, + hubbard_2012_ub, + hubbard_2017, + martin_nominal_aspect, + martin_ub_aspect, + martin_lb_aspect, + ] - # ========================================================================== - # Lackner-Gottardi stellarator scaling - elif i_confinement_time == 23: - t_electron_confinement = ( - confinement.lackner_gottardi_stellarator_confinement_time( - rmajor, - rminor, - dnla20, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - q95, - ) - ) +def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa, lhat): + """Function for calculating upstream temperature(keV) in Reinke model + author: H Lux, CCFE/UKAEA + b_plasma_toroidal_on_axis : input real : toroidal field on axis (T) + flh : input real : fraction of Psep/P_LH + qstar : input real : safety factor similar to q95 (see #707) + rmajor : input real : major radius (m) + eps : input real : inverse aspect ratio + fgw : input real : ratio of volume averaged density to n_GW + kappa : input real : elongation + lhat : input real : connection length factor + This function calculates the upstream temperature in the + divertor/SoL model used for the Reinke citerion. + Issue #707 + M.L. Reinke 2017 Nucl. Fusion 57 034004 + """ + kappa_0 = 2.0e3 # Stangeby W/m/eV^(7/2) - # ========================================================================== + return ( + ( + b_plasma_toroidal_on_axis**0.72 + * flh**0.29 + * fgw**0.21 + * qstar**0.08 + * rmajor**0.33 + ) + * (eps**0.15 * (1.0 + kappa**2.0) ** 0.34) + * (lhat**0.29 * kappa_0 ** (-0.29) * 0.285) + ) - # ITER_93 ELM-free H-mode scaling - elif i_confinement_time == 24: - t_electron_confinement = confinement.iter_93h_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - m_fuel_amu, - rmajor, - dnla20, - aspect, - kappa, - ) - # ========================================================================== - # Scaling removed - elif i_confinement_time == 25: - raise ProcessValueError("Scaling removed") - # ========================================================================== +class BetaNormMaxModel(IntEnum): + """Beta norm max (β_N_max) model types""" - # ELM-free: ITERH-97P - elif i_confinement_time == 26: - t_electron_confinement = confinement.iter_h97p_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - dnla19, - rmajor, - aspect, - kappa, - m_fuel_amu, - ) + USER_INPUT = 0 + WESSON = 1 + ORIGINAL_SCALING = 2 + MENARD = 3 + THLOREUS = 4 + STAMBAUGH = 5 - # ========================================================================== - # ELMy: ITERH-97P(y) - elif i_confinement_time == 27: - t_electron_confinement = confinement.iter_h97p_elmy_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - dnla19, - rmajor, - aspect, - kappa, - m_fuel_amu, - ) +class PlasmaBeta: + """Class to hold plasma beta calculations for plasma processing.""" - # ========================================================================== + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE - # ITER-96P (= ITER-97L) L-mode scaling - elif i_confinement_time == 28: - t_electron_confinement = confinement.iter_96p_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - kappa95, - rmajor, - aspect, - dnla19, - m_fuel_amu, - p_plasma_loss_mw, - ) + def get_beta_norm_max_value(self, model: BetaNormMaxModel) -> float: + """Get the beta norm max value (β_N_max) for the specified model.""" + model_map = { + BetaNormMaxModel.USER_INPUT: physics_variables.beta_norm_max, + BetaNormMaxModel.WESSON: physics_variables.beta_norm_max_wesson, + BetaNormMaxModel.ORIGINAL_SCALING: physics_variables.beta_norm_max_original_scaling, + BetaNormMaxModel.MENARD: physics_variables.beta_norm_max_menard, + BetaNormMaxModel.THLOREUS: physics_variables.beta_norm_max_thloreus, + BetaNormMaxModel.STAMBAUGH: physics_variables.beta_norm_max_stambaugh, + } + return model_map[model] - # ========================================================================== + def run(self) -> None: + """ + Calculate plasma beta values. + """ - # Valovic modified ELMy-H mode scaling - # WARNING: No reference found for this scaling. This may not be its real name - elif i_confinement_time == 29: - t_electron_confinement = confinement.valovic_elmy_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - m_fuel_amu, - rmajor, - rminor, - kappa, - p_plasma_loss_mw, - ) + # ----------------------------------------------------- + # Normalised Beta Limit + # ----------------------------------------------------- - # ========================================================================== + # Normalised beta from Troyon beta limit + physics_variables.beta_norm_total = self.calculate_normalised_beta( + beta=physics_variables.beta_total_vol_avg, + rminor=physics_variables.rminor, + c_plasma=physics_variables.plasma_current, + b_field=physics_variables.b_plasma_toroidal_on_axis, + ) - # Kaye PPPL Workshop April 1998 L-mode scaling - # WARNING: No reference found for this scaling. This may not be its real name - elif i_confinement_time == 30: - t_electron_confinement = confinement.kaye_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - kappa, - rmajor, - aspect, - dnla19, - m_fuel_amu, - p_plasma_loss_mw, - ) + # Define beta_norm_max calculations - # ========================================================================== + physics_variables.beta_norm_max_wesson = self.calculate_beta_norm_max_wesson( + ind_plasma_internal_norm=physics_variables.ind_plasma_internal_norm + ) - # ITERH-PB98P(y), ELMy H-mode scaling - # WARNING: No reference found for this scaling. This may not be its real name - elif i_confinement_time == 31: - t_electron_confinement = confinement.iter_pb98py_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, - ) + # Original scaling law + physics_variables.beta_norm_max_original_scaling = ( + self.calculate_beta_norm_max_original(eps=physics_variables.eps) + ) - # ========================================================================== + # J. Menard scaling law + physics_variables.beta_norm_max_menard = self.calculate_beta_norm_max_menard( + eps=physics_variables.eps + ) - # IPB98(y), ELMy H-mode scaling - elif i_confinement_time == 32: - t_electron_confinement = confinement.iter_ipb98y_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - kappa, - aspect, - m_fuel_amu, + # E. Tholerus scaling law + physics_variables.beta_norm_max_thloreus = self.calculate_beta_norm_max_thloreus( + c_beta=physics_variables.c_beta, + pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis, + pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg, + ) + + # R. D. Stambaugh scaling law + physics_variables.beta_norm_max_stambaugh = ( + self.calculate_beta_norm_max_stambaugh( + f_c_plasma_bootstrap=current_drive_variables.f_c_plasma_bootstrap, + kappa=physics_variables.kappa, + aspect=physics_variables.aspect, ) + ) - # ========================================================================== + # Calculate beta_norm_max based on i_beta_norm_max + try: + model = BetaNormMaxModel(int(physics_variables.i_beta_norm_max)) + physics_variables.beta_norm_max = self.get_beta_norm_max_value(model) + except ValueError: + raise ProcessValueError( + "Illegal value of i_beta_norm_max", + i_beta_norm_max=physics_variables.i_beta_norm_max, + ) from None - # IPB98(y,1), ELMy H-mode scaling - elif i_confinement_time == 33: - t_electron_confinement = confinement.iter_ipb98y1_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, - ) - - # ========================================================================== + # calculate_beta_limit() returns the beta_vol_avg_max for beta + physics_variables.beta_vol_avg_max = self.calculate_beta_limit_from_norm( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + beta_norm_max=physics_variables.beta_norm_max, + plasma_current=physics_variables.plasma_current, + rminor=physics_variables.rminor, + ) - # IPB98(y,2), ELMy H-mode scaling - elif i_confinement_time == 34: - t_electron_confinement = confinement.iter_ipb98y2_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, - ) + physics_variables.beta_toroidal_vol_avg = ( + physics_variables.beta_total_vol_avg + * physics_variables.b_plasma_total**2 + / physics_variables.b_plasma_toroidal_on_axis**2 + ) - # ========================================================================== + # Calculate physics_variables.beta poloidal [-] + physics_variables.beta_poloidal_vol_avg = self.calculate_poloidal_beta( + b_plasma_total=physics_variables.b_plasma_total, + b_plasma_poloidal_average=physics_variables.b_plasma_poloidal_average, + beta=physics_variables.beta_total_vol_avg, + ) - # IPB98(y,3), ELMy H-mode scaling - elif i_confinement_time == 35: - t_electron_confinement = confinement.iter_ipb98y3_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, - ) + physics_variables.beta_thermal_vol_avg = ( + physics_variables.beta_total_vol_avg + - physics_variables.beta_fast_alpha + - physics_variables.beta_beam + ) - # ========================================================================== + physics_variables.beta_poloidal_eps = ( + physics_variables.beta_poloidal_vol_avg * physics_variables.eps + ) - # IPB98(y,4), ELMy H-mode scaling - elif i_confinement_time == 36: - t_electron_confinement = confinement.iter_ipb98y4_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, + physics_variables.beta_thermal_poloidal_vol_avg = ( + physics_variables.beta_thermal_vol_avg + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_poloidal_average + ) + ** 2 + ) + physics_variables.beta_thermal_toroidal_vol_avg = ( + physics_variables.beta_thermal_vol_avg + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_toroidal_on_axis ) + ** 2 + ) - # ========================================================================== + # ======================================================= - # ISS95 stellarator scaling - elif i_confinement_time == 37: - # dummy argument q95 is actual argument iotabar for stellarators - iotabar = q95 - t_electron_confinement = confinement.iss95_stellarator_confinement_time( - rminor, - rmajor, - dnla19, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - iotabar, + # Mirror the pressure profiles to match the doubled toroidal field profile + pres_profile_total = np.concatenate([ + physics_variables.pres_plasma_thermal_total_profile[::-1], + physics_variables.pres_plasma_thermal_total_profile, + ]) + + physics_variables.beta_thermal_toroidal_profile = np.array([ + self.calculate_plasma_beta( + pres_plasma=pres_profile_total[i], + b_field=physics_variables.b_plasma_toroidal_profile[i], ) + for i in range(len(physics_variables.b_plasma_toroidal_profile)) + ]) - # ========================================================================== + # ======================================================= - # ISS04 stellarator scaling - elif i_confinement_time == 38: - # dummy argument q95 is actual argument iotabar for stellarators - iotabar = q95 - t_electron_confinement = confinement.iss04_stellarator_confinement_time( - rminor, - rmajor, - dnla19, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - iotabar, + physics_variables.beta_norm_thermal = self.calculate_normalised_beta( + beta=physics_variables.beta_thermal_vol_avg, + rminor=physics_variables.rminor, + c_plasma=physics_variables.plasma_current, + b_field=physics_variables.b_plasma_toroidal_on_axis, + ) + + physics_variables.beta_norm_toroidal = ( + physics_variables.beta_norm_total + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_toroidal_on_axis + ) + ** 2 + ) + physics_variables.beta_norm_poloidal = ( + physics_variables.beta_norm_total + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_poloidal_average ) + ** 2 + ) - # ========================================================================== + physics_variables.f_beta_alpha_beam_thermal = ( + physics_variables.beta_fast_alpha + physics_variables.beta_beam + ) / physics_variables.beta_thermal_vol_avg - # DS03 beta-independent H-mode scaling - elif i_confinement_time == 39: - t_electron_confinement = confinement.ds03_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - kappa95, - aspect, - m_fuel_amu, - ) + # Plasma thermal energy derived from the thermal beta + physics_variables.e_plasma_beta_thermal = self.calculate_plasma_energy_from_beta( + beta=physics_variables.beta_thermal_vol_avg, + b_field=physics_variables.b_plasma_total, + vol_plasma=physics_variables.vol_plasma, + ) - # ========================================================================== + # Plasma thermal energy derived from the total beta + physics_variables.e_plasma_beta = self.calculate_plasma_energy_from_beta( + beta=physics_variables.beta_total_vol_avg, + b_field=physics_variables.b_plasma_total, + vol_plasma=physics_variables.vol_plasma, + ) - # Murari "Non-power law" scaling - elif i_confinement_time == 40: - t_electron_confinement = confinement.murari_confinement_time( - pcur, - rmajor, - physics_variables.kappa_ipb, - dnla19, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - ) + @staticmethod + def calculate_plasma_beta( + pres_plasma: float | np.ndarray, b_field: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the plasma beta (β) for a given pressure and field. - # ========================================================================== + :param pres_plasma: Plasma pressure (in Pascals). + :type pres_plasma: float | np.ndarray + :param b_field: Magnetic field strength (in Tesla). + :type b_field: float | np.ndarray + :returns: The plasma beta (dimensionless). + :rtype: float | np.ndarray - # Petty08, beta independent dimensionless scaling - elif i_confinement_time == 41: - t_electron_confinement = confinement.petty08_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - ) + Plasma beta is the ratio of plasma pressure to magnetic pressure. + """ - # ========================================================================== + return 2 * constants.RMU0 * pres_plasma / (b_field**2) - # Lang high density relevant confinement scaling - elif i_confinement_time == 42: - t_electron_confinement = confinement.lang_high_density_confinement_time( - plasma_current, - b_plasma_toroidal_on_axis, - nd_plasma_electron_line, - p_plasma_loss_mw, - rmajor, - rminor, - q95, - qstar, - aspect, - m_fuel_amu, - physics_variables.kappa_ipb, - ) + @staticmethod + def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: + """ + Calculate the Wesson normalsied beta upper limit. - # ========================================================================== + :param ind_plasma_internal_norm: Plasma normalised internal inductance + :type ind_plasma_internal_norm: float - # Hubbard 2017 I-mode confinement time scaling - nominal - elif i_confinement_time == 43: - t_electron_confinement = confinement.hubbard_nominal_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla20, - p_plasma_loss_mw, - ) + :return: The Wesson normalised beta upper limit. + :rtype: float - # ========================================================================== + :Notes: + - It is recommended to use this method with the other Wesson relations for normalsied internal + inductance and current profile index. + - This fit is derived from the DIII-D database for β_N >= 2.5 - # Hubbard 2017 I-mode confinement time scaling - lower - elif i_confinement_time == 44: - t_electron_confinement = confinement.hubbard_lower_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla20, - p_plasma_loss_mw, - ) + :References: + - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, + International Series of Monographs on Physics, Volume 149. - # ========================================================================== + - T. T. S et al., “Profile Optimization and High Beta Discharges and Stability of High Elongation Plasmas in the DIII-D Tokamak,” + Osti.gov, Oct. 1990. https://www.osti.gov/biblio/6194284 (accessed Dec. 19, 2024). + """ + return 4 * ind_plasma_internal_norm - # Hubbard 2017 I-mode confinement time scaling - upper - elif i_confinement_time == 45: - t_electron_confinement = confinement.hubbard_upper_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla20, - p_plasma_loss_mw, - ) + @staticmethod + def calculate_beta_norm_max_original(eps: float) -> float: + """ + Calculate the original scaling law normalsied beta upper limit. - # ========================================================================== + :param eps: Plasma normalised internal inductance + :type eps: float - # Menard NSTX, ELMy H-mode scaling - elif i_confinement_time == 46: - t_electron_confinement = confinement.menard_nstx_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, - ) + :return: The original scaling law normalised beta upper limit. + :rtype: float - # ========================================================================== + :Notes: - # Menard NSTX-Petty08 Hybrid - elif i_confinement_time == 47: - t_electron_confinement = ( - confinement.menard_nstx_petty08_hybrid_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.kappa_ipb, - aspect, - m_fuel_amu, - ) - ) + :References: - # ========================================================================== + """ + return 2.7 * (1.0 + 5.0 * eps**3.5) - # NSTX gyro-Bohm (Buxton) - elif i_confinement_time == 48: - t_electron_confinement = confinement.nstx_gyro_bohm_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - rmajor, - dnla20, - ) + @staticmethod + def calculate_beta_norm_max_menard(eps: float) -> float: + """ + Calculate the Menard normalsied beta upper limit. - # ========================================================================== + :param eps: Plasma normalised internal inductance + :type eps: float - # ITPA20 H-mode scaling - elif i_confinement_time == 49: - t_electron_confinement = confinement.itpa20_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - dnla19, - p_plasma_loss_mw, - rmajor, - physics_variables.triang, - physics_variables.kappa_ipb, - eps, - physics_variables.m_ions_total_amu, - ) + :return: The Menard normalised beta upper limit. + :rtype: float - # ========================================================================== + :Notes: + - Found as a reasonable fit to the computed no wall limit at f_BS ≈ 50% + - Uses maximum κ data from NSTX at A = 1.45, A = 1.75. Along with record + β_T data from DIII-D at A = 2.9 and high κ. - # ITPA20-IL confinement time scaling - elif i_confinement_time == 50: - t_electron_confinement = confinement.itpa20_il_confinement_time( - pcur, - b_plasma_toroidal_on_axis, - p_plasma_loss_mw, - dnla19, - physics_variables.m_ions_total_amu, - rmajor, - physics_variables.triang, - physics_variables.kappa_ipb, - ) + :References: + - # J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” + Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, + doi: https://doi.org/10.1088/0029-5515/56/10/106023. - # ========================================================================== + """ + return 3.12 + 3.5 * eps**1.7 - else: - raise ProcessValueError( - "Illegal value for i_confinement_time", - i_confinement_time=i_confinement_time, - ) + @staticmethod + def calculate_beta_norm_max_thloreus( + c_beta: float, pres_plasma_on_axis: float, pres_plasma_vol_avg: float + ) -> float: + """ + Calculate the E. Tholerus normalized beta upper limit. - # Apply H-factor correction to chosen scaling - t_electron_energy_confinement = hfact * t_electron_confinement + :param c_beta: Pressure peaking factor coefficient. + :type c_beta: float + :param pres_plasma_on_axis: Central plasma pressure (Pa). + :type pres_plasma_on_axis: float + :param pres_plasma_vol_avg: Volume-averaged plasma pressure (Pa). + :type pres_plasma_vol_avg: float - # Ion energy confinement time - t_ion_energy_confinement = t_electron_energy_confinement + :return: The E. Tholerus normalized beta upper limit. + :rtype: float - # Calculation of the transport power loss terms - # Transport losses in Watts/m3 are 3/2 * n.e.T / tau , with T in eV - # (here, temp_plasma_ion_density_weighted_kev and temp_plasma_electron_density_weighted_kev are in keV, and pden_electron_transport_loss_mw and pden_ion_transport_loss_mw are in MW/m3) + :Notes: + - This method calculates the normalized beta upper limit based on the pressure peaking factor (Fp), + which is defined as the ratio of the peak pressure to the average pressure. + - The formula is derived from operational space studies of flat-top plasma in the STEP power plant. - # The transport losses is just the electron and ion thermal energies divided by the confinement time. - pden_ion_transport_loss_mw = ( - (3 / 2) - * (constants.ELECTRON_CHARGE / 1e3) - * nd_plasma_ions_total_vol_avg - * temp_plasma_ion_density_weighted_kev - / t_ion_energy_confinement - ) - pden_electron_transport_loss_mw = ( - (3 / 2) - * (constants.ELECTRON_CHARGE / 1e3) - * nd_plasma_electrons_vol_avg - * temp_plasma_electron_density_weighted_kev - / t_electron_energy_confinement + :References: + - E. Tholerus et al., “Flat-top plasma operational space of the STEP power plant,” + Nuclear Fusion, Aug. 2024, doi: https://doi.org/10.1088/1741-4326/ad6ea2. + """ + return 3.7 + ( + (c_beta / (pres_plasma_on_axis / pres_plasma_vol_avg)) + * (12.5 - 3.5 * (pres_plasma_on_axis / pres_plasma_vol_avg)) ) - ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * ( - temp_plasma_ion_density_weighted_kev - / temp_plasma_electron_density_weighted_kev - ) + @staticmethod + def calculate_beta_norm_max_stambaugh( + f_c_plasma_bootstrap: float, + kappa: float, + aspect: float, + ) -> float: + """ + Calculate the Stambaugh normalized beta upper limit. - # Global energy confinement time + :param f_c_plasma_bootstrap: Bootstrap current fraction. + :type f_c_plasma_bootstrap: float + :param kappa: Plasma separatrix elongation. + :type kappa: float + :param aspect: Plasma aspect ratio. + :type aspect: float - t_energy_confinement = (ratio + 1.0e0) / ( - ratio / t_ion_energy_confinement + 1.0e0 / t_electron_energy_confinement - ) + :return: The Stambaugh normalized beta upper limit. + :rtype: float - # For comparison directly calculate the confinement time from the stored energy calculated - # from the total plasma beta and the loss power used above. - physics_variables.t_energy_confinement_beta = ( - physics_variables.e_plasma_beta / 1e6 - ) / p_plasma_loss_mw + :Notes: + - This method calculates the normalized beta upper limit based on the Stambaugh scaling. + - The formula is derived from empirical fits to high-performance, steady-state tokamak equilibria. + + :References: + - R. D. Stambaugh et al., “Fusion Nuclear Science Facility Candidates,” + Fusion Science and Technology, vol. 59, no. 2, pp. 279-307, Feb. 2011, + doi: https://doi.org/10.13182/fst59-279. + - Y. R. Lin-Liu and R. D. Stambaugh, “Optimum equilibria for high performance, steady state tokamaks,” + Nuclear Fusion, vol. 44, no. 4, pp. 548-554, Mar. 2004, + doi: https://doi.org/10.1088/0029-5515/44/4/009. + """ return ( - pden_electron_transport_loss_mw, - pden_ion_transport_loss_mw, - t_electron_energy_confinement, - t_ion_energy_confinement, - t_energy_confinement, - p_plasma_loss_mw, + f_c_plasma_bootstrap + * 10 + * (-0.7748 + (1.2869 * kappa) - (0.2921 * kappa**2) + (0.0197 * kappa**3)) + / (aspect**0.5523 * np.tanh((1.8524 + (0.2319 * kappa)) / aspect**0.6163)) ) @staticmethod - def calculate_plasma_masses( - m_fuel_amu: float, - m_ions_total_amu: float, - nd_plasma_ions_total_vol_avg: float, - nd_plasma_fuel_ions_vol_avg: float, - nd_plasma_alphas_vol_avg: float, - vol_plasma: float, - nd_plasma_electrons_vol_avg: float, - ) -> tuple[float, float, float, float, float]: + def calculate_normalised_beta( + beta: float, rminor: float, c_plasma: float, b_field: float + ) -> float: + """Calculate normalised beta (β_N). + + :param beta: Plasma beta (fraction). + :type beta: float + :param rminor: Plasma minor radius (m). + :type rminor: float + :param c_plasma: Plasma current (A). + :type c_plasma: float + :param b_field: Magnetic field (T). + :type b_field: float + :return: Normalised beta. + :rtype: float + + :Notes: + - 1.0e8 is a conversion factor to get beta_N in standard units, as plasma current is normally in MA and + beta is in percentage instead of fraction. """ - Calculate the plasma masses. - :param m_fuel_amu: Average mass of fuel (amu). - :type m_fuel_amu: float - :param m_ions_total_amu: Average mass of all ions (amu). - :type m_ions_total_amu: float - :param nd_plasma_ions_total_vol_avg: Total ion density (/m3). - :type nd_plasma_ions_total_vol_avg: float - :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (/m3). - :type nd_plasma_fuel_ions_vol_avg: float - :param nd_plasma_alphas_vol_avg: Alpha ash density (/m3). - :type nd_plasma_alphas_vol_avg: float - :param vol_plasma: Plasma volume (m3). - :type vol_plasma: float - :param nd_plasma_electrons_vol_avg: Volume averaged electron density (/m3). - :type nd_plasma_electrons_vol_avg: float - - :returns: A tuple containing: - :rtype: tuple[float, float, float, float, float] - """ - - # Calculate mass of fuel ions - m_plasma_fuel_ions = (m_fuel_amu * constants.ATOMIC_MASS_UNIT) * ( - nd_plasma_fuel_ions_vol_avg * vol_plasma - ) - - m_plasma_ions_total = (m_ions_total_amu * constants.ATOMIC_MASS_UNIT) * ( - nd_plasma_ions_total_vol_avg * vol_plasma - ) - - m_plasma_alpha = (nd_plasma_alphas_vol_avg * vol_plasma) * constants.ALPHA_MASS + return 1.0e8 * (beta * rminor * b_field) / c_plasma - m_plasma_electron = constants.ELECTRON_MASS * ( - nd_plasma_electrons_vol_avg * vol_plasma - ) + @staticmethod + def calculate_plasma_energy_from_beta( + beta: float, b_field: float, vol_plasma: float + ) -> float: + """Calculate plasma thermal energy from beta. - m_plasma = m_plasma_electron + m_plasma_ions_total + E_plasma = 1.5 * β * B² / (2 * μ_0) * V - return ( - m_plasma_fuel_ions, - m_plasma_ions_total, - m_plasma_alpha, - m_plasma_electron, - m_plasma, - ) + :param beta: Plasma beta (fraction). + :type beta: float + :param b_field: Magnetic field (T). + :type b_field: float + :param vol_plasma: Plasma volume (m³). + :type vol_plasma: float + :return: Plasma energy (J). + :rtype: float + """ -def res_diff_time(rmajor, res_plasma, kappa95): - """Calculates resistive diffusion time + return (1.5e0 * beta * b_field**2) / (2.0e0 * constants.RMU0) * vol_plasma - Author: James Morris (UKAEA) + @staticmethod + def calculate_beta_limit_from_norm( + b_plasma_toroidal_on_axis: float, + beta_norm_max: float, + plasma_current: float, + rminor: float, + ) -> float: + """ + Calculate the maximum allowed beta (β) from a given normalised (β_N). - :param rmajor: plasma major radius (m) - :param res_plasma: plasma resistivity (Ohms) - :param kappa95: plasma elongation at 95% flux surface - """ + :param b_plasma_toroidal_on_axis: Toroidal B-field on plasma axis [T]. + :type b_plasma_toroidal_on_axis: float + :param beta_norm_max: Troyon-like g coefficient. + :type beta_norm_max: float + :param plasma_current: Plasma current [A]. + :type plasma_current: float + :param rminor: Plasma minor axis [m]. + :type rminor: float + :return: Beta limit as defined below. + :rtype: float - return 2 * constants.RMU0 * rmajor / (res_plasma * kappa95) + This subroutine calculates the beta limit using the algorithm documented in AEA FUS 172. + The limit applies to beta defined with respect to the total B-field. + Switch i_beta_component determines which components of beta to include. + Notes: + - If i_beta_component = 0, then the limit is applied to the total beta. + - If i_beta_component = 1, then the limit is applied to the thermal beta only. + - If i_beta_component = 2, then the limit is applied to the thermal + neutral beam beta components. + - If i_beta_component = 3, then the limit is applied to the toroidal beta. -def l_h_threshold_power( - nd_plasma_electron_line: float, - b_plasma_toroidal_on_axis: float, - rmajor: float, - rminor: float, - kappa: float, - a_plasma_surface: float, - m_ions_total_amu: float, - aspect: float, - plasma_current: float, -) -> list[float]: - """ - L-mode to H-mode power threshold calculation. + - The default value for the g coefficient is beta_norm_max = 3.5. - :param nd_plasma_electron_line: Line-averaged electron density (/m3) - :type nd_plasma_electron_line: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T) - :type b_plasma_toroidal_on_axis: float - :param rmajor: Plasma major radius (m) - :type rmajor: float - :param rminor: Plasma minor radius (m) - :type rminor: float - :param kappa: Plasma elongation - :type kappa: float - :param a_plasma_surface: Plasma surface area (m2) - :type a_plasma_surface: float - :param m_ions_total_amu: Average mass of all ions (amu) - :type m_ions_total_amu: float - :param aspect: Aspect ratio - :type aspect: float - :param plasma_current: Plasma current (A) - :type plasma_current: float + References: + - F. Troyon et.al, “Beta limit in tokamaks. Experimental and computational status,” + Plasma Physics and Controlled Fusion, vol. 30, no. 11, pp. 1597-1609, Oct. 1988, + doi: https://doi.org/10.1088/0741-3335/30/11/019. - :returns: Array of power thresholds - :rtype: list[float] - """ + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - dnla20 = 1e-20 * nd_plasma_electron_line + """ - # ======================================================================== + # Multiplied by 0.01 to convert from % to fraction + return ( + 0.01 + * beta_norm_max + * (plasma_current / 1.0e6) + / (rminor * b_plasma_toroidal_on_axis) + ) - # ITER-1996 H-mode power threshold database - # Fit to 1996 H-mode power threshold database: nominal + @staticmethod + def calculate_poloidal_beta( + b_plasma_total: float, b_plasma_poloidal_average: float, beta: float + ) -> float: + """Calculates total poloidal beta (β_p) - # i_l_h_threshold = 1 - iterdd = transition.calculate_iter1996_nominal( - dnla20, b_plasma_toroidal_on_axis, rmajor - ) + :type b_plasma_total: float + :param b_plasma_poloidal_average: The average poloidal magnetic field of the plasma (in Tesla). + :type b_plasma_poloidal_average: float + :param beta: The plasma beta, a dimensionless parameter representing the ratio of plasma pressure to magnetic pressure. + :type beta: float + :return: The calculated total poloidal beta. + :rtype: float - # Fit to 1996 H-mode power threshold database: upper bound - # i_l_h_threshold = 2 - iterdd_ub = transition.calculate_iter1996_upper( - dnla20, b_plasma_toroidal_on_axis, rmajor - ) + :references: + - J.P. Freidberg, "Plasma physics and fusion energy", Cambridge University Press (2007) + Page 270 ISBN 0521851076 - # Fit to 1996 H-mode power threshold database: lower bound - # i_l_h_threshold = 3 - iterdd_lb = transition.calculate_iter1996_lower( - dnla20, b_plasma_toroidal_on_axis, rmajor - ) + """ + return beta * (b_plasma_total / b_plasma_poloidal_average) ** 2 - # ======================================================================== + @staticmethod + def fast_alpha_beta( + b_plasma_poloidal_average: float, + b_plasma_toroidal_on_axis: float, + nd_plasma_electrons_vol_avg: float, + nd_plasma_fuel_ions_vol_avg: float, + nd_plasma_ions_total_vol_avg: float, + temp_plasma_electron_density_weighted_kev: float, + temp_plasma_ion_density_weighted_kev: float, + pden_alpha_total_mw: float, + pden_plasma_alpha_mw: float, + i_beta_fast_alpha: int, + ) -> float: + """ + Calculate the fast alpha beta (β_fast_alpha) component. - # Snipes 1997 ITER H-mode power threshold + This function computes the fast alpha beta contribution based on the provided plasma parameters. - # i_l_h_threshold = 4 - snipes_1997 = transition.calculate_snipes1997_iter( - dnla20, b_plasma_toroidal_on_axis, rmajor - ) + :param b_plasma_poloidal_average: Poloidal field (T). + :type b_plasma_poloidal_average: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param nd_plasma_electrons_vol_avg: Electron density (m⁻³). + :type nd_plasma_electrons_vol_avg: float + :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (m⁻³). + :type nd_plasma_fuel_ions_vol_avg: float + :param nd_plasma_ions_total_vol_avg: Total ion density (m⁻³). + :type nd_plasma_ions_total_vol_avg: float + :param temp_plasma_electron_density_weighted_kev: Density-weighted electron temperature (keV). + :type temp_plasma_electron_density_weighted_kev: float + :param temp_plasma_ion_density_weighted_kev: Density-weighted ion temperature (keV). + :type temp_plasma_ion_density_weighted_kev: float + :param pden_alpha_total_mw: Alpha power per unit volume, from beams and plasma (MW/m³). + :type pden_alpha_total_mw: float + :param pden_plasma_alpha_mw: Alpha power per unit volume just from plasma (MW/m³). + :type pden_plasma_alpha_mw: float + :param i_beta_fast_alpha: Switch for fast alpha pressure method. + :type i_beta_fast_alpha: int - # i_l_h_threshold = 5 - snipes_1997_kappa = transition.calculate_snipes1997_kappa( - dnla20, b_plasma_toroidal_on_axis, rmajor, kappa - ) + :return: Fast alpha beta component. + :rtype: float - # ======================================================================== + :Notes: + - For IPDG89 scaling applicability is Z_eff = 1.5, T_i/T_e = 1, 〈T〉 = 5-20 keV - # Martin et al (2008) for recent ITER scaling, with mass correction - # and 95% confidence limits - # i_l_h_threshold = 6 - martin_nominal = transition.calculate_martin08_nominal( - dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu - ) + :References: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', + https://inis.iaea.org/collection/NCLCollectionStore/_Public/21/068/21068960.pdf - # i_l_h_threshold = 7 - martin_ub = transition.calculate_martin08_upper( - dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu - ) + - Uckan, N. A., Tolliver, J. S., Houlberg, W. A., and Attenberger, S. E. + Influence of fast alpha diffusion and thermal alpha buildup on tokamak reactor performance. + United States: N. p., 1987. Web.https://www.osti.gov/servlets/purl/5611706 - # i_l_h_threshold = 8 - martin_lb = transition.calculate_martin08_lower( - dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu - ) + """ - # ======================================================================== + # Determine average fast alpha density + if physics_variables.f_plasma_fuel_deuterium < 1.0: + beta_thermal = ( + 2.0 + * constants.RMU0 + * constants.KILOELECTRON_VOLT + * ( + nd_plasma_electrons_vol_avg + * temp_plasma_electron_density_weighted_kev + + nd_plasma_ions_total_vol_avg * temp_plasma_ion_density_weighted_kev + ) + / (b_plasma_toroidal_on_axis**2 + b_plasma_poloidal_average**2) + ) - # Snipes et al (2000) scaling with mass correction - # Nominal, upper and lower - - # i_l_h_threshold = 9 - snipes_2000 = transition.calculate_snipes2000_nominal( - dnla20, b_plasma_toroidal_on_axis, rmajor, rminor, m_ions_total_amu - ) + # jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T) + # IPDG89 fast alpha scaling + if i_beta_fast_alpha == 0: + fact = min( + 0.3, + 0.29 + * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 + * ( + ( + temp_plasma_electron_density_weighted_kev + + temp_plasma_ion_density_weighted_kev + ) + / 20.0 + - 0.37 + ), + ) - # i_l_h_threshold = 10 - snipes_2000_ub = transition.calculate_snipes2000_upper( - dnla20, b_plasma_toroidal_on_axis, rmajor, rminor, m_ions_total_amu - ) + # Modified scaling, D J Ward + else: + fact = min( + 0.30, + 0.26 + * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 + * np.sqrt( + max( + 0.0, + ( + ( + temp_plasma_electron_density_weighted_kev + + temp_plasma_ion_density_weighted_kev + ) + / 20.0 + - 0.65 + ), + ) + ), + ) - # i_l_h_threshold = 11 - snipes_2000_lb = transition.calculate_snipes2000_lower( - dnla20, b_plasma_toroidal_on_axis, rmajor, rminor, m_ions_total_amu - ) + fact = max(fact, 0.0) + fact2 = pden_alpha_total_mw / pden_plasma_alpha_mw + beta_fast_alpha = beta_thermal * fact * fact2 - # ======================================================================== + else: # negligible alpha production, alpha_power_density = p_beam_alpha_mw = 0 + beta_fast_alpha = 0.0 - # Snipes et al (2000) scaling (closed divertor) with mass correction - # Nominal, upper and lower + return beta_fast_alpha - # i_l_h_threshold = 12 - snipes_2000_cd = transition.calculate_snipes2000_closed_divertor_nominal( - dnla20, b_plasma_toroidal_on_axis, rmajor, m_ions_total_amu - ) + def output_beta_information(self): + """Output beta information to file.""" - # i_l_h_threshold = 13 - snipes_2000_cd_ub = transition.calculate_snipes2000_closed_divertor_upper( - dnla20, b_plasma_toroidal_on_axis, rmajor, m_ions_total_amu - ) + po.osubhd(self.outfile, "Beta Information :") + if physics_variables.i_beta_component == 0: + po.ovarrf( + self.outfile, + "Upper limit on total beta", + "(beta_vol_avg_max)", + physics_variables.beta_vol_avg_max, + "OP ", + ) + elif physics_variables.i_beta_component == 1: + po.ovarrf( + self.outfile, + "Upper limit on thermal beta", + "(beta_vol_avg_max)", + physics_variables.beta_vol_avg_max, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Upper limit on thermal + NB beta", + "(beta_vol_avg_max)", + physics_variables.beta_vol_avg_max, + "OP ", + ) - # i_l_h_threshold = 14 - snipes_2000_cd_lb = transition.calculate_snipes2000_closed_divertor_lower( - dnla20, b_plasma_toroidal_on_axis, rmajor, m_ions_total_amu - ) + po.ovarre( + self.outfile, + "Total plasma beta", + "(beta_total_vol_avg)", + physics_variables.beta_total_vol_avg, + ) + if physics_variables.i_beta_component == 0: + po.ovarrf( + self.outfile, + "Lower limit on total beta", + "(beta_vol_avg_min)", + physics_variables.beta_vol_avg_min, + "IP", + ) + elif physics_variables.i_beta_component == 1: + po.ovarrf( + self.outfile, + "Lower limit on thermal beta", + "(beta_vol_avg_min)", + physics_variables.beta_vol_avg_min, + "IP", + ) + else: + po.ovarrf( + self.outfile, + "Lower limit on thermal + NB beta", + "(beta_vol_avg_min)", + physics_variables.beta_vol_avg_min, + "IP", + ) + po.ovarre( + self.outfile, + "Upper limit on poloidal beta", + "(beta_poloidal_max)", + constraint_variables.beta_poloidal_max, + "IP", + ) + po.ovarre( + self.outfile, + "Total poloidal beta", + "(beta_poloidal_vol_avg)", + physics_variables.beta_poloidal_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Volume averaged toroidal beta", + "(beta_toroidal_vol_avg)", + physics_variables.beta_toroidal_vol_avg, + "OP ", + ) + for i in range(len(physics_variables.beta_thermal_toroidal_profile)): + po.ovarre( + self.mfile, + f"Beta toroidal profile at point {i}", + f"beta_thermal_toroidal_profile{i}", + physics_variables.beta_thermal_toroidal_profile[i], + ) - # ======================================================================== + po.ovarre( + self.outfile, + "Fast alpha beta", + "(beta_fast_alpha)", + physics_variables.beta_fast_alpha, + "OP ", + ) + po.ovarre( + self.outfile, + "Neutral Beam ion beta", + "(beta_beam)", + physics_variables.beta_beam, + "OP ", + ) + po.ovarre( + self.outfile, + "Ratio of fast alpha and beam beta to thermal beta", + "(f_beta_alpha_beam_thermal)", + physics_variables.f_beta_alpha_beam_thermal, + "OP ", + ) - # Hubbard et al. 2012 L-I threshold scaling + po.ovarre( + self.outfile, + "Volume averaged thermal beta", + "(beta_thermal_vol_avg)", + physics_variables.beta_thermal_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermal poloidal beta", + "(beta_thermal_poloidal_vol_avg)", + physics_variables.beta_thermal_poloidal_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermal toroidal beta", + "(beta_thermal_toroidal_vol_avg)", + physics_variables.beta_thermal_toroidal_vol_avg, + "OP ", + ) - # i_l_h_threshold = 15 - hubbard_2012 = transition.calculate_hubbard2012_nominal(plasma_current, dnla20) + po.ovarrf( + self.outfile, + "Poloidal beta and inverse aspect ratio", + "(beta_poloidal_eps)", + physics_variables.beta_poloidal_eps, + "OP ", + ) + po.ovarrf( + self.outfile, + "Poloidal beta and inverse aspect ratio upper limit", + "(beta_poloidal_eps_max)", + physics_variables.beta_poloidal_eps_max, + ) + po.osubhd(self.outfile, "Normalised Beta Information :") + if stellarator_variables.istell == 0: + if physics_variables.i_beta_norm_max != 0: + po.ovarrf( + self.outfile, + "Beta g coefficient", + "(beta_norm_max)", + physics_variables.beta_norm_max, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Beta g coefficient", + "(beta_norm_max)", + physics_variables.beta_norm_max, + ) + po.ovarrf( + self.outfile, + "Normalised total beta", + "(beta_norm_total)", + physics_variables.beta_norm_total, + "OP ", + ) + po.ovarrf( + self.outfile, + "Normalised thermal beta", + "(beta_norm_thermal) ", + physics_variables.beta_norm_thermal, + "OP ", + ) - # i_l_h_threshold = 16 - hubbard_2012_lb = transition.calculate_hubbard2012_lower(plasma_current, dnla20) + po.ovarrf( + self.outfile, + "Normalised toroidal beta", + "(beta_norm_toroidal) ", + physics_variables.beta_norm_toroidal, + "OP ", + ) - # i_l_h_threshold = 17 - hubbard_2012_ub = transition.calculate_hubbard2012_upper(plasma_current, dnla20) - - # ======================================================================== - - # Hubbard et al. 2017 L-I threshold scaling - - # i_l_h_threshold = 18 - hubbard_2017 = transition.calculate_hubbard2017( - dnla20, a_plasma_surface, b_plasma_toroidal_on_axis - ) - - # ======================================================================== - - # Aspect ratio corrected Martin et al (2008) - - # i_l_h_threshold = 19 - martin_nominal_aspect = transition.calculate_martin08_aspect_nominal( - dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu, aspect - ) - - # i_l_h_threshold = 20 - martin_ub_aspect = transition.calculate_martin08_aspect_upper( - dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu, aspect - ) - - # i_l_h_threshold = 21 - martin_lb_aspect = transition.calculate_martin08_aspect_lower( - dnla20, b_plasma_toroidal_on_axis, a_plasma_surface, m_ions_total_amu, aspect - ) - - # ======================================================================== - - return [ - iterdd, - iterdd_ub, - iterdd_lb, - snipes_1997, - snipes_1997_kappa, - martin_nominal, - martin_ub, - martin_lb, - snipes_2000, - snipes_2000_ub, - snipes_2000_lb, - snipes_2000_cd, - snipes_2000_cd_ub, - snipes_2000_cd_lb, - hubbard_2012, - hubbard_2012_lb, - hubbard_2012_ub, - hubbard_2017, - martin_nominal_aspect, - martin_ub_aspect, - martin_lb_aspect, - ] + po.ovarrf( + self.outfile, + "Normalised poloidal beta", + "(beta_norm_poloidal) ", + physics_variables.beta_norm_poloidal, + "OP ", + ) + po.osubhd(self.outfile, "Maximum normalised beta scalings :") + po.ovarrf( + self.outfile, + "J. Wesson normalised beta upper limit", + "(beta_norm_max_wesson) ", + physics_variables.beta_norm_max_wesson, + "OP ", + ) + po.ovarrf( + self.outfile, + "Original normalsied beta upper limit", + "(beta_norm_max_original_scaling) ", + physics_variables.beta_norm_max_original_scaling, + "OP ", + ) + po.ovarrf( + self.outfile, + "J. Menard normalised beta upper limit", + "(beta_norm_max_menard) ", + physics_variables.beta_norm_max_menard, + "OP ", + ) + po.ovarrf( + self.outfile, + "E. Thloreus normalised beta upper limit", + "(beta_norm_max_thloreus) ", + physics_variables.beta_norm_max_thloreus, + "OP ", + ) + po.ovarrf( + self.outfile, + "R. Stambaugh normalised beta upper limit", + "(beta_norm_max_stambaugh) ", + physics_variables.beta_norm_max_stambaugh, + "OP ", + ) -def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa, lhat): - """Function for calculating upstream temperature(keV) in Reinke model - author: H Lux, CCFE/UKAEA - b_plasma_toroidal_on_axis : input real : toroidal field on axis (T) - flh : input real : fraction of Psep/P_LH - qstar : input real : safety factor similar to q95 (see #707) - rmajor : input real : major radius (m) - eps : input real : inverse aspect ratio - fgw : input real : ratio of volume averaged density to n_GW - kappa : input real : elongation - lhat : input real : connection length factor - This function calculates the upstream temperature in the - divertor/SoL model used for the Reinke citerion. - Issue #707 - M.L. Reinke 2017 Nucl. Fusion 57 034004 - """ - kappa_0 = 2.0e3 # Stangeby W/m/eV^(7/2) + po.osubhd(self.outfile, "Plasma energies derived from beta :") + po.ovarre( + self.outfile, + "Plasma thermal energy derived from thermal beta (J)", + "(e_plasma_beta_thermal) ", + physics_variables.e_plasma_beta_thermal, + "OP", + ) - return ( - ( - b_plasma_toroidal_on_axis**0.72 - * flh**0.29 - * fgw**0.21 - * qstar**0.08 - * rmajor**0.33 + po.ovarre( + self.outfile, + "Plasma thermal energy derived from the total beta (J)", + "(e_plasma_beta)", + physics_variables.e_plasma_beta, + "OP", ) - * (eps**0.15 * (1.0 + kappa**2.0) ** 0.34) - * (lhat**0.29 * kappa_0 ** (-0.29) * 0.285) - ) + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) -class BetaNormMaxModel(IntEnum): - """Beta norm max (β_N_max) model types""" + +class IndInternalNormModel(IntEnum): + """Normalised internal inductance (l_i) model types""" USER_INPUT = 0 WESSON = 1 - ORIGINAL_SCALING = 2 - MENARD = 3 - THLOREUS = 4 - STAMBAUGH = 5 + MENARD = 2 -class PlasmaBeta: - """Class to hold plasma beta calculations for plasma processing.""" +class PlasmaInductance: + """Class to hold plasma inductance calculations for plasma processing.""" def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE - def get_beta_norm_max_value(self, model: BetaNormMaxModel) -> float: - """Get the beta norm max value (β_N_max) for the specified model.""" - model_map = { - BetaNormMaxModel.USER_INPUT: physics_variables.beta_norm_max, - BetaNormMaxModel.WESSON: physics_variables.beta_norm_max_wesson, - BetaNormMaxModel.ORIGINAL_SCALING: physics_variables.beta_norm_max_original_scaling, - BetaNormMaxModel.MENARD: physics_variables.beta_norm_max_menard, - BetaNormMaxModel.THLOREUS: physics_variables.beta_norm_max_thloreus, - BetaNormMaxModel.STAMBAUGH: physics_variables.beta_norm_max_stambaugh, - } - return model_map[model] - - def run(self) -> None: - """ - Calculate plasma beta values. - """ - - # ----------------------------------------------------- - # Normalised Beta Limit - # ----------------------------------------------------- - - # Normalised beta from Troyon beta limit - physics_variables.beta_norm_total = self.calculate_normalised_beta( - beta=physics_variables.beta_total_vol_avg, - rminor=physics_variables.rminor, - c_plasma=physics_variables.plasma_current, - b_field=physics_variables.b_plasma_toroidal_on_axis, - ) - - # Define beta_norm_max calculations - - physics_variables.beta_norm_max_wesson = self.calculate_beta_norm_max_wesson( - ind_plasma_internal_norm=physics_variables.ind_plasma_internal_norm - ) - - # Original scaling law - physics_variables.beta_norm_max_original_scaling = ( - self.calculate_beta_norm_max_original(eps=physics_variables.eps) - ) - - # J. Menard scaling law - physics_variables.beta_norm_max_menard = self.calculate_beta_norm_max_menard( - eps=physics_variables.eps + def run(self): + physics_variables.ind_plasma_internal_norm_wesson = ( + self.calculate_internal_inductance_wesson(alphaj=physics_variables.alphaj) ) - # E. Tholerus scaling law - physics_variables.beta_norm_max_thloreus = self.calculate_beta_norm_max_thloreus( - c_beta=physics_variables.c_beta, - pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis, - pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg, + # Spherical Tokamak relation for internal inductance + # Menard et al. (2016), Nuclear Fusion, 56, 106023 + physics_variables.ind_plasma_internal_norm_menard = ( + self.calculate_internal_inductance_menard(kappa=physics_variables.kappa) ) - # R. D. Stambaugh scaling law - physics_variables.beta_norm_max_stambaugh = ( - self.calculate_beta_norm_max_stambaugh( - f_c_plasma_bootstrap=current_drive_variables.f_c_plasma_bootstrap, - kappa=physics_variables.kappa, - aspect=physics_variables.aspect, + physics_variables.ind_plasma_internal_norm_iter_3 = ( + self.calculate_normalised_internal_inductance_iter_3( + b_plasma_poloidal_vol_avg=physics_variables.b_plasma_poloidal_average, + c_plasma=physics_variables.plasma_current, + vol_plasma=physics_variables.vol_plasma, + rmajor=physics_variables.rmajor, ) ) - # Calculate beta_norm_max based on i_beta_norm_max + # Calculate ind_plasma_internal_norm based on i_ind_plasma_internal_norm try: - model = BetaNormMaxModel(int(physics_variables.i_beta_norm_max)) - physics_variables.beta_norm_max = self.get_beta_norm_max_value(model) + model = IndInternalNormModel( + int(physics_variables.i_ind_plasma_internal_norm) + ) + physics_variables.ind_plasma_internal_norm = ( + self.get_ind_internal_norm_value(model) + ) except ValueError: raise ProcessValueError( - "Illegal value of i_beta_norm_max", - i_beta_norm_max=physics_variables.i_beta_norm_max, + "Illegal value of i_ind_plasma_internal_norm", + i_ind_plasma_internal_norm=physics_variables.i_ind_plasma_internal_norm, ) from None - # calculate_beta_limit() returns the beta_vol_avg_max for beta - physics_variables.beta_vol_avg_max = self.calculate_beta_limit_from_norm( - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - beta_norm_max=physics_variables.beta_norm_max, - plasma_current=physics_variables.plasma_current, - rminor=physics_variables.rminor, - ) + def get_ind_internal_norm_value(self, model: IndInternalNormModel) -> float: + """Get the normalised internal inductance (l_i) for the specified model.""" + model_map = { + IndInternalNormModel.USER_INPUT: physics_variables.ind_plasma_internal_norm, + IndInternalNormModel.WESSON: physics_variables.ind_plasma_internal_norm_wesson, + IndInternalNormModel.MENARD: physics_variables.ind_plasma_internal_norm_menard, + } + return model_map[model] - physics_variables.beta_toroidal_vol_avg = ( - physics_variables.beta_total_vol_avg - * physics_variables.b_plasma_total**2 - / physics_variables.b_plasma_toroidal_on_axis**2 - ) + @staticmethod + def calculate_volt_second_requirements( + csawth: float, + eps: float, + f_c_plasma_inductive: float, + ejima_coeff: float, + kappa: float, + rmajor: float, + res_plasma: float, + plasma_current: float, + t_plant_pulse_fusion_ramp: float, + t_plant_pulse_burn: float, + ind_plasma_internal_norm: float, + ) -> tuple[float, float, float, float, float, float]: + """Calculate the volt-second requirements and related parameters for plasma physics. - # Calculate physics_variables.beta poloidal [-] - physics_variables.beta_poloidal_vol_avg = self.calculate_poloidal_beta( - b_plasma_total=physics_variables.b_plasma_total, - b_plasma_poloidal_average=physics_variables.b_plasma_poloidal_average, - beta=physics_variables.beta_total_vol_avg, - ) + :param csawth: Coefficient for sawteeth effects + :type csawth: float + :param eps: Inverse aspect ratio + :type eps: float + :param f_c_plasma_inductive: Fraction of plasma current produced inductively + :type f_c_plasma_inductive: float + :param ejima_coeff: Ejima coefficient for resistive start-up V-s component + :type ejima_coeff: float + :param kappa: Plasma elongation + :type kappa: float + :param rmajor: Plasma major radius (m) + :type rmajor: float + :param res_plasma: Plasma resistance (ohm) + :type res_plasma: float + :param plasma_current: Plasma current (A) + :type plasma_current: float + :param t_plant_pulse_fusion_ramp: Heating time (s) + :type t_plant_pulse_fusion_ramp: float + :param t_plant_pulse_burn: Burn time (s) + :type t_plant_pulse_burn: float + :param ind_plasma_internal_norm: Plasma normalized internal inductance + :type ind_plasma_internal_norm: float - physics_variables.beta_thermal_vol_avg = ( - physics_variables.beta_total_vol_avg - - physics_variables.beta_fast_alpha - - physics_variables.beta_beam - ) + :return: A tuple containing: + - vs_plasma_internal: Internal plasma volt-seconds (Wb) + - ind_plasma_internal: Plasma inductance (H) + - vs_plasma_burn_required: Volt-seconds needed during flat-top (heat+burn) (Wb) + - vs_plasma_ramp_required: Volt-seconds needed during ramp-up (Wb) + - ind_plasma_total,: Internal and external plasma inductance V-s (Wb) + - vs_res_ramp: Resistive losses in start-up volt-seconds (Wb) + - vs_plasma_total_required: Total volt-seconds needed (Wb) + :rtype: tuple[float, float, float, float, float, float] - physics_variables.beta_poloidal_eps = ( - physics_variables.beta_poloidal_vol_avg * physics_variables.eps - ) + :notes: - physics_variables.beta_thermal_poloidal_vol_avg = ( - physics_variables.beta_thermal_vol_avg - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_poloidal_average - ) - ** 2 - ) - physics_variables.beta_thermal_toroidal_vol_avg = ( - physics_variables.beta_thermal_vol_avg - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_toroidal_on_axis - ) - ** 2 - ) + :references: + - S. Ejima, R. W. Callis, J. L. Luxon, R. D. Stambaugh, T. S. Taylor, and J. C. Wesley, + “Volt-second analysis and consumption in Doublet III plasmas,” + Nuclear Fusion, vol. 22, no. 10, pp. 1313-1319, Oct. 1982, doi: + https://doi.org/10.1088/0029-5515/22/10/006. - # ======================================================= + - S. C. Jardin, C. E. Kessel, and N Pomphrey, + “Poloidal flux linkage requirements for the International Thermonuclear Experimental Reactor,” + Nuclear Fusion, vol. 34, no. 8, pp. 1145-1160, Aug. 1994, + doi: https://doi.org/10.1088/0029-5515/34/8/i07. - # Mirror the pressure profiles to match the doubled toroidal field profile - pres_profile_total = np.concatenate([ - physics_variables.pres_plasma_thermal_total_profile[::-1], - physics_variables.pres_plasma_thermal_total_profile, - ]) + - S. P. Hirshman and G. H. Neilson, “External inductance of an axisymmetric plasma,” + The Physics of Fluids, vol. 29, no. 3, pp. 790-793, Mar. 1986, + doi: https://doi.org/10.1063/1.865934. + """ + # Plasma internal inductance - physics_variables.beta_thermal_toroidal_profile = np.array([ - self.calculate_plasma_beta( - pres_plasma=pres_profile_total[i], - b_field=physics_variables.b_plasma_toroidal_profile[i], - ) - for i in range(len(physics_variables.b_plasma_toroidal_profile)) - ]) + ind_plasma_internal = constants.RMU0 * rmajor * ind_plasma_internal_norm / 2.0 - # ======================================================= + # Internal plasma flux (V-s) component + vs_plasma_internal = ind_plasma_internal * plasma_current - physics_variables.beta_norm_thermal = self.calculate_normalised_beta( - beta=physics_variables.beta_thermal_vol_avg, - rminor=physics_variables.rminor, - c_plasma=physics_variables.plasma_current, - b_field=physics_variables.b_plasma_toroidal_on_axis, - ) + # Start-up resistive component + # Uses ITER formula without the 10 V-s add-on - physics_variables.beta_norm_toroidal = ( - physics_variables.beta_norm_total - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_toroidal_on_axis - ) - ** 2 - ) - physics_variables.beta_norm_poloidal = ( - physics_variables.beta_norm_total - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_poloidal_average - ) - ** 2 - ) + vs_res_ramp = ejima_coeff * constants.RMU0 * plasma_current * rmajor - physics_variables.f_beta_alpha_beam_thermal = ( - physics_variables.beta_fast_alpha + physics_variables.beta_beam - ) / physics_variables.beta_thermal_vol_avg + # ====================================================================== - # Plasma thermal energy derived from the thermal beta - physics_variables.e_plasma_beta_thermal = self.calculate_plasma_energy_from_beta( - beta=physics_variables.beta_thermal_vol_avg, - b_field=physics_variables.b_plasma_total, - vol_plasma=physics_variables.vol_plasma, + # Hirshman and Neilson fit for external inductance + + aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( + 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps ) + beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) - # Plasma thermal energy derived from the total beta - physics_variables.e_plasma_beta = self.calculate_plasma_energy_from_beta( - beta=physics_variables.beta_total_vol_avg, - b_field=physics_variables.b_plasma_total, - vol_plasma=physics_variables.vol_plasma, + ind_plasma_external = ( + rmajor * constants.RMU0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) ) - @staticmethod - def calculate_plasma_beta( - pres_plasma: float | np.ndarray, b_field: float | np.ndarray - ) -> float | np.ndarray: - """ - Calculate the plasma beta (β) for a given pressure and field. + # ====================================================================== - :param pres_plasma: Plasma pressure (in Pascals). - :type pres_plasma: float | np.ndarray - :param b_field: Magnetic field strength (in Tesla). - :type b_field: float | np.ndarray - :returns: The plasma beta (dimensionless). - :rtype: float | np.ndarray + ind_plasma_total = ind_plasma_external + ind_plasma_internal - Plasma beta is the ratio of plasma pressure to magnetic pressure. - """ + # Inductive V-s component - return 2 * constants.RMU0 * pres_plasma / (b_field**2) + vs_self_ind_ramp = ind_plasma_total * plasma_current + vs_plasma_ramp_required = vs_res_ramp + vs_self_ind_ramp - @staticmethod - def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: - """ - Calculate the Wesson normalsied beta upper limit. + # Plasma loop voltage during flat-top + # Include enhancement factor in flattop V-s requirement + # to account for MHD sawtooth effects. - :param ind_plasma_internal_norm: Plasma normalised internal inductance - :type ind_plasma_internal_norm: float + v_plasma_loop_burn = plasma_current * res_plasma * f_c_plasma_inductive - :return: The Wesson normalised beta upper limit. - :rtype: float + v_burn_resistive = v_plasma_loop_burn * csawth - :Notes: - - It is recommended to use this method with the other Wesson relations for normalsied internal - inductance and current profile index. - - This fit is derived from the DIII-D database for β_N >= 2.5 + # N.B. t_plant_pulse_burn on first iteration will not be correct + # if the pulsed reactor option is used, but the value + # will be correct on subsequent calls. - :References: - - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, - International Series of Monographs on Physics, Volume 149. + vs_plasma_burn_required = v_burn_resistive * ( + t_plant_pulse_fusion_ramp + t_plant_pulse_burn + ) + vs_plasma_total_required = vs_plasma_ramp_required + vs_plasma_burn_required - - T. T. S et al., “Profile Optimization and High Beta Discharges and Stability of High Elongation Plasmas in the DIII-D Tokamak,” - Osti.gov, Oct. 1990. https://www.osti.gov/biblio/6194284 (accessed Dec. 19, 2024). - """ - return 4 * ind_plasma_internal_norm + return ( + vs_plasma_internal, + ind_plasma_total, + vs_plasma_burn_required, + vs_plasma_ramp_required, + vs_self_ind_ramp, + vs_res_ramp, + vs_plasma_total_required, + v_plasma_loop_burn, + ) @staticmethod - def calculate_beta_norm_max_original(eps: float) -> float: + def calculate_normalised_internal_inductance_iter_3( + b_plasma_poloidal_vol_avg: float, + c_plasma: float, + vol_plasma: float, + rmajor: float, + ) -> float: """ - Calculate the original scaling law normalsied beta upper limit. + Calculate the normalised internal inductance using ITER-3 scaling li(3). - :param eps: Plasma normalised internal inductance - :type eps: float + :param b_plasma_poloidal_vol_avg: Volume-averaged poloidal magnetic field (T). + :type b_plasma_poloidal_vol_avg: float + :param c_plasma: Plasma current (A). + :type c_plasma: float + :param vol_plasma: Plasma volume (m^3). + :type vol_plasma: float + :param rmajor: Plasma major radius (m). + :type rmajor: float - :return: The original scaling law normalised beta upper limit. + :returns: The li(3) normalised internal inductance. :rtype: float - :Notes: - - :References: - - """ - return 2.7 * (1.0 + 5.0 * eps**3.5) + :references: + - T. C. Luce, D. A. Humphreys, G. L. Jackson, and W. M. Solomon, + “Inductive flux usage and its optimization in tokamak operation,” + Nuclear Fusion, vol. 54, no. 9, p. 093005, Jul. 2014, + doi: https://doi.org/10.1088/0029-5515/54/9/093005. - @staticmethod - def calculate_beta_norm_max_menard(eps: float) -> float: + - G. L. Jackson et al., “ITER startup studies in the DIII-D tokamak,” + Nuclear Fusion, vol. 48, no. 12, p. 125002, Nov. 2008, + doi: https://doi.org/10.1088/0029-5515/48/12/125002. + ‌ """ - Calculate the Menard normalsied beta upper limit. - :param eps: Plasma normalised internal inductance - :type eps: float - - :return: The Menard normalised beta upper limit. - :rtype: float - - :Notes: - - Found as a reasonable fit to the computed no wall limit at f_BS ≈ 50% - - Uses maximum κ data from NSTX at A = 1.45, A = 1.75. Along with record - β_T data from DIII-D at A = 2.9 and high κ. - - :References: - - # J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” - Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, - doi: https://doi.org/10.1088/0029-5515/56/10/106023. - - """ - return 3.12 + 3.5 * eps**1.7 + return ( + 2 + * vol_plasma + * b_plasma_poloidal_vol_avg**2 + / (constants.RMU0**2 * c_plasma**2 * rmajor) + ) @staticmethod - def calculate_beta_norm_max_thloreus( - c_beta: float, pres_plasma_on_axis: float, pres_plasma_vol_avg: float - ) -> float: + def calculate_internal_inductance_menard(kappa: float) -> float: """ - Calculate the E. Tholerus normalized beta upper limit. + Calculate the Menard plasma normalized internal inductance. - :param c_beta: Pressure peaking factor coefficient. - :type c_beta: float - :param pres_plasma_on_axis: Central plasma pressure (Pa). - :type pres_plasma_on_axis: float - :param pres_plasma_vol_avg: Volume-averaged plasma pressure (Pa). - :type pres_plasma_vol_avg: float + :param kappa: Plasma separatrix elongation. + :type kappa: float - :return: The E. Tholerus normalized beta upper limit. + :return: The Menard plasma normalised internal inductance. :rtype: float :Notes: - - This method calculates the normalized beta upper limit based on the pressure peaking factor (Fp), - which is defined as the ratio of the peak pressure to the average pressure. - - The formula is derived from operational space studies of flat-top plasma in the STEP power plant. + - This relation is based off of data from NSTX for l_i in the range of 0.4-0.85 + - This model is only recommneded to be used for ST's with kappa > 2.5 :References: - - E. Tholerus et al., “Flat-top plasma operational space of the STEP power plant,” - Nuclear Fusion, Aug. 2024, doi: https://doi.org/10.1088/1741-4326/ad6ea2. + - J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” + Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, + doi: https://doi.org/10.1088/0029-5515/56/10/106023. """ - return 3.7 + ( - (c_beta / (pres_plasma_on_axis / pres_plasma_vol_avg)) - * (12.5 - 3.5 * (pres_plasma_on_axis / pres_plasma_vol_avg)) - ) + return 3.4 - kappa @staticmethod - def calculate_beta_norm_max_stambaugh( - f_c_plasma_bootstrap: float, - kappa: float, - aspect: float, - ) -> float: + def calculate_internal_inductance_wesson(alphaj: float) -> float: """ - Calculate the Stambaugh normalized beta upper limit. + Calculate the Wesson plasma normalized internal inductance. - :param f_c_plasma_bootstrap: Bootstrap current fraction. - :type f_c_plasma_bootstrap: float - :param kappa: Plasma separatrix elongation. - :type kappa: float - :param aspect: Plasma aspect ratio. - :type aspect: float + :param alphaj: Current profile index. + :type alphaj: float - :return: The Stambaugh normalized beta upper limit. + :return: The Wesson plasma normalised internal inductance. :rtype: float :Notes: - - This method calculates the normalized beta upper limit based on the Stambaugh scaling. - - The formula is derived from empirical fits to high-performance, steady-state tokamak equilibria. + - It is recommended to use this method with the other Wesson relations for normalised beta and + current profile index. + - This relation is only true for the cyclindrical plasma approximation with parabolic profiles. :References: - - R. D. Stambaugh et al., “Fusion Nuclear Science Facility Candidates,” - Fusion Science and Technology, vol. 59, no. 2, pp. 279-307, Feb. 2011, - doi: https://doi.org/10.13182/fst59-279. - - - Y. R. Lin-Liu and R. D. Stambaugh, “Optimum equilibria for high performance, steady state tokamaks,” - Nuclear Fusion, vol. 44, no. 4, pp. 548-554, Mar. 2004, - doi: https://doi.org/10.1088/0029-5515/44/4/009. + - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, + International Series of Monographs on Physics, Volume 149. """ - return ( - f_c_plasma_bootstrap - * 10 - * (-0.7748 + (1.2869 * kappa) - (0.2921 * kappa**2) + (0.0197 * kappa**3)) - / (aspect**0.5523 * np.tanh((1.8524 + (0.2319 * kappa)) / aspect**0.6163)) + return np.log(1.65 + 0.89 * alphaj) + + def output_volt_second_information(self): + """Output volt-second information to file.""" + + po.osubhd(self.outfile, "Plasma Volt-second Requirements :") + po.ovarre( + self.outfile, + "Total plasma volt-seconds required for pulse (Wb)", + "(vs_plasma_total_required)", + physics_variables.vs_plasma_total_required, + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Total plasma inductive flux consumption for plasma current ramp-up (Wb)", + "(vs_plasma_ind_ramp)", + physics_variables.vs_plasma_ind_ramp, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma resistive flux consumption for plasma current ramp-up (Wb)", + "(vs_plasma_res_ramp)", + physics_variables.vs_plasma_res_ramp, + "OP ", + ) + po.ovarre( + self.outfile, + "Total flux consumption for plasma current ramp-up (Wb)", + "(vs_plasma_ramp_required)", + physics_variables.vs_plasma_ramp_required, + "OP ", + ) + po.ovarrf( + self.outfile, + "Ejima coefficient", + "(ejima_coeff)", + physics_variables.ejima_coeff, + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Internal plasma V-s", + "(vs_plasma_internal)", + physics_variables.vs_plasma_internal, ) - @staticmethod - def calculate_normalised_beta( - beta: float, rminor: float, c_plasma: float, b_field: float - ) -> float: - """Calculate normalised beta (β_N). + po.ovarre( + self.outfile, + "Plasma volt-seconds needed for flat-top (heat + burn times) (Wb)", + "(vs_plasma_burn_required)", + physics_variables.vs_plasma_burn_required, + "OP ", + ) - :param beta: Plasma beta (fraction). - :type beta: float - :param rminor: Plasma minor radius (m). - :type rminor: float - :param c_plasma: Plasma current (A). - :type c_plasma: float - :param b_field: Magnetic field (T). - :type b_field: float - :return: Normalised beta. - :rtype: float + po.ovarre( + self.outfile, + "Plasma loop voltage during burn (V)", + "(v_plasma_loop_burn)", + physics_variables.v_plasma_loop_burn, + "OP ", + ) + po.ovarrf( + self.outfile, + "Coefficient for sawtooth effects on burn V-s requirement", + "(csawth)", + physics_variables.csawth, + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Plasma resistance (ohm)", + "(res_plasma)", + physics_variables.res_plasma, + "OP ", + ) - :Notes: - - 1.0e8 is a conversion factor to get beta_N in standard units, as plasma current is normally in MA and - beta is in percentage instead of fraction. - """ + po.ovarre( + self.outfile, + "Plasma resistive diffusion time (s)", + "(t_plasma_res_diffusion)", + physics_variables.t_plasma_res_diffusion, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma inductance (H)", + "(ind_plasma)", + physics_variables.ind_plasma, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma magnetic energy stored (J)", + "(e_plasma_magnetic_stored)", + physics_variables.e_plasma_magnetic_stored, + "OP ", + ) + po.ovarrf( + self.outfile, + "Plasma normalised internal inductance", + "(ind_plasma_internal_norm)", + physics_variables.ind_plasma_internal_norm, + "OP ", + ) - return 1.0e8 * (beta * rminor * b_field) / c_plasma + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) - @staticmethod - def calculate_plasma_energy_from_beta( - beta: float, b_field: float, vol_plasma: float - ) -> float: - """Calculate plasma thermal energy from beta. - E_plasma = 1.5 * β * B² / (2 * μ_0) * V +class BootstrapCurrentFractionModel(IntEnum): + """Bootstrap plasma current fraction (f_BS) model types""" - :param beta: Plasma beta (fraction). - :type beta: float - :param b_field: Magnetic field (T). - :type b_field: float - :param vol_plasma: Plasma volume (m³). - :type vol_plasma: float - :return: Plasma energy (J). - :rtype: float + USER_INPUT = 0 + ITER_89 = 1 + NEVINS = 2 + WILSON = 3 + SAUTER = 4 + SAKAI = 5 + ARIES = 6 + ANDRADE = 7 + HOANG = 8 + WONG = 9 + GI_1 = 10 + GI_2 = 11 + SUGIYAMA_L_MODE = 12 + SUGIYAMA_H_MODE = 13 + + +class PlasmaBootstrapCurrent: + """Class to hold plasma bootstrap current for plasma processing.""" - """ + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE - return (1.5e0 * beta * b_field**2) / (2.0e0 * constants.RMU0) * vol_plasma + def get_bootstrap_current_fraction_value( + self, model: BootstrapCurrentFractionModel + ) -> float: + """Get the plasma current bootstrap fraction (f_BS) for the specified model.""" + model_map = { + BootstrapCurrentFractionModel.USER_INPUT: current_drive_variables.f_c_plasma_bootstrap, + BootstrapCurrentFractionModel.ITER_89: current_drive_variables.f_c_plasma_bootstrap_iter89, + BootstrapCurrentFractionModel.NEVINS: current_drive_variables.f_c_plasma_bootstrap_nevins, + BootstrapCurrentFractionModel.WILSON: current_drive_variables.f_c_plasma_bootstrap_wilson, + BootstrapCurrentFractionModel.SAUTER: current_drive_variables.f_c_plasma_bootstrap_sauter, + BootstrapCurrentFractionModel.SAKAI: current_drive_variables.f_c_plasma_bootstrap_sakai, + BootstrapCurrentFractionModel.ARIES: current_drive_variables.f_c_plasma_bootstrap_aries, + BootstrapCurrentFractionModel.ANDRADE: current_drive_variables.f_c_plasma_bootstrap_andrade, + BootstrapCurrentFractionModel.HOANG: current_drive_variables.f_c_plasma_bootstrap_hoang, + BootstrapCurrentFractionModel.WONG: current_drive_variables.f_c_plasma_bootstrap_wong, + BootstrapCurrentFractionModel.GI_1: current_drive_variables.bscf_gi_i, + BootstrapCurrentFractionModel.GI_2: current_drive_variables.bscf_gi_ii, + BootstrapCurrentFractionModel.SUGIYAMA_L_MODE: current_drive_variables.f_c_plasma_bootstrap_sugiyama_l, + BootstrapCurrentFractionModel.SUGIYAMA_H_MODE: current_drive_variables.f_c_plasma_bootstrap_sugiyama_h, + } + return model_map[model] @staticmethod - def calculate_beta_limit_from_norm( + def bootstrap_fraction_iter89( + aspect: float, + beta: float, b_plasma_toroidal_on_axis: float, - beta_norm_max: float, plasma_current: float, - rminor: float, + q95: float, + q0: float, + rmajor: float, + vol_plasma: float, ) -> float: """ - Calculate the maximum allowed beta (β) from a given normalised (β_N). + Calculate the bootstrap-driven fraction of the plasma current. - :param b_plasma_toroidal_on_axis: Toroidal B-field on plasma axis [T]. - :type b_plasma_toroidal_on_axis: float - :param beta_norm_max: Troyon-like g coefficient. - :type beta_norm_max: float - :param plasma_current: Plasma current [A]. - :type plasma_current: float - :param rminor: Plasma minor axis [m]. - :type rminor: float - :return: Beta limit as defined below. - :rtype: float - - This subroutine calculates the beta limit using the algorithm documented in AEA FUS 172. - The limit applies to beta defined with respect to the total B-field. - Switch i_beta_component determines which components of beta to include. - - Notes: - - If i_beta_component = 0, then the limit is applied to the total beta. - - If i_beta_component = 1, then the limit is applied to the thermal beta only. - - If i_beta_component = 2, then the limit is applied to the thermal + neutral beam beta components. - - If i_beta_component = 3, then the limit is applied to the toroidal beta. - - - The default value for the g coefficient is beta_norm_max = 3.5. + Args: + aspect (float): Plasma aspect ratio. + beta (float): Plasma total beta. + b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). + plasma_current (float): Plasma current (A). + q95 (float): Safety factor at 95% surface. + q0 (float): Central safety factor. + rmajor (float): Plasma major radius (m). + vol_plasma (float): Plasma volume (m3). - References: - - F. Troyon et.al, “Beta limit in tokamaks. Experimental and computational status,” - Plasma Physics and Controlled Fusion, vol. 30, no. 11, pp. 1597-1609, Oct. 1988, - doi: https://doi.org/10.1088/0741-3335/30/11/019. + Returns: + float: The bootstrap-driven fraction of the plasma current. - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + This function performs the original ITER calculation of the plasma current bootstrap fraction. + Reference: + - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 """ - # Multiplied by 0.01 to convert from % to fraction - return ( - 0.01 - * beta_norm_max - * (plasma_current / 1.0e6) - / (rminor * b_plasma_toroidal_on_axis) - ) + # Calculate the bootstrap current coefficient + c_bs = 1.32 - 0.235 * (q95 / q0) + 0.0185 * (q95 / q0) ** 2 - @staticmethod - def calculate_poloidal_beta( - b_plasma_total: float, b_plasma_poloidal_average: float, beta: float - ) -> float: - """Calculates total poloidal beta (β_p) + # Calculate the average minor radius + average_a = np.sqrt(vol_plasma / (2 * np.pi**2 * rmajor)) - :type b_plasma_total: float - :param b_plasma_poloidal_average: The average poloidal magnetic field of the plasma (in Tesla). - :type b_plasma_poloidal_average: float - :param beta: The plasma beta, a dimensionless parameter representing the ratio of plasma pressure to magnetic pressure. - :type beta: float - :return: The calculated total poloidal beta. - :rtype: float + b_pa = (plasma_current / 1e6) / (5 * average_a) - :references: - - J.P. Freidberg, "Plasma physics and fusion energy", Cambridge University Press (2007) - Page 270 ISBN 0521851076 + # Calculate the poloidal beta for bootstrap current + betapbs = beta * (b_plasma_toroidal_on_axis / b_pa) ** 2 - """ - return beta * (b_plasma_total / b_plasma_poloidal_average) ** 2 + # Ensure betapbs is positive + if betapbs <= 0.0: + return 0.0 + + # Calculate and return the bootstrap current fraction + return c_bs * (betapbs / np.sqrt(aspect)) ** 1.3 @staticmethod - def fast_alpha_beta( - b_plasma_poloidal_average: float, - b_plasma_toroidal_on_axis: float, - nd_plasma_electrons_vol_avg: float, - nd_plasma_fuel_ions_vol_avg: float, - nd_plasma_ions_total_vol_avg: float, - temp_plasma_electron_density_weighted_kev: float, - temp_plasma_ion_density_weighted_kev: float, - pden_alpha_total_mw: float, - pden_plasma_alpha_mw: float, - i_beta_fast_alpha: int, + def bootstrap_fraction_wilson( + alphaj: float, + alphap: float, + alphat: float, + betpth: float, + q0: float, + q95: float, + rmajor: float, + rminor: float, ) -> float: """ - Calculate the fast alpha beta (β_fast_alpha) component. + Bootstrap current fraction from Wilson et al scaling - This function computes the fast alpha beta contribution based on the provided plasma parameters. + Args: + alphaj (float): Current profile index. + alphap (float): Pressure profile index. + alphat (float): Temperature profile index. + betpth (float): Thermal component of poloidal beta. + q0 (float): Safety factor on axis. + q95 (float): Edge safety factor. + rmajor (float): Major radius (m). + rminor (float): Minor radius (m). - :param b_plasma_poloidal_average: Poloidal field (T). - :type b_plasma_poloidal_average: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param nd_plasma_electrons_vol_avg: Electron density (m⁻³). - :type nd_plasma_electrons_vol_avg: float - :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (m⁻³). - :type nd_plasma_fuel_ions_vol_avg: float - :param nd_plasma_ions_total_vol_avg: Total ion density (m⁻³). - :type nd_plasma_ions_total_vol_avg: float - :param temp_plasma_electron_density_weighted_kev: Density-weighted electron temperature (keV). - :type temp_plasma_electron_density_weighted_kev: float - :param temp_plasma_ion_density_weighted_kev: Density-weighted ion temperature (keV). - :type temp_plasma_ion_density_weighted_kev: float - :param pden_alpha_total_mw: Alpha power per unit volume, from beams and plasma (MW/m³). - :type pden_alpha_total_mw: float - :param pden_plasma_alpha_mw: Alpha power per unit volume just from plasma (MW/m³). - :type pden_plasma_alpha_mw: float - :param i_beta_fast_alpha: Switch for fast alpha pressure method. - :type i_beta_fast_alpha: int + Returns: + float: The bootstrap current fraction. - :return: Fast alpha beta component. - :rtype: float + This function calculates the bootstrap current fraction using the numerically fitted algorithm written by Howard Wilson. - :Notes: - - For IPDG89 scaling applicability is Z_eff = 1.5, T_i/T_e = 1, 〈T〉 = 5-20 keV + Reference: + - AEA FUS 172: Physics Assessment for the European Reactor Study, 1989 + - H. R. Wilson, Nuclear Fusion 32 (1992) 257 + """ + term1 = np.log(0.5) + term2 = np.log(q0 / q95) - :References: - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', - https://inis.iaea.org/collection/NCLCollectionStore/_Public/21/068/21068960.pdf + # Re-arranging of parabolic profile to be equal to (r/a)^2 where the profile value is half of the core value - - Uckan, N. A., Tolliver, J. S., Houlberg, W. A., and Attenberger, S. E. - Influence of fast alpha diffusion and thermal alpha buildup on tokamak reactor performance. - United States: N. p., 1987. Web.https://www.osti.gov/servlets/purl/5611706 + termp = 1.0 - 0.5 ** (1.0 / alphap) + termt = 1.0 - 0.5 ** (1.0 / alphat) + termj = 1.0 - 0.5 ** (1.0 / alphaj) - """ + # Assuming a parabolic safety factor profile of the form q = q0 + (q95 - q0) * (r/a)^2 + # Substitute (r/a)^2 term from temperature,pressure and current profiles into q profile when values is 50% of core value + # Take natural log of q profile over q95 and q0 to get the profile index - # Determine average fast alpha density - if physics_variables.f_plasma_fuel_deuterium < 1.0: - beta_thermal = ( - 2.0 - * constants.RMU0 - * constants.KILOELECTRON_VOLT - * ( - nd_plasma_electrons_vol_avg - * temp_plasma_electron_density_weighted_kev - + nd_plasma_ions_total_vol_avg * temp_plasma_ion_density_weighted_kev - ) - / (b_plasma_toroidal_on_axis**2 + b_plasma_poloidal_average**2) + alfpnw = term1 / np.log(np.log((q0 + (q95 - q0) * termp) / q95) / term2) + alftnw = term1 / np.log(np.log((q0 + (q95 - q0) * termt) / q95) / term2) + aj = term1 / np.log(np.log((q0 + (q95 - q0) * termj) / q95) / term2) + + # Crude check for NaN errors or other illegal values. + if np.isnan(aj) or np.isnan(alfpnw) or np.isnan(alftnw) or aj < 0: + raise ProcessValueError( + "Illegal profile value found", aj=aj, alfpnw=alfpnw, alftnw=alftnw ) - # jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T) - # IPDG89 fast alpha scaling - if i_beta_fast_alpha == 0: - fact = min( - 0.3, - 0.29 - * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 - * ( - ( - temp_plasma_electron_density_weighted_kev - + temp_plasma_ion_density_weighted_kev - ) - / 20.0 - - 0.37 - ), - ) + # Ratio of ionic charge to electron charge - # Modified scaling, D J Ward - else: - fact = min( - 0.30, - 0.26 - * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 - * np.sqrt( - max( - 0.0, - ( - ( - temp_plasma_electron_density_weighted_kev - + temp_plasma_ion_density_weighted_kev - ) - / 20.0 - - 0.65 - ), - ) - ), - ) + z = 1.0 - fact = max(fact, 0.0) - fact2 = pden_alpha_total_mw / pden_plasma_alpha_mw - beta_fast_alpha = beta_thermal * fact * fact2 + # Inverse aspect ratio: r2 = maximum plasma radius, r1 = minimum + # This is the definition used in the paper + r2 = rmajor + rminor + r1 = rmajor - rminor + eps1 = (r2 - r1) / (r2 + r1) - else: # negligible alpha production, alpha_power_density = p_beam_alpha_mw = 0 - beta_fast_alpha = 0.0 + # Coefficients fitted using least squares techniques - return beta_fast_alpha + # Square root of current profile index term + saj = np.sqrt(aj) - def output_beta_information(self): - """Output beta information to file.""" + a = np.array([ + 1.41 * (1.0 - 0.28 * saj) * (1.0 + 0.12 / z), + 0.36 * (1.0 - 0.59 * saj) * (1.0 + 0.8 / z), + -0.27 * (1.0 - 0.47 * saj) * (1.0 + 3.0 / z), + 0.0053 * (1.0 + 5.0 / z), + -0.93 * (1.0 - 0.34 * saj) * (1.0 + 0.15 / z), + -0.26 * (1.0 - 0.57 * saj) * (1.0 - 0.27 * z), + 0.064 * (1.0 - 0.6 * aj + 0.15 * aj * aj) * (1.0 + 7.6 / z), + -0.0011 * (1.0 + 9.0 / z), + -0.33 * (1.0 - aj + 0.33 * aj * aj), + -0.26 * (1.0 - 0.87 / saj - 0.16 * aj), + -0.14 * (1.0 - 1.14 / saj - 0.45 * saj), + -0.0069, + ]) - po.osubhd(self.outfile, "Beta Information :") - if physics_variables.i_beta_component == 0: - po.ovarrf( - self.outfile, - "Upper limit on total beta", - "(beta_vol_avg_max)", - physics_variables.beta_vol_avg_max, - "OP ", - ) - elif physics_variables.i_beta_component == 1: - po.ovarrf( - self.outfile, - "Upper limit on thermal beta", - "(beta_vol_avg_max)", - physics_variables.beta_vol_avg_max, - "OP ", - ) - else: - po.ovarrf( - self.outfile, - "Upper limit on thermal + NB beta", - "(beta_vol_avg_max)", - physics_variables.beta_vol_avg_max, - "OP ", - ) + seps1 = np.sqrt(eps1) - po.ovarre( - self.outfile, - "Total plasma beta", - "(beta_total_vol_avg)", - physics_variables.beta_total_vol_avg, - ) - if physics_variables.i_beta_component == 0: - po.ovarrf( - self.outfile, - "Lower limit on total beta", - "(beta_vol_avg_min)", - physics_variables.beta_vol_avg_min, - "IP", - ) - elif physics_variables.i_beta_component == 1: - po.ovarrf( - self.outfile, - "Lower limit on thermal beta", - "(beta_vol_avg_min)", - physics_variables.beta_vol_avg_min, - "IP", - ) - else: - po.ovarrf( - self.outfile, - "Lower limit on thermal + NB beta", - "(beta_vol_avg_min)", - physics_variables.beta_vol_avg_min, - "IP", - ) - po.ovarre( - self.outfile, - "Upper limit on poloidal beta", - "(beta_poloidal_max)", - constraint_variables.beta_poloidal_max, - "IP", - ) - po.ovarre( - self.outfile, - "Total poloidal beta", - "(beta_poloidal_vol_avg)", - physics_variables.beta_poloidal_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Volume averaged toroidal beta", - "(beta_toroidal_vol_avg)", - physics_variables.beta_toroidal_vol_avg, - "OP ", - ) - for i in range(len(physics_variables.beta_thermal_toroidal_profile)): - po.ovarre( - self.mfile, - f"Beta toroidal profile at point {i}", - f"beta_thermal_toroidal_profile{i}", - physics_variables.beta_thermal_toroidal_profile[i], - ) + b = np.array([ + 1.0, + alfpnw, + alftnw, + alfpnw * alftnw, + seps1, + alfpnw * seps1, + alftnw * seps1, + alfpnw * alftnw * seps1, + eps1, + alfpnw * eps1, + alftnw * eps1, + alfpnw * alftnw * eps1, + ]) - po.ovarre( - self.outfile, - "Fast alpha beta", - "(beta_fast_alpha)", - physics_variables.beta_fast_alpha, - "OP ", - ) - po.ovarre( - self.outfile, - "Neutral Beam ion beta", - "(beta_beam)", - physics_variables.beta_beam, - "OP ", - ) - po.ovarre( - self.outfile, - "Ratio of fast alpha and beam beta to thermal beta", - "(f_beta_alpha_beam_thermal)", - physics_variables.f_beta_alpha_beam_thermal, - "OP ", - ) + # Empirical bootstrap current fraction + return seps1 * betpth * (a * b).sum() - po.ovarre( - self.outfile, - "Volume averaged thermal beta", - "(beta_thermal_vol_avg)", - physics_variables.beta_thermal_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Thermal poloidal beta", - "(beta_thermal_poloidal_vol_avg)", - physics_variables.beta_thermal_poloidal_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Thermal toroidal beta", - "(beta_thermal_toroidal_vol_avg)", - physics_variables.beta_thermal_toroidal_vol_avg, - "OP ", + @staticmethod + def bootstrap_fraction_nevins( + alphan: float, + alphat: float, + beta_toroidal: float, + b_plasma_toroidal_on_axis: float, + nd_plasma_electrons_vol_avg: float, + plasma_current: float, + q95: float, + q0: float, + rmajor: float, + rminor: float, + te: float, + zeff: float, + ) -> float: + """ + Calculate the bootstrap current fraction from Nevins et al scaling. + + Args: + alphan (float): Density profile index. + alphat (float): Temperature profile index. + beta_toroidal (float): Toroidal plasma beta. + b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). + nd_plasma_electrons_vol_avg (float): Electron density (/m3). + plasma_current (float): Plasma current (A). + q0 (float): Central safety factor. + q95 (float): Safety factor at 95% surface. + rmajor (float): Plasma major radius (m). + rminor (float): Plasma minor radius (m). + te (float): Volume averaged plasma temperature (keV). + zeff (float): Plasma effective charge. + + Returns: + float: The bootstrap current fraction. + + This function calculates the bootstrap current fraction using the Nevins et al method. + + Reference: See appendix of: + Keii Gi, Makoto Nakamura, Kenji Tobita, Yasushi Ono, + Bootstrap current fraction scaling for a tokamak reactor design study, + Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, + ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2014.07.009. + + Nevins, W. M. "Summary report: ITER specialists' meeting on heating and current drive." + ITER-TN-PH-8-4, June 1988. 1988. + + """ + # Calculate peak electron beta at plasma centre, this is not the form used in the paper + # The paper assumes parabolic profiles for calculating core values with the profile indexes. + # We instead use the directly calculated electron density and temperature values at the core. + # So that it is compatible with all profiles + + betae0 = ( + physics_variables.nd_plasma_electron_on_axis + * physics_variables.temp_plasma_electron_on_axis_kev + * 1.0e3 + * constants.ELECTRON_CHARGE + / (b_plasma_toroidal_on_axis**2 / (2.0 * constants.RMU0)) ) - po.ovarrf( - self.outfile, - "Poloidal beta and inverse aspect ratio", - "(beta_poloidal_eps)", - physics_variables.beta_poloidal_eps, - "OP ", + # Call integration routine using definite integral routine from scipy + + ainteg, _ = integrate.quad( + lambda y: _nevins_integral( + y, + nd_plasma_electrons_vol_avg, + te, + b_plasma_toroidal_on_axis, + rminor, + rmajor, + zeff, + alphat, + alphan, + q0, + q95, + beta_toroidal, + ), + 0, # Lower bound + 1.0, # Upper bound ) - po.ovarrf( - self.outfile, - "Poloidal beta and inverse aspect ratio upper limit", - "(beta_poloidal_eps_max)", - physics_variables.beta_poloidal_eps_max, + + # Calculate bootstrap current and fraction + + aibs = 2.5 * betae0 * rmajor * b_plasma_toroidal_on_axis * q95 * ainteg + return 1.0e6 * aibs / plasma_current + + @staticmethod + def bootstrap_fraction_sauter(plasma_profile: float) -> float: + """Calculate the bootstrap current fraction from the Sauter et al scaling. + + Args: + plasma_profile (PlasmaProfile): The plasma profile object. + + Returns: + float: The bootstrap current fraction. + + This function calculates the bootstrap current fraction using the Sauter, Angioni, and Lin-Liu scaling. + + Reference: + - O. Sauter, C. Angioni, Y. R. Lin-Liu; + Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. + Phys. Plasmas 1 July 1999; 6 (7): 2834-2839. https://doi.org/10.1063/1.873240 + - 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 + + Note: + The code was supplied by Emiliano Fable, IPP Garching (private communication). + """ + + # Radial points from 0 to 1 seperated by 1/profile_size + roa = plasma_profile.neprofile.profile_x + + # Local circularised minor radius + rho = np.sqrt(physics_variables.a_plasma_poloidal / np.pi) * roa + + # Square root of local aspect ratio + sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) + + # Calculate electron and ion density profiles + ne = plasma_profile.neprofile.profile_y * 1e-19 + ni = ( + physics_variables.nd_plasma_ions_total_vol_avg + / physics_variables.nd_plasma_electrons_vol_avg + ) * ne + + # Calculate electron and ion temperature profiles + tempe = plasma_profile.teprofile.profile_y + tempi = ( + physics_variables.temp_plasma_ion_vol_avg_kev + / physics_variables.temp_plasma_electron_vol_avg_kev + ) * tempe + + # Flat Zeff profile assumed + # Return tempi like array object filled with zeff + zeff = np.full_like(tempi, physics_variables.n_charge_plasma_effective_vol_avg) + + # inverse_q = 1/safety factor + # Parabolic q profile assumed + inverse_q = 1 / ( + physics_variables.q0 + + (physics_variables.q95 - physics_variables.q0) * roa**2 ) - po.osubhd(self.outfile, "Normalised Beta Information :") - if stellarator_variables.istell == 0: - if physics_variables.i_beta_norm_max != 0: - po.ovarrf( - self.outfile, - "Beta g coefficient", - "(beta_norm_max)", - physics_variables.beta_norm_max, - "OP ", + # Create new array of average mass of fuel portion of ions + amain = np.full_like(inverse_q, physics_variables.m_ions_total_amu) + + # Create new array of average main ion charge + zmain = np.full_like(inverse_q, 1.0 + physics_variables.f_plasma_fuel_helium3) + + # Calculate total bootstrap current (MA) by summing along profiles + # Looping from 2 because _calculate_l31_coefficient() etc should return 0 @ j == 1 + radial_elements = np.arange(2, plasma_profile.profile_size) + + # Change in localised minor radius to be used as delta term in derivative + drho = rho[radial_elements] - rho[radial_elements - 1] + + # Area of annulus, assuming circular plasma cross-section + da = 2 * np.pi * rho[radial_elements - 1] * drho # area of annulus + + # Create the partial derivatives using numpy gradient (central differences) + dlogte_drho = np.gradient(np.log(tempe), rho)[radial_elements - 1] + dlogti_drho = np.gradient(np.log(tempi), rho)[radial_elements - 1] + dlogne_drho = np.gradient(np.log(ne), rho)[radial_elements - 1] + + jboot = ( + 0.5 + * ( + _calculate_l31_coefficient( + radial_elements, + plasma_profile.profile_size, + physics_variables.rmajor, + physics_variables.b_plasma_toroidal_on_axis, + physics_variables.triang, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + zeff, + sqeps, ) - else: - po.ovarrf( - self.outfile, - "Beta g coefficient", - "(beta_norm_max)", - physics_variables.beta_norm_max, + * dlogne_drho + + _calculate_l31_32_coefficient( + radial_elements, + plasma_profile.profile_size, + physics_variables.rmajor, + physics_variables.b_plasma_toroidal_on_axis, + physics_variables.triang, + ne, + ni, + tempe, + tempi, + inverse_q, + rho, + zeff, + sqeps, ) - po.ovarrf( - self.outfile, - "Normalised total beta", - "(beta_norm_total)", - physics_variables.beta_norm_total, - "OP ", + * dlogte_drho + + _calculate_l34_alpha_31_coefficient( + radial_elements, + plasma_profile.profile_size, + physics_variables.rmajor, + physics_variables.b_plasma_toroidal_on_axis, + physics_variables.triang, + inverse_q, + sqeps, + tempi, + tempe, + amain, + zmain, + ni, + ne, + rho, + zeff, + ) + * dlogti_drho ) - po.ovarrf( - self.outfile, - "Normalised thermal beta", - "(beta_norm_thermal) ", - physics_variables.beta_norm_thermal, - "OP ", + * 1.0e6 + * ( + -physics_variables.b_plasma_toroidal_on_axis + / (0.2 * np.pi * physics_variables.rmajor) + * rho[radial_elements - 1] + * inverse_q[radial_elements - 1] ) + ) # A/m2 + + return (np.sum(da * jboot, axis=0) / physics_variables.plasma_current), jboot + + @staticmethod + def bootstrap_fraction_sakai( + beta_poloidal: float, + q95: float, + q0: float, + alphan: float, + alphat: float, + eps: float, + ind_plasma_internal_norm: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Sakai formula. + + Parameters: + beta_poloidal (float): Plasma poloidal beta. + q95 (float): Safety factor at 95% of the plasma radius. + q0 (float): Safety factor at the magnetic axis. + alphan (float): Density profile index + alphat (float): Temperature profile index + eps (float): Inverse aspect ratio. + ind_plasma_internal_norm (float): Plasma normalised internal inductance + + Returns: + float: The calculated bootstrap fraction. + + Notes: + The profile assumed for the alphan and alpat indexes is only a parabolic profile without a pedestal (L-mode). + The Root Mean Squared Error for the fitting database of this formula was 0.025 + Concentrating on the positive shear plasmas using the ACCOME code equilibria with the fully non-inductively driven + conditions with neutral beam (NB) injection only are calculated. + The electron temperature and the ion temperature were assumed to be equal + This can be used for all apsect ratios. + The diamagnetic fraction is included in this formula. + + References: + 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. + """ + # Sakai states that the ACCOME dataset used has the toridal diamagnetic current included in the bootstrap current + # So the diamganetic current should not be calculated with this. i_diamagnetic_current = 0 + return ( + 10 ** (0.951 * eps - 0.948) + * beta_poloidal ** (1.226 * eps + 1.584) + * ind_plasma_internal_norm ** (-0.184 * eps - 0.282) + * (q95 / q0) ** (-0.042 * eps - 0.02) + * alphan ** (0.13 * eps + 0.05) + * alphat ** (0.502 * eps - 0.273) + ) - po.ovarrf( - self.outfile, - "Normalised toroidal beta", - "(beta_norm_toroidal) ", - physics_variables.beta_norm_toroidal, - "OP ", - ) + @staticmethod + def bootstrap_fraction_aries( + beta_poloidal: float, + ind_plasma_internal_norm: float, + core_density: float, + average_density: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the ARIES formula. - po.ovarrf( - self.outfile, - "Normalised poloidal beta", - "(beta_norm_poloidal) ", - physics_variables.beta_norm_poloidal, - "OP ", - ) + Parameters: + beta_poloidal (float): Plasma poloidal beta. + ind_plasma_internal_norm (float): Plasma normalized internal inductance. + core_density (float): Core plasma density. + average_density (float): Average plasma density. + inverse_aspect (float): Inverse aspect ratio. - po.osubhd(self.outfile, "Maximum normalised beta scalings :") - po.ovarrf( - self.outfile, - "J. Wesson normalised beta upper limit", - "(beta_norm_max_wesson) ", - physics_variables.beta_norm_max_wesson, - "OP ", - ) - po.ovarrf( - self.outfile, - "Original normalsied beta upper limit", - "(beta_norm_max_original_scaling) ", - physics_variables.beta_norm_max_original_scaling, - "OP ", - ) - po.ovarrf( - self.outfile, - "J. Menard normalised beta upper limit", - "(beta_norm_max_menard) ", - physics_variables.beta_norm_max_menard, - "OP ", - ) - po.ovarrf( - self.outfile, - "E. Thloreus normalised beta upper limit", - "(beta_norm_max_thloreus) ", - physics_variables.beta_norm_max_thloreus, - "OP ", - ) - po.ovarrf( - self.outfile, - "R. Stambaugh normalised beta upper limit", - "(beta_norm_max_stambaugh) ", - physics_variables.beta_norm_max_stambaugh, - "OP ", - ) + Returns: + float: The calculated bootstrap fraction. - po.osubhd(self.outfile, "Plasma energies derived from beta :") - po.ovarre( - self.outfile, - "Plasma thermal energy derived from thermal beta (J)", - "(e_plasma_beta_thermal) ", - physics_variables.e_plasma_beta_thermal, - "OP", - ) + Notes: + - The source reference does not provide any info about the derivation of the formula. It is only stated - po.ovarre( - self.outfile, - "Plasma thermal energy derived from the total beta (J)", - "(e_plasma_beta)", - physics_variables.e_plasma_beta, - "OP", - ) + 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. - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) + """ + # Using the standard variable naming from the ARIES paper + a_1 = ( + 1.10 - 1.165 * ind_plasma_internal_norm + 0.47 * ind_plasma_internal_norm**2 + ) + b_1 = ( + 0.806 + - 0.885 * ind_plasma_internal_norm + + 0.297 * ind_plasma_internal_norm**2 + ) + c_bs = a_1 + b_1 * (core_density / average_density) -class IndInternalNormModel(IntEnum): - """Normalised internal inductance (l_i) model types""" + return c_bs * np.sqrt(inverse_aspect) * beta_poloidal - USER_INPUT = 0 - WESSON = 1 - MENARD = 2 + @staticmethod + def bootstrap_fraction_andrade( + beta_poloidal: float, + core_pressure: float, + average_pressure: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Andrade et al formula. + Parameters: + beta_poloidal (float): Plasma poloidal beta. + core_pressure (float): Core plasma pressure. + average_pressure (float): Average plasma pressure. + inverse_aspect (float): Inverse aspect ratio. -class PlasmaInductance: - """Class to hold plasma inductance calculations for plasma processing.""" + Returns: + float: The calculated bootstrap fraction. - def __init__(self): - self.outfile = constants.NOUT - self.mfile = constants.MFILE + 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 - def run(self): - physics_variables.ind_plasma_internal_norm_wesson = ( - self.calculate_internal_inductance_wesson(alphaj=physics_variables.alphaj) - ) + 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. - # Spherical Tokamak relation for internal inductance - # Menard et al. (2016), Nuclear Fusion, 56, 106023 - physics_variables.ind_plasma_internal_norm_menard = ( - self.calculate_internal_inductance_menard(kappa=physics_variables.kappa) - ) + """ - physics_variables.ind_plasma_internal_norm_iter_3 = ( - self.calculate_normalised_internal_inductance_iter_3( - b_plasma_poloidal_vol_avg=physics_variables.b_plasma_poloidal_average, - c_plasma=physics_variables.plasma_current, - vol_plasma=physics_variables.vol_plasma, - rmajor=physics_variables.rmajor, - ) - ) + # Using the standard variable naming from the Andrade et.al. paper + c_p = core_pressure / average_pressure - # Calculate ind_plasma_internal_norm based on i_ind_plasma_internal_norm - try: - model = IndInternalNormModel( - int(physics_variables.i_ind_plasma_internal_norm) - ) - physics_variables.ind_plasma_internal_norm = ( - self.get_ind_internal_norm_value(model) - ) - except ValueError: - raise ProcessValueError( - "Illegal value of i_ind_plasma_internal_norm", - i_ind_plasma_internal_norm=physics_variables.i_ind_plasma_internal_norm, - ) from None + # Error +- 0.0007 + c_bs = 0.2340 - def get_ind_internal_norm_value(self, model: IndInternalNormModel) -> float: - """Get the normalised internal inductance (l_i) for the specified model.""" - model_map = { - IndInternalNormModel.USER_INPUT: physics_variables.ind_plasma_internal_norm, - IndInternalNormModel.WESSON: physics_variables.ind_plasma_internal_norm_wesson, - IndInternalNormModel.MENARD: physics_variables.ind_plasma_internal_norm_menard, - } - return model_map[model] + return c_bs * np.sqrt(inverse_aspect) * beta_poloidal * c_p**0.8 @staticmethod - def calculate_volt_second_requirements( - csawth: float, - eps: float, - f_c_plasma_inductive: float, - ejima_coeff: float, - kappa: float, - rmajor: float, - res_plasma: float, - plasma_current: float, - t_plant_pulse_fusion_ramp: float, - t_plant_pulse_burn: float, - ind_plasma_internal_norm: float, - ) -> tuple[float, float, float, float, float, float]: - """Calculate the volt-second requirements and related parameters for plasma physics. + def bootstrap_fraction_hoang( + beta_poloidal: float, + pressure_index: float, + current_index: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Hoang et al formula. - :param csawth: Coefficient for sawteeth effects - :type csawth: float - :param eps: Inverse aspect ratio - :type eps: float - :param f_c_plasma_inductive: Fraction of plasma current produced inductively - :type f_c_plasma_inductive: float - :param ejima_coeff: Ejima coefficient for resistive start-up V-s component - :type ejima_coeff: float - :param kappa: Plasma elongation - :type kappa: float - :param rmajor: Plasma major radius (m) - :type rmajor: float - :param res_plasma: Plasma resistance (ohm) - :type res_plasma: float - :param plasma_current: Plasma current (A) - :type plasma_current: float - :param t_plant_pulse_fusion_ramp: Heating time (s) - :type t_plant_pulse_fusion_ramp: float - :param t_plant_pulse_burn: Burn time (s) - :type t_plant_pulse_burn: float - :param ind_plasma_internal_norm: Plasma normalized internal inductance - :type ind_plasma_internal_norm: float + Parameters: + beta_poloidal (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$ - :return: A tuple containing: - - vs_plasma_internal: Internal plasma volt-seconds (Wb) - - ind_plasma_internal: Plasma inductance (H) - - vs_plasma_burn_required: Volt-seconds needed during flat-top (heat+burn) (Wb) - - vs_plasma_ramp_required: Volt-seconds needed during ramp-up (Wb) - - ind_plasma_total,: Internal and external plasma inductance V-s (Wb) - - vs_res_ramp: Resistive losses in start-up volt-seconds (Wb) - - vs_plasma_total_required: Total volt-seconds needed (Wb) - :rtype: tuple[float, float, float, float, float, float] + # 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 - :notes: + # This could lead to large changes in the value depending on interpretation of the profile index - :references: - - S. Ejima, R. W. Callis, J. L. Luxon, R. D. Stambaugh, T. S. Taylor, and J. C. Wesley, - “Volt-second analysis and consumption in Doublet III plasmas,” - Nuclear Fusion, vol. 22, no. 10, pp. 1313-1319, Oct. 1982, doi: - https://doi.org/10.1088/0029-5515/22/10/006. + c_bs = np.sqrt((pressure_index + 1) / (current_index + 1)) - - S. C. Jardin, C. E. Kessel, and N Pomphrey, - “Poloidal flux linkage requirements for the International Thermonuclear Experimental Reactor,” - Nuclear Fusion, vol. 34, no. 8, pp. 1145-1160, Aug. 1994, - doi: https://doi.org/10.1088/0029-5515/34/8/i07. + return 0.4 * np.sqrt(inverse_aspect) * beta_poloidal**0.9 * c_bs - - S. P. Hirshman and G. H. Neilson, “External inductance of an axisymmetric plasma,” - The Physics of Fluids, vol. 29, no. 3, pp. 790-793, Mar. 1986, - doi: https://doi.org/10.1063/1.865934. + @staticmethod + def bootstrap_fraction_wong( + beta_poloidal: float, + density_index: float, + temperature_index: float, + inverse_aspect: float, + elongation: float, + ) -> float: """ - # Plasma internal inductance + Calculate the bootstrap fraction using the Wong et al formula. - ind_plasma_internal = constants.RMU0 * rmajor * ind_plasma_internal_norm / 2.0 + Parameters: + beta_poloidal (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. - # Internal plasma flux (V-s) component - vs_plasma_internal = ind_plasma_internal * plasma_current + Returns: + float: The calculated bootstrap fraction. - # Start-up resistive component - # Uses ITER formula without the 10 V-s add-on + 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 - vs_res_ramp = ejima_coeff * constants.RMU0 * plasma_current * rmajor + 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) - # Hirshman and Neilson fit for external inductance + c_bs = 0.773 + 0.019 * elongation - aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( - 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps - ) - beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) + return c_bs * f_peak**0.25 * beta_poloidal * np.sqrt(inverse_aspect) - ind_plasma_external = ( - rmajor * constants.RMU0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) - ) + @staticmethod + def bootstrap_fraction_gi_I( # noqa: N802 + beta_poloidal: 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: + beta_poloidal (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. - ind_plasma_total = ind_plasma_external + ind_plasma_internal + Returns: + float: The calculated bootstrap fraction. - # Inductive V-s component + 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. - vs_self_ind_ramp = ind_plasma_total * plasma_current - vs_plasma_ramp_required = vs_res_ramp + vs_self_ind_ramp + 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. + """ - # Plasma loop voltage during flat-top - # Include enhancement factor in flattop V-s requirement - # to account for MHD sawtooth effects. + # Using the standard variable naming from the Gi et.al. paper - v_plasma_loop_burn = plasma_current * res_plasma * f_c_plasma_inductive + c_bs = ( + 0.474 + * inverse_aspect**-0.1 + * pressure_index**0.974 + * temperature_index**-0.416 + * effective_charge**0.178 + * (q95 / q0) ** -0.133 + ) - v_burn_resistive = v_plasma_loop_burn * csawth + return c_bs * np.sqrt(inverse_aspect) * beta_poloidal - # N.B. t_plant_pulse_burn on first iteration will not be correct - # if the pulsed reactor option is used, but the value - # will be correct on subsequent calls. + @staticmethod + def bootstrap_fraction_gi_II( # noqa: N802 + beta_poloidal: 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. - vs_plasma_burn_required = v_burn_resistive * ( - t_plant_pulse_fusion_ramp + t_plant_pulse_burn - ) - vs_plasma_total_required = vs_plasma_ramp_required + vs_plasma_burn_required + Parameters: + beta_poloidal (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. - return ( - vs_plasma_internal, - ind_plasma_total, - vs_plasma_burn_required, - vs_plasma_ramp_required, - vs_self_ind_ramp, - vs_res_ramp, - vs_plasma_total_required, - v_plasma_loop_burn, + 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) * beta_poloidal + @staticmethod - def calculate_normalised_internal_inductance_iter_3( - b_plasma_poloidal_vol_avg: float, - c_plasma: float, - vol_plasma: float, - rmajor: float, + def bootstrap_fraction_sugiyama_l_mode( + eps: float, + beta_poloidal: float, + alphan: float, + alphat: float, + zeff: float, + q95: float, + q0: float, ) -> float: """ - Calculate the normalised internal inductance using ITER-3 scaling li(3). + Calculate the bootstrap fraction using the L-mode scaling from the Sugiyama et al formula. - :param b_plasma_poloidal_vol_avg: Volume-averaged poloidal magnetic field (T). - :type b_plasma_poloidal_vol_avg: float - :param c_plasma: Plasma current (A). - :type c_plasma: float - :param vol_plasma: Plasma volume (m^3). - :type vol_plasma: float - :param rmajor: Plasma major radius (m). - :type rmajor: float + :param eps: Inverse aspect ratio. + :type eps: float + :param beta_poloidal: Plasma poloidal beta. + :type beta_poloidal: float + :param alphan: Density profile index. + :type alphan: float + :param alphat: Temperature profile index. + :type alphat: float + :param zeff: Plasma effective charge. + :type zeff: float + :param q95: Safety factor at 95% of the plasma radius. + :type q95: float + :param q0: Safety factor at the magnetic axis. + :type q0: float - :returns: The li(3) normalised internal inductance. + :returns: The calculated bootstrap fraction. :rtype: float - :references: - - T. C. Luce, D. A. Humphreys, G. L. Jackson, and W. M. Solomon, - “Inductive flux usage and its optimization in tokamak operation,” - Nuclear Fusion, vol. 54, no. 9, p. 093005, Jul. 2014, - doi: https://doi.org/10.1088/0029-5515/54/9/093005. + :notes: + - This scaling is derived for L-mode plasmas. + - Ion and electron temperature are the same + - Z_eff has a uniform profile, with only fully stripped carbon impurity - - G. L. Jackson et al., “ITER startup studies in the DIII-D tokamak,” - Nuclear Fusion, vol. 48, no. 12, p. 125002, Nov. 2008, - doi: https://doi.org/10.1088/0029-5515/48/12/125002. - ‌ + :references: + - S. Sugiyama, T. Goto, H. Utoh, and Y. Sakamoto, “Improvement of core plasma power and + current balance models for tokamak systems code considering H-mode plasma profiles,” + Fusion Engineering and Design, vol. 216, p. 115022, Jul. 2025, doi: + https://doi.org/10.1016/j.fusengdes.2025.115022. """ return ( - 2 - * vol_plasma - * b_plasma_poloidal_vol_avg**2 - / (constants.RMU0**2 * c_plasma**2 * rmajor) + 0.740 + * eps**0.418 + * beta_poloidal**0.904 + * alphan**0.06 + * alphat**-0.138 + * zeff**0.230 + * (q95 / q0) ** -0.142 ) @staticmethod - def calculate_internal_inductance_menard(kappa: float) -> float: + def bootstrap_fraction_sugiyama_h_mode( + eps: float, + beta_poloidal: float, + alphan: float, + alphat: float, + tbeta: float, + zeff: float, + q95: float, + q0: float, + radius_plasma_pedestal_density_norm: float, + nd_plasma_pedestal_electron: float, + n_greenwald: float, + temp_plasma_pedestal_kev: float, + ) -> float: """ - Calculate the Menard plasma normalized internal inductance. + Calculate the bootstrap fraction using the H-mode scaling from the Sugiyama et al formula. - :param kappa: Plasma separatrix elongation. - :type kappa: float + :param eps: Inverse aspect ratio. + :type eps: float + :param beta_poloidal: Plasma poloidal beta. + :type beta_poloidal: float + :param alphan: Density profile index. + :type alphan: float + :param alphat: Temperature profile index. + :type alphat: float + :param tbeta: Second temperature profile index. + :type tbeta: float + :param zeff: Plasma effective charge. + :type zeff: float + :param q95: Safety factor at 95% of the plasma radius. + :type q95: float + :param q0: Safety factor at the magnetic axis. + :type q0: float + :param radius_plasma_pedestal_density_norm: Normalised plasma radius of density pedestal. + :type radius_plasma_pedestal_density_norm: float + :param nd_plasma_pedestal_electron: Electron number density at the pedestal [m^-3]. + :type nd_plasma_pedestal_electron: float + :param n_greenwald: Greenwald density limit [m^-3]. + :type n_greenwald: float + :param temp_plasma_pedestal_kev: Electron temperature at the pedestal [keV]. + :type temp_plasma_pedestal_kev: float - :return: The Menard plasma normalised internal inductance. + :returns: The calculated bootstrap fraction. :rtype: float - :Notes: - - This relation is based off of data from NSTX for l_i in the range of 0.4-0.85 - - This model is only recommneded to be used for ST's with kappa > 2.5 - - :References: - - J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” - Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, - doi: https://doi.org/10.1088/0029-5515/56/10/106023. - """ - return 3.4 - kappa + :notes: + - This scaling is derived for H-mode plasmas. + - The temperature and density pedestal positions are the same + - Separatrix temperature and density are zero + - Ion and electron temperature are the same + - Z_eff has a uniform profile, with only fully stripped carbon impurity - @staticmethod - def calculate_internal_inductance_wesson(alphaj: float) -> float: + :references: + - S. Sugiyama, T. Goto, H. Utoh, and Y. Sakamoto, “Improvement of core plasma power and + current balance models for tokamak systems code considering H-mode plasma profiles,” + Fusion Engineering and Design, vol. 216, p. 115022, Jul. 2025, doi: + https://doi.org/10.1016/j.fusengdes.2025.115022. """ - Calculate the Wesson plasma normalized internal inductance. - - :param alphaj: Current profile index. - :type alphaj: float - :return: The Wesson plasma normalised internal inductance. - :rtype: float + return ( + 0.789 + * eps**0.606 + * beta_poloidal**0.960 + * alphan**0.0319 + * alphat**0.00822 + * tbeta**-0.0783 + * zeff**0.241 + * (q95 / q0) ** -0.103 + * radius_plasma_pedestal_density_norm**0.367 + * (nd_plasma_pedestal_electron / n_greenwald) ** -0.174 + * temp_plasma_pedestal_kev**0.0552 + ) - :Notes: - - It is recommended to use this method with the other Wesson relations for normalised beta and - current profile index. - - This relation is only true for the cyclindrical plasma approximation with parabolic profiles. + def output_bootstrap_current_information(self): + """Output the calculated bootstrap current information to the output file.""" - :References: - - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, - International Series of Monographs on Physics, Volume 149. - """ - return np.log(1.65 + 0.89 * alphaj) + po.ovarrf( + self.outfile, + "Bootstrap current fraction multiplier", + "(cboot)", + current_drive_variables.cboot, + ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (ITER 1989)", + "(f_c_plasma_bootstrap_iter89)", + current_drive_variables.f_c_plasma_bootstrap_iter89, + "OP ", + ) - def output_volt_second_information(self): - """Output volt-second information to file.""" + po.ovarrf( + self.outfile, + "Bootstrap fraction (Sauter et al)", + "(f_c_plasma_bootstrap_sauter)", + current_drive_variables.f_c_plasma_bootstrap_sauter, + "OP ", + ) + for point in range(len(physics_variables.j_plasma_bootstrap_sauter_profile)): + po.ovarrf( + self.mfile, + f"Sauter et al bootstrap current density profile at point {point}", + f"(j_plasma_bootstrap_sauter_profile{point})", + physics_variables.j_plasma_bootstrap_sauter_profile[point], + "OP ", + ) - po.osubhd(self.outfile, "Plasma Volt-second Requirements :") - po.ovarre( + po.ovarrf( self.outfile, - "Total plasma volt-seconds required for pulse (Wb)", - "(vs_plasma_total_required)", - physics_variables.vs_plasma_total_required, + "Bootstrap fraction (Nevins et al)", + "(f_c_plasma_bootstrap_nevins)", + current_drive_variables.f_c_plasma_bootstrap_nevins, "OP ", ) - po.oblnkl(self.outfile) - po.ovarre( + po.ovarrf( self.outfile, - "Total plasma inductive flux consumption for plasma current ramp-up (Wb)", - "(vs_plasma_ind_ramp)", - physics_variables.vs_plasma_ind_ramp, + "Bootstrap fraction (Wilson)", + "(f_c_plasma_bootstrap_wilson)", + current_drive_variables.f_c_plasma_bootstrap_wilson, "OP ", ) - po.ovarre( + po.ovarrf( self.outfile, - "Plasma resistive flux consumption for plasma current ramp-up (Wb)", - "(vs_plasma_res_ramp)", - physics_variables.vs_plasma_res_ramp, + "Bootstrap fraction (Sakai)", + "(f_c_plasma_bootstrap_sakai)", + current_drive_variables.f_c_plasma_bootstrap_sakai, "OP ", ) - po.ovarre( + po.ovarrf( self.outfile, - "Total flux consumption for plasma current ramp-up (Wb)", - "(vs_plasma_ramp_required)", - physics_variables.vs_plasma_ramp_required, + "Bootstrap fraction (ARIES)", + "(f_c_plasma_bootstrap_aries)", + current_drive_variables.f_c_plasma_bootstrap_aries, "OP ", ) po.ovarrf( self.outfile, - "Ejima coefficient", - "(ejima_coeff)", - physics_variables.ejima_coeff, + "Bootstrap fraction (Andrade)", + "(f_c_plasma_bootstrap_andrade)", + current_drive_variables.f_c_plasma_bootstrap_andrade, + "OP ", ) - po.oblnkl(self.outfile) - po.ovarre( + po.ovarrf( self.outfile, - "Internal plasma V-s", - "(vs_plasma_internal)", - physics_variables.vs_plasma_internal, + "Bootstrap fraction (Hoang)", + "(f_c_plasma_bootstrap_hoang)", + current_drive_variables.f_c_plasma_bootstrap_hoang, + "OP ", ) - - po.ovarre( + po.ovarrf( self.outfile, - "Plasma volt-seconds needed for flat-top (heat + burn times) (Wb)", - "(vs_plasma_burn_required)", - physics_variables.vs_plasma_burn_required, + "Bootstrap fraction (Wong)", + "(f_c_plasma_bootstrap_wong)", + current_drive_variables.f_c_plasma_bootstrap_wong, "OP ", ) - - po.ovarre( + po.ovarrf( self.outfile, - "Plasma loop voltage during burn (V)", - "(v_plasma_loop_burn)", - physics_variables.v_plasma_loop_burn, + "Bootstrap fraction (Gi I)", + "(bscf_gi_i)", + current_drive_variables.bscf_gi_i, "OP ", ) po.ovarrf( self.outfile, - "Coefficient for sawtooth effects on burn V-s requirement", - "(csawth)", - physics_variables.csawth, + "Bootstrap fraction (Gi II)", + "(bscf_gi_ii)", + current_drive_variables.bscf_gi_ii, + "OP ", ) - po.oblnkl(self.outfile) - po.ovarre( + po.ovarrf( self.outfile, - "Plasma resistance (ohm)", - "(res_plasma)", - physics_variables.res_plasma, + "Bootstrap fraction (Sugiyama L-mode)", + "(f_c_plasma_bootstrap_sugiyama_l)", + current_drive_variables.f_c_plasma_bootstrap_sugiyama_l, "OP ", ) - - po.ovarre( + po.ovarrf( self.outfile, - "Plasma resistive diffusion time (s)", - "(t_plasma_res_diffusion)", - physics_variables.t_plasma_res_diffusion, + "Bootstrap fraction (Sugiyama H-mode)", + "(f_c_plasma_bootstrap_sugiyama_h)", + current_drive_variables.f_c_plasma_bootstrap_sugiyama_h, "OP ", ) - po.ovarre( + + po.ovarrf( self.outfile, - "Plasma inductance (H)", - "(ind_plasma)", - physics_variables.ind_plasma, + "Diamagnetic fraction (Hender)", + "(f_c_plasma_diamagnetic_hender)", + current_drive_variables.f_c_plasma_diamagnetic_hender, "OP ", ) - po.ovarre( + po.ovarrf( self.outfile, - "Plasma magnetic energy stored (J)", - "(e_plasma_magnetic_stored)", - physics_variables.e_plasma_magnetic_stored, + "Diamagnetic fraction (SCENE)", + "(f_c_plasma_diamagnetic_scene)", + current_drive_variables.f_c_plasma_diamagnetic_scene, "OP ", ) po.ovarrf( self.outfile, - "Plasma normalised internal inductance", - "(ind_plasma_internal_norm)", - physics_variables.ind_plasma_internal_norm, + "Pfirsch-Schlueter fraction (SCENE)", + "(f_c_plasma_pfirsch_schluter_scene)", + current_drive_variables.f_c_plasma_pfirsch_schluter_scene, "OP ", ) + # Error to catch if bootstap fraction limit has been enforced + if physics_variables.err242 == 1: + logger.error("Bootstrap fraction upper limit enforced") - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) + # Error to catch if self-driven current fraction limit has been enforced + if physics_variables.err243 == 1: + logger.error( + "Predicted plasma driven current is more than upper limit on non-inductive fraction" + ) + + if physics_variables.i_bootstrap_current == 0: + po.ocmmnt(self.outfile, " (User-specified bootstrap current fraction used)") + elif physics_variables.i_bootstrap_current == 1: + po.ocmmnt( + self.outfile, " (ITER 1989 bootstrap current fraction model used)" + ) + elif physics_variables.i_bootstrap_current == 2: + po.ocmmnt( + self.outfile, + " (Nevins et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 3: + po.ocmmnt(self.outfile, " (Wilson bootstrap current fraction model used)") + elif physics_variables.i_bootstrap_current == 4: + po.ocmmnt( + self.outfile, + " (Sauter et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 5: + po.ocmmnt( + 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(self.outfile, " (Diamagnetic current fraction not calculated)") + # Error to show if diamagnetic current is above 1% but not used + if current_drive_variables.f_c_plasma_diamagnetic_scene > 0.01e0: + logger.error( + "Diamagnetic fraction is more than 1%, but not calculated. " + "Consider using i_diamagnetic_current=2 and i_pfirsch_schluter_current=1" + ) + + elif physics_variables.i_diamagnetic_current == 1: + po.ocmmnt( + self.outfile, " (Hender diamagnetic current fraction scaling used)" + ) + elif physics_variables.i_diamagnetic_current == 2: + po.ocmmnt( + self.outfile, " (SCENE diamagnetic current fraction scaling used)" + ) + + if physics_variables.i_pfirsch_schluter_current == 0: + po.ocmmnt(self.outfile, " Pfirsch-Schluter current fraction not calculated") + elif physics_variables.i_pfirsch_schluter_current == 1: + po.ocmmnt( + self.outfile, + " (SCENE Pfirsch-Schluter current fraction scaling used)", + ) + + po.ovarrf( + self.outfile, + "Bootstrap fraction (enforced)", + "(f_c_plasma_bootstrap.)", + current_drive_variables.f_c_plasma_bootstrap, + "OP ", + ) + po.ovarrf( + self.outfile, + "Diamagnetic fraction (enforced)", + "(f_c_plasma_diamagnetic.)", + current_drive_variables.f_c_plasma_diamagnetic, + "OP ", + ) + po.ovarrf( + self.outfile, + "Pfirsch-Schlueter fraction (enforced)", + "(f_c_plasma_pfirsch_schluter.)", + current_drive_variables.f_c_plasma_pfirsch_schluter, + "OP ", + ) class DetailedPhysics: diff --git a/process/stellarator.py b/process/stellarator.py index ed5118b12..574f947bd 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -69,6 +69,7 @@ def __init__( neoclassics, plasma_beta, plasma_inductance, + plasma_bootstrap, ) -> None: """Initialises the Stellarator model's variables @@ -107,6 +108,7 @@ def __init__( self.neoclassics = neoclassics self.beta = plasma_beta self.inductance = plasma_inductance + self.bootstrap = plasma_bootstrap def run(self, output: bool): """Routine to call the physics and engineering modules diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 9992a24a7..48bfb8f7b 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -25,6 +25,7 @@ Physics, PlasmaBeta, PlasmaInductance, + PlasmaBootstrapCurrent, calculate_current_coefficient_hastie, calculate_plasma_current_peng, calculate_poloidal_field, @@ -56,6 +57,7 @@ def physics(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaBootstrapCurrent(), ) @@ -143,7 +145,7 @@ def test_bootstrap_fraction_iter89(bootstrapfractioniter89param, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - f_c_plasma_bootstrap = physics.bootstrap_fraction_iter89( + f_c_plasma_bootstrap = PlasmaBootstrapCurrent().bootstrap_fraction_iter89( aspect=bootstrapfractioniter89param.aspect, beta=bootstrapfractioniter89param.beta, b_plasma_toroidal_on_axis=bootstrapfractioniter89param.b_plasma_toroidal_on_axis, @@ -238,7 +240,7 @@ def test_bootstrap_fraction_nevins(bootstrapfractionnevinsparam, monkeypatch, ph bootstrapfractionnevinsparam.nd_plasma_electron_on_axis, ) - fibs = physics.bootstrap_fraction_nevins( + fibs = PlasmaBootstrapCurrent().bootstrap_fraction_nevins( alphan=bootstrapfractionnevinsparam.alphan, alphat=bootstrapfractionnevinsparam.alphat, beta_toroidal=bootstrapfractionnevinsparam.beta_toroidal, @@ -316,7 +318,7 @@ def test_bootstrap_fraction_wilson(bootstrapfractionwilsonparam, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - bfw = physics.bootstrap_fraction_wilson( + bfw = PlasmaBootstrapCurrent().bootstrap_fraction_wilson( alphaj=bootstrapfractionwilsonparam.alphaj, alphap=bootstrapfractionwilsonparam.alphap, alphat=bootstrapfractionwilsonparam.alphat, @@ -555,7 +557,7 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph monkeypatch.setattr(physics_variables, "alphat", bootstrapfractionsauterparam.alphat) physics.plasma_profile.run() - bfs, _ = physics.bootstrap_fraction_sauter(physics.plasma_profile) + bfs, _ = PlasmaBootstrapCurrent().bootstrap_fraction_sauter(physics.plasma_profile) assert bfs == pytest.approx(bootstrapfractionsauterparam.expected_bfs) @@ -639,7 +641,7 @@ def test_bootstrap_fraction_sakai(bootstrapfractionsakaiparam, monkeypatch, phys bootstrapfractionsakaiparam.ind_plasma_internal_norm, ) - bfs = physics.bootstrap_fraction_sakai( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_sakai( beta_poloidal=bootstrapfractionsakaiparam.beta_poloidal, q95=bootstrapfractionsakaiparam.q95, q0=bootstrapfractionsakaiparam.q0, @@ -689,7 +691,7 @@ def test_bootstrap_fraction_aries(bootstrapfractionariesparam, physics): :type bootstrapfractionsauterparam: bootstrapfractionsauterparam """ - bfs = physics.bootstrap_fraction_aries( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_aries( beta_poloidal=bootstrapfractionariesparam.beta_poloidal, ind_plasma_internal_norm=bootstrapfractionariesparam.ind_plasma_internal_norm, core_density=bootstrapfractionariesparam.core_density, @@ -734,7 +736,7 @@ def test_bootstrap_fraction_andrade(bootstrapfractionandradeparam, physics): :type bootstrapfractionsauterparam: bootstrapfractionsauterparam """ - bfs = physics.bootstrap_fraction_andrade( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_andrade( beta_poloidal=bootstrapfractionandradeparam.beta_poloidal, core_pressure=bootstrapfractionandradeparam.core_pressure, average_pressure=bootstrapfractionandradeparam.average_pressure, @@ -778,7 +780,7 @@ def test_bootstrap_fraction_hoang(bootstrapfractionhoangparam, physics): :type bootstrapfractionsauterparam: bootstrapfractionsauterparam """ - bfs = physics.bootstrap_fraction_hoang( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_hoang( beta_poloidal=bootstrapfractionhoangparam.beta_poloidal, pressure_index=bootstrapfractionhoangparam.pressure_index, current_index=bootstrapfractionhoangparam.current_index, @@ -825,7 +827,7 @@ def test_bootstrap_fraction_wong(bootstrapfractionwongparam, physics): :type bootstrapfractionsauterparam: bootstrapfractionsauterparam """ - bfs = physics.bootstrap_fraction_wong( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_wong( beta_poloidal=bootstrapfractionwongparam.beta_poloidal, density_index=bootstrapfractionwongparam.density_index, temperature_index=bootstrapfractionwongparam.temperature_index, @@ -879,7 +881,7 @@ def test_bootstrap_fraction_gi_I(bootstrapfractiongiiparam, physics): # noqa: N :type bootstrapfractionsauterparam: bootstrapfractionsauterparam """ - bfs = physics.bootstrap_fraction_gi_I( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_gi_I( beta_poloidal=bootstrapfractiongiiparam.beta_poloidal, pressure_index=bootstrapfractiongiiparam.pressure_index, temperature_index=bootstrapfractiongiiparam.temperature_index, @@ -929,7 +931,7 @@ def test_bootstrap_fraction_gi_II(bootstrapfractiongiiiparam, physics): # noqa: :type bootstrapfractionsauterparam: bootstrapfractionsauterparam """ - bfs = physics.bootstrap_fraction_gi_II( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_gi_II( beta_poloidal=bootstrapfractiongiiiparam.beta_poloidal, pressure_index=bootstrapfractiongiiiparam.pressure_index, temperature_index=bootstrapfractiongiiiparam.temperature_index, @@ -986,7 +988,7 @@ def test_bootstrap_fraction_sugiyama_l_mode(bootstrapfractionsugiyamalparam, phy :param bootstrapfractionsugiyamalparam: Parameters for the test case. :type bootstrapfractionsugiyamalparam: BootstrapFractionSugiyamaLModeParam """ - bfs = physics.bootstrap_fraction_sugiyama_l_mode( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_sugiyama_l_mode( eps=bootstrapfractionsugiyamalparam.eps, beta_poloidal=bootstrapfractionsugiyamalparam.beta_poloidal, alphan=bootstrapfractionsugiyamalparam.alphan, @@ -1092,7 +1094,7 @@ def test_bootstrap_fraction_sugiyama_h_mode(bootstrapfractionsugiyamahparam, phy :param bootstrapfractionsugiyamahparam: Parameters for the test case. :type bootstrapfractionsugiyamahparam: BootstrapFractionSugiyamaHModeParam """ - bfs = physics.bootstrap_fraction_sugiyama_h_mode( + bfs = PlasmaBootstrapCurrent().bootstrap_fraction_sugiyama_h_mode( eps=bootstrapfractionsugiyamahparam.eps, beta_poloidal=bootstrapfractionsugiyamahparam.beta_poloidal, alphan=bootstrapfractionsugiyamahparam.alphan, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index d3f38585e..57c5f0cea 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -29,7 +29,7 @@ ) from process.fw import FirstWall from process.hcpb import CCFE_HCPB -from process.physics import Physics, PlasmaBeta, PlasmaInductance +from process.physics import Physics, PlasmaBeta, PlasmaInductance, PlasmaBootstrapCurrent from process.plasma_profiles import PlasmaProfile from process.power import Power from process.stellarator import Neoclassics, Stellarator @@ -71,10 +71,12 @@ def stellarator(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaBootstrapCurrent(), ), Neoclassics(), plasma_beta=PlasmaBeta(), plasma_inductance=PlasmaInductance(), + plasma_bootstrap=PlasmaBootstrapCurrent(), )